fixed the rframe xforms.HEADmaster

also brought in the manpage.

2 files changed, 882 insertions, 16 deletions

diff --git a/geometry.2 b/geometry.2
new file mode 100644
index 0000000..9eeb57c
--- /dev/null
+++ b/geometry.2
@@ -0,0 +1,810 @@
+.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(\fIa\fP,\fIb\fP,\fIt\fP)
+Performs a linear interpolation by a factor of
+.I t
+between
+.I a
+and
+.IR b ,
+and returns the result.
+.TP
+.B fclamp(\fIn\fP,\fImin\fP,\fImax\fP)
+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(\fIx\fP,\fIy\fP,\fIw\fP)
+Constructor function for a Point2 point.
+.TP
+.B Vec2(\fIx\fP,\fIy\fP)
+Constructor function for a Point2 vector.
+.TP
+.B addpt2(\fIa\fP,\fIb\fP)
+Creates a new 2D point out of the sum of
+.IR a 's
+and
+.IR b 's
+components.
+.TP
+.B subpt2(\fIa\fP,\fIb\fP)
+Creates a new 2D point out of the substraction of
+.IR a 's
+by
+.IR b 's
+components.
+.TP
+.B mulpt2(\fIp\fP,\fIs\fP)
+Creates a new 2D point from multiplying
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+.B divpt2(\fIp\fP,\fIs\fP)
+Creates a new 2D point from dividing
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+.B lerp2(\fIa\fP,\fIb\fP,\fIt\fP)
+Performs a linear interpolation between the 2D points
+.I a
+and
+.I b
+by a factor of
+.IR t ,
+and returns the result.
+.TP
+.B dotvec2(\fIa\fP,\fIb\fP)
+Computes the dot product of vectors
+.I a
+and
+.IR b .
+.TP
+.B vec2len(\fIv\fP)
+Computes the length—magnitude—of vector
+.IR v .
+.TP
+.B normvec2(\fIv\fP)
+Normalizes the vector
+.I v
+and returns a new 2D point.
+.TP
+.B edgeptcmp(\fIe0\fP,\fIe1\fP,\fIp\fP)
+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
+.B ptinpoly(\fIp\fP,\fIpts\fP,\fInpts\fP)
+Returns 1 if the 2D point
+.I p
+lies within the
+.IR npts -vertex
+polygon defined by
+.IR pts ,
+0 otherwise.
+.TP
+.B Pt3(\fIx\fP,\fIy\fP,\fIz\fP,\fIw\fP)
+Constructor function for a Point3 point.
+.TP
+.B Vec3(\fIx\fP,\fIy\fP,\fIz\fP)
+Constructor function for a Point3 vector.
+.TP
+.B addpt3(\fIa\fP,\fIb\fP)
+Creates a new 3D point out of the sum of
+.IR a 's
+and
+.IR b 's
+components.
+.TP
+.B subpt3(\fIa\fP,\fIb\fP)
+Creates a new 3D point out of the substraction of
+.IR a 's
+by
+.IR b 's
+components.
+.TP
+.B mulpt3(\fIp\fP,\fIs\fP)
+Creates a new 3D point from multiplying
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+.B divpt3(\fIp\fP,\fIs\fP)
+Creates a new 3D point from dividing
+.IR p 's
+components by the scalar
+.IR s .
+.TP
+.B lerp3(\fIa\fP,\fIb\fP,\fIt\fP)
+Performs a linear interpolation between the 3D points
+.I a
+and
+.I b
+by a factor of
+.IR t ,
+and returns the result.
+.TP
+.B dotvec3(\fIa\fP,\fIb\fP)
+Computes the dot product of vectors
+.I a
+and
+.IR b .
+.TP
+.B crossvec3(\fIa\fP,\fIb\fP)
+Computes the cross product of vectors
+.I a
+and
+.IR b .
+.TP
+.B vec3len(\fIv\fP)
+Computes the length—magnitude—of vector
+.IR v .
+.TP
+.B normvec3(\fIv\fP)
+Normalizes the vector
+.I v
+and returns a new 3D point.
+.SS Matrices
+.TP
+Name
+Description
+.TP
+.B identity(\fIm\fP)
+Initializes
+.I m
+into an identity, 3x3 matrix.
+.TP
+.B addm(\fIa\fP,\fIb\fP)
+Sums the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+.B subm(\fIa\fP,\fIb\fP)
+Substracts the matrix
+.I a
+by
+.I b
+and stores the result back in
+.IR a .
+.TP
+.B mulm(\fIa\fP,\fIb\fP)
+Multiplies the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+.B smulm(\fIm\fP,\fIs\fP)
+Multiplies every element of
+.I m
+by the scalar
+.IR s ,
+storing the result in m.
+.TP
+.B transposem(\fIm\fP)
+Transforms the matrix
+.I m
+into its transpose.
+.TP
+.B detm(\fIm\fP)
+Computes the determinant of
+.I m
+and returns the result.
+.TP
+.B tracem(\fIm\fP)
+Computes the trace of
+.I m
+and returns the result.
+.TP
+.B adjm(\fIm\fP)
+Transforms the matrix
+.I m
+into its adjoint.
+.TP
+.B invm(\fIm\fP)
+Transforms the matrix
+.I m
+into its inverse.
+.TP
+.B xform(\fIp\fP,\fIm\fP)
+Transforms the point
+.I p
+by the matrix
+.I m
+and returns the new 2D point.
+.TP
+.B identity3(\fIm\fP)
+Initializes
+.I m
+into an identity, 4x4 matrix.
+.TP
+.B addm3(\fIa\fP,\fIb\fP)
+Sums the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+.B subm3(\fIa\fP,\fIb\fP)
+Substracts the matrix
+.I a
+by
+.I b
+and stores the result back in
+.IR a .
+.TP
+.B mulm3(\fIa\fP,\fIb\fP)
+Multiplies the matrices
+.I a
+and
+.I b
+and stores the result back in
+.IR a .
+.TP
+.B smulm3(\fIm\fP,\fIs\fP)
+Multiplies every element of
+.I m
+by the scalar
+.IR s ,
+storing the result in m.
+.TP
+.B transposem3(\fIm\fP)
+Transforms the matrix
+.I m
+into its transpose.
+.TP
+.B detm3(\fIm\fP)
+Computes the determinant of
+.I m
+and returns the result.
+.TP
+.B tracem3(\fIm\fP)
+Computes the trace of
+.I m
+and returns the result.
+.TP
+.B adjm3(\fIm\fP)
+Transforms the matrix
+.I m
+into its adjoint.
+.TP
+.B invm3(\fIm\fP)
+Transforms the matrix
+.I m
+into its inverse.
+.TP
+.B xform3(\fIp\fP,\fIm\fP)
+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
+.B Quat(\fIr\fP,\fIi\fP,\fIj\fP,\fIk\fP)
+Constructor function for a Quaternion.
+.TP
+.B Quatvec(\fIr\fP,\fIv\fP)
+Constructor function for a Quaternion that takes the imaginary part in
+the form of a vector
+.IR v .
+.TP
+.B addq(\fIa\fP,\fIb\fP)
+Creates a new quaternion out of the sum of
+.IR a 's
+and
+.IR b 's
+components.
+.TP
+.B subq(\fIa\fP,\fIb\fP)
+Creates a new quaternion from the substraction of
+.IR a 's
+by
+.IR b 's
+components.
+.TP
+.B mulq(\fIa\fP,\fIb\fP)
+Multiplies
+.I a
+and
+.I b
+and returns a new quaternion.
+.TP
+.B smulq(\fIq\fP,\fIs\fP)
+Multiplies each of the components of
+.I q
+by the scalar
+.IR s ,
+returning a new quaternion.
+.TP
+.B sdivq(\fIq\fP,\fIs\fP)
+Divides each of the components of
+.I q
+by the scalar
+.IR s ,
+returning a new quaternion.
+.TP
+.B dotq(\fIq\fP,\fIr\fP)
+Computes the dot-product of
+.I q
+and
+.IR r ,
+and returns the result.
+.TP
+.B invq(\fIq\fP)
+Computes the inverse of
+.I q
+and returns a new quaternion out of it.
+.TP
+.B qlen(\fIq\fP)
+Computes
+.IR q 's
+length—magnitude—and returns the result.
+.TP
+.B normq(\fIq\fP)
+Normalizes
+.I q
+and returns a new quaternion out of it.
+.TP
+.B slerp(\fIq\fP,\fIr\fP,\fIt\fP)
+Performs a spherical linear interpolation between the quaternions
+.I q
+and
+.I r
+by a factor of
+.IR t ,
+and returns the result.
+.TP
+.B qrotate(\fIp\fP,\fIaxis\fP,\fIθ\fP)
+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 ,
+and the remaining being the
+basis vectors
+.I b1,⋯,bn
+that define the metric within that frame.
+.PP
+The origin point and the bases are all defined in terms of an origin
+frame of reference O. Applying a forward transformation
+.RI ( rframexform
+and
+.IR rframexform3 )
+to a point relative to O will result in a point relative to the new frame.
+Applying an inverse transformation
+.RI ( invrframexform
+and
+.IR invrframexform3 )
+to that same point—now defined in terms of the new frame—will bring it back to O.
+.TP
+Name
+Description
+.TP
+.B rframexform(\fIp\fP,\fIrf\fP)
+Transforms the point
+.IR p ,
+relative to some origin frame of reference O, into the frame of reference
+.IR rf .
+It then returns the new 2D point.
+.TP
+.B rframexform3(\fIp\fP,\fIrf\fP)
+Transforms the point
+.IR p ,
+relative to some origin frame of reference O, into the frame of reference
+.IR rf .
+It then returns the new 3D point.
+.TP
+.B invrframexform(\fIp\fP,\fIrf\fP)
+Transforms the point
+.IR p ,
+relative to
+.IR rf ,
+into a point relative to the origin frame of reference O.
+It then returns the new 2D point.
+.TP
+.B invrframexform3(\fIp\fP,\fIrf\fP)
+Transforms the point
+.IR p ,
+relative to
+.IR rf ,
+into a point relative to the origin frame of reference O.
+It then returns the new 3D point.
+.SS Triangles
+.TP
+Name
+Description
+.TP
+.B centroid(\fIt\fP)
+Returns the geometric center of
+.B Triangle2
+.IR t .
+.TP
+.B barycoords(\fIt\fP,\fIp\fP)
+Returns a 3D point that represents the barycentric coordinates of the
+2D point
+.I p
+relative to the triangle
+.IR t .
+.TP
+.B centroid3(\fIt\fP)
+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 worldrf;
+
+/* from screen... */
+Point2
+toworld(Point p)
+{
+	return rframexform(p, worldrf);
+}
+
+/* ...to screen */
+Point
+fromworld(Point2 p)
+{
+	p = invrframexform(Pt2(p.x,p.y,1), worldrf);
+	return Pt(p.x,p.y);
+}
+
+void
+main(void)
+	⋯
+	worldrf.p = Pt2(screen->r.min.x+Dx(screen->r)/2,screen->r.max.y-Dy(screen->r)/2,1);
+	worldrf.bx = Vec2(1, 0);
+	worldrf.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] = fromworld(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] = fromworld(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/rframe.c b/rframe.c
index 2f7c682..aee3c03 100644
--- a/rframe.c
+++ b/rframe.c
@@ -2,50 +2,106 @@
 #include <libc.h>
 #include <geometry.h>
 
+/*
+ * implicit identity origin rframes
+ *
+ * 	static RFrame IRF2 = {
+ * 		.p  = {0,0,1},
+ * 		.bx = {1,0,0},
+ * 		.by = {0,1,0},
+ * 	};
+ * 	
+ * 	static RFrame3 IRF3 = {
+ * 		.p  = {0,0,0,1},
+ * 		.bx = {1,0,0,0},
+ * 		.by = {0,1,0,0},
+ * 		.bz = {0,0,1,0},
+ * 	};
+ *
+ * these rframes are used on every xform to keep the points in the
+ * correct plane (i.e. with proper w values); they are written here as a
+ * reference for future changes. the bases are ignored since they turn
+ * into an unnecessary identity xform.
+ *
+ * the implicitness comes from the fact that using the _irf* filters
+ * makes the rframexform equivalent to:
+ * 	rframexform(invrframexform(p, IRF), rf);
+ * and the invrframexform to:
+ * 	rframexform(invrframexform(p, rf), IRF);
+ */
+
+static Point2
+_irfxform(Point2 p)
+{
+	p.w--;
+	return p;
+}
+
+static Point2
+_irfxform⁻¹(Point2 p)
+{
+	p.w++;
+	return p;
+}
+
+static Point3
+_irfxform3(Point3 p)
+{
+	p.w--;
+	return p;
+}
+
+static Point3
+_irfxform3⁻¹(Point3 p)
+{
+	p.w++;
+	return p;
+}
+
 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),
+		rf.bx.x, rf.by.x, 0,
+		rf.bx.y, rf.by.y, 0,
 		0, 0, 1,
 	};
-	return xform(p, m);
+	invm(m);
+	return xform(subpt2(_irfxform⁻¹(p), rf.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),
+		rf.bx.x, rf.by.x, rf.bz.x, 0,
+		rf.bx.y, rf.by.y, rf.bz.y, 0,
+		rf.bx.z, rf.by.z, rf.bz.z, 0,
 		0, 0, 0, 1,
 	};
-	return xform3(p, m);
+	invm3(m);
+	return xform3(subpt3(_irfxform3⁻¹(p), rf.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),
+		rf.bx.x, rf.by.x, 0,
+		rf.bx.y, rf.by.y, 0,
 		0, 0, 1,
 	};
-	invm(m);
-	return xform(p, m);
+	return _irfxform(addpt2(xform(p, m), rf.p));
 }
 
 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),
+		rf.bx.x, rf.by.x, rf.bz.x, 0,
+		rf.bx.y, rf.by.y, rf.bz.y, 0,
+		rf.bx.z, rf.by.z, rf.bz.z, 0,
 		0, 0, 0, 1,
 	};
-	invm3(m);
-	return xform3(p, m);
+	return _irfxform3(addpt3(xform3(p, m), rf.p));
 }