diff 8188b4f4f0e07b6669e6ae3c6c1099af917eaab4 uncommitted --- a/sys/include/geometry.h +++ b/sys/include/geometry.h @@ -1,89 +1,137 @@ #pragma lib "libgeometry.a" #pragma src "/sys/src/libgeometry" -typedef double Matrix[4][4]; + +#define DEG 0.01745329251994330 /* π/180 */ + +typedef struct Point2 Point2; typedef struct Point3 Point3; +typedef double Matrix[3][3]; +typedef double Matrix3[4][4]; typedef struct Quaternion Quaternion; -typedef struct Space Space; -struct Point3{ +typedef struct RFrame RFrame; +typedef struct RFrame3 RFrame3; +typedef struct Triangle2 Triangle2; +typedef struct Triangle3 Triangle3; + +struct Point2 { + double x, y, w; +}; + +struct Point3 { double x, y, z, w; }; -struct Quaternion{ + +struct Quaternion { double r, i, j, k; }; -struct Space{ - Matrix t; - Matrix tinv; - Space *next; + +struct RFrame { + Point2 p; + Point2 bx, by; }; -/* - * 3-d point arithmetic - */ -Point3 add3(Point3 a, Point3 b); -Point3 sub3(Point3 a, Point3 b); -Point3 neg3(Point3 a); -Point3 div3(Point3 a, double b); -Point3 mul3(Point3 a, double b); -int eqpt3(Point3 p, Point3 q); -int closept3(Point3 p, Point3 q, double eps); -double dot3(Point3 p, Point3 q); -Point3 cross3(Point3 p, Point3 q); -double len3(Point3 p); -double dist3(Point3 p, Point3 q); -Point3 unit3(Point3 p); -Point3 midpt3(Point3 p, Point3 q); -Point3 lerp3(Point3 p, Point3 q, double alpha); -Point3 reflect3(Point3 p, Point3 p0, Point3 p1); -Point3 nearseg3(Point3 p0, Point3 p1, Point3 testp); -double pldist3(Point3 p, Point3 p0, Point3 p1); -double vdiv3(Point3 a, Point3 b); -Point3 vrem3(Point3 a, Point3 b); -Point3 pn2f3(Point3 p, Point3 n); -Point3 ppp2f3(Point3 p0, Point3 p1, Point3 p2); -Point3 fff2p3(Point3 f0, Point3 f1, Point3 f2); -Point3 pdiv4(Point3 a); -Point3 add4(Point3 a, Point3 b); -Point3 sub4(Point3 a, Point3 b); -/* - * Quaternion arithmetic - */ -void qtom(Matrix, Quaternion); -Quaternion mtoq(Matrix); -Quaternion qadd(Quaternion, Quaternion); -Quaternion qsub(Quaternion, Quaternion); -Quaternion qneg(Quaternion); -Quaternion qmul(Quaternion, Quaternion); -Quaternion qdiv(Quaternion, Quaternion); -Quaternion qunit(Quaternion); -Quaternion qinv(Quaternion); + +struct RFrame3 { + Point3 p; + Point3 bx, by, bz; +}; + +struct Triangle2 +{ + Point2 p0, p1, p2; +}; + +struct Triangle3 { + Point3 p0, p1, p2; +}; + +/* utils */ +double flerp(double, double, double); +double fclamp(double, double, double); + +/* Point2 */ +Point2 Pt2(double, double, double); +Point2 Vec2(double, double); +Point2 addpt2(Point2, Point2); +Point2 subpt2(Point2, Point2); +Point2 mulpt2(Point2, double); +Point2 divpt2(Point2, double); +Point2 lerp2(Point2, Point2, double); +double dotvec2(Point2, Point2); +double vec2len(Point2); +Point2 normvec2(Point2); +int edgeptcmp(Point2, Point2, Point2); +int ptinpoly(Point2, Point2*, ulong); + +/* Point3 */ +Point3 Pt3(double, double, double, double); +Point3 Vec3(double, double, double); +Point3 addpt3(Point3, Point3); +Point3 subpt3(Point3, Point3); +Point3 mulpt3(Point3, double); +Point3 divpt3(Point3, double); +Point3 lerp3(Point3, Point3, double); +double dotvec3(Point3, Point3); +Point3 crossvec3(Point3, Point3); +double vec3len(Point3); +Point3 normvec3(Point3); + +/* Matrix */ +void identity(Matrix); +void addm(Matrix, Matrix); +void subm(Matrix, Matrix); +void mulm(Matrix, Matrix); +void smulm(Matrix, double); +void transposem(Matrix); +double detm(Matrix); +double tracem(Matrix); +void adjm(Matrix); +void invm(Matrix); +Point2 xform(Point2, Matrix); + +/* Matrix3 */ +void identity3(Matrix3); +void addm3(Matrix3, Matrix3); +void subm3(Matrix3, Matrix3); +void mulm3(Matrix3, Matrix3); +void smulm3(Matrix3, double); +void transposem3(Matrix3); +double detm3(Matrix3); +double tracem3(Matrix3); +void adjm3(Matrix3); +void invm3(Matrix3); +Point3 xform3(Point3, Matrix3); + +/* Quaternion */ +Quaternion Quat(double, double, double, double); +Quaternion Quatvec(double, Point3); +Quaternion addq(Quaternion, Quaternion); +Quaternion subq(Quaternion, Quaternion); +Quaternion mulq(Quaternion, Quaternion); +Quaternion smulq(Quaternion, double); +Quaternion sdivq(Quaternion, double); +double dotq(Quaternion, Quaternion); +Quaternion invq(Quaternion); double qlen(Quaternion); +Quaternion normq(Quaternion); Quaternion slerp(Quaternion, Quaternion, double); -Quaternion qmid(Quaternion, Quaternion); -Quaternion qsqrt(Quaternion); -void qball(Rectangle, Mouse *, Quaternion *, void (*)(void), Quaternion *); -/* - * Matrix arithmetic - */ -void ident(Matrix); -void matmul(Matrix, Matrix); -void matmulr(Matrix, Matrix); -double determinant(Matrix); -void adjoint(Matrix, Matrix); -double invertmat(Matrix, Matrix); -/* - * Space stack routines - */ -Space *pushmat(Space *); -Space *popmat(Space *); -void rot(Space *, double, int); -void qrot(Space *, Quaternion); -void scale(Space *, double, double, double); -void move(Space *, double, double, double); -void xform(Space *, Matrix); -void ixform(Space *, Matrix, Matrix); -void look(Space *, Point3, Point3, Point3); -int persp(Space *, double, double, double); -void viewport(Space *, Rectangle, double); -Point3 xformpoint(Point3, Space *, Space *); -Point3 xformpointd(Point3, Space *, Space *); -Point3 xformplane(Point3, Space *, Space *); -#define radians(d) ((d)*.01745329251994329572) +Point3 qrotate(Point3, Point3, double); + +/* RFrame */ +Point2 rframexform(Point2, RFrame); +Point3 rframexform3(Point3, RFrame3); +Point2 invrframexform(Point2, RFrame); +Point3 invrframexform3(Point3, RFrame3); + +/* Triangle2 */ +Point2 centroid(Triangle2); +Point3 barycoords(Triangle2, Point2); + +/* Triangle3 */ +Point3 centroid3(Triangle3); + +/* Fmt */ +#pragma varargck type "v" Point2 +#pragma varargck type "V" Point3 +int vfmt(Fmt*); +int Vfmt(Fmt*); +void GEOMfmtinstall(void); --- a/sys/man/2/arith3 +++ /dev/null @@ -1,268 +1,0 @@ -.TH ARITH3 2 -.SH NAME -add3, sub3, neg3, div3, mul3, eqpt3, closept3, dot3, cross3, len3, dist3, unit3, midpt3, lerp3, reflect3, nearseg3, pldist3, vdiv3, vrem3, pn2f3, ppp2f3, fff2p3, pdiv4, add4, sub4 \- operations on 3-d points and planes -.SH SYNOPSIS -.B -#include -.br -.B -#include -.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 +#include +#include +.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 +#include +#include +#include + +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 -.PP -.B -#include -.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 -.br -.B -#include -.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 -#include -#include -#include -/* - * 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)=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 +#include +#include + +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 #include -#include #include -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 @@ +#include +#include + +/* 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 -#include -#include -#include -#include -#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.xxy)); - 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 #include -#include #include -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 +#include +#include + +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 -#include -#include -#include -/* - * 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 +#include +#include + +/* 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 -#include -#include -#include -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 +#include +#include + +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; +}