ref: a5c6374b77610cb2bcb794551475e092d990ef8b
parent: 08a080e8c2c775eda149d3e830bd4fad2c35f249
author: rodri <[email protected]>
date: Sun Jan 29 18:11:05 EST 2023
libgeometry revamp
--- 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;
+}