libgeometry revamp

4 files changed, 804 insertions, 769 deletions

diff --git a/sys/man/2/arith3 b/sys/man/2/arith3
deleted file mode 100644
index d2c4acb08..000000000
--- a/sys/man/2/arith3
+++ /dev/null
@@ -1,268 +0,0 @@
-.TH ARITH3 2
-.SH NAME
-add3, sub3, neg3, div3, mul3, eqpt3, closept3, dot3, cross3, len3, dist3, unit3, midpt3, lerp3, reflect3, nearseg3, pldist3, vdiv3, vrem3, pn2f3, ppp2f3, fff2p3, pdiv4, add4, sub4 \- operations on 3-d points and planes
-.SH SYNOPSIS
-.B
-#include <draw.h>
-.br
-.B
-#include <geometry.h>
-.PP
-.B
-Point3 add3(Point3 a, Point3 b)
-.PP
-.B
-Point3 sub3(Point3 a, Point3 b)
-.PP
-.B
-Point3 neg3(Point3 a)
-.PP
-.B
-Point3 div3(Point3 a, double b)
-.PP
-.B
-Point3 mul3(Point3 a, double b)
-.PP
-.B
-int eqpt3(Point3 p, Point3 q)
-.PP
-.B
-int closept3(Point3 p, Point3 q, double eps)
-.PP
-.B
-double dot3(Point3 p, Point3 q)
-.PP
-.B
-Point3 cross3(Point3 p, Point3 q)
-.PP
-.B
-double len3(Point3 p)
-.PP
-.B
-double dist3(Point3 p, Point3 q)
-.PP
-.B
-Point3 unit3(Point3 p)
-.PP
-.B
-Point3 midpt3(Point3 p, Point3 q)
-.PP
-.B
-Point3 lerp3(Point3 p, Point3 q, double alpha)
-.PP
-.B
-Point3 reflect3(Point3 p, Point3 p0, Point3 p1)
-.PP
-.B
-Point3 nearseg3(Point3 p0, Point3 p1, Point3 testp)
-.PP
-.B
-double pldist3(Point3 p, Point3 p0, Point3 p1)
-.PP
-.B
-double vdiv3(Point3 a, Point3 b)
-.PP
-.B
-Point3 vrem3(Point3 a, Point3 b)
-.PP
-.B
-Point3 pn2f3(Point3 p, Point3 n)
-.PP
-.B
-Point3 ppp2f3(Point3 p0, Point3 p1, Point3 p2)
-.PP
-.B
-Point3 fff2p3(Point3 f0, Point3 f1, Point3 f2)
-.PP
-.B
-Point3 pdiv4(Point3 a)
-.PP
-.B
-Point3 add4(Point3 a, Point3 b)
-.PP
-.B
-Point3 sub4(Point3 a, Point3 b)
-.SH DESCRIPTION
-These routines do arithmetic on points and planes in affine or projective 3-space.
-Type
-.B Point3
-is
-.IP
-.EX
-.ta 6n
-typedef struct Point3 Point3;
-struct Point3{
-	double x, y, z, w;
-};
-.EE
-.PP
-Routines whose names end in
-.B 3
-operate on vectors or ordinary points in affine 3-space, represented by their Euclidean
-.B (x,y,z)
-coordinates.
-(They assume
-.B w=1
-in their arguments, and set
-.B w=1
-in their results.)
-.TF reflect3
-.TP
-Name
-Description
-.TP
-.B add3
-Add the coordinates of two points.
-.TP
-.B sub3
-Subtract coordinates of two points.
-.TP
-.B neg3
-Negate the coordinates of a point.
-.TP
-.B mul3
-Multiply coordinates by a scalar.
-.TP
-.B div3
-Divide coordinates by a scalar.
-.TP
-.B eqpt3
-Test two points for exact equality.
-.TP
-.B closept3
-Is the distance between two points smaller than 
-.IR eps ?
-.TP
-.B dot3
-Dot product.
-.TP
-.B cross3
-Cross product.
-.TP
-.B len3
-Distance to the origin.
-.TP
-.B dist3
-Distance between two points.
-.TP
-.B unit3
-A unit vector parallel to
-.IR p .
-.TP
-.B midpt3
-The midpoint of line segment 
-.IR pq .
-.TP
-.B lerp3
-Linear interpolation between 
-.I p
-and
-.IR q .
-.TP
-.B reflect3
-The reflection of point
-.I p
-in the segment joining 
-.I p0
-and
-.IR p1 .
-.TP
-.B nearseg3
-The closest point to 
-.I testp
-on segment
-.IR "p0 p1" .
-.TP
-.B pldist3
-The distance from 
-.I p
-to segment
-.IR "p0 p1" .
-.TP
-.B vdiv3
-Vector divide \(em the length of the component of 
-.I a
-parallel to
-.IR b ,
-in units of the length of
-.IR b .
-.TP
-.B vrem3
-Vector remainder \(em the component of 
-.I a
-perpendicular to
-.IR b .
-Ignoring roundoff, we have 
-.BR "eqpt3(add3(mul3(b, vdiv3(a, b)), vrem3(a, b)), a)" .
-.PD
-.PP
-The following routines convert amongst various representations of points
-and planes.  Planes are represented identically to points, by duality;
-a point
-.B p
-is on a plane
-.B q
-whenever
-.BR p.x*q.x+p.y*q.y+p.z*q.z+p.w*q.w=0 .
-Although when dealing with affine points we assume
-.BR p.w=1 ,
-we can't make the same assumption for planes.
-The names of these routines are extra-cryptic.  They contain an
-.B f
-(for `face') to indicate a plane,
-.B p
-for a point and
-.B n
-for a normal vector.
-The number
-.B 2
-abbreviates the word `to.'
-The number
-.B 3
-reminds us, as before, that we're dealing with affine points.
-Thus
-.B pn2f3
-takes a point and a normal vector and returns the corresponding plane.
-.TF reflect3
-.TP
-Name
-Description
-.TP
-.B pn2f3
-Compute the plane passing through
-.I p
-with normal
-.IR n .
-.TP
-.B ppp2f3
-Compute the plane passing through three points.
-.TP
-.B fff2p3
-Compute the intersection point of three planes.
-.PD
-.PP
-The names of the following routines end in
-.B 4
-because they operate on points in projective 4-space,
-represented by their homogeneous coordinates.
-.TP
-pdiv4
-Perspective division.  Divide
-.B p.w
-into
-.IR p 's
-coordinates, converting to affine coordinates.
-If
-.B p.w
-is zero, the result is the same as the argument.
-.TP
-add4
-Add the coordinates of two points.
-.PD
-.TP
-sub4
-Subtract the coordinates of two points.
-.SH SOURCE
-.B /sys/src/libgeometry
-.SH "SEE ALSO
-.IR matrix (2)
diff --git a/sys/man/2/geometry b/sys/man/2/geometry
new file mode 100644
index 000000000..a8a44247e
--- /dev/null
+++ b/sys/man/2/geometry
@@ -0,0 +1,804 @@
+.TH GEOMETRY 2
+.SH NAME
+Flerp, fclamp, Pt2, Vec2, addpt2, subpt2, mulpt2, divpt2, lerp2, dotvec2, vec2len, normvec2, edgeptcmp, ptinpoly, Pt3, Vec3, addpt3, subpt3, mulpt3, divpt3, lerp3, dotvec3, crossvec3, vec3len, normvec3, identity, addm, subm, mulm, smulm, transposem, detm, tracem, adjm, invm, xform, identity3, addm3, subm3, mulm3, smulm3, transposem3, detm3, tracem3, adjm3, invm3, xform3, Quat, Quatvec, addq, subq, mulq, smulq, sdivq, dotq, invq, qlen, normq, slerp, qrotate, rframexform, rframexform3, invrframexform, invrframexform3, centroid, barycoords, centroid3, vfmt, Vfmt, GEOMfmtinstall \- computational geometry library
+.SH SYNOPSIS
+.de PB
+.PP
+.ft L
+.nf
+..
+.PB
+#include <u.h>
+#include <libc.h>
+#include <geometry.h>
+.PB
+#define DEG 0.01745329251994330	/* π/180 */
+.PB
+typedef struct Point2 Point2;
+typedef struct Point3 Point3;
+typedef double Matrix[3][3];
+typedef double Matrix3[4][4];
+typedef struct Quaternion Quaternion;
+typedef struct RFrame RFrame;
+typedef struct RFrame3 RFrame3;
+typedef struct Triangle2 Triangle2;
+typedef struct Triangle3 Triangle3;
+.PB
+struct Point2 {
+	double x, y, w;
+};
+.PB
+struct Point3 {
+	double x, y, z, w;
+};
+.PB
+struct Quaternion {
+	double r, i, j, k;
+};
+.PB
+struct RFrame {
+	Point2 p;
+	Point2 bx, by;
+};
+.PB
+struct RFrame3 {
+	Point3 p;
+	Point3 bx, by, bz;
+};
+.PB
+struct Triangle2
+{
+	Point2 p0, p1, p2;
+};
+.PB
+struct Triangle3 {
+	Point3 p0, p1, p2;
+};
+.PB
+/* utils */
+double flerp(double a, double b, double t);
+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);
+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);
+double dotvec3(Point3 a, Point3 b);
+Point3 crossvec3(Point3 a, Point3 b);
+double vec3len(Point3 v);
+Point3 normvec3(Point3 v);
+.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);
+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);
+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 θ);
+.PB
+/* RFrame */
+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);
+.PB
+/* Triangle3 */
+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);
+.SH DESCRIPTION
+This library provides routines to operate with homogeneous coordinates
+in 2D and 3D projective spaces by means of points, matrices,
+quaternions and frames of reference.
+.PP
+Besides their many mathematical properties and applications, the data
+structures and algorithms used here to represent these abstractions
+are specifically tailored to the world of computer graphics and
+simulators, and so it uses the conventions associated with these
+fields, such as the right-hand rule for coordinate systems and column
+vectors for matrix operations.
+.SS UTILS
+These utility functions provide extra floating-point operations that
+are not available in the standard libc.
+.TP
+Name
+Description
+.TP
+.B flerp
+Performs a linear interpolation by a factor of
+.I t
+between
+.I a
+and
+.IR b ,
+and returns the result.
+.TP
+.B fclamp
+Constrains
+.I n
+to a value between
+.I min
+and
+.IR max ,
+and returns the result.
+.SS Points
+A point
+.B (x,y,w)
+in projective space results in the point
+.B (x/w,y/w)
+in Euclidean space. Vectors are represented by setting
+.B w
+to zero, since they don't belong to any projective plane themselves.
+.TP
+Name
+Description
+.TP
+.B Pt2
+Constructor function for a Point2 point.
+.TP
+.B Vec2
+Constructor function for a Point2 vector.
+.TP
+.B addpt2
+Creates a new 2D point out of the sum of
+.IR a 's
+and
+.IR b 's
+components.
+.TP
+subpt2
+Creates a new 2D point out of the substraction of
+.IR a 's
+by
+.IR b 's
+components.
+.TP
+mulpt2
+Creates a new 2D point from multiplying
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+divpt2
+Creates a new 2D point from dividing
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+lerp2
+Performs a linear interpolation between the 2D points
+.I a
+and
+.I b
+by a factor of
+.IR t ,
+and returns the result.
+.TP
+dotvec2
+Computes the dot product of vectors
+.I a
+and
+.IR b .
+.TP
+vec2len
+Computes the length—magnitude—of vector
+.IR v .
+.TP
+normvec2
+Normalizes the vector
+.I v
+and returns a new 2D point.
+.TP
+edgeptcmp
+Performs a comparison between an edge, defined by a directed line from
+.I e0
+to
+.IR e1 ,
+and the point
+.IR p .
+If the point is to the right of the line, the result is >0; if it's to
+the left, the result is <0; otherwise—when the point is on the line—,
+it returns 0.
+.TP
+ptinpoly
+Returns 1 if the 2D point
+.I p
+lies within the
+.IR npts -vertex
+polygon defined by
+.IR pts ,
+0 otherwise.
+.TP
+Pt3
+Constructor function for a Point3 point.
+.TP
+Vec3
+Constructor function for a Point3 vector.
+.TP
+addpt3
+Creates a new 3D point out of the sum of
+.IR a 's
+and
+.IR b 's
+components.
+.TP
+subpt3
+Creates a new 3D point out of the substraction of
+.IR a 's
+by
+.IR b 's
+components.
+.TP
+mulpt3
+Creates a new 3D point from multiplying
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+divpt3
+Creates a new 3D point from dividing
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+lerp3
+Performs a linear interpolation between the 3D points
+.I a
+and
+.I b
+by a factor of
+.IR t ,
+and returns the result.
+.TP
+dotvec3
+Computes the dot/inner product of vectors
+.I a
+and
+.IR b .
+.TP
+crossvec3
+Computes the cross/outer product of vectors
+.I a
+and
+.IR b .
+.TP
+vec3len
+Computes the length—magnitude—of vector
+.IR v .
+.TP
+normvec3
+Normalizes the vector
+.I v
+and returns a new 3D point.
+.SS Matrices
+.TP
+Name
+Description
+.TP
+identity
+Initializes
+.I m
+into an identity, 3x3 matrix.
+.TP
+addm
+Sums the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+subm
+Substracts the matrix
+.I a
+by
+.I b
+and stores the result back in
+.IR a .
+.TP
+mulm
+Multiplies the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+smulm
+Multiplies every element of
+.I m
+by the scalar
+.IR s ,
+storing the result in m.
+.TP
+transposem
+Transforms the matrix
+.I m
+into its transpose.
+.TP
+detm
+Computes the determinant of
+.I m
+and returns the result.
+.TP
+tracem
+Computes the trace of
+.I m
+and returns the result.
+.TP
+adjm
+Transforms the matrix
+.I m
+into its adjoint.
+.TP
+invm
+Transforms the matrix
+.I m
+into its inverse.
+.TP
+xform
+Transforms the point
+.I p
+by the matrix
+.I m
+and returns the new 2D point.
+.TP
+identity3
+Initializes
+.I m
+into an identity, 4x4 matrix.
+.TP
+addm3
+Sums the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+subm3
+Substracts the matrix
+.I a
+by
+.I b
+and stores the result back in
+.IR a .
+.TP
+mulm3
+Multiplies the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+smulm3
+Multiplies every element of
+.I m
+by the scalar
+.IR s ,
+storing the result in m.
+.TP
+transposem3
+Transforms the matrix
+.I m
+into its transpose.
+.TP
+detm3
+Computes the determinant of
+.I m
+and returns the result.
+.TP
+tracem3
+Computes the trace of
+.I m
+and returns the result.
+.TP
+adjm3
+Transforms the matrix
+.I m
+into its adjoint.
+.TP
+invm3
+Transforms the matrix
+.I m
+into its inverse.
+.TP
+xform3
+Transforms the point
+.I p
+by the matrix
+.I m
+and returns the new 3D point.
+.SS Quaternions
+Quaternions are an extension of the complex numbers conceived as a
+tool to analyze 3-dimensional points.  They are most commonly used to
+orient and rotate objects in 3D space.
+.TP
+Name
+Description
+.TP
+Quat
+Constructor function for a Quaternion.
+.TP
+Quatvec
+Constructor function for a Quaternion that takes the imaginary part in
+the form of a vector
+.IR v .
+.TP
+addq
+Creates a new quaternion out of the sum of
+.IR a 's
+and
+.IR b 's
+components.
+.TP
+subq
+Creates a new quaternion from the substraction of
+.IR a 's
+by
+.IR b 's
+components.
+.TP
+mulq
+Multiplies
+.I a
+and
+.I b
+and returns a new quaternion.
+.TP
+smulq
+Multiplies each of the components of
+.I q
+by the scalar
+.IR s ,
+returning a new quaternion.
+.TP
+sdivq
+Divides each of the components of
+.I q
+by the scalar
+.IR s ,
+returning a new quaternion.
+.TP
+dotq
+Computes the dot-product of
+.I q
+and
+.IR r ,
+and returns the result.
+.TP
+invq
+Computes the inverse of
+.I q
+and returns a new quaternion out of it.
+.TP
+qlen
+Computes
+.IR q 's
+length—magnitude—and returns the result.
+.TP
+normq
+Normalizes
+.I q
+and returns a new quaternion out of it.
+.TP
+slerp
+Performs a spherical linear interpolation between the quaternions
+.I q
+and
+.I r
+by a factor of
+.IR t ,
+and returns the result.
+.TP
+qrotate
+Returns the result of rotating the point
+.I p
+around the vector
+.I axis
+by
+.I θ
+radians.
+.SS Frames of reference
+A frame of reference in a
+.IR n -dimensional
+space is made out of n+1 points, one being the origin
+.IR p ,
+relative to some other frame of reference, and the remaining being the
+basis vectors
+.I b1,⋯,bn
+that define the metric within that frame.
+.PP
+Every one of these routines assumes the origin reference frame
+.B O
+has an orthonormal basis when performing an inverse transformation;
+it's up to the user to apply a forward transformation to the resulting
+point with the proper reference frame if that's not the case.
+.TP
+Name
+Description
+.TP
+rframexform
+Transforms the point
+.IR p ,
+relative to origin O, into the frame of reference
+.I rf
+with origin in
+.BR rf.p ,
+which is itself also relative to O. It then returns the new 2D point.
+.TP
+rframexform3
+Transforms the point
+.IR p ,
+relative to origin O, into the frame of reference
+.I rf
+with origin in
+.BR rf.p ,
+which is itself also relative to O. It then returns the new 3D point.
+.TP
+invrframexform
+Transforms the point
+.IR p ,
+relative to
+.BR rf.p ,
+into the frame of reference O, assumed to have an orthonormal basis.
+.TP
+invrframexform3
+Transforms the point
+.IR p ,
+relative to
+.BR rf.p ,
+into the frame of reference O, assumed to have an orthonormal basis.
+.SS Triangles
+.TP
+Name
+Description
+.TP
+centroid
+Returns the geometric center of
+.B Triangle2
+.IR t .
+.TP
+barycoords
+Returns a 3D point that represents the barycentric coordinates of the
+2D point
+.I p
+relative to the triangle
+.IR t .
+.TP
+centroid3
+Returns the geometric center of
+.B Triangle3
+.IR t .
+.SH EXAMPLE
+The following is a common example of an
+.B RFrame
+being used to define the coordinate system of a
+.IR rio (3)
+window.  It places the origin at the center of the window and sets up
+an orthonormal basis with the
+.IR y -axis
+pointing upwards, to contrast with the window system where
+.IR y -values
+grow downwards (see
+.IR graphics (2)).
+.PP
+.EX
+#include <u.h>
+#include <libc.h>
+#include <draw.h>
+#include <geometry.h>
+
+RFrame screenrf;
+
+Point
+toscreen(Point2 p)
+{
+	p = invrframexform(p, screenrf);
+	return Pt(p.x,p.y);
+}
+
+Point2
+fromscreen(Point p)
+{
+	return rframexform(Pt2(p.x,p.y,1), screenrf);
+}
+
+void
+main(void)
+	⋯
+	screenrf.p = Pt2(screen->r.min.x+Dx(screen->r)/2,screen->r.max.y-Dy(screen->r)/2,1);
+	screenrf.bx = Vec2(1, 0);
+	screenrf.by = Vec2(0,-1);
+	⋯
+.EE
+.PP
+The following snippet shows how to use the
+.B RFrame
+declared earlier to locate and draw a ship based on its orientation,
+for which we use matrix translation
+.B T
+and rotation
+.BR R
+transformations.
+.PP
+.EX
+⋯
+typedef struct Ship Ship;
+typedef struct Shipmdl Shipmdl;
+
+struct Ship
+{
+	RFrame;
+	double θ; /* orientation (yaw) */
+	Shipmdl mdl;
+};
+
+struct Shipmdl
+{
+	Point2 pts[3]; /* a free-form triangle */
+};
+
+Ship *ship;
+
+void
+redraw(void)
+{
+	int i;
+	Point pts[3+1];
+	Point2 *p;
+	Matrix T = {
+		1, 0, ship->p.x,
+		0, 1, ship->p.y,
+		0, 0, 1,
+	}, R = {
+		cos(ship->θ), -sin(ship->θ), 0,
+		sin(ship->θ),  cos(ship->θ), 0,
+		0, 0, 1,
+	};
+
+	mulm(T, R); /* rotate, then translate */
+	p = ship->mdl.pts;
+	for(i = 0; i < nelem(pts)-1; i++)
+		pts[i] = toscreen(xform(p[i], T));
+	pts[i] = pts[0];
+	draw(screen, screen->r, display->white, nil, ZP);
+	poly(screen, pts, nelem(pts), 0, 0, 0, display->black, ZP);
+}
+⋯
+main(void)
+	⋯
+	ship = malloc(sizeof(Ship));
+	ship->p = Pt2(0,0,1); /* place it at the origin */
+	ship->θ = 45*DEG; /* counter-clockwise */
+	ship->mdl.pts[0] = Pt2( 10, 0,1);
+	ship->mdl.pts[1] = Pt2(-10, 5,1);
+	ship->mdl.pts[2] = Pt2(-10,-5,1);
+	⋯
+	redraw();
+⋯
+.EE
+.PP
+Notice how we could've used the
+.B RFrame
+embedded in the
+.B ship
+to transform the
+.B Shipmdl
+into the window.  Instead of applying the matrices to every point, the
+ship's local frame of reference can be rotated, effectively changing
+the model coordinates after an
+.IR invrframexform .
+We are also getting rid of the
+.B θ
+variable, since it's no longer needed.
+.PP
+.EX
+⋯
+struct Ship
+{
+	RFrame;
+	Shipmdl mdl;
+};
+⋯
+redraw(void)
+	⋯
+		pts[i] = toscreen(invrframexform(p[i], *ship));
+⋯
+main(void)
+	⋯
+	Matrix R = {
+		cos(45*DEG), -sin(45*DEG), 0,
+		sin(45*DEG),  cos(45*DEG), 0,
+		0, 0, 1,
+	};
+	⋯
+	//ship->θ = 45*DEG; /* counter-clockwise */
+	ship->bx = xform(ship->bx, R);
+	ship->by = xform(ship->by, R);
+⋯
+.EE
+.SH SOURCE
+.B /sys/src/libgeometry
+.SH SEE ALSO
+.IR sin (2),
+.IR floor (2),
+.IR graphics (2)
+.br
+Philip J. Schneider, David H. Eberly,
+“Geometric Tools for Computer Graphics”,
+.I
+Morgan Kaufmann Publishers, 2003.
+.br
+Jonathan Blow,
+“Understanding Slerp, Then Not Using it”,
+.I
+The Inner Product, April 2004.
+.br
+https://www.3dgep.com/understanding-quaternions/
+.SH BUGS
+No care is taken to avoid numeric overflows.
+.SH HISTORY
+Libgeometry first appeared in Plan 9 from Bell Labs.  It was revamped
+for 9front in January of 2023.
diff --git a/sys/man/2/matrix b/sys/man/2/matrix
deleted file mode 100644
index 6291aa6de..000000000
--- a/sys/man/2/matrix
+++ /dev/null
@@ -1,350 +0,0 @@
-.TH MATRIX 2
-.SH NAME
-ident, matmul, matmulr, determinant, adjoint, invertmat, xformpoint, xformpointd, xformplane, pushmat, popmat, rot, qrot, scale, move, xform, ixform, persp, look, viewport \- Geometric transformations
-.SH SYNOPSIS
-.PP
-.B
-#include <draw.h>
-.PP
-.B
-#include <geometry.h>
-.PP
-.B
-void ident(Matrix m)
-.PP
-.B
-void matmul(Matrix a, Matrix b)
-.PP
-.B
-void matmulr(Matrix a, Matrix b)
-.PP
-.B
-double determinant(Matrix m)
-.PP
-.B
-void adjoint(Matrix m, Matrix madj)
-.PP
-.B
-double invertmat(Matrix m, Matrix inv)
-.PP
-.B
-Point3 xformpoint(Point3 p, Space *to, Space *from)
-.PP
-.B
-Point3 xformpointd(Point3 p, Space *to, Space *from)
-.PP
-.B
-Point3 xformplane(Point3 p, Space *to, Space *from)
-.PP
-.B
-Space *pushmat(Space *t)
-.PP
-.B
-Space *popmat(Space *t)
-.PP
-.B
-void rot(Space *t, double theta, int axis)
-.PP
-.B
-void qrot(Space *t, Quaternion q)
-.PP
-.B
-void scale(Space *t, double x, double y, double z)
-.PP
-.B
-void move(Space *t, double x, double y, double z)
-.PP
-.B
-void xform(Space *t, Matrix m)
-.PP
-.B
-void ixform(Space *t, Matrix m, Matrix inv)
-.PP
-.B
-int persp(Space *t, double fov, double n, double f)
-.PP
-.B
-void look(Space *t, Point3 eye, Point3 look, Point3 up)
-.PP
-.B
-void viewport(Space *t, Rectangle r, double aspect)
-.SH DESCRIPTION
-These routines manipulate 3-space affine and projective transformations,
-represented as 4\(mu4 matrices, thus:
-.IP
-.EX
-.ta 6n
-typedef double Matrix[4][4];
-.EE
-.PP
-.I Ident
-stores an identity matrix in its argument.
-.I Matmul
-stores
-.I a\(mub
-in
-.IR a .
-.I Matmulr
-stores
-.I b\(mua
-in
-.IR b .
-.I Determinant
-returns the determinant of matrix
-.IR m .
-.I Adjoint
-stores the adjoint (matrix of cofactors) of
-.I m
-in
-.IR madj .
-.I Invertmat
-stores the inverse of matrix
-.I m
-in
-.IR minv ,
-returning
-.IR m 's
-determinant.
-Should
-.I m
-be singular (determinant zero),
-.I invertmat
-stores its
-adjoint in
-.IR minv .
-.PP
-The rest of the routines described here
-manipulate
-.I Spaces
-and transform
-.IR Point3s .
-A
-.I Point3
-is a point in three-space, represented by its
-homogeneous coordinates:
-.IP
-.EX
-typedef struct Point3 Point3;
-struct Point3{
-	double x, y, z, w;
-};
-.EE
-.PP
-The homogeneous coordinates
-.RI ( x ,
-.IR y ,
-.IR z ,
-.IR w )
-represent the Euclidean point
-.RI ( x / w ,
-.IR y / w ,
-.IR z / w )
-if
-.IR w ≠0,
-and a ``point at infinity'' if
-.IR w =0.
-.PP
-A
-.I Space
-is just a data structure describing a coordinate system:
-.IP
-.EX
-typedef struct Space Space;
-struct Space{
-	Matrix t;
-	Matrix tinv;
-	Space *next;
-};
-.EE
-.PP
-It contains a pair of transformation matrices and a pointer
-to the
-.IR Space 's
-parent.  The matrices transform points to and from the ``root
-coordinate system,'' which is represented by a null
-.I Space
-pointer.
-.PP
-.I Pushmat
-creates a new
-.IR Space .
-Its argument is a pointer to the parent space.  Its result
-is a newly allocated copy of the parent, but with its
-.B next
-pointer pointing at the parent.
-.I Popmat
-discards the
-.B Space
-that is its argument, returning a pointer to the stack.
-Nominally, these two functions define a stack of transformations,
-but
-.B pushmat
-can be called multiple times
-on the same
-.B Space
-multiple times, creating a transformation tree.
-.PP
-.I Xformpoint
-and
-.I Xformpointd
-both transform points from the
-.B Space
-pointed to by
-.I from
-to the space pointed to by
-.IR to .
-Either pointer may be null, indicating the root coordinate system.
-The difference between the two functions is that
-.B xformpointd
-divides
-.IR x ,
-.IR y ,
-.IR z ,
-and
-.I w
-by
-.IR w ,
-if
-.IR w ≠0,
-making
-.RI ( x ,
-.IR y ,
-.IR z )
-the Euclidean coordinates of the point.
-.PP
-.I Xformplane
-transforms planes or normal vectors.  A plane is specified by the
-coefficients
-.RI ( a ,
-.IR b ,
-.IR c ,
-.IR d )
-of its implicit equation
-.IR ax+by+cz+d =0.
-Since this representation is dual to the homogeneous representation of points,
-.B libgeometry
-represents planes by
-.B Point3
-structures, with
-.RI ( a ,
-.IR b ,
-.IR c ,
-.IR d )
-stored in
-.RI ( x ,
-.IR y ,
-.IR z ,
-.IR w ).
-.PP
-The remaining functions transform the coordinate system represented
-by a
-.BR Space .
-Their
-.B Space *
-argument must be non-null \(em you can't modify the root
-.BR Space .
-.I Rot
-rotates by angle
-.I theta
-(in radians) about the given
-.IR axis ,
-which must be one of
-.BR XAXIS ,
-.B YAXIS
-or
-.BR ZAXIS .
-.I Qrot
-transforms by a rotation about an arbitrary axis, specified by
-.B Quaternion
-.IR q .
-.PP
-.I Scale
-scales the coordinate system by the given scale factors in the directions of the three axes.
-.IB Move
-translates by the given displacement in the three axial directions.
-.PP
-.I Xform
-transforms the coordinate system by the given
-.BR Matrix .
-If the matrix's inverse is known
-.I a
-.IR priori ,
-calling
-.I ixform
-will save the work of recomputing it.
-.PP
-.I Persp
-does a perspective transformation.
-The transformation maps the frustum with apex at the origin,
-central axis down the positive
-.I y
-axis, and apex angle
-.I fov
-and clipping planes
-.IR y = n
-and
-.IR y = f
-into the double-unit cube.
-The plane
-.IR y = n
-maps to
-.IR y '=-1,
-.IR y = f
-maps to
-.IR y '=1.
-.PP
-.I Look
-does a view-pointing transformation.  The
-.B eye
-point is moved to the origin.
-The line through the
-.I eye
-and
-.I look
-points is aligned with the y axis,
-and the plane containing the
-.BR eye ,
-.B look
-and
-.B up
-points is rotated into the
-.IR x - y
-plane.
-.PP
-.I Viewport
-maps the unit-cube window into the given screen viewport.
-The viewport rectangle
-.I r
-has
-.IB r .min
-at the top left-hand corner, and
-.IB r .max
-just outside the lower right-hand corner.
-Argument
-.I aspect
-is the aspect ratio
-.RI ( dx / dy )
-of the viewport's pixels (not of the whole viewport).
-The whole window is transformed to fit centered inside the viewport with equal
-slop on either top and bottom or left and right, depending on the viewport's
-aspect ratio.
-The window is viewed down the
-.I y
-axis, with
-.I x
-to the left and
-.I z
-up.  The viewport
-has
-.I x
-increasing to the right and
-.I y
-increasing down.  The window's
-.I y
-coordinates are mapped, unchanged, into the viewport's
-.I z
-coordinates.
-.SH SOURCE
-.B /sys/src/libgeometry/matrix.c
-.SH "SEE ALSO
-.IR arith3 (2)
diff --git a/sys/man/2/quaternion b/sys/man/2/quaternion
deleted file mode 100644
index 3b637b184..000000000
--- a/sys/man/2/quaternion
+++ /dev/null
@@ -1,151 +0,0 @@
-.TH QUATERNION 2
-.SH NAME
-qtom, mtoq, qadd, qsub, qneg, qmul, qdiv, qunit, qinv, qlen, slerp, qmid, qsqrt \- Quaternion arithmetic
-.SH SYNOPSIS
-.B
-#include <draw.h>
-.br
-.B
-#include <geometry.h>
-.PP
-.B
-Quaternion qadd(Quaternion q, Quaternion r)
-.PP
-.B
-Quaternion qsub(Quaternion q, Quaternion r)
-.PP
-.B
-Quaternion qneg(Quaternion q)
-.PP
-.B
-Quaternion qmul(Quaternion q, Quaternion r)
-.PP
-.B
-Quaternion qdiv(Quaternion q, Quaternion r)
-.PP
-.B
-Quaternion qinv(Quaternion q)
-.PP
-.B
-double qlen(Quaternion p)
-.PP
-.B
-Quaternion qunit(Quaternion q)
-.PP
-.B
-void qtom(Matrix m, Quaternion q)
-.PP
-.B
-Quaternion mtoq(Matrix mat)
-.PP
-.B
-Quaternion slerp(Quaternion q, Quaternion r, double a)
-.PP
-.B
-Quaternion qmid(Quaternion q, Quaternion r)
-.PP
-.B
-Quaternion qsqrt(Quaternion q)
-.SH DESCRIPTION
-The Quaternions are a non-commutative extension field of the Real numbers, designed
-to do for rotations in 3-space what the complex numbers do for rotations in 2-space.
-Quaternions have a real component
-.I r
-and an imaginary vector component \fIv\fP=(\fIi\fP,\fIj\fP,\fIk\fP).
-Quaternions add componentwise and multiply according to the rule
-(\fIr\fP,\fIv\fP)(\fIs\fP,\fIw\fP)=(\fIrs\fP-\fIv\fP\v'-.3m'.\v'.3m'\fIw\fP, \fIrw\fP+\fIvs\fP+\fIv\fP×\fIw\fP),
-where \v'-.3m'.\v'.3m' and × are the ordinary vector dot and cross products.
-The multiplicative inverse of a non-zero quaternion (\fIr\fP,\fIv\fP)
-is (\fIr\fP,\fI-v\fP)/(\fIr\^\fP\u\s-22\s+2\d-\fIv\fP\v'-.3m'.\v'.3m'\fIv\fP).
-.PP
-The following routines do arithmetic on quaternions, represented as
-.IP
-.EX
-.ta 6n
-typedef struct Quaternion Quaternion;
-struct Quaternion{
-	double r, i, j, k;
-};
-.EE
-.TF qunit
-.TP
-Name
-Description
-.TP
-.B qadd
-Add two quaternions.
-.TP
-.B qsub
-Subtract two quaternions.
-.TP
-.B qneg
-Negate a quaternion.
-.TP
-.B qmul
-Multiply two quaternions.
-.TP
-.B qdiv
-Divide two quaternions.
-.TP
-.B qinv
-Return the multiplicative inverse of a quaternion.
-.TP
-.B qlen
-Return
-.BR sqrt(q.r*q.r+q.i*q.i+q.j*q.j+q.k*q.k) ,
-the length of a quaternion.
-.TP
-.B qunit
-Return a unit quaternion 
-.RI ( length=1 )
-with components proportional to
-.IR q 's.
-.PD
-.PP
-A rotation by angle \fIθ\fP about axis
-.I A
-(where
-.I A
-is a unit vector) can be represented by
-the unit quaternion \fIq\fP=(cos \fIθ\fP/2, \fIA\fPsin \fIθ\fP/2).
-The same rotation is represented by \(mi\fIq\fP; a rotation by \(mi\fIθ\fP about \(mi\fIA\fP is the same as a rotation by \fIθ\fP about \fIA\fP.
-The quaternion \fIq\fP transforms points by
-(0,\fIx',y',z'\fP) = \%\fIq\fP\u\s-2-1\s+2\d(0,\fIx,y,z\fP)\fIq\fP.
-Quaternion multiplication composes rotations.
-The orientation of an object in 3-space can be represented by a quaternion
-giving its rotation relative to some `standard' orientation.
-.PP
-The following routines operate on rotations or orientations represented as unit quaternions:
-.TF slerp
-.TP
-.B mtoq
-Convert a rotation matrix (see
-.IR matrix (2))
-to a unit quaternion.
-.TP
-.B qtom
-Convert a unit quaternion to a rotation matrix.
-.TP
-.B slerp
-Spherical lerp.  Interpolate between two orientations.
-The rotation that carries
-.I q
-to
-.I r
-is \%\fIq\fP\u\s-2-1\s+2\d\fIr\fP, so
-.B slerp(q, r, t)
-is \fIq\fP(\fIq\fP\u\s-2-1\s+2\d\fIr\fP)\u\s-2\fIt\fP\s+2\d.
-.TP
-.B qmid
-.B slerp(q, r, .5)
-.TP
-.B qsqrt
-The square root of
-.IR q .
-This is just a rotation about the same axis by half the angle.
-.PD
-.SH SOURCE
-.B /sys/src/libgeometry/quaternion.c
-.SH SEE ALSO
-.IR matrix (2),
-.IR qball (2)