GEOMETRY(2) GEOMETRY(2) 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 SYNOPSIS #include <u.h> #include <libc.h> #include <geometry.h> #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 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 { double r, i, j, k; }; struct RFrame { Point2 p; Point2 bx, by; }; struct RFrame3 { Point3 p; Point3 bx, by, bz; GEOMETRY(2) GEOMETRY(2) }; struct Triangle2 { Point2 p0, p1, p2; }; struct Triangle3 { Point3 p0, p1, p2; }; /* utils */ double flerp(double a, double b, double t); double fclamp(double n, double min, double max); /* 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) /* 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); /* 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); GEOMETRY(2) GEOMETRY(2) void invm(Matrix m); Point2 xform(Point2 p, Matrix m); /* 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); /* 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 θ); /* 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); /* Triangle2 */ Point2 centroid(Triangle2 t); Point3 barycoords(Triangle2 t, Point2 p); /* Triangle3 */ Point3 centroid3(Triangle3 t); /* Fmt */ #pragma varargck type "v" Point2 #pragma varargck type "V" Point3 int vfmt(Fmt*); int Vfmt(Fmt*); void GEOMfmtinstall(void); DESCRIPTION GEOMETRY(2) GEOMETRY(2) 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. 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. UTILS These utility functions provide extra floating-point operations that are not available in the standard libc. Name Description flerp Performs a linear interpolation by a factor of t between a and b, and returns the result. fclamp Constrains n to a value between min and max, and returns the result. Points A point (x,y,w) in projective space results in the point (x/w,y/w) in Euclidean space. Vectors are represented by setting w to zero, since they don't belong to any projective plane themselves. Name Description Pt2 Constructor function for a Point2 point. Vec2 Constructor function for a Point2 vector. addpt2 Creates a new 2D point out of the sum of a's and b's components. subpt2 Creates a new 2D point out of the substraction of a's by b's components. mulpt2 Creates a new 2D point from multiplying p's components by the scalar s. divpt2 Creates a new 2D point from dividing p's components by GEOMETRY(2) GEOMETRY(2) the scalar s. lerp2 Performs a linear interpolation between the 2D points a and b by a factor of t, and returns the result. dotvec2 Computes the dot product of vectors a and b. vec2len Computes the lengthmagnitudeof vector v. normvec2 Normalizes the vector v and returns a new 2D point. edgeptcmp Performs a comparison between an edge, defined by a directed line from e0 to e1, and the point 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; otherwisewhen the point is on the line, it returns 0. ptinpoly Returns 1 if the 2D point p lies within the npts-vertex polygon defined by pts, 0 otherwise. Pt3 Constructor function for a Point3 point. Vec3 Constructor function for a Point3 vector. addpt3 Creates a new 3D point out of the sum of a's and b's components. subpt3 Creates a new 3D point out of the substraction of a's by b's components. mulpt3 Creates a new 3D point from multiplying p's components by the scalar s. divpt3 Creates a new 3D point from dividing p's components by the scalar s. lerp3 Performs a linear interpolation between the 3D points a and b by a factor of t, and returns the result. dotvec3 Computes the dot/inner product of vectors a and b. GEOMETRY(2) GEOMETRY(2) crossvec3 Computes the cross/outer product of vectors a and b. vec3len Computes the lengthmagnitudeof vector v. normvec3 Normalizes the vector v and returns a new 3D point. Matrices Name Description identity Initializes m into an identity, 3x3 matrix. addm Sums the matrices a and b and stores the result back in a. subm Substracts the matrix a by b and stores the result back in a. mulm Multiplies the matrices a and b and stores the result back in a. smulm Multiplies every element of m by the scalar s, storing the result in m. transposem Transforms the matrix m into its transpose. detm Computes the determinant of m and returns the result. tracem Computes the trace of m and returns the result. adjm Transforms the matrix m into its adjoint. invm Transforms the matrix m into its inverse. xform Transforms the point p by the matrix m and returns the new 2D point. identity3 Initializes m into an identity, 4x4 matrix. addm3 Sums the matrices a and b and stores the result back in a. subm3 GEOMETRY(2) GEOMETRY(2) Substracts the matrix a by b and stores the result back in a. mulm3 Multiplies the matrices a and b and stores the result back in a. smulm3 Multiplies every element of m by the scalar s, storing the result in m. transposem3 Transforms the matrix m into its transpose. detm3 Computes the determinant of m and returns the result. tracem3 Computes the trace of m and returns the result. adjm3 Transforms the matrix m into its adjoint. invm3 Transforms the matrix m into its inverse. xform3 Transforms the point p by the matrix m and returns the new 3D point. Quaternions Quaternions are an extension of the complex numbers con- ceived as a tool to analyze 3-dimensional points. They are most commonly used to orient and rotate objects in 3D space. Name Description Quat Constructor function for a Quaternion. Quatvec Constructor function for a Quaternion that takes the imaginary part in the form of a vector v. addq Creates a new quaternion out of the sum of a's and b's components. subq Creates a new quaternion from the substraction of a's by b's components. mulq Multiplies a and b and returns a new quaternion. smulq GEOMETRY(2) GEOMETRY(2) Multiplies each of the components of q by the scalar s, returning a new quaternion. sdivq Divides each of the components of q by the scalar s, returning a new quaternion. dotq Computes the dot-product of q and r, and returns the result. invq Computes the inverse of q and returns a new quaternion out of it. qlen Computes q's lengthmagnitudeand returns the result. normq Normalizes q and returns a new quaternion out of it. slerp Performs a spherical linear interpolation between the quaternions q and r by a factor of t, and returns the result. qrotate Returns the result of rotating the point p around the vector axis by θ radians. Frames of reference A frame of reference in a n-dimensional space is made out of n+1 points, one being the origin p, relative to some other frame of reference, and the remaining being the basis vec- tors b1,⋯,bn that define the metric within that frame. Every one of these routines assumes the origin reference frame 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 refer- ence frame if that's not the case. Name Description rframexform Transforms the point p, relative to origin O, into the frame of reference rf with origin in rf.p, which is itself also relative to O. It then returns the new 2D point. rframexform3 Transforms the point p, relative to origin O, into the frame of reference rf with origin in rf.p, which is itself also relative to O. It then returns the new 3D point. GEOMETRY(2) GEOMETRY(2) invrframexform Transforms the point p, relative to rf.p, into the frame of reference O, assumed to have an orthonormal basis. invrframexform3 Transforms the point p, relative to rf.p, into the frame of reference O, assumed to have an orthonormal basis. Triangles Name Description centroid Returns the geometric center of Triangle2 t. barycoords Returns a 3D point that represents the barycentric coordinates of the 2D point p relative to the triangle t. centroid3 Returns the geometric center of Triangle3 t. EXAMPLE The following is a common example of an RFrame being used to define the coordinate system of a rio(3) window. It places the origin at the center of the window and sets up an orthonormal basis with the y-axis pointing upwards, to con- trast with the window system where y-values grow downwards (see graphics(2)). #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); } GEOMETRY(2) GEOMETRY(2) 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); ⋯ The following snippet shows how to use the RFrame declared earlier to locate and draw a ship based on its orientation, for which we use matrix translation T and rotation R trans- formations. ⋯ 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]; GEOMETRY(2) GEOMETRY(2) 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(); ⋯ Notice how we could've used the RFrame embedded in the ship to transform the Shipmdl into the window. Instead of apply- ing the matrices to every point, the ship's local frame of reference can be rotated, effectively changing the model coordinates after an invrframexform. We are also getting rid of the θ variable, since it's no longer needed. ⋯ 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); ⋯ SOURCE /sys/src/libgeometry SEE ALSO sin(2), floor(2), graphics(2) GEOMETRY(2) GEOMETRY(2) Philip J. Schneider, David H. Eberly, Geometric Tools for Computer Graphics, Morgan Kaufmann Publishers, 2003. Jonathan Blow, Understanding Slerp, Then Not Using it, The Inner Product, April 2004. https://www.3dgep.com/understanding-quaternions/ BUGS No care is taken to avoid numeric overflows. HISTORY Libgeometry first appeared in Plan 9 from Bell Labs. It was revamped for 9front in January of 2023.