shithub: front

Download patch

ref: 65ea63e52174e262884b751b4d9a86b54092e559
parent: 84170d89b70bf328df3d6be192390614d2938993
author: rodri <rgl@antares-labs.eu>
date: Mon Oct 14 07:14:50 EDT 2024

libgeometry: add quaternion sandwich product functions

also took the chance to remove the semicolons from
the function signatures in the manpage synopsis, as
is custom in the other manuals.

--- a/sys/include/geometry.h
+++ b/sys/include/geometry.h
@@ -122,6 +122,8 @@
 double qlen(Quaternion);
 Quaternion normq(Quaternion);
 Quaternion slerp(Quaternion, Quaternion, double);
+Quaternion qsandwich(Quaternion, Quaternion);
+Point3 qsandwichpt3(Quaternion, Point3);
 Point3 qrotate(Point3, Point3, double);
 
 /* RFrame */
--- a/sys/man/2/geometry
+++ b/sys/man/2/geometry
@@ -56,107 +56,109 @@
 };
 .PB
 /* utils */
-double flerp(double a, double b, double t);
-double fberp(double a, double b, double c, Point3 bc);
-double fclamp(double n, double min, double max);
+double flerp(double a, double b, double t)
+double fberp(double a, double b, double c, Point3 bc)
+double fclamp(double n, double min, double max)
 .PB
 /* Point2 */
-Point2 Pt2(double x, double y, double w);
-Point2 Vec2(double x, double y);
-Point2 addpt2(Point2 a, Point2 b);
-Point2 subpt2(Point2 a, Point2 b);
-Point2 mulpt2(Point2 p, double s);
-Point2 divpt2(Point2 p, double s);
-Point2 lerp2(Point2 a, Point2 b, double t);
-Point2 berp2(Point2 a, Point2 b, Point2 c, Point3 bc);
-double dotvec2(Point2 a, Point2 b);
-double vec2len(Point2 v);
-Point2 normvec2(Point2 v);
-int edgeptcmp(Point2 e0, Point2 e1, Point2 p);
+Point2 Pt2(double x, double y, double w)
+Point2 Vec2(double x, double y)
+Point2 addpt2(Point2 a, Point2 b)
+Point2 subpt2(Point2 a, Point2 b)
+Point2 mulpt2(Point2 p, double s)
+Point2 divpt2(Point2 p, double s)
+Point2 lerp2(Point2 a, Point2 b, double t)
+Point2 berp2(Point2 a, Point2 b, Point2 c, Point3 bc)
+double dotvec2(Point2 a, Point2 b)
+double vec2len(Point2 v)
+Point2 normvec2(Point2 v)
+int edgeptcmp(Point2 e0, Point2 e1, Point2 p)
 int ptinpoly(Point2 p, Point2 *pts, ulong npts)
 .PB
 /* Point3 */
-Point3 Pt3(double x, double y, double z, double w);
-Point3 Vec3(double x, double y, double z);
-Point3 addpt3(Point3 a, Point3 b);
-Point3 subpt3(Point3 a, Point3 b);
-Point3 mulpt3(Point3 p, double s);
-Point3 divpt3(Point3 p, double s);
-Point3 lerp3(Point3 a, Point3 b, double t);
-Point3 berp3(Point3 a, Point3 b, Point3 c, Point3 bc);
-double dotvec3(Point3 a, Point3 b);
-Point3 crossvec3(Point3 a, Point3 b);
-double vec3len(Point3 v);
-Point3 normvec3(Point3 v);
+Point3 Pt3(double x, double y, double z, double w)
+Point3 Vec3(double x, double y, double z)
+Point3 addpt3(Point3 a, Point3 b)
+Point3 subpt3(Point3 a, Point3 b)
+Point3 mulpt3(Point3 p, double s)
+Point3 divpt3(Point3 p, double s)
+Point3 lerp3(Point3 a, Point3 b, double t)
+Point3 berp3(Point3 a, Point3 b, Point3 c, Point3 bc)
+double dotvec3(Point3 a, Point3 b)
+Point3 crossvec3(Point3 a, Point3 b)
+double vec3len(Point3 v)
+Point3 normvec3(Point3 v)
 int lineXsphere(Point3 *rp, Point3 p0, Point3 p1,
-		Point3 c, double rad, int isaray);
+		Point3 c, double rad, int isaray)
 .PB
 /* Matrix */
-void identity(Matrix m);
-void addm(Matrix a, Matrix b);
-void subm(Matrix a, Matrix b);
-void mulm(Matrix a, Matrix b);
-void smulm(Matrix m, double s);
-void transposem(Matrix m);
-double detm(Matrix m);
-double tracem(Matrix m);
-double minorm(Matrix m, int row, int col);
-double cofactorm(Matrix m, int row, int col);
-void adjm(Matrix m);
-void invm(Matrix m);
-Point2 xform(Point2 p, Matrix m);
+void identity(Matrix m)
+void addm(Matrix a, Matrix b)
+void subm(Matrix a, Matrix b)
+void mulm(Matrix a, Matrix b)
+void smulm(Matrix m, double s)
+void transposem(Matrix m)
+double detm(Matrix m)
+double tracem(Matrix m)
+double minorm(Matrix m, int row, int col)
+double cofactorm(Matrix m, int row, int col)
+void adjm(Matrix m)
+void invm(Matrix m)
+Point2 xform(Point2 p, Matrix m)
 .PB
 /* Matrix3 */
-void identity3(Matrix3 m);
-void addm3(Matrix3 a, Matrix3 b);
-void subm3(Matrix3 a, Matrix3 b);
-void mulm3(Matrix3 a, Matrix3 b);
-void smulm3(Matrix3 m, double s);
-void transposem3(Matrix3 m);
-double detm3(Matrix3 m);
-double tracem3(Matrix3 m);
-double minorm3(Matrix3 m, int row, int col);
-double cofactorm3(Matrix3 m, int row, int col);
-void adjm3(Matrix3 m);
-void invm3(Matrix3 m);
-Point3 xform3(Point3 p, Matrix3 m);
+void identity3(Matrix3 m)
+void addm3(Matrix3 a, Matrix3 b)
+void subm3(Matrix3 a, Matrix3 b)
+void mulm3(Matrix3 a, Matrix3 b)
+void smulm3(Matrix3 m, double s)
+void transposem3(Matrix3 m)
+double detm3(Matrix3 m)
+double tracem3(Matrix3 m)
+double minorm3(Matrix3 m, int row, int col)
+double cofactorm3(Matrix3 m, int row, int col)
+void adjm3(Matrix3 m)
+void invm3(Matrix3 m)
+Point3 xform3(Point3 p, Matrix3 m)
 .PB
 /* Quaternion */
-Quaternion Quat(double r, double i, double j, double k);
-Quaternion Quatvec(double r, Point3 v);
-Quaternion addq(Quaternion a, Quaternion b);
-Quaternion subq(Quaternion a, Quaternion b);
-Quaternion mulq(Quaternion q, Quaternion r);
-Quaternion smulq(Quaternion q, double s);
-Quaternion sdivq(Quaternion q, double s);
-double dotq(Quaternion q, Quaternion r);
-Quaternion invq(Quaternion q);
-double qlen(Quaternion q);
-Quaternion normq(Quaternion q);
-Quaternion slerp(Quaternion q, Quaternion r, double t);
-Point3 qrotate(Point3 p, Point3 axis, double θ);
+Quaternion Quat(double r, double i, double j, double k)
+Quaternion Quatvec(double r, Point3 v)
+Quaternion addq(Quaternion a, Quaternion b)
+Quaternion subq(Quaternion a, Quaternion b)
+Quaternion mulq(Quaternion q, Quaternion r)
+Quaternion smulq(Quaternion q, double s)
+Quaternion sdivq(Quaternion q, double s)
+double dotq(Quaternion q, Quaternion r)
+Quaternion invq(Quaternion q)
+double qlen(Quaternion q)
+Quaternion normq(Quaternion q)
+Quaternion slerp(Quaternion q, Quaternion r, double t)
+Quaternion qsandwich(Quaternion q, Quaternion r)
+Point3 qsandwichpt3(Quaternion q, Point3 p)
+Point3 qrotate(Point3 p, Point3 axis, double θ)
 .PB
 /* RFrame */
-void rframematrix(Matrix, RFrame);
-void rframematrix3(Matrix3, RFrame3);
-Point2 rframexform(Point2 p, RFrame rf);
-Point3 rframexform3(Point3 p, RFrame3 rf);
-Point2 invrframexform(Point2 p, RFrame rf);
-Point3 invrframexform3(Point3 p, RFrame3 rf);
+void rframematrix(Matrix, RFrame)
+void rframematrix3(Matrix3, RFrame3)
+Point2 rframexform(Point2 p, RFrame rf)
+Point3 rframexform3(Point3 p, RFrame3 rf)
+Point2 invrframexform(Point2 p, RFrame rf)
+Point3 invrframexform3(Point3 p, RFrame3 rf)
 .PB
 /* Triangle2 */
-Point2 centroid(Triangle2 t);
-Point3 barycoords(Triangle2 t, Point2 p);
+Point2 centroid(Triangle2 t)
+Point3 barycoords(Triangle2 t, Point2 p)
 .PB
 /* Triangle3 */
-Point3 centroid3(Triangle3 t);
+Point3 centroid3(Triangle3 t)
 .PB
 /* Fmt */
 #pragma varargck type "v" Point2
 #pragma varargck type "V" Point3
-int vfmt(Fmt*);
-int Vfmt(Fmt*);
-void GEOMfmtinstall(void);
+int vfmt(Fmt*)
+int Vfmt(Fmt*)
+void GEOMfmtinstall(void)
 .SH DESCRIPTION
 This library provides routines to operate with homogeneous coordinates
 in 2D and 3D projective spaces by means of points, matrices,
@@ -628,6 +630,21 @@
 by a factor of
 .IR t ,
 and returns the result.
+.TP
+.B qsandwich(\fIq\fP,\fIr\fP)
+Returns the result of "sandwiching"
+.I q
+with
+.IR r ,
+i.e. r' = qrq⁻¹.
+.TP
+.B qsandwichpt3(\fIq\fP,\fIp\fP)
+Same as
+.IR qsandwich (2),
+but treating
+.I p
+as a pure quaternion when computing the product, and keeping its
+projective coordinate (w) intact.
 .TP
 .B qrotate(\fIp\fP,\fIaxis\fP,\fIθ\fP)
 Returns the result of rotating the point
--- a/sys/src/libgeometry/quaternion.c
+++ b/sys/src/libgeometry/quaternion.c
@@ -96,13 +96,27 @@
 	return addq(smulq(q, cos(θ)), smulq(v, sin(θ))); /* q cos(θ) + v sin(θ) */
 }
 
+Quaternion
+qsandwich(Quaternion q, Quaternion p)
+{
+	return mulq(mulq(q, p), invq(q)); /* qpq⁻¹ */
+}
+
 Point3
+qsandwichpt3(Quaternion q, Point3 p)
+{
+	Quaternion r;
+
+	r = qsandwich(q, Quatvec(0, p));
+	return Pt3(r.i, r.j, r.k, p.w);
+}
+
+Point3
 qrotate(Point3 p, Point3 axis, double θ)
 {
-	Quaternion qaxis, qr;
+	Quaternion qaxis;
 
 	θ /= 2;
 	qaxis = Quatvec(cos(θ), mulpt3(axis, sin(θ)));
-	qr = mulq(mulq(qaxis, Quatvec(0, p)), invq(qaxis)); /* qpq⁻¹ */
-	return Pt3(qr.i, qr.j, qr.k, p.w);
+	return qsandwichpt3(qaxis, p);
 }
--