[9front] [patchset] libgeometry revamp

[9front] [patchset] libgeometry revamp

[rgl]

[-- Attachment #1: Type: text/plain, Size: 1317 bytes --]

hello everyone,

this is a proposal to renovate the current libgeometry with a more
modern implementation that i've been using in my personal projects for
the last four years.  it provides almost all of the functionality
available in the old one, with a cleaner and more consistent interface
(similar naming conventions as those used in libdraw's geometric
routines), plus the addition of 2D transformations (not just 3D), a
frame of reference data type for 2- and 3-spaces, and utilities such
as point-in-polygon testing and barycentric coordinate computations.

i documented everything in a single man page—geometry(2)—, so i got
rid of those related to the old library, with the exception of
qball(2), on which i'm still working and whose implementation i
believe to be broken; i'm using a version by BurnZeZ with some minor
changes that actually works, but it doesn't perform rotations the way
it should (demo: https://www.youtube.com/watch?v=g3pha1EXhIk) and i
will need more time to come up with a proper Arcball controller.

some examples of projects where i use the library:

 - http://git.antares-labs.eu/boids/
 - http://git.antares-labs.eu/etoys/
 - http://git.antares-labs.eu/gamephysics/
 - http://git.antares-labs.eu/3dee/


if it looks reasonable to you, i can push it upstream.



cheers!

-rodri

[-- Attachment #2: Type: text/plain, Size: 76088 bytes --]

diff 8188b4f4f0e07b6669e6ae3c6c1099af917eaab4 uncommitted
--- a/sys/include/geometry.h
+++ b/sys/include/geometry.h
@@ -1,89 +1,137 @@
 #pragma lib "libgeometry.a"
 #pragma src "/sys/src/libgeometry"
-typedef double Matrix[4][4];
+
+#define DEG 0.01745329251994330	/* π/180 */
+
+typedef struct Point2 Point2;
 typedef struct Point3 Point3;
+typedef double Matrix[3][3];
+typedef double Matrix3[4][4];
 typedef struct Quaternion Quaternion;
-typedef struct Space Space;
-struct Point3{
+typedef struct RFrame RFrame;
+typedef struct RFrame3 RFrame3;
+typedef struct Triangle2 Triangle2;
+typedef struct Triangle3 Triangle3;
+
+struct Point2 {
+	double x, y, w;
+};
+
+struct Point3 {
 	double x, y, z, w;
 };
-struct Quaternion{
+
+struct Quaternion {
 	double r, i, j, k;
 };
-struct Space{
-	Matrix t;
-	Matrix tinv;
-	Space *next;
+
+struct RFrame {
+	Point2 p;
+	Point2 bx, by;
 };
-/*
- * 3-d point arithmetic
- */
-Point3 add3(Point3 a, Point3 b);
-Point3 sub3(Point3 a, Point3 b);
-Point3 neg3(Point3 a);
-Point3 div3(Point3 a, double b);
-Point3 mul3(Point3 a, double b);
-int eqpt3(Point3 p, Point3 q);
-int closept3(Point3 p, Point3 q, double eps);
-double dot3(Point3 p, Point3 q);
-Point3 cross3(Point3 p, Point3 q);
-double len3(Point3 p);
-double dist3(Point3 p, Point3 q);
-Point3 unit3(Point3 p);
-Point3 midpt3(Point3 p, Point3 q);
-Point3 lerp3(Point3 p, Point3 q, double alpha);
-Point3 reflect3(Point3 p, Point3 p0, Point3 p1);
-Point3 nearseg3(Point3 p0, Point3 p1, Point3 testp);
-double pldist3(Point3 p, Point3 p0, Point3 p1);
-double vdiv3(Point3 a, Point3 b);
-Point3 vrem3(Point3 a, Point3 b);
-Point3 pn2f3(Point3 p, Point3 n);
-Point3 ppp2f3(Point3 p0, Point3 p1, Point3 p2);
-Point3 fff2p3(Point3 f0, Point3 f1, Point3 f2);
-Point3 pdiv4(Point3 a);
-Point3 add4(Point3 a, Point3 b);
-Point3 sub4(Point3 a, Point3 b);
-/*
- * Quaternion arithmetic
- */
-void qtom(Matrix, Quaternion);
-Quaternion mtoq(Matrix);
-Quaternion qadd(Quaternion, Quaternion);
-Quaternion qsub(Quaternion, Quaternion);
-Quaternion qneg(Quaternion);
-Quaternion qmul(Quaternion, Quaternion);
-Quaternion qdiv(Quaternion, Quaternion);
-Quaternion qunit(Quaternion);
-Quaternion qinv(Quaternion);
+
+struct RFrame3 {
+	Point3 p;
+	Point3 bx, by, bz;
+};
+
+struct Triangle2
+{
+	Point2 p0, p1, p2;
+};
+
+struct Triangle3 {
+	Point3 p0, p1, p2;
+};
+
+/* utils */
+double flerp(double, double, double);
+double fclamp(double, double, double);
+
+/* Point2 */
+Point2 Pt2(double, double, double);
+Point2 Vec2(double, double);
+Point2 addpt2(Point2, Point2);
+Point2 subpt2(Point2, Point2);
+Point2 mulpt2(Point2, double);
+Point2 divpt2(Point2, double);
+Point2 lerp2(Point2, Point2, double);
+double dotvec2(Point2, Point2);
+double vec2len(Point2);
+Point2 normvec2(Point2);
+int edgeptcmp(Point2, Point2, Point2);
+int ptinpoly(Point2, Point2*, ulong);
+
+/* Point3 */
+Point3 Pt3(double, double, double, double);
+Point3 Vec3(double, double, double);
+Point3 addpt3(Point3, Point3);
+Point3 subpt3(Point3, Point3);
+Point3 mulpt3(Point3, double);
+Point3 divpt3(Point3, double);
+Point3 lerp3(Point3, Point3, double);
+double dotvec3(Point3, Point3);
+Point3 crossvec3(Point3, Point3);
+double vec3len(Point3);
+Point3 normvec3(Point3);
+
+/* Matrix */
+void identity(Matrix);
+void addm(Matrix, Matrix);
+void subm(Matrix, Matrix);
+void mulm(Matrix, Matrix);
+void smulm(Matrix, double);
+void transposem(Matrix);
+double detm(Matrix);
+double tracem(Matrix);
+void adjm(Matrix);
+void invm(Matrix);
+Point2 xform(Point2, Matrix);
+
+/* Matrix3 */
+void identity3(Matrix3);
+void addm3(Matrix3, Matrix3);
+void subm3(Matrix3, Matrix3);
+void mulm3(Matrix3, Matrix3);
+void smulm3(Matrix3, double);
+void transposem3(Matrix3);
+double detm3(Matrix3);
+double tracem3(Matrix3);
+void adjm3(Matrix3);
+void invm3(Matrix3);
+Point3 xform3(Point3, Matrix3);
+
+/* Quaternion */
+Quaternion Quat(double, double, double, double);
+Quaternion Quatvec(double, Point3);
+Quaternion addq(Quaternion, Quaternion);
+Quaternion subq(Quaternion, Quaternion);
+Quaternion mulq(Quaternion, Quaternion);
+Quaternion smulq(Quaternion, double);
+Quaternion sdivq(Quaternion, double);
+double dotq(Quaternion, Quaternion);
+Quaternion invq(Quaternion);
 double qlen(Quaternion);
+Quaternion normq(Quaternion);
 Quaternion slerp(Quaternion, Quaternion, double);
-Quaternion qmid(Quaternion, Quaternion);
-Quaternion qsqrt(Quaternion);
-void qball(Rectangle, Mouse *, Quaternion *, void (*)(void), Quaternion *);
-/*
- * Matrix arithmetic
- */
-void ident(Matrix);
-void matmul(Matrix, Matrix);
-void matmulr(Matrix, Matrix);
-double determinant(Matrix);
-void adjoint(Matrix, Matrix);
-double invertmat(Matrix, Matrix);
-/*
- * Space stack routines
- */
-Space *pushmat(Space *);
-Space *popmat(Space *);
-void rot(Space *, double, int);
-void qrot(Space *, Quaternion);
-void scale(Space *, double, double, double);
-void move(Space *, double, double, double);
-void xform(Space *, Matrix);
-void ixform(Space *, Matrix, Matrix);
-void look(Space *, Point3, Point3, Point3);
-int persp(Space *, double, double, double);
-void viewport(Space *, Rectangle, double);
-Point3 xformpoint(Point3, Space *, Space *);
-Point3 xformpointd(Point3, Space *, Space *);
-Point3 xformplane(Point3, Space *, Space *);
-#define	radians(d)	((d)*.01745329251994329572)
+Point3 qrotate(Point3, Point3, double);
+
+/* RFrame */
+Point2 rframexform(Point2, RFrame);
+Point3 rframexform3(Point3, RFrame3);
+Point2 invrframexform(Point2, RFrame);
+Point3 invrframexform3(Point3, RFrame3);
+
+/* Triangle2 */
+Point2 centroid(Triangle2);
+Point3 barycoords(Triangle2, Point2);
+
+/* Triangle3 */
+Point3 centroid3(Triangle3);
+
+/* Fmt */
+#pragma varargck type "v" Point2
+#pragma varargck type "V" Point3
+int vfmt(Fmt*);
+int Vfmt(Fmt*);
+void GEOMfmtinstall(void);
--- a/sys/man/2/arith3
+++ /dev/null
@@ -1,268 +1,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)
--- /dev/null
+++ b/sys/man/2/geometry
@@ -1,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.
--- a/sys/man/2/matrix
+++ /dev/null
@@ -1,350 +1,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)
--- a/sys/man/2/quaternion
+++ /dev/null
@@ -1,151 +1,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)
--- a/sys/src/libgeometry/arith3.c
+++ /dev/null
@@ -1,215 +1,0 @@
-#include <u.h>
-#include <libc.h>
-#include <draw.h>
-#include <geometry.h>
-/*
- * Routines whose names end in 3 work on points in Affine 3-space.
- * They ignore w in all arguments and produce w=1 in all results.
- * Routines whose names end in 4 work on points in Projective 3-space.
- */
-Point3 add3(Point3 a, Point3 b){
-	a.x+=b.x;
-	a.y+=b.y;
-	a.z+=b.z;
-	a.w=1.;
-	return a;
-}
-Point3 sub3(Point3 a, Point3 b){
-	a.x-=b.x;
-	a.y-=b.y;
-	a.z-=b.z;
-	a.w=1.;
-	return a;
-}
-Point3 neg3(Point3 a){
-	a.x=-a.x;
-	a.y=-a.y;
-	a.z=-a.z;
-	a.w=1.;
-	return a;
-}
-Point3 div3(Point3 a, double b){
-	a.x/=b;
-	a.y/=b;
-	a.z/=b;
-	a.w=1.;
-	return a;
-}
-Point3 mul3(Point3 a, double b){
-	a.x*=b;
-	a.y*=b;
-	a.z*=b;
-	a.w=1.;
-	return a;
-}
-int eqpt3(Point3 p, Point3 q){
-	return p.x==q.x && p.y==q.y && p.z==q.z;
-}
-/*
- * Are these points closer than eps, in a relative sense
- */
-int closept3(Point3 p, Point3 q, double eps){
-	return 2.*dist3(p, q)<eps*(len3(p)+len3(q));
-}
-double dot3(Point3 p, Point3 q){
-	return p.x*q.x+p.y*q.y+p.z*q.z;
-}
-Point3 cross3(Point3 p, Point3 q){
-	Point3 r;
-	r.x=p.y*q.z-p.z*q.y;
-	r.y=p.z*q.x-p.x*q.z;
-	r.z=p.x*q.y-p.y*q.x;
-	r.w=1.;
-	return r;
-}
-double len3(Point3 p){
-	return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
-}
-double dist3(Point3 p, Point3 q){
-	p.x-=q.x;
-	p.y-=q.y;
-	p.z-=q.z;
-	return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
-}
-Point3 unit3(Point3 p){
-	double len=sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
-	p.x/=len;
-	p.y/=len;
-	p.z/=len;
-	p.w=1.;
-	return p;
-}
-Point3 midpt3(Point3 p, Point3 q){
-	p.x=.5*(p.x+q.x);
-	p.y=.5*(p.y+q.y);
-	p.z=.5*(p.z+q.z);
-	p.w=1.;
-	return p;
-}
-Point3 lerp3(Point3 p, Point3 q, double alpha){
-	p.x+=(q.x-p.x)*alpha;
-	p.y+=(q.y-p.y)*alpha;
-	p.z+=(q.z-p.z)*alpha;
-	p.w=1.;
-	return p;
-}
-/*
- * Reflect point p in the line joining p0 and p1
- */
-Point3 reflect3(Point3 p, Point3 p0, Point3 p1){
-	Point3 a, b;
-	a=sub3(p, p0);
-	b=sub3(p1, p0);
-	return add3(a, mul3(b, 2*dot3(a, b)/dot3(b, b)));
-}
-/*
- * Return the nearest point on segment [p0,p1] to point testp
- */
-Point3 nearseg3(Point3 p0, Point3 p1, Point3 testp){
-	double num, den;
-	Point3 q, r;
-	q=sub3(p1, p0);
-	r=sub3(testp, p0);
-	num=dot3(q, r);;
-	if(num<=0) return p0;
-	den=dot3(q, q);
-	if(num>=den) return p1;
-	return add3(p0, mul3(q, num/den));
-}
-/*
- * distance from point p to segment [p0,p1]
- */
-#define	SMALL	1e-8	/* what should this value be? */
-double pldist3(Point3 p, Point3 p0, Point3 p1){
-	Point3 d, e;
-	double dd, de, dsq;
-	d=sub3(p1, p0);
-	e=sub3(p, p0);
-	dd=dot3(d, d);
-	de=dot3(d, e);
-	if(dd<SMALL*SMALL) return len3(e);
-	dsq=dot3(e, e)-de*de/dd;
-	if(dsq<SMALL*SMALL) return 0;
-	return sqrt(dsq);
-}
-/*
- * vdiv3(a, b) is the magnitude of the projection of a onto b
- * measured in units of the length of b.
- * vrem3(a, b) is the component of a perpendicular to b.
- */
-double vdiv3(Point3 a, Point3 b){
-	return (a.x*b.x+a.y*b.y+a.z*b.z)/(b.x*b.x+b.y*b.y+b.z*b.z);
-}
-Point3 vrem3(Point3 a, Point3 b){
-	double quo=(a.x*b.x+a.y*b.y+a.z*b.z)/(b.x*b.x+b.y*b.y+b.z*b.z);
-	a.x-=b.x*quo;
-	a.y-=b.y*quo;
-	a.z-=b.z*quo;
-	a.w=1.;
-	return a;
-}
-/*
- * Compute face (plane) with given normal, containing a given point
- */
-Point3 pn2f3(Point3 p, Point3 n){
-	n.w=-dot3(p, n);
-	return n;
-}
-/*
- * Compute face containing three points
- */
-Point3 ppp2f3(Point3 p0, Point3 p1, Point3 p2){
-	Point3 p01, p02;
-	p01=sub3(p1, p0);
-	p02=sub3(p2, p0);
-	return pn2f3(p0, cross3(p01, p02));
-}
-/*
- * Compute point common to three faces.
- * Cramer's rule, yuk.
- */
-Point3 fff2p3(Point3 f0, Point3 f1, Point3 f2){
-	double det;
-	Point3 p;
-	det=dot3(f0, cross3(f1, f2));
-	if(fabs(det)<SMALL){	/* parallel planes, bogus answer */
-		p.x=0.;
-		p.y=0.;
-		p.z=0.;
-		p.w=0.;
-		return p;
-	}
-	p.x=(f0.w*(f2.y*f1.z-f1.y*f2.z)
-		+f1.w*(f0.y*f2.z-f2.y*f0.z)+f2.w*(f1.y*f0.z-f0.y*f1.z))/det;
-	p.y=(f0.w*(f2.z*f1.x-f1.z*f2.x)
-		+f1.w*(f0.z*f2.x-f2.z*f0.x)+f2.w*(f1.z*f0.x-f0.z*f1.x))/det;
-	p.z=(f0.w*(f2.x*f1.y-f1.x*f2.y)
-		+f1.w*(f0.x*f2.y-f2.x*f0.y)+f2.w*(f1.x*f0.y-f0.x*f1.y))/det;
-	p.w=1.;
-	return p;
-}
-/*
- * pdiv4 does perspective division to convert a projective point to affine coordinates.
- */
-Point3 pdiv4(Point3 a){
-	if(a.w==0) return a;
-	a.x/=a.w;
-	a.y/=a.w;
-	a.z/=a.w;
-	a.w=1.;
-	return a;
-}
-Point3 add4(Point3 a, Point3 b){
-	a.x+=b.x;
-	a.y+=b.y;
-	a.z+=b.z;
-	a.w+=b.w;
-	return a;
-}
-Point3 sub4(Point3 a, Point3 b){
-	a.x-=b.x;
-	a.y-=b.y;
-	a.z-=b.z;
-	a.w-=b.w;
-	return a;
-}
--- /dev/null
+++ b/sys/src/libgeometry/fmt.c
@@ -1,0 +1,28 @@
+#include <u.h>
+#include <libc.h>
+#include <geometry.h>
+
+int
+vfmt(Fmt *f)
+{
+	Point2 p;
+
+	p = va_arg(f->args, Point2);
+	return fmtprint(f, "[%g %g %g]", p.x, p.y, p.w);
+}
+
+int
+Vfmt(Fmt *f)
+{
+	Point3 p;
+
+	p = va_arg(f->args, Point3);
+	return fmtprint(f, "[%g %g %g %g]", p.x, p.y, p.z, p.w);
+}
+
+void
+GEOMfmtinstall(void)
+{
+	fmtinstall('v', vfmt);
+	fmtinstall('V', Vfmt);
+}
--- a/sys/src/libgeometry/matrix.c
+++ b/sys/src/libgeometry/matrix.c
@@ -1,106 +1,348 @@
-/*
- * ident(m)		store identity matrix in m
- * matmul(a, b)		matrix multiply a*=b
- * matmulr(a, b)	matrix multiply a=b*a
- * determinant(m)	returns det(m)
- * adjoint(m, minv)	minv=adj(m)
- * invertmat(m, minv)	invert matrix m, result in minv, returns det(m)
- *			if m is singular, minv=adj(m)
- */
 #include <u.h>
 #include <libc.h>
-#include <draw.h>
 #include <geometry.h>
-void ident(Matrix m){
-	register double *s=&m[0][0];
-	*s++=1;*s++=0;*s++=0;*s++=0;
-	*s++=0;*s++=1;*s++=0;*s++=0;
-	*s++=0;*s++=0;*s++=1;*s++=0;
-	*s++=0;*s++=0;*s++=0;*s=1;
+
+/* 2D */
+
+void
+identity(Matrix m)
+{
+	memset(m, 0, 3*3*sizeof(double));
+	m[0][0] = m[1][1] = m[2][2] = 1;
 }
-void matmul(Matrix a, Matrix b){
-	register i, j, k;
-	double sum;
+
+void
+addm(Matrix a, Matrix b)
+{
+	int i, j;
+
+	for(i = 0; i < 3; i++)
+		for(j = 0; j < 3; j++)
+			a[i][j] += b[i][j];
+}
+
+void
+subm(Matrix a, Matrix b)
+{
+	int i, j;
+
+	for(i = 0; i < 3; i++)
+		for(j = 0; j < 3; j++)
+			a[i][j] -= b[i][j];
+}
+
+void
+mulm(Matrix a, Matrix b)
+{
+	int i, j, k;
 	Matrix tmp;
-	for(i=0;i!=4;i++) for(j=0;j!=4;j++){
-		sum=0;
-		for(k=0;k!=4;k++)
-			sum+=a[i][k]*b[k][j];
-		tmp[i][j]=sum;
-	}
-	for(i=0;i!=4;i++) for(j=0;j!=4;j++)
-		a[i][j]=tmp[i][j];
+
+	for(i = 0; i < 3; i++)
+		for(j = 0; j < 3; j++){
+			tmp[i][j] = 0;
+			for(k = 0; k < 3; k++)
+				tmp[i][j] += a[i][k]*b[k][j];
+		}
+	memmove(a, tmp, 3*3*sizeof(double));
 }
-void matmulr(Matrix a, Matrix b){
-	register i, j, k;
-	double sum;
+
+void
+smulm(Matrix m, double s)
+{
+	int i, j;
+
+	for(i = 0; i < 3; i++)
+		for(j = 0; j < 3; j++)
+			m[i][j] *= s;
+}
+
+void
+transposem(Matrix m)
+{
+	int i, j;
+	double tmp;
+
+	for(i = 0; i < 3; i++)
+		for(j = i; j < 3; j++){
+			tmp = m[i][j];
+			m[i][j] = m[j][i];
+			m[j][i] = tmp;
+		}
+}
+
+double
+detm(Matrix m)
+{
+	return m[0][0]*(m[1][1]*m[2][2] - m[1][2]*m[2][1])+
+	       m[0][1]*(m[1][2]*m[2][0] - m[1][0]*m[2][2])+
+	       m[0][2]*(m[1][0]*m[2][1] - m[1][1]*m[2][0]);
+}
+
+double
+tracem(Matrix m)
+{
+	return m[0][0] + m[1][1] + m[2][2];
+}
+
+void
+adjm(Matrix m)
+{
 	Matrix tmp;
-	for(i=0;i!=4;i++) for(j=0;j!=4;j++){
-		sum=0;
-		for(k=0;k!=4;k++)
-			sum+=b[i][k]*a[k][j];
-		tmp[i][j]=sum;
-	}
-	for(i=0;i!=4;i++) for(j=0;j!=4;j++)
-		a[i][j]=tmp[i][j];
+
+	tmp[0][0] =  m[1][1]*m[2][2] - m[1][2]*m[2][1];
+	tmp[0][1] = -m[0][1]*m[2][2] + m[0][2]*m[2][1];
+	tmp[0][2] =  m[0][1]*m[1][2] - m[0][2]*m[1][1];
+	tmp[1][0] = -m[1][0]*m[2][2] + m[1][2]*m[2][0];
+	tmp[1][1] =  m[0][0]*m[2][2] - m[0][2]*m[2][0];
+	tmp[1][2] = -m[0][0]*m[1][2] + m[0][2]*m[1][0];
+	tmp[2][0] =  m[1][0]*m[2][1] - m[1][1]*m[2][0];
+	tmp[2][1] = -m[0][0]*m[2][1] + m[0][1]*m[2][0];
+	tmp[2][2] =  m[0][0]*m[1][1] - m[0][1]*m[1][0];
+	memmove(m, tmp, 3*3*sizeof(double));
 }
-/*
- * Return det(m)
- */
-double determinant(Matrix m){
-	return m[0][0]*(m[1][1]*(m[2][2]*m[3][3]-m[2][3]*m[3][2])+
-			m[1][2]*(m[2][3]*m[3][1]-m[2][1]*m[3][3])+
-			m[1][3]*(m[2][1]*m[3][2]-m[2][2]*m[3][1]))
-	      -m[0][1]*(m[1][0]*(m[2][2]*m[3][3]-m[2][3]*m[3][2])+
-			m[1][2]*(m[2][3]*m[3][0]-m[2][0]*m[3][3])+
-			m[1][3]*(m[2][0]*m[3][2]-m[2][2]*m[3][0]))
-	      +m[0][2]*(m[1][0]*(m[2][1]*m[3][3]-m[2][3]*m[3][1])+
-			m[1][1]*(m[2][3]*m[3][0]-m[2][0]*m[3][3])+
-			m[1][3]*(m[2][0]*m[3][1]-m[2][1]*m[3][0]))
-	      -m[0][3]*(m[1][0]*(m[2][1]*m[3][2]-m[2][2]*m[3][1])+
-			m[1][1]*(m[2][2]*m[3][0]-m[2][0]*m[3][2])+
-			m[1][2]*(m[2][0]*m[3][1]-m[2][1]*m[3][0]));
+
+/* Cayley-Hamilton */
+//void
+//invertm(Matrix m)
+//{
+//	Matrix m², r;
+//	double det, trm, trm²;
+//
+//	det = detm(m);
+//	if(det == 0)
+//		return;
+//	trm = tracem(m);
+//	memmove(m², m, 3*3*sizeof(double));
+//	mulm(m², m²);
+//	trm² = tracem(m²);
+//	identity(r);
+//	smulm(r, (trm*trm - trm²)/2);
+//	smulm(m, trm);
+//	subm(r, m);
+//	addm(r, m²);
+//	smulm(r, 1/det);
+//	memmove(m, r, 3*3*sizeof(double));
+//}
+
+/* Cramer's */
+void
+invm(Matrix m)
+{
+	double det;
+
+	det = detm(m);
+	if(det == 0)
+		return; /* singular matrices are not invertible */
+	adjm(m);
+	smulm(m, 1/det);
 }
-/*
- * Store the adjoint (matrix of cofactors) of m in madj.
- * Works fine even if m and madj are the same matrix.
- */
-void adjoint(Matrix m, Matrix madj){
-	double m00=m[0][0], m01=m[0][1], m02=m[0][2], m03=m[0][3];
-	double m10=m[1][0], m11=m[1][1], m12=m[1][2], m13=m[1][3];
-	double m20=m[2][0], m21=m[2][1], m22=m[2][2], m23=m[2][3];
-	double m30=m[3][0], m31=m[3][1], m32=m[3][2], m33=m[3][3];
-	madj[0][0]=m11*(m22*m33-m23*m32)+m21*(m13*m32-m12*m33)+m31*(m12*m23-m13*m22);
-	madj[0][1]=m01*(m23*m32-m22*m33)+m21*(m02*m33-m03*m32)+m31*(m03*m22-m02*m23);
-	madj[0][2]=m01*(m12*m33-m13*m32)+m11*(m03*m32-m02*m33)+m31*(m02*m13-m03*m12);
-	madj[0][3]=m01*(m13*m22-m12*m23)+m11*(m02*m23-m03*m22)+m21*(m03*m12-m02*m13);
-	madj[1][0]=m10*(m23*m32-m22*m33)+m20*(m12*m33-m13*m32)+m30*(m13*m22-m12*m23);
-	madj[1][1]=m00*(m22*m33-m23*m32)+m20*(m03*m32-m02*m33)+m30*(m02*m23-m03*m22);
-	madj[1][2]=m00*(m13*m32-m12*m33)+m10*(m02*m33-m03*m32)+m30*(m03*m12-m02*m13);
-	madj[1][3]=m00*(m12*m23-m13*m22)+m10*(m03*m22-m02*m23)+m20*(m02*m13-m03*m12);
-	madj[2][0]=m10*(m21*m33-m23*m31)+m20*(m13*m31-m11*m33)+m30*(m11*m23-m13*m21);
-	madj[2][1]=m00*(m23*m31-m21*m33)+m20*(m01*m33-m03*m31)+m30*(m03*m21-m01*m23);
-	madj[2][2]=m00*(m11*m33-m13*m31)+m10*(m03*m31-m01*m33)+m30*(m01*m13-m03*m11);
-	madj[2][3]=m00*(m13*m21-m11*m23)+m10*(m01*m23-m03*m21)+m20*(m03*m11-m01*m13);
-	madj[3][0]=m10*(m22*m31-m21*m32)+m20*(m11*m32-m12*m31)+m30*(m12*m21-m11*m22);
-	madj[3][1]=m00*(m21*m32-m22*m31)+m20*(m02*m31-m01*m32)+m30*(m01*m22-m02*m21);
-	madj[3][2]=m00*(m12*m31-m11*m32)+m10*(m01*m32-m02*m31)+m30*(m02*m11-m01*m12);
-	madj[3][3]=m00*(m11*m22-m12*m21)+m10*(m02*m21-m01*m22)+m20*(m01*m12-m02*m11);
+
+Point2
+xform(Point2 p, Matrix m)
+{
+	return (Point2){
+		p.x*m[0][0] + p.y*m[0][1] + p.w*m[0][2],
+		p.x*m[1][0] + p.y*m[1][1] + p.w*m[1][2],
+		p.x*m[2][0] + p.y*m[2][1] + p.w*m[2][2]
+	};
 }
-/*
- * Store the inverse of m in minv.
- * If m is singular, minv is instead its adjoint.
- * Returns det(m).
- * Works fine even if m and minv are the same matrix.
- */
-double invertmat(Matrix m, Matrix minv){
-	double d, dinv;
+
+/* 3D */
+
+void
+identity3(Matrix3 m)
+{
+	memset(m, 0, 4*4*sizeof(double));
+	m[0][0] = m[1][1] = m[2][2] = m[3][3] = 1;
+}
+
+void
+addm3(Matrix3 a, Matrix3 b)
+{
 	int i, j;
-	d=determinant(m);
-	adjoint(m, minv);
-	if(d!=0.){
-		dinv=1./d;
-		for(i=0;i!=4;i++) for(j=0;j!=4;j++) minv[i][j]*=dinv;
-	}
-	return d;
+
+	for(i = 0; i < 4; i++)
+		for(j = 0; j < 4; j++)
+			a[i][j] += b[i][j];
+}
+
+void
+subm3(Matrix3 a, Matrix3 b)
+{
+	int i, j;
+
+	for(i = 0; i < 4; i++)
+		for(j = 0; j < 4; j++)
+			a[i][j] -= b[i][j];
+}
+
+void
+mulm3(Matrix3 a, Matrix3 b)
+{
+	int i, j, k;
+	Matrix3 tmp;
+
+	for(i = 0; i < 4; i++)
+		for(j = 0; j < 4; j++){
+			tmp[i][j] = 0;
+			for(k = 0; k < 4; k++)
+				tmp[i][j] += a[i][k]*b[k][j];
+		}
+	memmove(a, tmp, 4*4*sizeof(double));
+}
+
+void
+smulm3(Matrix3 m, double s)
+{
+	int i, j;
+
+	for(i = 0; i < 4; i++)
+		for(j = 0; j < 4; j++)
+			m[i][j] *= s;
+}
+
+void
+transposem3(Matrix3 m)
+{
+	int i, j;
+	double tmp;
+
+	for(i = 0; i < 4; i++)
+		for(j = i; j < 4; j++){
+			tmp = m[i][j];
+			m[i][j] = m[j][i];
+			m[j][i] = tmp;
+		}
+}
+
+double
+detm3(Matrix3 m)
+{
+	return m[0][0]*(m[1][1]*(m[2][2]*m[3][3] - m[2][3]*m[3][2])+
+			m[1][2]*(m[2][3]*m[3][1] - m[2][1]*m[3][3])+
+			m[1][3]*(m[2][1]*m[3][2] - m[2][2]*m[3][1]))
+	      -m[0][1]*(m[1][0]*(m[2][2]*m[3][3] - m[2][3]*m[3][2])+
+			m[1][2]*(m[2][3]*m[3][0] - m[2][0]*m[3][3])+
+			m[1][3]*(m[2][0]*m[3][2] - m[2][2]*m[3][0]))
+	      +m[0][2]*(m[1][0]*(m[2][1]*m[3][3] - m[2][3]*m[3][1])+
+			m[1][1]*(m[2][3]*m[3][0] - m[2][0]*m[3][3])+
+			m[1][3]*(m[2][0]*m[3][1] - m[2][1]*m[3][0]))
+	      -m[0][3]*(m[1][0]*(m[2][1]*m[3][2] - m[2][2]*m[3][1])+
+			m[1][1]*(m[2][2]*m[3][0] - m[2][0]*m[3][2])+
+			m[1][2]*(m[2][0]*m[3][1] - m[2][1]*m[3][0]));
+}
+
+double
+tracem3(Matrix3 m)
+{
+	return m[0][0] + m[1][1] + m[2][2] + m[3][3];
+}
+
+void
+adjm3(Matrix3 m)
+{
+	Matrix3 tmp;
+
+	tmp[0][0]=m[1][1]*(m[2][2]*m[3][3] - m[2][3]*m[3][2])+
+		  m[2][1]*(m[1][3]*m[3][2] - m[1][2]*m[3][3])+
+		  m[3][1]*(m[1][2]*m[2][3] - m[1][3]*m[2][2]);
+	tmp[0][1]=m[0][1]*(m[2][3]*m[3][2] - m[2][2]*m[3][3])+
+		  m[2][1]*(m[0][2]*m[3][3] - m[0][3]*m[3][2])+
+		  m[3][1]*(m[0][3]*m[2][2] - m[0][2]*m[2][3]);
+	tmp[0][2]=m[0][1]*(m[1][2]*m[3][3] - m[1][3]*m[3][2])+
+		  m[1][1]*(m[0][3]*m[3][2] - m[0][2]*m[3][3])+
+		  m[3][1]*(m[0][2]*m[1][3] - m[0][3]*m[1][2]);
+	tmp[0][3]=m[0][1]*(m[1][3]*m[2][2] - m[1][2]*m[2][3])+
+		  m[1][1]*(m[0][2]*m[2][3] - m[0][3]*m[2][2])+
+		  m[2][1]*(m[0][3]*m[1][2] - m[0][2]*m[1][3]);
+	tmp[1][0]=m[1][0]*(m[2][3]*m[3][2] - m[2][2]*m[3][3])+
+		  m[2][0]*(m[1][2]*m[3][3] - m[1][3]*m[3][2])+
+		  m[3][0]*(m[1][3]*m[2][2] - m[1][2]*m[2][3]);
+	tmp[1][1]=m[0][0]*(m[2][2]*m[3][3] - m[2][3]*m[3][2])+
+		  m[2][0]*(m[0][3]*m[3][2] - m[0][2]*m[3][3])+
+		  m[3][0]*(m[0][2]*m[2][3] - m[0][3]*m[2][2]);
+	tmp[1][2]=m[0][0]*(m[1][3]*m[3][2] - m[1][2]*m[3][3])+
+		  m[1][0]*(m[0][2]*m[3][3] - m[0][3]*m[3][2])+
+		  m[3][0]*(m[0][3]*m[1][2] - m[0][2]*m[1][3]);
+	tmp[1][3]=m[0][0]*(m[1][2]*m[2][3] - m[1][3]*m[2][2])+
+		  m[1][0]*(m[0][3]*m[2][2] - m[0][2]*m[2][3])+
+		  m[2][0]*(m[0][2]*m[1][3] - m[0][3]*m[1][2]);
+	tmp[2][0]=m[1][0]*(m[2][1]*m[3][3] - m[2][3]*m[3][1])+
+		  m[2][0]*(m[1][3]*m[3][1] - m[1][1]*m[3][3])+
+		  m[3][0]*(m[1][1]*m[2][3] - m[1][3]*m[2][1]);
+	tmp[2][1]=m[0][0]*(m[2][3]*m[3][1] - m[2][1]*m[3][3])+
+		  m[2][0]*(m[0][1]*m[3][3] - m[0][3]*m[3][1])+
+		  m[3][0]*(m[0][3]*m[2][1] - m[0][1]*m[2][3]);
+	tmp[2][2]=m[0][0]*(m[1][1]*m[3][3] - m[1][3]*m[3][1])+
+		  m[1][0]*(m[0][3]*m[3][1] - m[0][1]*m[3][3])+
+		  m[3][0]*(m[0][1]*m[1][3] - m[0][3]*m[1][1]);
+	tmp[2][3]=m[0][0]*(m[1][3]*m[2][1] - m[1][1]*m[2][3])+
+		  m[1][0]*(m[0][1]*m[2][3] - m[0][3]*m[2][1])+
+		  m[2][0]*(m[0][3]*m[1][1] - m[0][1]*m[1][3]);
+	tmp[3][0]=m[1][0]*(m[2][2]*m[3][1] - m[2][1]*m[3][2])+
+		  m[2][0]*(m[1][1]*m[3][2] - m[1][2]*m[3][1])+
+		  m[3][0]*(m[1][2]*m[2][1] - m[1][1]*m[2][2]);
+	tmp[3][1]=m[0][0]*(m[2][1]*m[3][2] - m[2][2]*m[3][1])+
+		  m[2][0]*(m[0][2]*m[3][1] - m[0][1]*m[3][2])+
+		  m[3][0]*(m[0][1]*m[2][2] - m[0][2]*m[2][1]);
+	tmp[3][2]=m[0][0]*(m[1][2]*m[3][1] - m[1][1]*m[3][2])+
+		  m[1][0]*(m[0][1]*m[3][2] - m[0][2]*m[3][1])+
+		  m[3][0]*(m[0][2]*m[1][1] - m[0][1]*m[1][2]);
+	tmp[3][3]=m[0][0]*(m[1][1]*m[2][2] - m[1][2]*m[2][1])+
+		  m[1][0]*(m[0][2]*m[2][1] - m[0][1]*m[2][2])+
+		  m[2][0]*(m[0][1]*m[1][2] - m[0][2]*m[1][1]);
+	memmove(m, tmp, 4*4*sizeof(double));
+}
+
+/* Cayley-Hamilton */
+//void
+//invertm3(Matrix3 m)
+//{
+//	Matrix3 m², m³, r;
+//	double det, trm, trm², trm³;
+//
+//	det = detm3(m);
+//	if(det == 0)
+//		return;
+//	trm = tracem3(m);
+//	memmove(m³, m, 4*4*sizeof(double));
+//	mulm(m³, m³);
+//	mulm(m³, m);
+//	trm³ = tracem3(m³);
+//	memmove(m², m, 4*4*sizeof(double));
+//	mulm(m², m²);
+//	trm² = tracem3(m²);
+//	identity3(r);
+//	smulm3(r, (trm*trm*trm - 3*trm*trm² + 2*trm³)/6);
+//	smulm3(m, (trm*trm - trm²)/2);
+//	smulm3(m², trm);
+//	subm(r, m);
+//	addm(r, m²);
+//	subm(r, m³);
+//	smulm(r, 1/det);
+//	memmove(m, r, 4*4*sizeof(double));
+//}
+
+/* Cramer's */
+void
+invm3(Matrix3 m)
+{
+	double det;
+
+	det = detm3(m);
+	if(det == 0)
+		return; /* singular matrices are not invertible */
+	adjm3(m);
+	smulm3(m, 1/det);
+}
+
+Point3
+xform3(Point3 p, Matrix3 m)
+{
+	return (Point3){
+		p.x*m[0][0] + p.y*m[0][1] + p.z*m[0][2] + p.w*m[0][3],
+		p.x*m[1][0] + p.y*m[1][1] + p.z*m[1][2] + p.w*m[1][3],
+		p.x*m[2][0] + p.y*m[2][1] + p.z*m[2][2] + p.w*m[2][3],
+		p.x*m[3][0] + p.y*m[3][1] + p.z*m[3][2] + p.w*m[3][3],
+	};
 }
--- a/sys/src/libgeometry/mkfile
+++ b/sys/src/libgeometry/mkfile
@@ -1,23 +1,22 @@
 </$objtype/mkfile
 
 LIB=/$objtype/lib/libgeometry.a
+
 OFILES=\
-	arith3.$O\
+	point.$O\
 	matrix.$O\
-	qball.$O\
 	quaternion.$O\
-	transform.$O\
-	tstack.$O\
+	rframe.$O\
+	triangle.$O\
+	utils.$O\
+	fmt.$O\
 
-HFILES=/sys/include/geometry.h
+HFILES=\
+	/sys/include/geometry.h
 
-</sys/src/cmd/mksyslib
-
 UPDATE=\
 	mkfile\
 	$HFILES\
 	${OFILES:%.$O=%.c}\
-	${LIB:/$objtype/%=/386/%}\
 
-listing:V:
-	pr mkfile $HFILES $CFILES|lp -du
+</sys/src/cmd/mksyslib
--- /dev/null
+++ b/sys/src/libgeometry/point.c
@@ -1,0 +1,200 @@
+#include <u.h>
+#include <libc.h>
+#include <geometry.h>
+
+/* 2D */
+
+Point2
+Pt2(double x, double y, double w)
+{
+	return (Point2){x, y, w};
+}
+
+Point2
+Vec2(double x, double y)
+{
+	return (Point2){x, y, 0};
+}
+
+Point2
+addpt2(Point2 a, Point2 b)
+{
+	return Pt2(a.x+b.x, a.y+b.y, a.w+b.w);
+}
+
+Point2
+subpt2(Point2 a, Point2 b)
+{
+	return Pt2(a.x-b.x, a.y-b.y, a.w-b.w);
+}
+
+Point2
+mulpt2(Point2 p, double s)
+{
+	return Pt2(p.x*s, p.y*s, p.w*s);
+}
+
+Point2
+divpt2(Point2 p, double s)
+{
+	return Pt2(p.x/s, p.y/s, p.w/s);
+}
+
+Point2
+lerp2(Point2 a, Point2 b, double t)
+{
+	t = fclamp(t, 0, 1);
+	return Pt2(
+		flerp(a.x, b.x, t),
+		flerp(a.y, b.y, t),
+		flerp(a.w, b.w, t)
+	);
+}
+
+double
+dotvec2(Point2 a, Point2 b)
+{
+	return a.x*b.x + a.y*b.y;
+}
+
+double
+vec2len(Point2 v)
+{
+	return sqrt(dotvec2(v, v));
+}
+
+Point2
+normvec2(Point2 v)
+{
+	double len;
+
+	len = vec2len(v);
+	if(len == 0)
+		return Pt2(0,0,0);
+	return Pt2(v.x/len, v.y/len, 0);
+}
+
+/*
+ * the edge function, from:
+ *
+ * Juan Pineda, “A Parallel Algorithm for Polygon Rasterization”,
+ * Computer Graphics, Vol. 22, No. 8, August 1988
+ *
+ * comparison of a point p with an edge [e0 e1]
+ * p to the right: +
+ * p to the left: -
+ * p on the edge: 0
+ */
+int
+edgeptcmp(Point2 e0, Point2 e1, Point2 p)
+{
+	Point3 e0p, e01, r;
+
+	p = subpt2(p, e0);
+	e1 = subpt2(e1, e0);
+	e0p = Vec3(p.x,p.y,0);
+	e01 = Vec3(e1.x,e1.y,0);
+	r = crossvec3(e0p, e01);
+
+	/* clamp to avoid overflow */
+	return fclamp(r.z, -1, 1); /* e0.x*e1.y - e0.y*e1.x */
+}
+
+/*
+ * (PNPOLY) algorithm by W. Randolph Franklin
+ */
+int
+ptinpoly(Point2 p, Point2 *pts, ulong npts)
+{
+	int i, j, c;
+
+	for(i = c = 0, j = npts-1; i < npts; j = i++)
+		if(p.y < pts[i].y != p.y < pts[j].y &&
+			p.x < (pts[j].x - pts[i].x) * (p.y - pts[i].y)/(pts[j].y - pts[i].y) + pts[i].x)
+		c ^= 1;
+	return c;
+}
+
+/* 3D */
+
+Point3
+Pt3(double x, double y, double z, double w)
+{
+	return (Point3){x, y, z, w};
+}
+
+Point3
+Vec3(double x, double y, double z)
+{
+	return (Point3){x, y, z, 0};
+}
+
+Point3
+addpt3(Point3 a, Point3 b)
+{
+	return Pt3(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
+}
+
+Point3
+subpt3(Point3 a, Point3 b)
+{
+	return Pt3(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w);
+}
+
+Point3
+mulpt3(Point3 p, double s)
+{
+	return Pt3(p.x*s, p.y*s, p.z*s, p.w*s);
+}
+
+Point3
+divpt3(Point3 p, double s)
+{
+	return Pt3(p.x/s, p.y/s, p.z/s, p.w/s);
+}
+
+Point3
+lerp3(Point3 a, Point3 b, double t)
+{
+	t = fclamp(t, 0, 1);
+	return Pt3(
+		flerp(a.x, b.x, t),
+		flerp(a.y, b.y, t),
+		flerp(a.z, b.z, t),
+		flerp(a.w, b.w, t)
+	);
+}
+
+double
+dotvec3(Point3 a, Point3 b)
+{
+	return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+Point3
+crossvec3(Point3 a, Point3 b)
+{
+	return Pt3(
+		a.y*b.z - a.z*b.y,
+		a.z*b.x - a.x*b.z,
+		a.x*b.y - a.y*b.x,
+		0
+	);
+}
+
+double
+vec3len(Point3 v)
+{
+	return sqrt(dotvec3(v, v));
+}
+
+Point3
+normvec3(Point3 v)
+{
+	double len;
+
+	len = vec3len(v);
+	if(len == 0)
+		return Pt3(0,0,0,0);
+	return Pt3(v.x/len, v.y/len, v.z/len, 0);
+}
--- a/sys/src/libgeometry/qball.c
+++ /dev/null
@@ -1,65 +1,0 @@
-/*
- * Ken Shoemake's Quaternion rotation controller
- */
-#include <u.h>
-#include <libc.h>
-#include <draw.h>
-#include <event.h>
-#include <geometry.h>
-#define	BORDER	4
-static Point ctlcen;		/* center of qball */
-static int ctlrad;		/* radius of qball */
-static Quaternion *axis;	/* constraint plane orientation, 0 if none */
-/*
- * Convert a mouse point into a unit quaternion, flattening if
- * constrained to a particular plane.
- */
-static Quaternion mouseq(Point p){
-	double qx=(double)(p.x-ctlcen.x)/ctlrad;
-	double qy=(double)(p.y-ctlcen.y)/ctlrad;
-	double rsq=qx*qx+qy*qy;
-	double l;
-	Quaternion q;
-	if(rsq>1){
-		rsq=sqrt(rsq);
-		q.r=0.;
-		q.i=qx/rsq;
-		q.j=qy/rsq;
-		q.k=0.;
-	}
-	else{
-		q.r=0.;
-		q.i=qx;
-		q.j=qy;
-		q.k=sqrt(1.-rsq);
-	}
-	if(axis){
-		l=q.i*axis->i+q.j*axis->j+q.k*axis->k;
-		q.i-=l*axis->i;
-		q.j-=l*axis->j;
-		q.k-=l*axis->k;
-		l=sqrt(q.i*q.i+q.j*q.j+q.k*q.k);
-		if(l!=0.){
-			q.i/=l;
-			q.j/=l;
-			q.k/=l;
-		}
-	}
-	return q;
-}
-void qball(Rectangle r, Mouse *m, Quaternion *result, void (*redraw)(void), Quaternion *ap){
-	Quaternion q, down;
-	Point rad;
-	axis=ap;
-	ctlcen=divpt(addpt(r.min, r.max), 2);
-	rad=divpt(subpt(r.max, r.min), 2);
-	ctlrad=(rad.x<rad.y?rad.x:rad.y)-BORDER;
-	down=qinv(mouseq(m->xy));
-	q=*result;
-	for(;;){
-		*m=emouse();
-		if(!m->buttons) break;
-		*result=qmul(q, qmul(down, mouseq(m->xy)));
-		(*redraw)();
-	}
-}
--- a/sys/src/libgeometry/quaternion.c
+++ b/sys/src/libgeometry/quaternion.c
@@ -1,242 +1,108 @@
-/*
- * Quaternion arithmetic:
- *	qadd(q, r)	returns q+r
- *	qsub(q, r)	returns q-r
- *	qneg(q)		returns -q
- *	qmul(q, r)	returns q*r
- *	qdiv(q, r)	returns q/r, can divide check.
- *	qinv(q)		returns 1/q, can divide check.
- *	double qlen(p)	returns modulus of p
- *	qunit(q)	returns a unit quaternion parallel to q
- * The following only work on unit quaternions and rotation matrices:
- *	slerp(q, r, a)	returns q*(r*q^-1)^a
- *	qmid(q, r)	slerp(q, r, .5) 
- *	qsqrt(q)	qmid(q, (Quaternion){1,0,0,0})
- *	qtom(m, q)	converts a unit quaternion q into a rotation matrix m
- *	mtoq(m)		returns a quaternion equivalent to a rotation matrix m
- */
 #include <u.h>
 #include <libc.h>
-#include <draw.h>
 #include <geometry.h>
-void qtom(Matrix m, Quaternion q){
-#ifndef new
-	m[0][0]=1-2*(q.j*q.j+q.k*q.k);
-	m[0][1]=2*(q.i*q.j+q.r*q.k);
-	m[0][2]=2*(q.i*q.k-q.r*q.j);
-	m[0][3]=0;
-	m[1][0]=2*(q.i*q.j-q.r*q.k);
-	m[1][1]=1-2*(q.i*q.i+q.k*q.k);
-	m[1][2]=2*(q.j*q.k+q.r*q.i);
-	m[1][3]=0;
-	m[2][0]=2*(q.i*q.k+q.r*q.j);
-	m[2][1]=2*(q.j*q.k-q.r*q.i);
-	m[2][2]=1-2*(q.i*q.i+q.j*q.j);
-	m[2][3]=0;
-	m[3][0]=0;
-	m[3][1]=0;
-	m[3][2]=0;
-	m[3][3]=1;
-#else
-	/*
-	 * Transcribed from Ken Shoemake's new code -- not known to work
-	 */
-	double Nq = q.r*q.r+q.i*q.i+q.j*q.j+q.k*q.k;
-	double s = (Nq > 0.0) ? (2.0 / Nq) : 0.0;
-	double xs = q.i*s,		ys = q.j*s,		zs = q.k*s;
-	double wx = q.r*xs,		wy = q.r*ys,		wz = q.r*zs;
-	double xx = q.i*xs,		xy = q.i*ys,		xz = q.i*zs;
-	double yy = q.j*ys,		yz = q.j*zs,		zz = q.k*zs;
-	m[0][0] = 1.0 - (yy + zz); m[1][0] = xy + wz;         m[2][0] = xz - wy;
-	m[0][1] = xy - wz;         m[1][1] = 1.0 - (xx + zz); m[2][1] = yz + wx;
-	m[0][2] = xz + wy;         m[1][2] = yz - wx;         m[2][2] = 1.0 - (xx + yy);
-	m[0][3] = m[1][3] = m[2][3] = m[3][0] = m[3][1] = m[3][2] = 0.0;
-	m[3][3] = 1.0;
-#endif
+
+Quaternion
+Quat(double r, double i, double j, double k)
+{
+	return (Quaternion){r, i, j, k};
 }
-Quaternion mtoq(Matrix mat){
-#ifndef new
-#define	EPS	1.387778780781445675529539585113525e-17	/* 2^-56 */
-	double t;
-	Quaternion q;
-	q.r=0.;
-	q.i=0.;
-	q.j=0.;
-	q.k=1.;
-	if((t=.25*(1+mat[0][0]+mat[1][1]+mat[2][2]))>EPS){
-		q.r=sqrt(t);
-		t=4*q.r;
-		q.i=(mat[1][2]-mat[2][1])/t;
-		q.j=(mat[2][0]-mat[0][2])/t;
-		q.k=(mat[0][1]-mat[1][0])/t;
-	}
-	else if((t=-.5*(mat[1][1]+mat[2][2]))>EPS){
-		q.i=sqrt(t);
-		t=2*q.i;
-		q.j=mat[0][1]/t;
-		q.k=mat[0][2]/t;
-	}
-	else if((t=.5*(1-mat[2][2]))>EPS){
-		q.j=sqrt(t);
-		q.k=mat[1][2]/(2*q.j);
-	}
-	return q;
-#else
-	/*
-	 * Transcribed from Ken Shoemake's new code -- not known to work
-	 */
-	/* This algorithm avoids near-zero divides by looking for a large
-	 * component -- first r, then i, j, or k.  When the trace is greater than zero,
-	 * |r| is greater than 1/2, which is as small as a largest component can be.
-	 * Otherwise, the largest diagonal entry corresponds to the largest of |i|,
-	 * |j|, or |k|, one of which must be larger than |r|, and at least 1/2.
-	 */
-	Quaternion qu;
-	double tr, s;
-	
-	tr = mat[0][0] + mat[1][1] + mat[2][2];
-	if (tr >= 0.0) {
-		s = sqrt(tr + mat[3][3]);
-		qu.r = s*0.5;
-		s = 0.5 / s;
-		qu.i = (mat[2][1] - mat[1][2]) * s;
-		qu.j = (mat[0][2] - mat[2][0]) * s;
-		qu.k = (mat[1][0] - mat[0][1]) * s;
-	}
-	else {
-		int i = 0;
-		if (mat[1][1] > mat[0][0]) i = 1;
-		if (mat[2][2] > mat[i][i]) i = 2;
-		switch(i){
-		case 0:
-			s = sqrt( (mat[0][0] - (mat[1][1]+mat[2][2])) + mat[3][3] );
-			qu.i = s*0.5;
-			s = 0.5 / s;
-			qu.j = (mat[0][1] + mat[1][0]) * s;
-			qu.k = (mat[2][0] + mat[0][2]) * s;
-			qu.r = (mat[2][1] - mat[1][2]) * s;
-			break;
-		case 1:
-			s = sqrt( (mat[1][1] - (mat[2][2]+mat[0][0])) + mat[3][3] );
-			qu.j = s*0.5;
-			s = 0.5 / s;
-			qu.k = (mat[1][2] + mat[2][1]) * s;
-			qu.i = (mat[0][1] + mat[1][0]) * s;
-			qu.r = (mat[0][2] - mat[2][0]) * s;
-			break;
-		case 2:
-			s = sqrt( (mat[2][2] - (mat[0][0]+mat[1][1])) + mat[3][3] );
-			qu.k = s*0.5;
-			s = 0.5 / s;
-			qu.i = (mat[2][0] + mat[0][2]) * s;
-			qu.j = (mat[1][2] + mat[2][1]) * s;
-			qu.r = (mat[1][0] - mat[0][1]) * s;
-			break;
-		}
-	}
-	if (mat[3][3] != 1.0){
-		s=1/sqrt(mat[3][3]);
-		qu.r*=s;
-		qu.i*=s;
-		qu.j*=s;
-		qu.k*=s;
-	}
-	return (qu);
-#endif
+
+Quaternion
+Quatvec(double s, Point3 v)
+{
+	return (Quaternion){s, v.x, v.y, v.z};
 }
-Quaternion qadd(Quaternion q, Quaternion r){
-	q.r+=r.r;
-	q.i+=r.i;
-	q.j+=r.j;
-	q.k+=r.k;
-	return q;
+
+Quaternion
+addq(Quaternion a, Quaternion b)
+{
+	return Quat(a.r+b.r, a.i+b.i, a.j+b.j, a.k+b.k);
 }
-Quaternion qsub(Quaternion q, Quaternion r){
-	q.r-=r.r;
-	q.i-=r.i;
-	q.j-=r.j;
-	q.k-=r.k;
-	return q;
+
+Quaternion
+subq(Quaternion a, Quaternion b)
+{
+	return Quat(a.r-b.r, a.i-b.i, a.j-b.j, a.k-b.k);
 }
-Quaternion qneg(Quaternion q){
-	q.r=-q.r;
-	q.i=-q.i;
-	q.j=-q.j;
-	q.k=-q.k;
-	return q;
+
+Quaternion
+mulq(Quaternion q, Quaternion r)
+{
+	Point3 qv, rv, tmp;
+
+	qv = Vec3(q.i, q.j, q.k);
+	rv = Vec3(r.i, r.j, r.k);
+	tmp = addpt3(addpt3(mulpt3(rv, q.r), mulpt3(qv, r.r)), crossvec3(qv, rv));
+	return Quatvec(q.r*r.r - dotvec3(qv, rv), tmp);
 }
-Quaternion qmul(Quaternion q, Quaternion r){
-	Quaternion s;
-	s.r=q.r*r.r-q.i*r.i-q.j*r.j-q.k*r.k;
-	s.i=q.r*r.i+r.r*q.i+q.j*r.k-q.k*r.j;
-	s.j=q.r*r.j+r.r*q.j+q.k*r.i-q.i*r.k;
-	s.k=q.r*r.k+r.r*q.k+q.i*r.j-q.j*r.i;
-	return s;
+
+Quaternion
+smulq(Quaternion q, double s)
+{
+	return Quat(q.r*s, q.i*s, q.j*s, q.k*s);
 }
-Quaternion qdiv(Quaternion q, Quaternion r){
-	return qmul(q, qinv(r));
+
+Quaternion
+sdivq(Quaternion q, double s)
+{
+	return Quat(q.r/s, q.i/s, q.j/s, q.k/s);
 }
-Quaternion qunit(Quaternion q){
-	double l=qlen(q);
-	q.r/=l;
-	q.i/=l;
-	q.j/=l;
-	q.k/=l;
-	return q;
+
+double
+dotq(Quaternion q, Quaternion r)
+{
+	return q.r*r.r + q.i*r.i + q.j*r.j + q.k*r.k;
 }
-/*
- * Bug?: takes no action on divide check
- */
-Quaternion qinv(Quaternion q){
-	double l=q.r*q.r+q.i*q.i+q.j*q.j+q.k*q.k;
-	q.r/=l;
-	q.i=-q.i/l;
-	q.j=-q.j/l;
-	q.k=-q.k/l;
-	return q;
+
+Quaternion
+invq(Quaternion q)
+{
+	double len²;
+
+	len² = dotq(q, q);
+	if(len² == 0)
+		return Quat(0,0,0,0);
+	return Quat(q.r/len², -q.i/len², -q.j/len², -q.k/len²);
 }
-double qlen(Quaternion p){
-	return sqrt(p.r*p.r+p.i*p.i+p.j*p.j+p.k*p.k);
+
+double
+qlen(Quaternion q)
+{
+	return sqrt(dotq(q, q));
 }
-Quaternion slerp(Quaternion q, Quaternion r, double a){
-	double u, v, ang, s;
-	double dot=q.r*r.r+q.i*r.i+q.j*r.j+q.k*r.k;
-	ang=dot<-1?PI:dot>1?0:acos(dot); /* acos gives NaN for dot slightly out of range */
-	s=sin(ang);
-	if(s==0) return ang<PI/2?q:r;
-	u=sin((1-a)*ang)/s;
-	v=sin(a*ang)/s;
-	q.r=u*q.r+v*r.r;
-	q.i=u*q.i+v*r.i;
-	q.j=u*q.j+v*r.j;
-	q.k=u*q.k+v*r.k;
-	return q;
+
+Quaternion
+normq(Quaternion q)
+{
+	return sdivq(q, qlen(q));
 }
+
 /*
- * Only works if qlen(q)==qlen(r)==1
+ * based on the implementation from:
+ *
+ * Jonathan Blow, “Understanding Slerp, Then Not Using it”,
+ * The Inner Product, April 2004.
  */
-Quaternion qmid(Quaternion q, Quaternion r){
-	double l;
-	q=qadd(q, r);
-	l=qlen(q);
-	if(l<1e-12){
-		q.r=r.i;
-		q.i=-r.r;
-		q.j=r.k;
-		q.k=-r.j;
-	}
-	else{
-		q.r/=l;
-		q.i/=l;
-		q.j/=l;
-		q.k/=l;
-	}
-	return q;
+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(θ) */
 }
-/*
- * Only works if qlen(q)==1
- */
-static Quaternion qident={1,0,0,0};
-Quaternion qsqrt(Quaternion q){
-	return qmid(q, qident);
+
+Point3
+qrotate(Point3 p, Point3 axis, double θ)
+{
+	Quaternion qaxis, qr;
+
+	θ /= 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);
 }
--- /dev/null
+++ b/sys/src/libgeometry/rframe.c
@@ -1,0 +1,51 @@
+#include <u.h>
+#include <libc.h>
+#include <geometry.h>
+
+Point2
+rframexform(Point2 p, RFrame rf)
+{
+	Matrix m = {
+		rf.bx.x, rf.bx.y, -dotvec2(rf.bx, rf.p),
+		rf.by.x, rf.by.y, -dotvec2(rf.by, rf.p),
+		0, 0, 1,
+	};
+	return xform(p, m);
+}
+
+Point3
+rframexform3(Point3 p, RFrame3 rf)
+{
+	Matrix3 m = {
+		rf.bx.x, rf.bx.y, rf.bx.z, -dotvec3(rf.bx, rf.p),
+		rf.by.x, rf.by.y, rf.by.z, -dotvec3(rf.by, rf.p),
+		rf.bz.x, rf.bz.y, rf.bz.z, -dotvec3(rf.bz, rf.p),
+		0, 0, 0, 1,
+	};
+	return xform3(p, m);
+}
+
+Point2
+invrframexform(Point2 p, RFrame rf)
+{
+	Matrix m = {
+		rf.bx.x, rf.bx.y, -dotvec2(rf.bx, rf.p),
+		rf.by.x, rf.by.y, -dotvec2(rf.by, rf.p),
+		0, 0, 1,
+	};
+	invm(m);
+	return xform(p, m);
+}
+
+Point3
+invrframexform3(Point3 p, RFrame3 rf)
+{
+	Matrix3 m = {
+		rf.bx.x, rf.bx.y, rf.bx.z, -dotvec3(rf.bx, rf.p),
+		rf.by.x, rf.by.y, rf.by.z, -dotvec3(rf.by, rf.p),
+		rf.bz.x, rf.bz.y, rf.bz.z, -dotvec3(rf.bz, rf.p),
+		0, 0, 0, 1,
+	};
+	invm3(m);
+	return xform3(p, m);
+}
--- a/sys/src/libgeometry/transform.c
+++ /dev/null
@@ -1,75 +1,0 @@
-/*
- * The following routines transform points and planes from one space
- * to another.  Points and planes are represented by their
- * homogeneous coordinates, stored in variables of type Point3.
- */
-#include <u.h>
-#include <libc.h>
-#include <draw.h>
-#include <geometry.h>
-/*
- * Transform point p.
- */
-Point3 xformpoint(Point3 p, Space *to, Space *from){
-	Point3 q, r;
-	register double *m;
-	if(from){
-		m=&from->t[0][0];
-		q.x=*m++*p.x; q.x+=*m++*p.y; q.x+=*m++*p.z; q.x+=*m++*p.w;
-		q.y=*m++*p.x; q.y+=*m++*p.y; q.y+=*m++*p.z; q.y+=*m++*p.w;
-		q.z=*m++*p.x; q.z+=*m++*p.y; q.z+=*m++*p.z; q.z+=*m++*p.w;
-		q.w=*m++*p.x; q.w+=*m++*p.y; q.w+=*m++*p.z; q.w+=*m  *p.w;
-	}
-	else
-		q=p;
-	if(to){
-		m=&to->tinv[0][0];
-		r.x=*m++*q.x; r.x+=*m++*q.y; r.x+=*m++*q.z; r.x+=*m++*q.w;
-		r.y=*m++*q.x; r.y+=*m++*q.y; r.y+=*m++*q.z; r.y+=*m++*q.w;
-		r.z=*m++*q.x; r.z+=*m++*q.y; r.z+=*m++*q.z; r.z+=*m++*q.w;
-		r.w=*m++*q.x; r.w+=*m++*q.y; r.w+=*m++*q.z; r.w+=*m  *q.w;
-	}
-	else
-		r=q;
-	return r;
-}
-/*
- * Transform point p with perspective division.
- */
-Point3 xformpointd(Point3 p, Space *to, Space *from){
-	p=xformpoint(p, to, from);
-	if(p.w!=0){
-		p.x/=p.w;
-		p.y/=p.w;
-		p.z/=p.w;
-		p.w=1;
-	}
-	return p;
-}
-/*
- * Transform plane p -- same as xformpoint, except multiply on the
- * other side by the inverse matrix.
- */
-Point3 xformplane(Point3 p, Space *to, Space *from){
-	Point3 q, r;
-	register double *m;
-	if(from){
-		m=&from->tinv[0][0];
-		q.x =*m++*p.x; q.y =*m++*p.x; q.z =*m++*p.x; q.w =*m++*p.x;
-		q.x+=*m++*p.y; q.y+=*m++*p.y; q.z+=*m++*p.y; q.w+=*m++*p.y;
-		q.x+=*m++*p.z; q.y+=*m++*p.z; q.z+=*m++*p.z; q.w+=*m++*p.z;
-		q.x+=*m++*p.w; q.y+=*m++*p.w; q.z+=*m++*p.w; q.w+=*m  *p.w;
-	}
-	else
-		q=p;
-	if(to){
-		m=&to->t[0][0];
-		r.x =*m++*q.x; r.y =*m++*q.x; r.z =*m++*q.x; r.w =*m++*q.x;
-		r.x+=*m++*q.y; r.y+=*m++*q.y; r.z+=*m++*q.y; r.w+=*m++*q.y;
-		r.x+=*m++*q.z; r.y+=*m++*q.z; r.z+=*m++*q.z; r.w+=*m++*q.z;
-		r.x+=*m++*q.w; r.y+=*m++*q.w; r.z+=*m++*q.w; r.w+=*m  *q.w;
-	}
-	else
-		r=q;
-	return r;
-}
--- /dev/null
+++ b/sys/src/libgeometry/triangle.c
@@ -1,0 +1,40 @@
+#include <u.h>
+#include <libc.h>
+#include <geometry.h>
+
+/* 2D */
+
+Point2
+centroid(Triangle2 t)
+{
+	return divpt2(addpt2(t.p0, addpt2(t.p1, t.p2)), 3);
+}
+
+/*
+ * based on the implementation from:
+ *
+ * Dmitry V. Sokolov, “Tiny Renderer: Lesson 2”,
+ * https://github.com/ssloy/tinyrenderer/wiki/Lesson-2:-Triangle-rasterization-and-back-face-culling
+ */
+Point3
+barycoords(Triangle2 t, Point2 p)
+{
+	Point2 p0p1 = subpt2(t.p1, t.p0);
+	Point2 p0p2 = subpt2(t.p2, t.p0);
+	Point2 pp0  = subpt2(t.p0, p);
+
+	Point3 v = crossvec3(Vec3(p0p2.x, p0p1.x, pp0.x), Vec3(p0p2.y, p0p1.y, pp0.y));
+
+	/* handle degenerate triangles—i.e. the ones where every point lies on the same line */
+	if(fabs(v.z) < 1)
+		return Pt3(-1,-1,-1,1);
+	return Pt3(1 - (v.x + v.y)/v.z, v.y/v.z, v.x/v.z, 1);
+}
+
+/* 3D */
+
+Point3
+centroid3(Triangle3 t)
+{
+	return divpt3(addpt3(t.p0, addpt3(t.p1, t.p2)), 3);
+}
--- a/sys/src/libgeometry/tstack.c
+++ /dev/null
@@ -1,169 +1,0 @@
-/*% cc -gpc %
- * These transformation routines maintain stacks of transformations
- * and their inverses.  
- * t=pushmat(t)		push matrix stack
- * t=popmat(t)		pop matrix stack
- * rot(t, a, axis)	multiply stack top by rotation
- * qrot(t, q)		multiply stack top by rotation, q is unit quaternion
- * scale(t, x, y, z)	multiply stack top by scale
- * move(t, x, y, z)	multiply stack top by translation
- * xform(t, m)		multiply stack top by m
- * ixform(t, m, inv)	multiply stack top by m.  inv is the inverse of m.
- * look(t, e, l, u)	multiply stack top by viewing transformation
- * persp(t, fov, n, f)	multiply stack top by perspective transformation
- * viewport(t, r, aspect)
- *			multiply stack top by window->viewport transformation.
- */
-#include <u.h>
-#include <libc.h>
-#include <draw.h>
-#include <geometry.h>
-Space *pushmat(Space *t){
-	Space *v;
-	v=malloc(sizeof(Space));
-	if(t==0){
-		ident(v->t);
-		ident(v->tinv);
-	}
-	else
-		*v=*t;
-	v->next=t;
-	return v;
-}
-Space *popmat(Space *t){
-	Space *v;
-	if(t==0) return 0;
-	v=t->next;
-	free(t);
-	return v;
-}
-void rot(Space *t, double theta, int axis){
-	double s=sin(radians(theta)), c=cos(radians(theta));
-	Matrix m, inv;
-	register i=(axis+1)%3, j=(axis+2)%3;
-	ident(m);
-	m[i][i] = c;
-	m[i][j] = -s;
-	m[j][i] = s;
-	m[j][j] = c;
-	ident(inv);
-	inv[i][i] = c;
-	inv[i][j] = s;
-	inv[j][i] = -s;
-	inv[j][j] = c;
-	ixform(t, m, inv);
-}
-void qrot(Space *t, Quaternion q){
-	Matrix m, inv;
-	int i, j;
-	qtom(m, q);
-	for(i=0;i!=4;i++) for(j=0;j!=4;j++) inv[i][j]=m[j][i];
-	ixform(t, m, inv);
-}
-void scale(Space *t, double x, double y, double z){
-	Matrix m, inv;
-	ident(m);
-	m[0][0]=x;
-	m[1][1]=y;
-	m[2][2]=z;
-	ident(inv);
-	inv[0][0]=1/x;
-	inv[1][1]=1/y;
-	inv[2][2]=1/z;
-	ixform(t, m, inv);
-}
-void move(Space *t, double x, double y, double z){
-	Matrix m, inv;
-	ident(m);
-	m[0][3]=x;
-	m[1][3]=y;
-	m[2][3]=z;
-	ident(inv);
-	inv[0][3]=-x;
-	inv[1][3]=-y;
-	inv[2][3]=-z;
-	ixform(t, m, inv);
-}
-void xform(Space *t, Matrix m){
-	Matrix inv;
-	if(invertmat(m, inv)==0) return;
-	ixform(t, m, inv);
-}
-void ixform(Space *t, Matrix m, Matrix inv){
-	matmul(t->t, m);
-	matmulr(t->tinv, inv);
-}
-/*
- * multiply the top of the matrix stack by a view-pointing transformation
- * with the eyepoint at e, looking at point l, with u at the top of the screen.
- * The coordinate system is deemed to be right-handed.
- * The generated transformation transforms this view into a view from
- * the origin, looking in the positive y direction, with the z axis pointing up,
- * and x to the right.
- */
-void look(Space *t, Point3 e, Point3 l, Point3 u){
-	Matrix m, inv;
-	Point3 r;
-	l=unit3(sub3(l, e));
-	u=unit3(vrem3(sub3(u, e), l));
-	r=cross3(l, u);
-	/* make the matrix to transform from (rlu) space to (xyz) space */
-	ident(m);
-	m[0][0]=r.x; m[0][1]=r.y; m[0][2]=r.z;
-	m[1][0]=l.x; m[1][1]=l.y; m[1][2]=l.z;
-	m[2][0]=u.x; m[2][1]=u.y; m[2][2]=u.z;
-	ident(inv);
-	inv[0][0]=r.x; inv[0][1]=l.x; inv[0][2]=u.x;
-	inv[1][0]=r.y; inv[1][1]=l.y; inv[1][2]=u.y;
-	inv[2][0]=r.z; inv[2][1]=l.z; inv[2][2]=u.z;
-	ixform(t, m, inv);
-	move(t, -e.x, -e.y, -e.z);
-}
-/*
- * generate a transformation that maps the frustum with apex at the origin,
- * apex angle=fov and clipping planes y=n and y=f into the double-unit cube.
- * plane y=n maps to y'=-1, y=f maps to y'=1
- */
-int persp(Space *t, double fov, double n, double f){
-	Matrix m;
-	double z;
-	if(n<=0 || f<=n || fov<=0 || 180<=fov) /* really need f!=n && sin(v)!=0 */
-		return -1;
-	z=1/tan(radians(fov)/2);
-	m[0][0]=z; m[0][1]=0;           m[0][2]=0; m[0][3]=0;
-	m[1][0]=0; m[1][1]=(f+n)/(f-n); m[1][2]=0; m[1][3]=f*(1-m[1][1]);
-	m[2][0]=0; m[2][1]=0;           m[2][2]=z; m[2][3]=0;
-	m[3][0]=0; m[3][1]=1;           m[3][2]=0; m[3][3]=0;
-	xform(t, m);
-	return 0;
-}
-/*
- * Map the unit-cube window into the given screen viewport.
- * r has min at the top left, max just outside the lower right.  Aspect is the
- * aspect ratio (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 y axis, with x to the left and z up.  The viewport
- * has x increasing to the right and y increasing down.  The window's y coordinates
- * are mapped, unchanged, into the viewport's z coordinates.
- */
-void viewport(Space *t, Rectangle r, double aspect){
-	Matrix m;
-	double xc, yc, wid, hgt, scale;
-	xc=.5*(r.min.x+r.max.x);
-	yc=.5*(r.min.y+r.max.y);
-	wid=(r.max.x-r.min.x)*aspect;
-	hgt=r.max.y-r.min.y;
-	scale=.5*(wid<hgt?wid:hgt);
-	ident(m);
-	m[0][0]=scale;
-	m[0][3]=xc;
-	m[1][1]=0;
-	m[1][2]=-scale;
-	m[1][3]=yc;
-	m[2][1]=1;
-	m[2][2]=0;
-	/* should get inverse by hand */
-	xform(t, m);
-}
--- /dev/null
+++ b/sys/src/libgeometry/utils.c
@@ -1,0 +1,15 @@
+#include <u.h>
+#include <libc.h>
+#include <geometry.h>
+
+double
+flerp(double a, double b, double t)
+{
+	return a + (b - a)*t;
+}
+
+double
+fclamp(double n, double min, double max)
+{
+	return n < min? min: n > max? max: n;
+}

Re: [9front] [patchset] libgeometry revamp

[cinap_lenrek]

wow, thats impressive.

afaik, there where no changes made to libgeometry since
the start of 9front with one excetion where someone realized
<stdio.h> wasnt needed.

so you are pretty much the maintainer now :)

make sure you got commit access (if you havnt already).

lgtm.

--
cinap

Re: [9front] [patchset] libgeometry revamp

[rgl]

thanks cinap!

yeah, i was expecting that in a way since i've
never seen any program that actually uses it
—in- nor outside of the system. i'll take care
of it :)

it will serve as the basis for some more
advanced projects i have in mind.


tty soon!

-rodri

Re: [9front] [patchset] libgeometry revamp

[Anthony Martin]

rgl@antares-labs.eu once said:
> yeah, i was expecting that in a way since i've
> never seen any program that actually uses it
> —in- nor outside of the system.

The only one I know of is cpr(9) from 2ed.
I think Tom Duff wrote it. It's an interesting
little program.

https://man.cat-v.org/plan_9_2nd_ed/9/cpr

Cheers,
  Anthony

Re: [9front] [patchset] libgeometry revamp

[rgl]

hi anthony,

> The only one I know of is cpr(9) from 2ed.
> I think Tom Duff wrote it. It's an interesting
> little program.

i was completely unaware of this program. very
interesting indeed, and it makes sense it was
td's because the source code looks very much his.
although rob wrote tight as well, but not *that*
tight i would think.


thanks for sharing.

-rodri

Re: [9front] [patchset] libgeometry revamp

[sirjofri]

Good evening,

as a professional game developer who also specializes in 3d rendering with only experience with GPUs and shaders I have to ask: How do you handle sampling textures via coordinates? I mean, you get the barycentric coordinates, calculate the UV for that pixel, and then sample the pixel from the texture, I guess, but isn't it quite slow on a CPU to read that texture like that, especially considering caching etc.?

I mean, obviously you wouldn't draw the pixels directly (each pixel a draw call), but instead generate the image in-memory and draw it completely within one draw call.

For me it's kinda hard to learn these things since I don't know what exactly to search for (and again, most of my experience is from shaders), so some hint would be helpful, even if it's just a response like "no, this is the way to do it usually and it's fast enough".

sirjofri

Re: [9front] [patchset] libgeometry revamp

[william]

Brilliant!

Quoth rgl@antares-labs.eu:
> hello everyone,
> 
> this is a proposal to renovate the current libgeometry with a more
> modern implementation that i've been using in my personal projects for
> the last four years.  it provides almost all of the functionality
> available in the old one, with a cleaner and more consistent interface
> (similar naming conventions as those used in libdraw's geometric
> routines), plus the addition of 2D transformations (not just 3D), a
> frame of reference data type for 2- and 3-spaces, and utilities such
> as point-in-polygon testing and barycentric coordinate computations.
> 
> i documented everything in a single man page—geometry(2)—, so i got
> rid of those related to the old library, with the exception of
> qball(2), on which i'm still working and whose implementation i
> believe to be broken; i'm using a version by BurnZeZ with some minor
> changes that actually works, but it doesn't perform rotations the way
> it should (demo: https://www.youtube.com/watch?v=g3pha1EXhIk) and i
> will need more time to come up with a proper Arcball controller.
> 
> some examples of projects where i use the library:
> 
>  - http://git.antares-labs.eu/boids/
>  - http://git.antares-labs.eu/etoys/
>  - http://git.antares-labs.eu/gamephysics/
>  - http://git.antares-labs.eu/3dee/
> 
> 
> if it looks reasonable to you, i can push it upstream.
> 
> 
> 
> cheers!
> 
> -rodri
>