# Homogeneous Geometric Algebra

Place a unit cube so that one corner lies at the origin and the opposite corner lies at (1, 1, 1). Consider an axis through the points (0, 1/2, 0) and (1, 1/2, 1), that is, the midpoints of two diagonally opposed edges.

Imagine a plane rotating around this axis. It cuts our cube. What does the cross-section look like?

Geometric algebra’s rotors excel at describing rotations around axes through the origin. But what about arbitrary axes? We could try sandwiching a rotor between two translations, but translation is additive, while rotation is multiplicative, making them awkward to combine.

3D vector algebra suffers from the same affliction. The cure was discovered long ago: add one more dimension, and work with homogeneous coordinates.

The same cure applies to geometric algebra. In fact, homogeneous geometric algebra turns out to be cleaner than its vector algebra counterpart. We gain shockingly concise equations for subspaces of any dimension (lines, planes, and beyond), for linearly transforming subspaces, and for finding the intersection of two subspaces.

On the downside, we wind up representing linear transformations with 4x4 matrices. We must forgo our graceful algebraic formulas for rotations.

We shall see we can avoid matrices by adding two dimensions instead of one to get conformal geometric algebra. For now, we’ll put up with them.

We have deliberately chosen an elementary problem for clarity. We must resist shortcuts: we could avoid translations altogether for our specific problem by centering the cube at the origin. Instead, we’ll insist on general techniques.

We’ll start our code from scratch, though we’ll copy over and tweak some definitions from before. This way, there are no dependencies so we can download the literate Haskell source of this page and try it interactively right away.

{-# LANGUAGE OverloadedLists #-}
module HomogeneousGeometricAlgebra where
import Data.List
import Data.Map (Map)
import qualified Data.Map as M
import Data.Maybe

type Multivector = Map String Double

m!s = M.findWithDefault 0 s m

a ~== b = abs (a - b) < 10**(-6)

gMul :: Multivector -> Multivector -> Multivector
gMul x y = M.filter (/= 0) $M.fromListWith (+)$ basisMul <$> M.assocs x <*> M.assocs y where basisMul (s, m) (t, n) = gnome [] (s ++ t) (m * n) where gnome pre (a:b:rest) c | a < b = gnome (pre ++ [a]) (b:rest) c | a == b = back pre rest c | a > b = back pre (b:a:rest) (-c) where back [] rest c = gnome [] rest c back pre rest c = gnome (init pre) (last pre:rest) c gnome pre rest c = (pre ++ rest, c) canon = gMul [("", 1)] rotor :: [Double] -> Double -> Multivector rotor axis theta = canon [("", -c), ("yz", nx*s), ("zx", ny*s), ("xy", nz*s)] where c = cos (theta / 2) s = sin (theta / 2) [nx, ny, nz] = map (/ sqrt(sum$ map (^2) axis)) axis


Haskell’s list monad succinctly generates the coordinates of the vertices of our unit cube. Then for each occurrence of 0 in a vertex coordinate, we create an edge to the vertex where that 0 is replaced with 1.

cubeVs :: [[Double]]
cubeVs = replicateM 3 [0, 1]

cubeEs = concatMap adj cubeVs where
adj v = map (v:) [[as ++ 1:cs] | v@(as, b:cs) <-
M.fromListWith (+) $catMaybes$ f <$> M.assocs m <*> M.assocs n where f a@(s, _) b@(t, _) | length u == length s + length t = Just c | otherwise = Nothing where c@(u, _) = basisMul a b prop_exampleOuter = wedge [("e", 3), ("f", 4)] [("e", -4), ("f", 3)] == [("ef", 25)]  We follow Macdonald (who cites Dorst) and choose the left contraction for the inner product: $A \cdot B = \langle AB \rangle_{k-j}$ dot :: Multivector -> Multivector -> Multivector dot m n = M.filter (/= 0)$
M.fromListWith (+) $catMaybes$ f <$> M.assocs m <*> M.assocs n where f a@(s, _) b@(t, _) | length u == length t - length s = Just c | otherwise = Nothing where c@(u, _) = basisMul a b prop_exampleInner = dot [("e", 1), ("ef", 2)] [("efg", 3), ("fgh", 4)] == [("fg", 3), ("g", -6)]  For any multivectors $$A$$ and $$B$$, the inner and outer products are dual: $\begin{array}{lll} A \cdot B^* &=& (A \wedge B)^* \\ A \wedge B^* &=& (A \cdot B)^* \end{array}$ We also have the extended fundamental identity: for any vector $$\vv{a}$$ and any blade $$\B$$: $\vv{a} \B = \vv{a} \cdot \B + \vv{a} \wedge \B$ Now we justify our boast. Given a vector $$\vv{a}$$ and a blade $$\B = \vv{b}_1 ... \vv{b}_k$$, suppose we wish to find the unique solution to: $\vv{a} = \vv{a}_\parallel + \vv{a}_\perp$ where $$\vv{a}_\parallel \in \B$$ is the projection and $$\vv{a}_\perp \perp \B$$ is the rejection. Then if $$\vv{a}_\parallel \cdot \B = \vv{a}_\parallel \B$$ is nonzero, then it is a $$(k-1)$$-blade in $$\B$$, and $$\vv{a}_\parallel \wedge \B = 0$$. Dually, if $$\vv{a}_\perp \wedge \B = \vv{a}_\perp \B$$ is nonzero, then it is a $$(k+1)$$-blade in $$\B$$ representing the span of $$\vv{a}_\perp$$ and $$\B$$, and $$\vv{a}_\perp \cdot \B = 0$$. Since every non-zero vector has an inverse, the blade $$\B$$ has an inverse (it is some multiple of $$\vv{b}_k...\vv{b}_1$$). We have: $\begin{array}{lll} \vv{a}_\parallel &=& (\vv{a} \cdot \B) \B^{-1} \\ \vv{a}_\perp &=& (\vv{a} \wedge \B) \B^{-1} \end{array}$ Thus the projection and rejection of a vector with respect to a given blade is indeed readily computed. We also have a test for subspace membership given its direct or dual representation: $\vv{a} \in \B \iff \vv{a} \wedge \B = 0 \iff \vv{a} \cdot \B^* = 0$ ## Wedged blades Recall the geometric product of two vectors is a scalar (the inner product) and a bivector (the outer product) which represents the plane the two vectors lie in. Using our new terminology, we can say if $$\vv{u}_1, \vv{u}_2$$ are linearly independent (and not necessarily orthogonal) then $$\vv{u}_1 \wedge \vv{u}_2$$ is the 2-blade they generate. Gram-Schmidt orthogonalization generalizes this. Given linearly independent $$\vv{u}_1,...,\vv{u}_k$$, the $$k$$-blade they generate is $$\vv{u}_1\wedge...\wedge\vv{u}_k$$. We sketch a geometric algebra proof. The $$k = 1$$ case is vacuously true. We’ll show how the $$k = 2$$ case implies the $$k = 3$$ case and leave the rest to the reader. By inductive assumption, $$\vv{u}_1\wedge\vv{u}_2 = \vv{b}_1 \vv{b}_2$$ for some $$\vv{b}_1 \perp \vv{b}_2$$. Let $$\vv{b}_3 = \vv{u}_3 \wedge \vv{u}_1 \wedge \vv{u}_2 / (\vv{u}_1 \wedge \vv{u}_2)$$, that is, the rejection of $$\vv{u}_3$$ by the blade $$\vv{b}_1 \vv{b}_2$$. Then $$\vv{b}_3$$ is nonzero because the $$\vv{u}_i$$ are linearly independent, and: $\vv{u}_1\wedge\vv{u}_2\wedge\vv{u}_3 = (\vv{u}_1\wedge\vv{u}_2)\wedge\vv{b}_3 = \vv{b}_1\vv{b}_2\wedge\vv{b}_3 = \vv{b}_1\vv{b}_2\vv{b}_3$ Thus we may view a blade as a geometric product of orthogonal vectors, or as the outer product of linearly independent vectors. ## Enter a new dimension The geometric algebra homogeneous model for $$\mathbb{R}^n$$ begins in lockstep with its vector algebra cousin. Define a unit vector $$e$$ orthogonal to $$\mathbb{R}^n$$. (We write $$e$$ with italics to emphasize its special nature; vectors in $$\mathbb{R}^n$$ remain in bold.) We use $$\mathbb{G}^{n+1}$$ to represent $$\mathbb{R}^n$$. A point $$P \in \mathbb{R}^n$$ is represented by the vector $$\newcommand{\p}{\vv{p}} e + \p$$, where $$\p$$ is the vector from the origin to $$P$$. The representation $$e + \p$$ is normalized; for any nonzero scalar $$\lambda$$, $$\lambda(e + \p)$$ also represents $$P$$. The following functions convert back and forth from coordinates of a 3D point and its representation in the geometric algebra homogeneous model. The return journey is hazardous because we must normalize, and also because no point is represented by a vector with no $$e$$ component. (Actually, nonzero vectors with a zero $$e$$ component represent points at infinity, but these lie outside $$\mathbb{R}^3$$.) toHGA :: [Double] -> Multivector toHGA xyz = canon$ M.fromList $("e", 1) : zip (pure <$> "xyz") xyz

fromHGA :: Multivector -> Maybe [Double]
fromHGA m
| d ~== 0 = Nothing
| otherwise = Just $map (get . pure) "xyz" where d = m!"e" get s = m!s / d pt = map (get . pure) "xyz"  How about lines and planes? Here, vector algebra takes a turn for the worse. In computer graphics, typically $$n = 3$$, so the homogeneous model has 4 dimensions. Unfortunately, vector algebra is tied to three dimensions because of the cross product, so either we convert back and forth from 3D to deal with normal vectors, or we have to learn Plücker coordinates. [And what of higher dimensions?] Meanwhile, geometric algebra seems impossibly simple. Our definitions work in any dimension, so four dimensions presents no particular difficulties. In fact, the extra dimension required by the homogeneous model works in our favour. Two distinct points define a line. In the homogeneous model, a point $$\p$$ is represented by scalar multiples of $$p = e + \p$$, that is, a line going through $$p$$ and the origin. Thus the line defined by the two points $$\p, \vv{q}$$ is represented by the plane going through $$p$$, $$q$$, and the origin, that is, the subspace spanned by $$p$$ and $$q$$. Therefore, in geometric algebra, the oriented line passing through $$p$$ and $$q$$ is: $p \wedge q$ Similarly, three distinct points define a plane. In geometric algebra an oriented plane through $$p$$, $$q$$ and $$r$$ is: $p \wedge q \wedge r$ These turn out to be exactly the Plücker coordinates of lines and planes. As claimed, Plücker coordinates arise naturally in geometric algebra. -- https://en.wikipedia.org/wiki/Pl%C3%BCcker_coordinates -- Wikipedia has a worked example showing the line between (2, 3, 7) and -- (2, 1, 0) has Plücker coordinates (0:-2:-7:-7:14:-4). subspace :: [[Double]] -> Multivector subspace = foldl1' wedge . map toHGA prop_examplePlucker = subspace [[2, 3, 7], [2, 1, 0]] == canon [("ex", 0), ("ey",-2), ("ez", -7), ("yz", -7), ("zx", 14), ("xy", -4)]  We can extend this construction to higher dimensions: we represent a $$k$$-blade simply by taking the outer product of the $$k$$ vectors representing linearly independent points defining the subspace. We might ask: why is there no analogous construction the geometric algbra vector model $$\mathbb{G}^n$$? What stops us taking the outer product of two points to represent a line in the vector model? Answer: in the vector model, points are represented by 1-blades (vectors), but for $$k \ge 1$$, $$k$$-blades represent $$k$$-blades. If we take the outer product of two 1-blades representing two points, we wind up with a 2-blade, which represents a 2-blade (plane), and not a 1-blade (line). In the homogeneous model, $$k$$-blades are represented by $$(k+1)$$-blades for all $$k$$. The outer product of two linearly independent 1-blades is a 2-blade, so from two representations of points we get a representation of a line. The outer product of three linearly independent 1-blades is a 3-blade, so three points yield a plane. And so on. In this way, the extra dimension of the homogeneous model works to our advantage. ## Outermorphisms Let $$f$$ be a linear transformation. Then the line defined by $$p$$ and $$q$$ is mapped to the line defined by $$f(p)$$ and $$f(q)$$, thus: $f(p\wedge q) = f(p) \wedge f(q)$ Similarly, we know where $$f$$ maps a given plane if we know where it maps three of its points: $f(p\wedge q\wedge r) = f(p) \wedge f(q) \wedge f(r)$ This extends to blades of any dimension, that is, $$f$$ commutes with the outer product. For this reason, in geometric algebra, a linear transformation is termed an outermorphism. Hence with geometric algbra, we can readily compute the effect of any linear transformation on any subspace. Contrast this with vector algebra, where linearly transforming a plane does not in general transform its normal vector in the same way. For example, a shear may preserve a plane but distort its normal vector, or vice versa. ## Intersections The join of blades $$\vv{A}$$ and $$\vv{B}$$ is the smallest subspace containing both blades, and we denote it by $$\vv{A} \cup \vv{B}$$. The meet of blades $$\vv{A}$$ and $$\vv{B}$$ is the largest subspace contained by both blades, and we denote it by $$\vv{A} \cap \vv{B}$$. For this section we follow the notation in Geometric Algebra for Subspace Operations by Bouma, Dorst, and Pijls. Macdonald writes the vee operator for the meet, but I’m worried this will make it seem it’s closely related to the outer product. Commandeering the union and intersection symbols from set theory does deprive them of their usual meaning, but this is seldom missed in geometric algebra. While efficiently computable (see above paper), neither the join or meet have a simple general geometric algebra formula. However, one can be easily written in terms of the other. If $$\vv{J} = \vv{A} \cup \vv{B}$$ then: $\vv{A} \cap \vv{B} = \vv{A} \vv{J}^{-1} \cdot \vv{B}$ In practice, we often encounter special cases where the join is trivially determined. Frequently, two blades only have the origin in common, in which case their join is simply their outer product. Another common case is when two blades span the whole space: $\vv{A} \cap \vv{B} = \vv{A}^* \cdot \vv{B}$ There’s an issue involving the scalar multiplier possessed by a meet or join, which we blissfully ignore. It means that the sign of our result might be flipped, but this is automatically corrected when normalizing homogeneous coordinates. This equation is reminiscent of De Morgan’s law: a point is a member of $$\vv{A} \cap \vv{B}$$ if and only if it is orthogonal to both $$\vv{A}^*$$ and $$\vv{B}^*$$, that is, it lies in $$(\vv{A}^* \wedge \vv{B}^*)^*$$. Duality of the inner and outer products implies we may replace this with $$\vv{A}^* \cdot \vv{B}$$. In Haskell, for the homogeneous model of $$\mathbb{R}^3$$: dual = (gMul M.singleton "zyxe" 1) -- Only works when inputs span the whole space. meet a b = dual a dot b  We’re working in a projective space, so if given a line parallel to a plane in 3D, meet will return the point at infinity for the direction of the line. This feature is unneeded so we simply use our fromHGA function to remove these answers. ## Cutting the cube Despite their elegance and simplicity, the above equations are practical. We use them to solve our problem. To find the cross-section produced by slicing a plane through a cube, for each edge, we compute its point of intersection with the given plane by finding the meet of the plane and the line defined by the points of the edge. We ignore intersection points that lie outside the cube, and the remaining points determine the polygon we see. clip = filter$ all (\x -> x >= 0 && x <= 1)

cutCube p = clip $mapMaybe (fromHGA . meet p . subspace) cubeEs  Let’s try this for one particular plane containing the axis through (0, 1/2, 0) and (1, 1/2, 1): hexample = cutCube$ subspace [[0, 0.5, 0], [1, 0.5, 1], [0, -0.5, 1]]

prop_hexample = let h = 0.5 in hexample ==
[[0,h,0], [0,0,h], [h,0,1], [h,1,0], [1,h,1], [1,1,h]]


These are the vertices of a regular hexagon.

Next, let’s find the cutting plane for any given angle. We’ll use matrices so we define functions familiar to anyone with experience in computer graphics: embedding a 3x3 matrix into a 4x4 matrix for the homoegeneous model, matrix multiplication by another matrix or by a vector, and building a translation matrix:

mat3To4 :: [[Double]] -> [[Double]]
mat3To4 m = map (++ [0]) m ++ [[0, 0, 0, 1]]

matMul :: [[Double]] -> [[Double]] -> [[Double]]
matMul m n = map (zipWith f (transpose n) . repeat) m where
f bs as = sum $zipWith (*) bs as matVec :: [[Double]] -> [Double] -> [Double] matVec m v = zipWith ((sum .) . zipWith (*)) m$ repeat v

traMat4 :: [Double] -> [[Double]]
traMat4 [x, y, z] = [ [1, 0, 0, x]
, [0, 1, 0, y]
, [0, 0, 1, z]
, [0, 0, 0, 1]
]


We still cling to geometric algebra for the rotations. After computing the rotor from the axis and angle of rotation, we derive the rotation matrix by observing where the rotor sends our basis vectors.

rotMat3 :: [Double] -> Double -> [[Double]]
rotMat3 axis theta = transpose $f . pure <$> "xyz"
where
f s  = flip (M.findWithDefault 0) (g s) . pure <$> "xyz" g s = conj (rotor axis theta) [(s, 1)] conj r u = gMul r$ gMul u $rotorInv r rotMat4 = (mat3To4 .) . rotMat3 rotorInv :: Multivector -> Multivector rotorInv = M.mapWithKey f where f [] v = v f _ v = -v  We start with the plane defined by (0, 0, 0), (0, 0, 1), and (1, 0, 1). Rotating this around the axis going through the origin and (1, 0, 1), then translating it by 0.5 along the $$y$$-axis is equivalent to the rotating plane from our problem. We use an outermorphism to compute where a matrix maps a plane, that is, we apply the matrix to the three points that define the plane, then compute their outer product. Thus the following gives the coordinates of the polygon resulting from cutting the cube with the plane rotated to the given angle: angleCutCube theta = cutCube p where m = matMul (traMat4 [0, 0.5, 0]) (rotMat4 [1, 0, 1] theta) p = foldl1' wedge$ M.fromList . zip (pure <$> "xyze") . matVec m . hom <$> [[0, 0, 0], [0, 0, 1], [1, 0, 1]]


With the help of Haste and a little more code we animate our rotating plane cutting a cube at the top of this very webpage.

Ben Lynn blynn@cs.stanford.edu 💡