plan9front/sys/man/2/geometry

.TH GEOMETRY 2
.SH NAME
Flerp, fberp, fclamp, Pt2, Vec2, addpt2, subpt2, mulpt2, divpt2, lerp2, berp2, dotvec2, vec2len, normvec2, edgeptcmp, ptinpoly, Pt3, Vec3, addpt3, subpt3, mulpt3, divpt3, lerp3, berp3, dotvec3, crossvec3, vec3len, normvec3, identity, addm, subm, mulm, smulm, transposem, detm, tracem, minorm, cofactorm, adjm, invm, xform, identity3, addm3, subm3, mulm3, smulm3, transposem3, detm3, tracem3, minorm3, cofactorm3, adjm3, invm3, xform3, Quat, Quatvec, addq, subq, mulq, smulq, sdivq, dotq, invq, qlen, normq, slerp, qrotate, rframematrix, rframematrix3, 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 fberp(double a, double b, double c, Point3 bc)
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)
Point2 berp2(Point2 a, Point2 b, Point2 c, Point3 bc)
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)
Point3 berp3(Point3 a, Point3 b, Point3 c, Point3 bc)
double dotvec3(Point3 a, Point3 b)
Point3 crossvec3(Point3 a, Point3 b)
double vec3len(Point3 v)
Point3 normvec3(Point3 v)
int lineXsphere(Point3 *rp, Point3 p0, Point3 p1,
		Point3 c, double rad, int isaray)
int ptincylinder(Point3 p, Point3 p0, Point3 p1, double r)
int ptincone(Point3 p, Point3 p0, Point3 p1, double br)
.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)
double minorm(Matrix m, int row, int col)
double cofactorm(Matrix m, int row, int col)
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)
double minorm3(Matrix3 m, int row, int col)
double cofactorm3(Matrix3 m, int row, int col)
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)
Quaternion qsandwich(Quaternion q, Quaternion r)
Point3 qsandwichpt3(Quaternion q, Point3 p)
Point3 qrotate(Point3 p, Point3 axis, double θ)
.PB
/* RFrame */
void rframematrix(Matrix, RFrame)
void rframematrix3(Matrix3, RFrame3)
Point2 rframexform(Point2 p, RFrame rf)
Point3 rframexform3(Point3 p, RFrame3 rf)
Point2 invrframexform(Point2 p, RFrame rf)
Point3 invrframexform3(Point3 p, RFrame3 rf)
.PB
/* Triangle2 */
Point2 centroid(Triangle2 t)
Point3 barycoords(Triangle2 t, Point2 p)
.PB
/* Triangle3 */
Point3 centroid3(Triangle3 t)
.PB
/* Fmt */
#pragma varargck type "v" Point2
#pragma varargck type "V" Point3
int vfmt(Fmt*)
int Vfmt(Fmt*)
void GEOMfmtinstall(void)
.SH DESCRIPTION
This library provides routines to operate with homogeneous coordinates
in 2D and 3D projective spaces by means of points, matrices,
quaternions and frames of reference.
.PP
Besides their many mathematical properties and applications, the data
structures and algorithms used here to represent these abstractions
are specifically tailored to the world of computer graphics and
simulators, and so it uses the conventions associated with these
fields, such as the right-hand rule for coordinate systems and column
vectors for matrix operations.
.SS UTILS
These utility functions provide extra floating-point operations that
are not available in the standard libc.
.TP
Name
Description
.TP
.B flerp(\fIa\fP,\fIb\fP,\fIt\fP)
Performs a linear interpolation by a factor of
.I t
between
.I a
and
.IR b ,
and returns the result.
.TP
.B fberp(\fIa\fP,\fIb\fP,\fIc\fP,\fIbc\fP)
Performs a barycentric interpolation of the values in
.IR a ,
.I b
and
.IR c ,
based on the weights provided by
.IR bc .
.TP
.B fclamp(\fIn\fP,\fImin\fP,\fImax\fP)
Constrains
.I n
to a value between
.I min
and
.IR max ,
and returns the result.
.SS Points
A point
.B (x,y,w)
in projective space results in the point
.B (x/w,y/w)
in Euclidean space. Vectors are represented by setting
.B w
to zero, since they don't belong to any projective plane themselves.
.TP
Name
Description
.TP
.B Pt2(\fIx\fP,\fIy\fP,\fIw\fP)
Constructor function for a Point2 point.
.TP
.B Vec2(\fIx\fP,\fIy\fP)
Constructor function for a Point2 vector.
.TP
.B addpt2(\fIa\fP,\fIb\fP)
Creates a new 2D point out of the sum of
.IR a 's
and
.IR b 's
components.
.TP
.B subpt2(\fIa\fP,\fIb\fP)
Creates a new 2D point out of the substraction of
.IR a 's
by
.IR b 's
components.
.TP
.B mulpt2(\fIp\fP,\fIs\fP)
Creates a new 2D point from multiplying
.IR p 's
components by the scalar
.IR s .
.TP
.B divpt2(\fIp\fP,\fIs\fP)
Creates a new 2D point from dividing
.IR p 's
components by the scalar
.IR s .
.TP
.B lerp2(\fIa\fP,\fIb\fP,\fIt\fP)
Performs a linear interpolation between the 2D points
.I a
and
.I b
by a factor of
.IR t ,
and returns the result.
.TP
.B berp2(\fIa\fP,\fIb\fP,\fIc\fP,\fIbc\fP)
Performs a barycentric interpolation between the 2D points
.IR a ,
.I b
and
.IR c ,
using the weights in
.IR bc .
.TP
.B dotvec2(\fIa\fP,\fIb\fP)
Computes the dot product of vectors
.I a
and
.IR b .
.TP
.B vec2len(\fIv\fP)
Computes the length—magnitude—of vector
.IR v .
.TP
.B normvec2(\fIv\fP)
Normalizes the vector
.I v
and returns a new 2D point.
.TP
.B edgeptcmp(\fIe0\fP,\fIe1\fP,\fIp\fP)
Performs a comparison between an edge, defined by a directed line from
.I e0
to
.IR e1 ,
and the point
.IR p .
If the point is to the right of the line, the result is >0; if it's to
the left, the result is <0; otherwise—when the point is on the line—,
it returns 0.
.TP
.B ptinpoly(\fIp\fP,\fIpts\fP,\fInpts\fP)
Returns 1 if the 2D point
.I p
lies within the
.IR npts -vertex
polygon defined by
.IR pts ,
0 otherwise.
.TP
.B Pt3(\fIx\fP,\fIy\fP,\fIz\fP,\fIw\fP)
Constructor function for a Point3 point.
.TP
.B Vec3(\fIx\fP,\fIy\fP,\fIz\fP)
Constructor function for a Point3 vector.
.TP
.B addpt3(\fIa\fP,\fIb\fP)
Creates a new 3D point out of the sum of
.IR a 's
and
.IR b 's
components.
.TP
.B subpt3(\fIa\fP,\fIb\fP)
Creates a new 3D point out of the substraction of
.IR a 's
by
.IR b 's
components.
.TP
.B mulpt3(\fIp\fP,\fIs\fP)
Creates a new 3D point from multiplying
.IR p 's
components by the scalar
.IR s .
.TP
.B divpt3(\fIp\fP,\fIs\fP)
Creates a new 3D point from dividing
.IR p 's
components by the scalar
.IR s .
.TP
.B lerp3(\fIa\fP,\fIb\fP,\fIt\fP)
Performs a linear interpolation between the 3D points
.I a
and
.I b
by a factor of
.IR t ,
and returns the result.
.TP
.B berp3(\fIa\fP,\fIb\fP,\fIc\fP,\fIbc\fP)
Performs a barycentric interpolation between the 3D points
.IR a ,
.I b
and
.IR c ,
using the weights in
.IR bc .
.TP
.B dotvec3(\fIa\fP,\fIb\fP)
Computes the dot product of vectors
.I a
and
.IR b .
.TP
.B crossvec3(\fIa\fP,\fIb\fP)
Computes the cross product of vectors
.I a
and
.IR b .
.TP
.B vec3len(\fIv\fP)
Computes the length—magnitude—of vector
.IR v .
.TP
.B normvec3(\fIv\fP)
Normalizes the vector
.I v
and returns a new 3D point.
.TP
.B lineXsphere(\fIrp\fP,\fIp0\fP,\fIp1\fP,\fIc\fP,\fIrad\fP,\fIisaray\fP)
Finds the intersection between the line defined by
.I p0
and
.I p1
and the sphere with center at
.I c
and a radius of
.IR rad .
If
.I isaray
is not zero, the line will be interpreted as a ray directed from
.I p0
to
.IR p1 .
The function returns the number of intersections (up to two) and, if
.I rp
is not nil, it is filled with the resulting points.
.TP
.B ptincylinder(\fIp\fP,\fIp0\fP,\fIp1\fP,\fIr\fP)
Tests if the point
.I p
is inside a cylinder with an axis defined by the line between
.I p0
and
.IR p1 ,
and a radius r.
.TP
.B ptincone(\fIp\fP,\fIp0\fP,\fIp1\fP,\fIbr\fP)
Tests if the point
.I p
is inside a cone with its apex at
.I p0
and its base centered at
.IR p1 ,
with a radius of br.
.SS Matrices
.TP
Name
Description
.TP
.B identity(\fIm\fP)
Initializes
.I m
into an identity, 3x3 matrix.
.TP
.B addm(\fIa\fP,\fIb\fP)
Sums the matrices
.I a
and
.I b
and stores the result back in
.IR a .
.TP
.B subm(\fIa\fP,\fIb\fP)
Substracts the matrix
.I a
by
.I b
and stores the result back in
.IR a .
.TP
.B mulm(\fIa\fP,\fIb\fP)
Multiplies the matrices
.I a
and
.I b
and stores the result back in
.IR a .
.TP
.B smulm(\fIm\fP,\fIs\fP)
Multiplies every element of
.I m
by the scalar
.IR s ,
storing the result in m.
.TP
.B transposem(\fIm\fP)
Transforms the matrix
.I m
into its transpose.
.TP
.B detm(\fIm\fP)
Computes the determinant of
.I m
and returns the result.
.TP
.B tracem(\fIm\fP)
Computes the trace of
.I m
and returns the result.
.TP
.B minorm(\fIm\fP,\fIrow\fP,\fIcol\fP)
Computes the minor of
.I m
and returns the result.
.TP
.B cofactorm(\fIm\fP,\fIrow\fP,\fIcol\fP)
Computes the cofactor of
.I m
and returns the result.
.TP
.B adjm(\fIm\fP)
Transforms the matrix
.I m
into its adjoint.
.TP
.B invm(\fIm\fP)
Transforms the matrix
.I m
into its inverse.
.TP
.B xform(\fIp\fP,\fIm\fP)
Transforms the point
.I p
by the matrix
.I m
and returns the new 2D point.
.TP
.B identity3(\fIm\fP)
Initializes
.I m
into an identity, 4x4 matrix.
.TP
.B addm3(\fIa\fP,\fIb\fP)
Sums the matrices
.I a
and
.I b
and stores the result back in
.IR a .
.TP
.B subm3(\fIa\fP,\fIb\fP)
Substracts the matrix
.I a
by
.I b
and stores the result back in
.IR a .
.TP
.B mulm3(\fIa\fP,\fIb\fP)
Multiplies the matrices
.I a
and
.I b
and stores the result back in
.IR a .
.TP
.B smulm3(\fIm\fP,\fIs\fP)
Multiplies every element of
.I m
by the scalar
.IR s ,
storing the result in m.
.TP
.B transposem3(\fIm\fP)
Transforms the matrix
.I m
into its transpose.
.TP
.B detm3(\fIm\fP)
Computes the determinant of
.I m
and returns the result.
.TP
.B tracem3(\fIm\fP)
Computes the trace of
.I m
and returns the result.
.TP
.B minorm3(\fIm\fP,\fIrow\fP,\fIcol\fP)
Computes the minor of
.I m
and returns the result.
.TP
.B cofactorm3(\fIm\fP,\fIrow\fP,\fIcol\fP)
Computes the cofactor of
.I m
and returns the result.
.TP
.B adjm3(\fIm\fP)
Transforms the matrix
.I m
into its adjoint.
.TP
.B invm3(\fIm\fP)
Transforms the matrix
.I m
into its inverse.
.TP
.B xform3(\fIp\fP,\fIm\fP)
Transforms the point
.I p
by the matrix
.I m
and returns the new 3D point.
.SS Quaternions
Quaternions are an extension of the complex numbers conceived as a
tool to analyze 3-dimensional points.  They are most commonly used to
orient and rotate objects in 3D space.
.TP
Name
Description
.TP
.B Quat(\fIr\fP,\fIi\fP,\fIj\fP,\fIk\fP)
Constructor function for a Quaternion.
.TP
.B Quatvec(\fIr\fP,\fIv\fP)
Constructor function for a Quaternion that takes the imaginary part in
the form of a vector
.IR v .
.TP
.B addq(\fIa\fP,\fIb\fP)
Creates a new quaternion out of the sum of
.IR a 's
and
.IR b 's
components.
.TP
.B subq(\fIa\fP,\fIb\fP)
Creates a new quaternion from the substraction of
.IR a 's
by
.IR b 's
components.
.TP
.B mulq(\fIa\fP,\fIb\fP)
Multiplies
.I a
and
.I b
and returns a new quaternion.
.TP
.B smulq(\fIq\fP,\fIs\fP)
Multiplies each of the components of
.I q
by the scalar
.IR s ,
returning a new quaternion.
.TP
.B sdivq(\fIq\fP,\fIs\fP)
Divides each of the components of
.I q
by the scalar
.IR s ,
returning a new quaternion.
.TP
.B dotq(\fIq\fP,\fIr\fP)
Computes the dot-product of
.I q
and
.IR r ,
and returns the result.
.TP
.B invq(\fIq\fP)
Computes the inverse of
.I q
and returns a new quaternion out of it.
.TP
.B qlen(\fIq\fP)
Computes
.IR q 's
length—magnitude—and returns the result.
.TP
.B normq(\fIq\fP)
Normalizes
.I q
and returns a new quaternion out of it.
.TP
.B slerp(\fIq\fP,\fIr\fP,\fIt\fP)
Performs a spherical linear interpolation between the quaternions
.I q
and
.I r
by a factor of
.IR t ,
and returns the result.
.TP
.B qsandwich(\fIq\fP,\fIr\fP)
Returns the result of "sandwiching"
.I q
with
.IR r ,
i.e. r' = qrq⁻¹.
.TP
.B qsandwichpt3(\fIq\fP,\fIp\fP)
Same as
.IR qsandwich (2),
but treating
.I p
as a pure quaternion when computing the product, and keeping its
projective coordinate (w) intact.
.TP
.B qrotate(\fIp\fP,\fIaxis\fP,\fIθ\fP)
Returns the result of rotating the point
.I p
around the vector
.I axis
by
.I θ
radians.
.SS Frames of reference
A frame of reference (or rframe) in a
.IR n -dimensional
space is made out of n+1 points, one being the origin
.IR p ,
and the remaining being the
basis vectors
.I b1,⋯,bn
that define the metric within that frame.
.TP
Name
Description
.TP
.B rframematrix(\fIm\fP,\fIrf\fP)
Initializes
.I m
as an rframe transformation matrix based on
.IR rf .
.TP
.B rframematrix3(\fIm\fP,\fIrf\fP)
Initializes
.I m
as an rframe transformation matrix based on
.IR rf .
.TP
.B rframexform(\fIp\fP,\fIrf\fP)
Transforms the point
.IR p ,
relative to some origin frame of reference O, into the rframe
.IR rf .
It then returns the new 2D point.
.TP
.B rframexform3(\fIp\fP,\fIrf\fP)
Transforms the point
.IR p ,
relative to some origin frame of reference O, into the rframe
.IR rf .
It then returns the new 3D point.
.TP
.B invrframexform(\fIp\fP,\fIrf\fP)
Transforms the point
.IR p ,
relative to
.IR rf ,
into a point relative to the origin rframe O.
It then returns the new 2D point.
.TP
.B invrframexform3(\fIp\fP,\fIrf\fP)
Transforms the point
.IR p ,
relative to
.IR rf ,
into a point relative to the origin rframe O.
It then returns the new 3D point.
.SS Triangles
.TP
Name
Description
.TP
.B centroid(\fIt\fP)
Returns the geometric center of
.B Triangle2
.IR t .
.TP
.B barycoords(\fIt\fP,\fIp\fP)
Returns a 3D point that represents the barycentric coordinates of the
2D point
.I p
relative to the triangle
.IR t .
.TP
.B centroid3(\fIt\fP)
Returns the geometric center of
.B Triangle3
.IR t .
.SH EXAMPLES
.PP
Rotate a point p by θ, scale it 2x, then translate it by vector v:
.IP
.EX
Matrix R = {
	cos(θ), -sin(θ), 0,
	sin(θ),  cos(θ), 0,
	0, 0, 1,
}, S = {
	2, 0, 0,
	0, 2, 0,
	0, 0, 1,
}, T = {
	1, 0, v.x,
	0, 1, v.y,
	0, 0, 1,
};

mulm(S, R);
mulm(T, S);
p = xform(p, T); /* p' = T·S·R·p */
.EE

.PP
Given a space with two observers, A and B, and a point p relative to
A, find its location relative to B:
.IP
.EX
pb = rframexform(invrframexform(p, A), B);
.EE
.PP
Now let's say observer C is located relative to A; find the point's
location in C:
.IP
.EX
pc = rframexform(p, C);
.EE
.PP
Finally, to obtain its location relative to the space itself (its
global position):
.IP
.EX
pg = invrframexform(p, A);
.EE

.PP
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)).
.IP
.EX
RFrame worldrf;

/* from screen... */
Point2
toworld(Point p)
{
	return rframexform(Pt2(p.x,p.y,1), worldrf);
}

/* ...to screen */
Point
fromworld(Point2 p)
{
	p = invrframexform(p, worldrf);
	return Pt(p.x,p.y);
}

void
main(void)
	⋯
	worldrf.p = Pt2(screen->r.min.x+Dx(screen->r)/2,
			screen->r.min.y+Dy(screen->r)/2, 1);
	worldrf.bx = Vec2(1, 0);
	worldrf.by = Vec2(0,-1);
	⋯
.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/
.br
https://motion.cs.illinois.edu/RoboticSystems/CoordinateTransformations.html
.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.