From fdeaa2423c582df8d9e9ae6b0e26d386b5531262 Mon Sep 17 00:00:00 2001
From: rodri <rgl@antares-labs.eu>
Date: Fri, 27 Jan 2023 22:57:21 +0000
Subject: bring the slerp function, make ptinpoly public and clean some of the
 code.

---
 geometry.h   |  2 ++
 point.c      | 10 ++++------
 quaternion.c | 18 ++++++++++++++++++
 3 files changed, 24 insertions(+), 6 deletions(-)

diff --git a/geometry.h b/geometry.h
index 7b4402d..efe1abe 100644
--- a/geometry.h
+++ b/geometry.h
@@ -56,6 +56,7 @@ Point2 lerp2(Point2, Point2, double);
 double dotvec2(Point2, Point2);
 double vec2len(Point2);
 Point2 normvec2(Point2);
+int ptinpoly(Point2, Point2*, ulong);
 
 /* Point3 */
 Point3 Pt3(double, double, double, double);
@@ -108,6 +109,7 @@ double dotq(Quaternion, Quaternion);
 Quaternion invq(Quaternion);
 double qlen(Quaternion);
 Quaternion normq(Quaternion);
+Quaternion slerp(Quaternion, Quaternion, double);
 Point3 qrotate(Point3, Point3, double);
 
 /* RFrame */
diff --git a/point.c b/point.c
index e891b02..91634f0 100644
--- a/point.c
+++ b/point.c
@@ -43,8 +43,7 @@ divpt2(Point2 p, double s)
 Point2
 lerp2(Point2 a, Point2 b, double t)
 {
-	if(t < 0) t = 0;
-	if(t > 1) t = 1;
+	t = fclamp(t, 0, 1);
 	return Pt2(
 		flerp(a.x, b.x, t),
 		flerp(a.y, b.y, t),
@@ -78,7 +77,7 @@ normvec2(Point2 v)
 /*
  * the edge function, from:
  *
- * Juan Pineda, A Parallel Algorithm for Polygon Rasterization,
+ * Juan Pineda, “A Parallel Algorithm for Polygon Rasterization”,
  * Computer Graphics, Vol. 22, No. 8, August 1988
  */
 int
@@ -97,7 +96,7 @@ edgeptcmp(Point2 e0, Point2 e1, Point2 p)
 }
 
 /*
- * algorithm by W. Randolph Franklin
+ * (PNPOLY) algorithm by W. Randolph Franklin
  */
 int
 ptinpoly(Point2 p, Point2 *pts, ulong npts)
@@ -152,8 +151,7 @@ divpt3(Point3 p, double s)
 Point3
 lerp3(Point3 a, Point3 b, double t)
 {
-	if(t < 0) t = 0;
-	if(t > 1) t = 1;
+	t = fclamp(t, 0, 1);
 	return Pt3(
 		flerp(a.x, b.x, t),
 		flerp(a.y, b.y, t),
diff --git a/quaternion.c b/quaternion.c
index 686c04e..ca44747 100644
--- a/quaternion.c
+++ b/quaternion.c
@@ -78,6 +78,24 @@ normq(Quaternion q)
 	return sdivq(q, qlen(q));
 }
 
+/*
+ * based on the implementation from:
+ *
+ * Jonathan Blow, “Understanding Slerp, Then Not Using it”,
+ * The Inner Product, April 2004.
+ */
+Quaternion
+slerp(Quaternion q, Quaternion r, double t)
+{
+	Quaternion v;
+	double θ, q·r;
+
+	q·r = fclamp(dotq(q, r), -1, 1); /* stay within the domain of acos(2) */
+	θ = acos(q·r)*t;
+	v = normq(subq(r, smulq(q, q·r))); /* v = r - (q·r)q / |v| */
+	return addq(smulq(q, cos(θ)), smulq(v, sin(θ))); /* q cos(θ) + v sin(θ) */
+}
+
 Point3
 qrotate(Point3 p, Point3 axis, double θ)
 {
-- 
cgit v1.2.3