diff options
| -rw-r--r-- | geometry.h | 2 | ||||
| -rw-r--r-- | point.c | 10 | ||||
| -rw-r--r-- | quaternion.c | 18 |
3 files changed, 24 insertions, 6 deletions
@@ -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 */ @@ -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 θ) { |