In "Clean up your Mesh" we propose the use of PGA as a geometric representation for $k$-simplices and $k$-complexes,
and demonstrate how this leads to elegant formulas for $k$-magnitude, center of mass and moments of inertia.
On this page, we provide a brief overview and example implementations.
In the above interactive example, the volume under the plane is calculated without the explicit construction of caps, using only a single scalar for runtime storage. The source code is available here
After embedding our vertices as dual vectors, joining $k+1$ points produces a $k$-simplex carrier $S_k$. $$ v_i = (\mathbf e_0 + \vec v)^*,\qquad S_k = v_0 \vee \dots \vee v_k,\qquad C_k = \sum S_k $$ While the $S_k$ does not contain the vertex positions, it represents the carrier line, plane, etc as well as the length, area, volume, etc of the segment, triangle, tetra, ... respectively. Similarly, while summing the individual $S_k$ into a complex carrier $C_k$ retains no information of the original simplices, the resulting $C_k$ does still encode the magnitude of the complex.
Boundary and Gap
A $k$-complex (e.g. a 3D mesh) $M$ can be represented by the $k$-simplices in its interior (a set of tetrahedrons),
or by the $k-1$-simplices on its boundary $\partial M$ (a set of triangles).
Such a boundary $\partial M$ is said to be closed when each $k-2$-simplex it contains is connected to exactly two $k-1$ simplices.
(each point on the boundary of a polygon belongs to two edges, or each edge of a 3D mesh belongs to two triangles, etc).
When the boundary is not closed, the missing $k-1$ simplices are called the gap and written $\xcancel \partial M$.
The $k$-magnitude (amount, length, area, volume, ...) of any $k$-simplex is now given by $$ {1 \over k!} \lVert S_k \rVert = {1 \over k!} \lVert v_0 \vee \dots \vee v_k \rVert $$ The $k$-magnitude of a $k$-simplex can also be calculated from its boundary $k-1$ complex $$ {1 \over k!} \lVert C_{\partial k-1} \rVert_\infty $$
The $k-1$-magnitude of the sum of the missing $k-1$ simplices from a boundary complex is given by $${1 \over (k-1)!} \lVert C_{\xcancel \partial k-1} \rVert = {1 \over (k-1)!} \lVert C_{\partial k-1} \rVert$$After calculating the 0th (volume) and 1st (com) moments, it is natural to wonder if the calculation of the 2nd moments of inertia also can be peformed within the PGA framework. The example below (available here) demonstrates the method from the paper for the known case of an arbitrary box.
PGA expressions are not more expensive than their vector or coordinate equivalents, in fact, they compile to them. A modern GA library indeed allows you to use geometrically meaningful, dimension agnostic and expressive formulas without imposing an extra computational cost. Let us examine how that works out for the basic construction of $k$-simplices. The results below were produced with GAmphetamine.js, but other libraries such as kingdon (py) provide the same functionality.
| Vector algebra | muls/adds | PGA | |
|---|---|---|---|
| Length | $\small \lVert \vec v_1 - \vec v_0 \rVert$ | 3/5 | $\small \lVert v_0 \vee v_1 \rVert$ |
| Area | $\small \tfrac 1 2 \lVert (\vec v_1 - \vec v_0) \times (\vec v_2 - \vec v_0) \rVert$ | 9/11 | $\small \tfrac 1 2 \lVert v_0 \vee v_1 \vee v_2 \rVert$ |
| Volume | $\small \tfrac 1 6 \lvert \big((\vec v_1 - \vec v_0) \times (\vec v_2 - \vec v_0)\big)\cdot(\vec v_3 - \vec v_0) \rvert$ | 9/14 | $\small \tfrac 1 6 \lvert v_0 \vee v_1 \vee v_2 \vee v_3 \rvert$ |
For brevity, let us ignore the square root and scalar prefactor that is equal on both sides and reconstruct the expressions
for the squared norm. Does the PGA expression $(v_0 \vee v_1)\widetilde{(v_0 \vee v_1)}$ for the (square of) the length of a line segment, which involves a regressive and a geometric product, really come at just the cost of that single vector subtraction and dot
we get with $(\vec v_1 - \vec v_0)\cdot(\vec v_1 - \vec v_0)$ in the Vector algebra case? One might be implied to think not, after all the regressive product alone is already far more expensive.
Let us begin by examining the code generated for the join of two points into a line.
GAmphetamine("3DPGA",()=>{
var join = (a, b) => a & b;
var func = Element.compile( join, [Element.point(), Element.point()]);
console.log(func.toString());
});
// output:
// a & b
// 6 muls / 6 adds
function join_point_point (a,b,res=new classes.bivector()) {
const a0=a[0],a1=a[1],a2=a[2],b0=b[0],b1=b[1],b2=b[2];
res[0]=b0-a0;
res[1]=b1-a1;
res[2]=b2-a2;
res[3]=a1*b2-a2*b1;
res[4]=a2*b0-a0*b2;
res[5]=a0*b1-a1*b0;
return res;
}GAmphetamine("3DPGA",()=>{
var squared_norm = (a, b) => (a & b) * ~(a & b);
var func = Element.compile( squared_norm, [Element.point(), Element.point()]);
console.log(func.toString());
});
// output:
// (a & b) * ~(a & b)
// 3 muls / 5 adds
function squared_norm_point_point (a,b) {
const t0=a[0]-b[0],t1=a[1]-b[1],t2=a[2]-b[2];
return t0*t0+t1*t1+t2*t2;
}// (a & b & c) * ~(a & b & c)
// 9 muls / 11 adds
function squared_norm_point_point_point (a,b,c) {
const a0=a[0],a1=a[1],a2=a[2],
_d0=a1-b[1],_d1=a0-c[0],_d2=a0-b[0],_d3=a1-c[1],_d4=a2-b[2],_d5=a2-c[2],
u0=_d2*_d3-_d0*_d1,u1=_d2*_d5-_d4*_d1,u2=_d0*_d5-_d4*_d3;
return u0*u0+u1*u1+u2*u2;
}// (a & b & c & d)
// 9 muls / 14 adds
function join_point_point_point_point (a,b,c,d) {
const a0=a[0],a1=a[1],a2=a[2],d0=d[0],d1=d[1],d2=d[2],
_e0=b[0]-a0,_e1=b[1]-a1,_e2=b[2]-a2,
_d0=a2-c[2],_d1=a1-d1,_d2=a1-c[1],_d3=a2-d2,_d4=a0-c[0],_d5=a0-d0;
return _e0*(_d2*_d3-_d0*_d1)+_e1*(_d5*_d0-_d3*_d4)+_e2*(_d4*_d1-_d2*_d5);
}