= 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?
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
link:ga.html[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 link:cga.html[conformal geometric algebra], or by allowing degeneracy
to get link:pga.html[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
link:ga.html[some definitions from before]. This way, there are no dependencies
so we can download link:hga.lhs[the literate Haskell source of this page] and
try it interactively right away.
\begin{code}
{-# 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
\end{code}
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.
\begin{code}
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]
\end{code}
== 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:
\begin{code}
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]]
\end{code}
== 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
link:ga.html[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.]
[cols="1,3",options="header",width="80%",align="center"]
|=============================================================================
| 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:
\begin{code}
angleBrackets :: Multivector -> Int -> Multivector
angleBrackets m k = M.filterWithKey f m where f s _ = length s == k
\end{code}
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:
\begin{code}
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)]
\end{code}
We follow Macdonald (who cites Dorst) and choose the 'left contraction'
for the inner product:
\[ A \cdot B = \langle AB \rangle_{k-j} \]
\begin{code}
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)]
\end{code}
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$.)
\begin{code}
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"
\end{code}
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.
\begin{code}
-- 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)]
\end{code}
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
http://lanl.arxiv.org/pdf/math/0104102v1['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$:
\begin{code}
dual = (`mul` M.singleton "zyxe" 1)
-- Only works when inputs span the whole space.
meet a b = dual a `dot` b
\end{code}
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.
\begin{code}
clip01 = filter $ all (\x -> x >= 0 && x <= 1)
cutCube p = clip01 $ mapMaybe (fromHGA . meet p . subspace) cubeEs
\end{code}
Let's try this for one particular plane containing the axis through
(0, 1/2, 0) and (1, 1/2, 1):
\begin{code}
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]]
\end{code}
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:
\begin{code}
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]
]
\end{code}
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.
\begin{code}
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
\end{code}
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:
\begin{code}
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]]
\end{code}
We throw together some code toanimate our rotating plane cutting a cube at the
top of this very webpage:
\begin{code}
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
\end{code}
link:cga.html[With one more dimension, we can do much more].