Projective Geometric Algebra

We copy over definitions from older articles:

The bluntest approach might be to define one null vector for each dimension in our space. Points will certainly be represented by null vectors, but this takes degeneracy too far. In particular, if we want dot products to give useful information, then we need \(n\) non-null basis vectors.

Therefore we start with the standard inner product of \(\mathbb{R}^n\), then add a null vector \(\e_0\) giving us an \((n+1)\)-dimensional space known as \(\mathbb{R}_{n,0,1}\). In our code, we write 0 for \(\e_0\) and handle this case specially, defaulting to the standard behaviour for all other basis vectors.

Next, we move to the dual space \(\mathbb{R}^*_{n,0,1}\) to spread \(\e_0\) around. Recall the dual of a \(k\)-blade is an \((n-k)\)-blade representing its orthogonal complement.

For \(\mathbb{R}^*_{2,0,1}\), if \(\e_1\) and \(\e_2\) are non-null basis vectors (which show up as 1 and 2 in our code), then we represent 1-vectors with the duals of \(\e_0, \e_1, \e_2\).

Already we are hampered by degeneracy. We usually multiply by \(\I\) to compute a dual, but this fails for deriving \(\E_0\) from \(\e_0\), since \(\e_0\) is a null vector. Furthermore, without this portal gun to teleport us between a space and its dual, we are forced to think "backwards" all the time: wedge products describe intersections and meet products describe spans. (If an orthogonal complement grows, then the original subspace shrinks, and vice versa.)

It so happens we can temporarily pretend to have a standard inner product and multiply \(\e_0, \e_1, \e_2\) by \(\I\) to get their duals:

\[ \E_0 = \e_1 \e_2, \E_1 = \e_2 \e_0, \E_2 = \e_0 \e_1 \]

The \(\E_1\) and \(\E_2\) represent the \(x\) and \(y\) directions, and \(\E_0\) is the origin. All but \(\E_0\) are null.

We represent euclidean points and ideal points as follows:

With these representations, for a finite point \(\P\) we have \(\P^2 = -1\), and for an ideal point \(V\) we have \(V^2 = 0\):

As we’re in a projective space, \(\lambda \P\) also represents \(\P\) for any nonzero scalar \(\lambda\); a finite point is normalized when \(\P^2 = -1\). Ideal points always square to zero.

Like before, we avoid geometrically meaningless squares thanks to null vectors. The square of a (normalized) point is a constant, so is independent of the origin.

We test our code by reproducing the multiplication table of the geometric algebra basis vectors from one of Gunn’s papers:

The line \(a x + b y + c = 0\) is represented by:

A euclidean line \(\a\) is normalized when \(a^2 = 1\).

We provide similar routines for 3D geometry:

Our demo needs to retrieve the coordinates of a euclidean point for the 3D case:

Define the squared norm of a euclidean \(k\)-vector \(\vv{X}\) by:

\[ || \vv{X} ||^2 = \vv{X}^2 \]

(This formula does not apply to all \(k\)-vectors, but rather only those representing points, lines, planes, and so on.)

The square of an ideal element is always zero, so we also define the square of the ideal norm \(||\vv{X}||^2_\infty\) to be the sum of the squares of the coefficients, and say an ideal \(k\)-vector \(\vv{X}\) is normalized when \(||\vv{X}||^2_\infty = 1\).

This can be formulated without coordinates:

\[ ||\vv{X}||^2_\infty = ||\vv{X}\vee\P||^2 \]

where \(\P\) is any normalized finite point. However, our implementation simply reads off the coefficients:

Recall in the conformal model we computed wedge products to find the span of subspaces, such as the line defined by two points, or the plane defined by three points. Because PGA uses the dual space, we need meet products instead.

Usually we take advantage of \(A \vee B = (A^* \wedge B^*)^*\). But without an easy dual operation, we resort to a shuffle product to compute the meet, a sort of geometric version of the inclusion-exclusion principle.

We need a helper function that returns all partitions of a sequence into two subsequences, along with the sign of the permutation obtained by concatening the subsequences.

To compute meets with a shuffle, we need the dimension of the space, which in a cleaner library might be held somewhere by each Multivector value. (Such a library could also prevent mixing up multivectors from different spaces.) Instead, we require the caller to supply the dimension.

Let’s start with the geometric product in 2D on normalized inputs.

Multiplication by \(\I\) is commutative, and computes the orthogonal complement in the Euclidean 2D space. (Not the geometric algebra itself, otherwise we would have easy access to the dual!)

For a euclidean line \(\a\), the product \(\a\I\) is the ideal point in the direction perpendicular to the line.

For a euclidean point \(\P\), the product \(\P\I\) is \(\e_0\), the ideal line.

The product of two euclidean lines is:

\[ \a\b = \a \cdot \b + \a \wedge \b \]

The dot product \(\a \cdot \b\) is the cosine of the angle \(\theta\) between them.

Let \(\P\) denote their normalized point of intersection. If \(\P\) is finite, that is, if the lines intersect, then \(\a \wedge \b = (\sin \theta) \P\). If \(\P\) is ideal, that is, the lines are parallel, then \(\a \wedge \b = d_{ab} \P\) where \(d_{ab}\) denotes the distance between the lines.

The product of two euclidean points is:

\[ \P \Q = -1 + d_{PQ} \vv{V} \]

where \(d_{PQ}\) is the distance between them and \(\vv{V}\) is an ideal line representing the direction perpendicular to the line joining them. We define

\[ \P \times \Q = d_{PQ} \vv{V} \]

hence \(d_{PQ} = || \P \times \Q ||_\infty\).

The product of a euclidean line and a euclidean point is:

\[ \a \P = \a \cdot \P + \a \wedge \P = \a^\perp_\P + d_{\a\P}\I \]

where \(\a^\perp_\P = \a \cdot \P\) is the line through \(\P\) perpendicular to \(\a\) and \(d_{\a\P}\) is the distance between the point and the line.

It turns out we also have meaningful products when one or both of the inputs are ideal.

There are more geometric interpretations to explore and we have yet to touch the 3-dimensional case. But let’s fast-forward to the part we’ve all been waiting for: rotors!

The reflection of \(\vv{X}\) in the \((n-1)\)-hyperplane \(\a\) is:

\[ \a\vv{X}\a \]

This formula reflects points, lines, planes, and so on.

We could finish here, as all Euclidean isometries can be generated from reflections. Composing reflections through \(\a\) then \(\b\) is simply:

\[ \b\a\vv{X}\a\b \]

If \(\vv{R} = \b \a\) is normalized then we call it a rotor, and transform a subspace \(\vv{X}\) by computing \(\vv{R} \vv{X} \tilde{\vv{R}}\), where \(\tilde{\vv{R}}\) is the reversal of \(\vv{R}\), meaning we multiply the factors in reverse order. Calculation shows \(\vv{R}\tilde{\vv{R}} = 1\).

Rotations have an exponential form:

\[ e^{\frac{\theta}{2} \P} \]

where \(\P\) is the normalized point (in 2D) or axis (in 3D) at the center of a rotation by the angle \(\theta\). For euclidean points, \(\P^2 = -1\) so \(\P\) takes the role played by \(\vv{i}\) in other models.

For example, a rotation around the \(y\)-axis by 120 degrees is given by:

Observe we computed the meet product of two points to find the multivector describing the \(y\)-axis. It turns out to have a simple form, though in general lines in 3D are less straightforward in PGA:

In 2D, translation by \(d\) units can be viewed as a rotation of size \(d\) with the center at an ideal point \(\P\) in the perpendicular direction. By considering the Taylor series expansion, and recalling \(\P^2 = 0\), we find:

\[ e^{\frac{d}{2} \P} = 1 + \frac{d}{2} \P \]

For a 3D translation, the axis is an ideal line, which Gunn states can be found by computing \((\E_0 \vee \vv{V})\I\) where \(\vv{V}\) is the direction of travel. It seems we find the line joining the origin and \(\vv{V}\) then take its orthogonal complement to yield the axis.

\[ e^{\frac{d}{2} (\E_0 \vee \vv{V})\I} \]

For example, the following rotor translates 12 units along the \(z\)-axis:

Rotations and translations involving the coordinate axes have a simple form so are great for demos.

Above, our definition for \(V\) uses notation from our adventure in the conformal model: the symbols \(\infty, \x, \y, \z\) correspond to \(\e_0,\e_1,\e_2,\e_3\) which are 0 1 2 3 in our code.

Lastly, we slap together some code to animate a cube on this webpage.