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, or by allowing degeneracy to get projective geometric algebra. For now, we’ll put up with matrices.

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 CPP #-}
{-# LANGUAGE OverloadedLists #-}
import Control.Monad
import Data.List
import Data.Map (Map)
import qualified Data.Map as M
import Data.Maybe
#ifdef __HASTE__
import Data.IORef
import Haste.DOM
import Haste.Events
import Haste.Graphics.Canvas
import Haste.Graphics.AnimationFrame
#endif

type Multivector = Map String Double

m ! s = M.findWithDefault 0 s m

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

mul :: Multivector -> Multivector -> Multivector
mul x y = M.filter (/= 0) $
  M.fromListWith (+) $ basisMul <$> M.assocs x <*> M.assocs y

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 = mul [("", 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) <-
    init $ zip (inits v) (tails v), b == 0]

Much ado about 0

Suppose we wish to manipulate a rectangular image. We might use Cartesian coordinates, defining, say, the bottom left corner to be the origin, and the \(\newcommand{\vv}[1]{\mathbf{#1}} \newcommand{\x}{\vv{x}} \newcommand{\y}{\vv{y}} \newcommand{\z}{\vv{z}} x\) and \(y\) axes to be the horizontal and vertical directions. We call this the vector space model of 2D space; the vector \((x, y)\) represents a point.

We can scale, shear, reflect, or rotate with a linear transformation, but never translate. A linear transformation must fix the origin, and a nonzero translation must not.

The trick is to work in a plane not containing the origin, which requires one more dimension. The equations are simplest if we move to a plane that is one unit away from the origin, that is, points of the form \((x, y, 1)\). We can trivially rewrite 2D linear transformations in 3D so they preserve the \(z\) coordinate, and furthermore, with appropriate 3D shear transformations, we can translate our picture.

What about points outside this plane? For any nonzero \(\lambda\), we define \((\lambda x, \lambda y, \lambda)\) to represent the same 2D point as \((x, y, 1)\). With this interpretation, perspective projection is also a linear transformation. We call this the homogeneous model of 2D space.

This trick generalizes to any dimension. For example, in the homogeneous model, every vertex of our unit cube gains a 4th coordinate of 1:

hom = (++ [1])

prop_homCubeCoords = map hom cubeVs ==
  [[0,0,0,1],[0,0,1,1],[0,1,0,1],[0,1,1,1],
   [1,0,0,1],[1,0,1,1],[1,1,0,1],[1,1,1,1]]

The meaning of a multivector

Before we introduce the geometric algebra homogeneous model, we fulfill some promises. We’ve been praising the multivectors of \(\mathbb{G}^n\) because they have useful geometric interpretations. We now reveal the details.

We have the usual interpretations of vectors (1-vectors) and scalars (0-vectors), but that’s unimpressive.

Recall from our discussion of 3D rotations that a plane is represented by the product of two orthonormal vectors it contains. Thus we can geometrically interpret a bivector (2-vector) if it happens to be the product of two orthonormal vectors.

More generally, a k-blade \(\newcommand{\B}{\vv{B}} \B\) or blade of grade \(k\) is the product of \(k\) orthogonal vectors \(\vv{b}_1 …​ \vv{b}_k\), and it represents a \(k\)-dimensional subspace, namely the subspace generated by \(\vv{b}_1,…​,\vv{b}_k\). For example, pseudoscalars like \(\vv{I}\) represent the entire space. Positive multiples of a blade represent the same subspace, while negative multiples represent the same subspace with the opposite orientation, analogous to the anticommutativity of the cross product.

In vector algebra, the norm of the cross product of two vectors is the area of the parallelogram they define. More generally, The norm of \(\B\) is \(|\B| = |\vv{b}_1| …​ |\vv{b}_2|\), and is the \(k\)-volume of a parallelepiped defined by the \(\vv{b}_i\) vectors. We will later build blades from any number of linearly independent vectors, justifying our use of the word "parallelepiped" instead of merely "box".

Recall also with 3D rotations, we multiplied a normal of a plane by \(\vv{I}\) to find the equivalent bivector. More generally, for any dimension, define the dual of a multivector \(A\) to be:

\[ A^* = A \vv{I}^{-1} \]

In \(n\) dimensions, the dual of a \(k\)-blade \(\B\) is the \((n - k)\)-blade \(\B^*\), which represents the orthogonal complement of \(\B\). This, along with many other results, can be seen by choosing the orthonormal basis of \(\mathbb{R}^n\) to include the vectors of \(\B\). We say \(\B\) is the direct representation of the subspace, and \(\B^*\) is the dual representation of the subspace. The dual operation is self-inverse up to sign.

In 3D, a normal vector indirectly represents a plane, that is, a 1-blade is a dual representation of a 2-blade. More generally, in \(n\)-dimensions, all \((n-1)\)-vectors are \((n-1)\)-blades, which we prove with a "double negative": the dual of an \((n-1)\)-vector is a sum of vectors, which is a single vector, whose dual is an \((n-1)\)-blade. We generalize this further to subspaces later.

Thus in 3D, we have \((\vv{u} \wedge \vv{v})^* = \vv{u} \times \vv{v}\), hence as promised, the cross product arises naturally in geometric algebra.

[Multiplication by any unit pseudoscalar results in the orthogonal complement. I’m not sure why right multiplication by \(\vv{I}^{-1}\) is the official dual operation.]

Algebra Geometry

0-vector

scalar

1-vector

line

\((n-1)\)-vector

\((n-1)\)-dimensional subspace

\(n\)-vector

entire space

\(k\)-blade

\(k\)-dimensional subspace

\(e^{\vv{i}\theta}\)

rotation in the plane \(\vv{i}\)

The above table shows some possible geometric interpretations of multivectors. Depending on context, we may attach different meanings, just as vectors in vector algebra are variously interpreted as points, directions, or line segments.

Unlike vector algebra, geometric algebra possesses straightforward notation for subspaces and their orthogonal complements.

Decomposition made easy

In 3D, we had claimed the decomposition of a vector with respect to a plane is "readily computed". Let us make an even bolder boast: in any dimension, the decomposition of a vector with respect to a blade is readily computed. To see how, we generalize the inner and outer products.

For a multivector \(A\), write \(\langle A \rangle_k\) for the \(k\)-vector part of \(A\), that is:

angleBrackets :: Multivector -> Int -> Multivector
angleBrackets m k = M.filterWithKey f m where f s _ = length s == k

There are several ways to generalize the inner and outer products. The following is a natural choice for the outer product. For a \(j\)-vector \(A\) and a \(k\)-vector \(B\), we define:

\[ A \wedge B = \langle AB \rangle_{j+k} \]

We extend this definition to arbitrary multivectors by linearity, that is:

wedge :: Multivector -> Multivector -> Multivector
wedge m n = M.filter (/= 0) $
  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.

Why was there no analogous construction before? 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 = (`mul` 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.

clip01 = filter $ all (\x -> x >= 0 && x <= 1)

cutCube p = clip01 $ 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 = mul r $ mul 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]]

We throw together some code toanimate our rotating plane cutting a cube at the top of this very webpage:

signedArea (x0, y0) (x1, y1) (x2, y2) =
  (x1 - x0)*(y2 - y0) - (x2 - x0)*(y1 - y0)
chain hull          []                 = tail hull
chain hull@(a:b:hs) (p:ps) | s > 0     = chain (b:hs) (p:ps)
                           | otherwise = chain (p:hull) ps
                           where s = signedArea a b p
chain hull          (p:ps)             = chain (p:hull) ps
convexHull = ((++) <$> chain [] <*> (chain [] . reverse)) . sort

#ifdef __HASTE__
wireframe canvas = sequence_ $ renderOnTop canvas . edge <$> cubeEs where
  edge :: [[Double]] -> Picture ()
  edge [p, q] = color (RGB 191 191 191) $ stroke $
    line (project p) (project q)

project :: [Double] -> (Double, Double)
project = scr . zipWith (+) [2, 2, -6] where
  -- Arrange signs so that x-axis goes right, y-axis goes up,
  -- and z-axis goes out of the screen.
  scr :: [Double] -> (Double, Double)
  scr [x, y, z] = ((x/z + 0.5) * 3 * 320 + 160, (-y/z - 0.5) * 3 * 320 + 160)

-- "Sorting" with convex hull avoids self-intersecting polygons.
paint canvas t = renderOnTop canvas $ color (RGB 0 0 0) . fill . path $
  convexHull $ project <$> angleCutCube theta where
  theta = fromIntegral (round t `mod` dur) * 2 * pi / fromIntegral dur
  dur = 20000

main = withElems ["canvas"] $ \[cElem] -> do
  Just canvas <- fromElem cElem
  paused <- newIORef False
  let
    anim t = do
      render canvas $ pure ()
      wireframe canvas
      paint canvas t
      renderOnTop canvas $ color (RGB 255 0 0) $ stroke $
        line (project [0, 0.5, 0]) (project [1, 0.5, 1])
      b <- readIORef paused
      unless b $ void $ requestAnimationFrame anim
    pause = do
      b <- readIORef paused
      writeIORef paused $ not b
      when b $ void $ requestAnimationFrame anim
  _ <- cElem `onEvent` MouseDown $ const pause
  requestAnimationFrame anim
#endif

Ben Lynn blynn@cs.stanford.edu 💡