Geometric Algebra (GA) has found success in various areas of the physical sciences and engineering over the last decade but remains relatively underutilised in industry and several key topics in the field remain unexplored. This thesis focuses on the practical applications of Geometric Algebra in various interconnected areas of mathematical engineering. In Part I we explore the properties of the objects resulting from the addition of blades in Conformal Geometric Algebra (CGA) and how we might use these objects in computer graphics and robotics algorithms. In Part II we explore how Screw Theory embeds into CGA, how to use this embedding for simulation of the dynamics of rigid bodies, and how practitioners can leverage the geometric primitives built into CGA to represent and solve constraints in multi-body robotic systems.

Computer Graphics, Computer Vision and Visualisation

In this Part of the Thesis we develop techniques that may find application in the fields of graphics, vision and visualisation. Graphics, vision and visualisation have long been a promising application area for geometric algebra and indeed much of the present interest in the field is driven from the perspective of computer graphics. Here we link the generation of rotors directly from geometric objects to the idea of direct addition of geometric primitives in CGA. We explore applications of addition of CGA objects in 3D computer vision scenarios and then extend the idea of addition of objects to interpolation and extrapolation and construct and visualise evolved splines and surfaces.

Calculating the rotor between conformal
objects

We shape our tools, and thereafter our tools shape us.John M. Culkin

Introduction

In this chapter we will address the problem of recovering covariant transformations between objects – specifically; lines, planes, circles, spheres and point pairs. Using the covariant language of conformal geometric algebra (CGA), we will derive such transformations in a very simple manner. In CGA, rotations, translations, dilations and inversions can be written as a single rotor, which is itself an element of the algebra. We will show that the rotor which takes a line to a line (or plane to a plane etc.) can easily be formed and we will investigate the nature of the rotors formed in this way.

If we can recover the rotor between one object and another of the same type, a useable metric which tells us how close one line (plane etc.) is to another, can be a function of how close this rotor is to the identity. Using these ideas, we find that we can define metrics for a number of common problems, specifically recovering the transformation between sets of noisy objects.

Related Work

Our primary aim in this chapter is to simultaneously estimate the rotation and translation that takes one object (line to line/circle to circle/plane to plane/sphere to sphere/point-pair to point-pair) to another. There are many methods that estimate rigid body transformations with points Eggert et al. [1997]Valkenburg and Dorst [2011]Tingelstad and Egeland [2017]De Keninck and Dorst [2019]. In Hitzer et al. [2009] the authors estimate a general rotor between arbitrary objects using the idea of carriers – while interesting, this method lacks simplicity and does not deal directly with the objects themselves.

Conformal Geometric Algebra

, plus two additional basis vectors, $e$

and $\bar{e}$ with signatures, $e^2=1$

, $\bar{e}^2=-1$ . Two null vectors can then be defined as: $n_{\infty}=e + \bar{e}$ and $n_0 = \frac{e - \bar{e}}{2}$ . The mapping of a 3D vector $x$

to its conformal representation $X$

is given by $X = F(x) = \frac{1}{2}(x^2 n_{\infty} + 2x - 2n_0)$ . Many of our target applications will be in computer vision, and in investigating algorithms which use more than just points, which is the case with most conventional computer vision algorithms.

A Rotor between Objects

Suppose we wish to find the rotor (rotation, translation, dilation) which takes an object $X_1$ to an object $X_2$ (where $X_1$ and $X_2$ are conformal $n$ -blades representing the lines/circles/planes/ spheres/point pairs). If we firstly take lines as an example, conventionally we would translate along the common perpendicular and then rotate about the intersection point – which requires a series of non-trivial geometric operations for two arbitrary lines in space. Here we seek a method which will not require reverting to the geometric properties of the lines, but which will give the transformation in terms of the lines themselves – and we wish this method to be valid for all objects. In CGA, let the rotor which takes $X_1$ to $X_2$ be $R_x$ , where this comprises both rotation, translation and dilation rotors. We assume both objects are normalised such that $X_1^2 = X_2^2 = \gamma$ , where $\gamma= 1$ for lines, circles and point pairs, and $\gamma = -1$ for planes and spheres:

$\displaystyle X_2 = R_x X_1 \tilde{R}_x$

Note, that $\tilde{X}= -\gamma X$ . We motivate our approach by considering the quantity $(X_1 + X_2)$ which is in some sense the `average' object; ie, if we reflect $X_1$ in , we should get some function of $X_2$ (we assume for convenience that $X^2=1$ , ie $\gamma= 1$ ):

$\displaystyle (X_1 + X_2)X_1(X_1 + X_2)$	$\displaystyle =$	$\displaystyle (1 + X_2X_1)(X_1 + X_2)$
	$\displaystyle =$	$\displaystyle [2 + (X_1X_2 + X_2X_1)]X_2 \equiv KX_2$	(1)

So the reflection does indeed produce a multiple, though the multiple is a scalar plus 4-vector, of $X_2$ . Since we can write the LHS of equation 2.1 as

$\displaystyle (X_1 + X_2)X_1(X_1 + X_2) = [(X_1 + X_2)X_1]X_1[X_1(X_1 + X_2)] = (1 + X_2X_1)X_1(1 + X_2X_1)\,\tilde{} \vspace{0.3cm}$

we propose to use the spinor quantity $Z = 1 + \gamma X_2X_1$ to form $R_x$

. As above (but now with $\gamma$ included) $Z X_1 \tilde{Z}$ gives ;

$\displaystyle Y = Z X_1 \tilde{Z} = 2 X_2 + \gamma(X_1 X_2 + X_2 X_1)X_2 = (2 + \gamma M_{12})X_2 = K X_2 \vspace{0.3cm}$

(2)

where $M_{12} = X_1 X_2 + X_2 X_1$ is the anticommutator of $X_1$

and

. Thus, we see that takes to a multiple of , where this multiple involves the anticommutator of the objects. In general this anticommutator will have scalar and 4-vector parts (the bivector part of $X_1X_2$

cancels with the bivector part of $X_2X_1$

Since all 4-vectors square to give a scalar, we can take $K^* = \langle K \rangle_0 - \langle K \rangle_4$ , such that $K K^* = \langle K \rangle^2_0 - \langle K \rangle^2_4$ , is a scalar, which we call $\mu$ . We show later that $\mu$ is always positive. We now multiply both sides of equation 2.2 by $K^*$ to give:

$\displaystyle \frac{1}{\mu}K^*Z X_1 \tilde{Z} = X_2$

We now look to split up $K^*$ such that $S^2 = K^*$ , where $S = \alpha + \beta M_{12} \equiv (\alpha +\beta\langle M_{12}\rangle_0) + \beta \langle M_{12}\rangle_4$ and $\alpha$ and $\beta$ are scalars. If $S$ takes this form, it is clear that it is both self-reverse and commutes with $Z$ and $X_1$ ; we can therefore write

$\displaystyle \left(\frac{1}{\sqrt{\mu}}SZ\right) X_1 {\left(\frac{1}{\sqrt{\mu}}SZ\right)}^{\sim} = X_2 \vspace{0.3cm}$

so that $\frac{1}{\sqrt{\mu}}SZ$ is our required rotor and $\mu = K^*K$ . To find such an $S$

we can use the square root formula given in Dorst and Valkenburg [2011] or simply equate scalar and 4-vector parts of the equation $S^2 = K^*$

. We do the latter first in order to see how the particular form of our scalar plus 4-vector behaves and then confirm that it agrees with the formula in Dorst and Valkenburg [2011]:

$\displaystyle (\alpha +\beta\langle M_{12}\rangle_0)^2 + 2\beta (\alpha +\beta\... ... \beta^2 \langle M_{12}\rangle^2_4 = \langle K \rangle_0 - \langle K \rangle_4$

Since $\langle K \rangle_0=2 + \gamma\langle M_{12}\rangle_0$ and $\langle K \rangle_4 = \gamma\langle M_{12} \rangle_4$ , we have:

$\displaystyle (\alpha +\beta\langle M_{12}\rangle_0)^2 + \beta^2 \langle M_{12}\rangle^2_4$	$\displaystyle =$	$\displaystyle \langle K \rangle_0$
$\displaystyle 2\beta (\alpha +\beta\langle M_{12}\rangle_0)\langle M_{12} \rangle_4$	$\displaystyle =$	$\displaystyle - \gamma\langle M_{12} \rangle_4$

From equating 4-vector parts we see that $2\beta(\alpha +\beta\langle M_{12}\rangle_0) =-\gamma$ so that, provided $\langle M_{12} \rangle_4 \ne 0$ ;

$\displaystyle S = -\frac{\gamma}{2\beta} +\beta \langle M_{12} \rangle_4$

If $\langle M_{12} \rangle_4 =0$ we simply have $S = \sqrt{ \langle K\rangle_0}$ if $\langle K\rangle_0$ is positive, which it is for lines, planes, circles and point pairs. $\langle K\rangle_0$ can take negative values for some sphere cases. If $\langle M_{12} \rangle_4 \ne 0$ we then find $\beta$ from the equation which equates scalar parts:

$\displaystyle \frac{1}{4\beta^2} - \beta^2 \lambda = \langle K \rangle_0 \vspace{0.3cm}$

where $\langle M_{12}\rangle^2_4 \equiv \langle K \rangle_4^2 = -\lambda$ , since the 4-vectors always square to give zero or a negative scalar. This is a quadratic in $u = \beta^2$ :

$\displaystyle 4\lambda u^2 + 4\langle K \rangle_0 u -1 = 0 \vspace{0.3cm}$

(3)

with solutions given by:

$\displaystyle u = \frac{- 4\langle K \rangle_0 \pm 4\sqrt{ \langle K \rangle^2_0 + \lambda}}{8\lambda}$

As $\beta^2 = u$ we need the solution which is guaranteed to be positive:

$\displaystyle \beta^2 = \frac{1}{2\lambda}\left(\sqrt{ \langle K \rangle^2_0 + ... ..._0 \right) = \frac{1}{2\lambda}\left(\sqrt{ \mu} - \langle K \rangle_0 \right)$

Recall $K = 2 + \gamma (X_1X_2 + X_2 X_1) = \langle K \rangle_0 + \langle K \rangle_4$ , $K^* = \langle K \rangle_0 - \langle K \rangle_4$ , $\lambda = -\langle K \rangle^2_4$ , $\mu = K^*K = \langle K \rangle^2_0 + \lambda$ , so is always positive (as $\lambda \ge 0$ ). We can now write the explicit form of the rotor as:

If $\langle M_{12} \rangle_4 \ne 0$ :

$\displaystyle \boxed{ R_x = \frac{1}{\sqrt{\mu}}\left( -\frac{1}{2\beta} +\beta \langle K \rangle_4 \right)(1 + \gamma X_2 X_1)}$ (4)

$\displaystyle \boxed{ \beta^2 = \frac{1}{2\left(\sqrt{ \mu} +\langle K \rangle_0 \right)} }$ (5)
If $\langle M_{12} \rangle_4 =0$ and $\langle K \rangle_0 > 0$

$\displaystyle \boxed{ R_x = \frac{1}{\sqrt{\langle K \rangle_0}}(1 + \gamma X_2 X_1)}$ (6)
If $\langle M_{12} \rangle_4 =0$ and $\langle K \rangle_0 < 0$ ,

$\displaystyle \boxed{ R_x = \frac{1}{\sqrt{\vert\langle K' \rangle_0\vert}}(1 + \gamma \bar{X}_2 X_1)}$ (7)

where $\bar{X}_2 = -X_2$ and $K' = 2 + \gamma (X_1\bar{X}_2 + \bar{X}_2 X_1)$ .

Taking the positive or negative square root for $\beta$ simply changes the sign of the rotor, which makes no difference to the transformation. These expressions hold for all CGA objects: lines, planes, circles, spheres, point pairs. The following subsection will give the explicit forms for each of these objects and will discuss the third case which can occur for spheres.

Before looking in more detail at the nature of the rotors formed by the process outlined here, we return to equation 2.1 and note that we can now take $X_1$ to $X_2$ via a reflection in the quantity $X_m$ where

$\displaystyle X_m = \frac{S}{\sqrt{\mu}}(X_1 + X_2) \vspace{0.3cm}$

where

and $\mu$ are as given previously, ie $\mu = K^*K$ and $S$

takes the form in equations 2.4,2.6,2.7 depending on the nature of $M_{12}$ . We will see in Chapter 3 that the quantity $\frac{S}{\sqrt{\mu}}$ projects the $m$

-vector obtained from the addition of the two blades $X_1$

and

onto an

-blade and therefore an object – the object being that in which we reflect $X_1$

in to get

We can also confirm the solutions in equations 2.4,2.6,2.7 using the result in Dorst and Valkenburg [2011], where the square root of the scalar plus 4-vector, $\Sigma$ , is given by

$\displaystyle \sqrt{\Sigma} = \frac{\Sigma \pm [[\Sigma]]}{\sqrt{2}\sqrt{\langl... ...rangle_4}{\sqrt{2}\sqrt{\langle{\Sigma}\rangle \pm [[\Sigma]]}} \vspace{0.3cm}$

where $[[\Sigma]] = \sqrt{\langle{\Sigma}\rangle^2 - \langle{\Sigma}\rangle_4^2}$ . Here, our $\Sigma = K^* = \langle K\rangle_0 -\langle K\rangle_4$ , so that (taking the solution corresponding to the $+$

sign):

$\displaystyle [[\Sigma]] = \sqrt{\mu} \;\;\;$ and $\displaystyle \;\;\; \sqrt{2}\sqrt{\langle{\Sigma}\rangle + [[\Sigma]]} = \pm \frac{1}{\beta} \vspace{0.3cm}$

giving $\sqrt{\Sigma} = -\frac{1}{2\beta} + \beta \langle K\rangle_4$ , as required (taking $-\frac{1}{\beta}$ ).

Lines

Conformal lines take the form $L = A\wedge B\wedge n_{\infty}$ , with $A$

being the conformal representations of two points lying on the line, and $n_{\infty}$ the point at infinity. $\tilde{L} = -L$ and we normalise such that $L^2=1$

, therefore $\gamma= 1$ . For lines, the 4 vector part of the anticommutator takes the form $\beta I_5 n_{\infty}\equiv \beta I_3 n_{\infty}$ , thus the square of this is always zero, which means $\lambda = 0$ and $\mu = \langle K \rangle^2_0$ , which reduces equation 2.3 to $u = 1/(4\langle K \rangle_0)$ and $\beta = \pm 1/(2\sqrt{\langle K \rangle_0}$ [note that it does not matter which sign we take], giving us the simpler form of the rotor as:

$\displaystyle \boxed{ R = \frac{1}{{\langle K \rangle_0}}\left( \sqrt{\langle K... ...eft(1 - \frac{\langle K \rangle_4}{2\langle K \rangle_0} \right) (1 + L_2 L_1)}$

(8)

Planes

With planes, as with lines, there is no issue of scaling as the objects are infinite. A plane $\Pi$ is taken to be the conformal 4-blade of the form $A\wedge B\wedge C\wedge n_{\infty}$ , with $A,B,C$ any 3 conformal points lying on the plane. Conformal planes square to a negative number, so we assume that planes are normalised such that $\Pi^2=-1$ , therefore $\gamma = -1$ . Note that $\tilde{\Pi} = \Pi$ .

For planes the anticommutator is a scalar and it is not hard to show that (for normalised planes) $\langle K\rangle_0$ is always positive. Thus, the form for the rotor in the plane-to-plane case is particularly simple as the $\langle K \rangle_4$ term vanishes:

$\displaystyle \boxed{ R_{\Pi} = \frac{1}{\sqrt{\langle K \rangle_0}} (1 - \Pi_2\Pi_1) } \vspace{0.3cm}$

(9)

where $K = 2 - (\Pi_1\Pi_2 + \Pi_2\Pi_1)$ .

Circles

One might think that the case of circles-to-circles would be more complex, as a transformation which takes one arbitrary circle to another involves a dilation as well as a rotation and translation. However, nothing in the above derivation assumed anything specific about the rotor, and we find that we can use precisely the same formula to move between arbitrary circles.

Let us start with two conformal circles, $C_1$ and $C_2$ not necessarily of the same radius. A conformal circle is a 3-blade of the form $P\wedge Q\wedge R$ , where $P,Q,R$ lie on the circle. Circles square to a positive scalar, so we will assume that our circles are normalised such that $C^2 = 1$ and therefore $\gamma= 1$ . Note that $\tilde{C}=-C$ .

The anticommutator, $M_{12}$ , is in general a scalar plus 4-vector, so we must use the form given in equations 2.4,2.5 and little simplification is possible:

$\displaystyle \boxed{ R_c = \frac{1}{\sqrt{\mu}}\left( -\frac{1}{2\beta} +\beta \langle K \rangle_4 \right)(1 + C_2 C_1)}$

(10)

$\displaystyle \boxed{ \beta^2 = \frac{1}{2\lambda}\left(\sqrt{\mu} - \langle K \rangle_0 \right)} \vspace{0.3cm}$

(11)

with $K = 2 + (C_1C_2 + C_2C_1), \;\; \mu = K^*K$ .

Spheres

We start with two conformal spheres, $S_1$ and $S_2$ not necessarily of the same radius. A conformal sphere is a 4-blade of the form $N\wedge P\wedge Q\wedge R$ , where $N,P,Q,R$ lie on the sphere. Circles square to a negative scalar, so we will assume that our spheres are normalised such that $S^2 = -1$ and therefore $\gamma = -1$ . Note that $\tilde{S}=S$ .

As for planes, $\langle K \rangle^2_4$ is zero, so the rotor takes a very simple form:

$\displaystyle \boxed{ R_{s} = \frac{1}{\sqrt{\vert\langle K \rangle_0\vert}} (1 - \bar{S}_2S_1) } \vspace{0.3cm}$

(12)

where $K = 2 - (S_1\bar{S}_2 + \bar{S}_2S_1)$ , $\bar{S_2} = S_2$ if $\langle K \rangle_0 > 0$ and $\bar{S_2} = -S_2$ if $\langle K \rangle_0 < 0$ . $-S_2$

is the same sphere as $S_2$

, so in a sense it does not matter whether we take $S_1$

or to

– this additional complexity occurs with spheres as they lack any intrinsic orientation, which is not the case for lines, planes, circles and point pairs.

Point Pairs

In the conformal setting, point pairs take the form $A\wedge B$ where $A,B$

are conformal points – we can think of a point pair as a line segment. For a point pair, $P$

, clearly $\tilde{P} = -P$ and $P^2$

gives a positive scalar. We will therefore assume that point pairs are normalised so that $P^2 = 1$

Since the anticommutator will generally have both scalar and 4-vector parts, we again have the general form taken from equations 2.4,2.5:

$\displaystyle \boxed{ R_p = \frac{1}{\sqrt{\mu}}\left( -\frac{1}{2\beta} +\beta \langle K \rangle_4 \right)(1 + P_2 P_1)}$

(13)

$\displaystyle \boxed{ \beta^2 = \frac{1}{2\lambda}\left(\sqrt{\mu} - \langle K \rangle_0 \right)} \vspace{0.3cm}$

(14)

with $K = 2 + (P_1P_2 + P_2P_1), \;\; \mu = K^*K$ .

Lines to Circles: Planes to Spheres

Note that in the previous rotor derivation we assumed $X_1$ and $X_2$ were blades of the same grade, but nothing further. Therefore, we should, and indeed do, find that the rotor formulae in equations 2.4-2.7 work for moving between lines and circles and between planes and spheres.

The Non-Uniqueness of the Recovered Rotors

Although we have recovered rotors for each case of lines, planes, circles, spheres and point pairs, it is clear that these rotors are not unique. For example, if we transform one line into another, we can then translate along the second line without altering the result. So, a natural question to ask is exactly what is the transformation we are recovering with the expression.

To investigate this further we extract the bivector, $B$ , for each recovered rotor, with $R=$ e $^B$ , and plot the interpolated objects for each of $\lambda_i$ , $i=1,..n$ , with $X_i =$ e $^{\lambda_iB}X_1$ e $^{-\lambda_iB}M_0$ , where $X_2 =RX_1\tilde{R}$ and $\lambda_i = i/n$ . Figure 2.1 shows these interpolations for each class of object.

Conclusion

In this chapter we have presented a general framework for extracting the conformal rotor that takes a conformal object of a given grade to another conformal object of the same grade. The technique works for point pairs, lines, circles, planes and spheres. In the process of investigating these rotors we have touched on the form of the object required to reflect one object into another and by visualising intermediate objects we have verified that the rotors take the objects smoothly to each other. Code that implements this rotor extraction algorithm is available in the clifford Hadfield et al. [sent] python package and novel applications of this technique are additionally presented in Eide and Lasenby [2018] Hadfield and Lasenby [2019] and Hadfield et al. [2019]. It is also interesting to note that the nature of the quantity $X_2X_1$ was investigated first in Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004], and then in Dorst et al. [2007], and noted to produce a quantity which was $R^2$ , where $R$ is the rotor taking $X_1$ to $X_2$ . This has also been used for interpolations between objects in Colapinto [2011]. Here we have given explicit expressions for the rotor itself and investigated the a range of use cases.

**Fig. 2.1:** The red objects in each of these images show the interpolations of the rotors formed from pairs of black objects. Here we see that a range of geometric primitives including circles, lines and point pairs are all handled elegantly by the same framework.

Direct linear interpolation of geometric objects in conformal geometric algebra

I write rhymes with addition and algebra, mental geometry.Ice T, Mind Over Matter

Abstract

Typically we do not add objects in conformal geometric algebra (CGA), rather we apply operations that preserve grade, usually via rotors, such as rotation, translation, dilation, or via reflection and inversion. However, here we show that direct linear interpolation of conformal geometric objects can be both intuitive and of practical use. We present a method that generates useful interpolations of point pairs, lines, circles, planes and spheres and describe algorithms and proofs of interest for computer vision applications that use this direct averaging of geometric objects.

Introduction

In this chapter we will look at adding CGA objects and adjusting the resulting multivectors to produce useful interpolations of the objects. We will present a general technique that is valid for all geometric objects of grade 2 or above. This technique uses the decompositions presented in Dorst and Valkenburg [2011].

The objects we work with here will be CGA objects unless explicitly stated otherwise. We will use the standard extension of the 3d geometric algebra, where our 5D CGA space is made up of the standard spatial basis vectors $\{e_i\}$ $i=1,2,3$ , plus two additional basis vectors, $e$ and $\bar{e}$ with signatures, $e^2=1$ , $\bar{e}^2=-1$ . Two null vectors can therefore be defined as: $n_{\infty}=e + \bar{e}$ and $n_0 = \frac{e - \bar{e}}{2}$ . The mapping of a 3d vector $x$ to its conformal representation $X$ is given by $X = F(x) = \frac{1}{2}(x^2 n_{\infty} + 2x - 2n_0)$ .

Motivation

In our conformal representation of space the blades of a given grade $n$

represent specific classes of object and lie on a manifold within the overall subspace of grade $n$

. Typically in geometric algebra we traverse this manifold using rotors and reflections. These transformations are incredibly useful and make up the vast majority of operations used in the field of applied geometric algebra. Unfortunately, while being of geometric significance, rotors and reflections are often unintuitive ways of thinking about a problem and traditional algorithms often require significant rehashing to fit within this framework.

For example: Given a cluster of geometric objects we would like to be able to create an `average' object that lies in some sense in the middle of the bundle. Most geometric algebra approaches to this problem would likely require the explicit design of a geometrically motivated cost function followed by constrained optimisation on the blade manifold either directly or via a parametrisation of rotors over the space.

While this approach has been very effective for a variety of problems, it requires the careful crafting of clever cost functions, consideration of the convexity of the underlying space, efficient implementation of the given optimisation scheme etc. The question we aim to answer here is: what if we just decided to add all the objects together?

Linearly interpolating conformal points

The result of linear combinations of conformal points is well known Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]. Consider two arbitrary points in 3d space $a$ and $b$ represented as $A$ and $B$ in our conformal model. Linear interpolation of these points followed by our conformal embedding produces a linear interpolation of our conformal points with an additional term in $n_{\infty}$ :

$\displaystyle F( \alpha a + (1-\alpha) b ) = \alpha A + (1-\alpha) B + \alpha(1-\alpha)(A\cdot B)n_{\infty}$

(15)

We can therefore get a useful interpolation of points by taking a direct linear interpolation and simply adding the final $\alpha(1-\alpha)(A\cdot B)n_{\infty}$ term to the result. If $Y= \alpha A + (1-\alpha) B$ , we can recover $Y' = F( \alpha a + (1-\alpha) b )$ via the following formula (assuming $Y'\cdot n_{\infty}=-1$ ):

$\displaystyle Y' = \frac{-Yn_{\infty}Y}{2(Y\cdot n_{\infty})^2}$

Linearly interpolating higher grade conformal objects

Objects of grade 2 and above are more difficult to interpolate in a sensible and computationally efficient way. Typically, schemes that have been found are either only valid for certain objects in specific cases Doran [2003], or the problem is attacked indirectly via carriers Hitzer et al. [2009] or by forming the rotor between the objects, extracting the corresponding bivector, which is then interpolated Wareham and Lasenby [2008] and applied to the first object.

It was shown in Chapter 2 that we can represent the mirror object $X_m$ that reflects one object $X_1$ , into another $X_2$ , as the left multiplication of the summation of the blades by a scalar + 4-vector factor $S$ :

$\displaystyle X_m = S(X_1 + X_2) = (\beta + \gamma(X_1X_2 + X_2X_1))(X_1 + X_2)$

(16)

where $\beta$ and $\gamma$ are scalars and the 4-vector part of $S$

is proportional to the anticommutator of $X_1$

and

For the previously known cases in which the linear interpolation of higher grade objects gives a blade, such as with circles Doran [2003] and point pairs Dorst et al. [2007] both with common points, the factor $S$ is a scalar and the object $X_m$ is simply `half-way' between the objects. We can extend this notion to the cases where the addition of objects is not a blade by using our object $X_m$ , which has been corrected to being a blade, as the half-way object. We can use this idea of the half-way object to recursively subdivide the space between $X_1$ and $X_2$ allowing us to create objects that are any fraction of $X_1$ and $X_2$ . While this technique allows us to generate interpolant objects from any two objects (of the same type), it is nevertheless clumsy to represent fractional interpolant objects via repeated subdivision. This subdivision technique also provides no obvious way of performing an average of many objects. What we would really like is some way of directly dealing with the linear interpolation $\alpha X_1 + (1-\alpha) X_2$ .

Creating a blade from a pure grade multivector

Consider the general interpolant, $X_{\alpha}' = \alpha X_1 + (1-\alpha)X_2$ where $X_1$ and $X_2$ are blades of the same grade. We claim that we can project $X_{\alpha}'$ into object space in a simple and general way. First we will generalise equation (3.2) to the interpolation case:

$\displaystyle X_\alpha = S(\alpha X_1 + (1-\alpha)X_2) = (\beta_\alpha + \gamma_\alpha(X_1X_2 + X_2X_1))(\alpha X_1 + (1-\alpha)X_2)$

(17)

where $\beta_\alpha$ and $\gamma_\alpha$ are once again scalars.

Since $S$ is of the form (scalar + 4-vector) it is self reverse. Defining $S^- = \langle S \rangle _0 - \langle S \rangle _4$ we get the result that $S^-S$ is a scalar, and can therefore write $X'_\alpha = k S^- X_\alpha$ where $k$ is a scalar and $k = \frac{1}{S^-S}$ .

To use this decomposition we need to extract $S$ from $X'_\alpha$ . To do this we can use the methods of Chapter 2, or as follows using the square root operator of Dorst and Valkenburg Dorst and Valkenburg [2011].

Let $SX'_\alpha = X_\alpha$ , where $X_\alpha$ is a valid object (squaring to $\pm$ 1). Now define $\Sigma \in Cl^{0,4}_{4,1}$ , ie. it only contains 0 and 4 grade coefficients and is an element of the conformal algebra. Then, defining $[[\Sigma]] = \sqrt{\langle{\Sigma}\rangle_0^2 - \langle{\Sigma}\rangle_4^2}$ , the square root can be found as:

$\displaystyle \sqrt{\Sigma} = \frac{\Sigma \pm [[\Sigma]]}{\sqrt{2}\sqrt{\langl... ...angle{\Sigma}\rangle_4}{\sqrt{2}\sqrt{\langle{\Sigma}\rangle_0 \pm [[\Sigma]]}}$

(18)

To use this method to find $S^-$ we multiply our non-blade object by its own reverse:

$\displaystyle X'_\alpha\tilde{X}'_\alpha = (kS^-X_\alpha )(kS^-X_\alpha )\,\til... ...X_\alpha \tilde{X}_\alpha(S^-)\,\tilde{} = -k^2S^-(S^-)\,\tilde{} = -k^2(S^-)^2$

(19)

This is now in a form where we can apply the above square root formula:

$\displaystyle kS^- = \sqrt{-X'_\alpha\tilde{X}'_\alpha}$

(20)

It now simply remains to isolate $X_\alpha$ via multiplication by $kS$

where $kS = \langle kS^- \rangle _0 - \langle kS^- \rangle _4$ . Since $(kS)(kS^-)$

is a scalar, we have

$\displaystyle X_\alpha = \frac{kS}{(kS)(kS^-)} X'_\alpha \equiv SX'_\alpha$

(21)

This result is particularly important as we have identified a way of projecting any pure grade object of the form $S^-X$ (with $X$ a blade) back to the blade manifold. An immediate application of this is that we can now deal with arbitrary linear combinations of objects, allowing us to smoothly interpolate as well as to average and cluster geometric primitives. Additionally we can correct numerical errors that result from arithmetic operations to give true blades again. Figure 3.1 shows examples of interpolating various geometric objects.

As shown in Chapter 2 this method holds for all the standard normalised conformal objects of grade 2 or above (point pairs, lines, circles, planes, spheres). The direct interpolation method is potentially more computationally efficient than the bivector interpolation method, and its form indicates that it is covariant, ie, for a rotor transformation given by $R$ ,

$\displaystyle R[\alpha X_1 + (1-\alpha)X_2]\tilde{R} = \alpha RX_1\tilde{R} + (1-\alpha)RX_2\tilde{R}$

Unfortunately, we cannot apply this form of interpolation to points as we encounter a problem due to the fact that for a conformal point $P$ , $P\tilde P = 0$ . However, we saw in equation (3.1) that points can be interpolated very easily using known explicit formulae.

**Fig. 3.1:** Linear interpolation between different geometric objects.
(a) circles, (b) point pairs, (c) planes, (d) spheres. The pure red and green objects here represent and and the intermediate colours show the interpolations between them. Here we are stepping linearly through $\alpha$ between 0 and 1 with the number of steps chosen to show the interpolations as clearly as possible.
$\begin{subfigure}[b]{0.4\textwidth} \includegraphics[width=\textwidth]{circles.png} \end{subfigure}$ $\begin{subfigure}[b]{0.4\textwidth} \includegraphics[width=\textwidth]{pairPair.png} \end{subfigure}$ $\begin{subfigure}[b]{0.4\textwidth} \includegraphics[width=\textwidth]{planes.png} \end{subfigure}$ $\begin{subfigure}[b]{0.4\textwidth} \includegraphics[width=\textwidth]{spheres.png} \end{subfigure}$

Techniques for understanding interpolant properties

In order to use the interpolant blades it is useful to get a handle on some of their properties. In several cases it is possible to get good insight into how the interpolant behaves by looking at the interpolant of the dual of the blades, but in others we need to consider the form of the (scalar + 4-vector) required to project the interpolant back to the blade manifold. As before we write our blades as:

$\displaystyle X_\alpha = (\langle S\rangle _0 + \langle S\rangle _4)(\alpha X_1 + (1-\alpha) X_2)$

From this we immediately see that for the multiplication to be grade preserving we require $\langle S\rangle _4X_1$ and $\langle S\rangle _4 X_2$ to give only objects of grade $n$

where

is the grade of $X_1$

and

. Table 3.1 shows the resultant grades from the geometric product of pure grade objects and Table 3.2 shows the resultant grades from the inner product. These tables are presented here for reference and will be returned to when dealing with individual grade blades.

Table 3.1: Resulting grades from the geometric product of pure grade objects. The grade 4 row is highlighted here as it is of specific interest for the form of the $S$

in equation (3.3).

	0	1	2	3	4	5
0	0	1	2	3	4	5
1	1	0,2	1,3	2,4	3,5	4
2	2	1,3	0,2,4	1,3,5	2,4	3
3	3	2,4	1,3,5	0,2,4	1,3	2
Gray 4	4	3,5	2,4	1,3	0,2	1
5	5	4	3	2	1	0

Table 3.2: Resulting grades from the inner product of pure grade objects. In the case of the inner product of a multivector and a scalar by definition the result is always 0 rather than some other scalar valued function of the scalar and multivector. Here again the grade 4 row is highlighted here as it is of specific interest for the form of the $S$

in equation (3.3).

	1	2	3	4	5
0	0	0	0	0	0
1	0	1	2	3	4
2	1	0	1	2	3
3	2	1	0	1	2
Gray 4	3	2	1	0	1
5	4	3	2	1	0

Point pairs

We start with point pairs. Previous work Dorst et al. [2007] has shown that when an end point is shared between point pairs $A$ and $B$ the interpolant point pairs are also blades and their end points trace out the circumference of the circle formed by the shared point and the additional separate end points. Three points $X,Y,Z$ define a circle $% latex2html id marker 16714 $ C \propto X\wedge Y \wedge Z$$ and a fourth point $V$ lying on the circle will satisfy $V\wedge C = 0$ , this allows us to define a check to see if two point pairs are chords of the same circle. Point pairs $A = V\wedge X$ , $B = Y\wedge Z$ will satisfy $A\wedge B = 0$ if they are both chords and any additional chord, $W$ , of the same circle will satisfy $W\wedge A = 0, W\wedge B = 0$ and thus $W\wedge (A + B) = 0$ . This leads to:

Theorem 1 If point pairs and are both chords of a common circle the interpolant point pairs $% latex2html id marker 16743 $ T \propto \alpha A + (1-\alpha)B$$ are blades and also have end points lying on as $\langle AB\rangle _4=A\wedge B=0$ and $(\alpha A + (1-\alpha)B)\wedge(A + B) = 0$ .

Note, since $\langle AB\rangle _4=A\wedge B=0$ , the projector $S$ is a scalar. The common circle itself is the `join' of the two original point pairs and can be computed with the algorithms supplied in Chapter 21 of Dorst, Fontijne and Mann Dorst et al. [2007]. Figure 3.2 shows two cases of the interpolation of co-planar point pairs that lie on the same circle.

**Fig. 3.2:** The interpolation of point pairs with endpoints lying in the same plane and on a common circle is a blade and also lies on the same common circle, even in cases in which there are no shared endpoints or intersections. In this figure the red lines are the interpolation of the black lines.

**Fig. 3.3:** The interpolation in red of point pairs and , here shown in black, lie on the surface of the sphere $% latex2html id marker 16764 $ \propto A\wedge B$$ , shown in blue with black equator circle.

Turning to the more general case of two point pairs in arbitrary positions in space we can get insight into the form of the interpolant by considering the components of the scalar + 4-vector projection factor. In the case of the geometric product between grade 4 and grade 2 objects we see from Table 3.1 that we produce both 2 and 4-vector grades. The 2-vector part of the geometric product comes from the inner product between the point pairs and the 4-vector. ie. for point pairs $A$ and $B$ , $\langle S\rangle _4(\alpha A + (1- \alpha)B) = \langle S\rangle _4\cdot(\alpha A + (1- \alpha)B)$ . For the general case of two point pairs not lying in plane ie. $A \wedge B \ne 0$ , we can show that there is only one object that behaves in this way, the sphere $% latex2html id marker 16781 $ \Sigma \propto A\wedge B$$ , as it passes through both end points of both point pairs. This is illustrated in Figure 3.3 and suggests that the sphere $\Sigma$ is intrinsically tied to the form of the interpolant objects. Indeed we can see from the same visualisation that the interpolant $C$ of point pairs $A$ and $B$ always has endpoints lying on the surface of the sphere $\Sigma$ .

We can prove this by showing that $C\wedge B$ or $C\wedge A$ also produces the sphere:

$\displaystyle \Sigma$

$% latex2html id marker 16797 $\displaystyle \propto A\wedge B$$

First consider an interpolant object C and its outer product with one of the original objects, $B$

$\displaystyle C$	$\displaystyle = (\langle S\rangle _0 + \langle S\rangle _4)(\alpha A + (1-\alpha)B)$
$\displaystyle C\wedge B$	$\displaystyle = ((\langle S\rangle _0 + \langle S\rangle _4)(\alpha A + (1-\alpha)B))\wedge B$
	$\displaystyle = \langle S\rangle _0(\alpha A + (1-\alpha)B)\wedge B + (\langle S\rangle _4(\alpha A + (1-\alpha)B))\wedge B$

as $B\wedge B = 0$ we see that

$\displaystyle \alpha\langle S\rangle _0 A\wedge B = \langle S\rangle _0(\alpha A + (1-\alpha)B)\wedge B$

Now we just need to prove that $\langle S\rangle _4((\alpha A + (1-\alpha)B)\wedge B)$ is a scalar multiple of $A\wedge B$ . From equation (3.3) we know that

$% latex2html id marker 16813 $\displaystyle \langle S\rangle _4 \propto \langle AB + BA \rangle _4$$

it is therefore sufficient to prove that:

$\displaystyle T$

$% latex2html id marker 16815 $\displaystyle = (\langle AB + BA \rangle _4(\alpha A + (1-\alpha)B))\wedge B \propto A\wedge B$$

We can convert the outer product into a geometric product followed by a projection and thus can write:

$\displaystyle T$	$\displaystyle = \langle \langle AB + BA \rangle _4(\alpha A + (1-\alpha)B)B\rangle _4$
	$\displaystyle = \langle \langle AB + BA \rangle _4(\alpha AB + (1-\alpha))\rangle _4$
	$\displaystyle = (1-\alpha)(\langle AB\rangle _4 + \langle BA \rangle _4) + \alpha\langle (\langle AB\rangle _4 + \langle BA \rangle _4)AB\rangle _4$

As $\langle AB\rangle _4 \equiv \langle BA\rangle _4$ we can write this as:

$\displaystyle T = 2 (1-\alpha)\langle AB\rangle _4 + 2\alpha\langle \langle AB\rangle _4 AB\rangle _4$

and this can further simplified using the fact that:

$\displaystyle \langle \langle AB\rangle _4 AB\rangle _4 = \langle \langle AB\rangle _4 \langle AB\rangle _0 \rangle _4$

As $\langle AB\rangle _0 = A\cdot B$ is a scalar:

$\displaystyle T$	$\displaystyle = 2(1-\alpha)\langle AB\rangle _4 + 2\alpha\langle AB\rangle _4 (A \cdot B)$
	$\displaystyle = 2((1-\alpha) + \alpha A \cdot B)\langle AB\rangle _4$

As $2((1-\alpha) + \alpha A \cdot B)$ is a scalar we see that the proof is completed

$% latex2html id marker 16832 $\displaystyle 2((1-\alpha) + \alpha A \cdot B)\langle AB\rangle _4 \propto \langle AB\rangle _4 \equiv A\wedge B$$

Figure 3.3 shows a graphical representation of the interpolant point pairs lying on the surface of the sphere. To summarise:

Theorem 2 For non-coplanar point pairs and , all interpolant point pairs lie on the surface of the sphere $% latex2html id marker 16843 $ \Sigma \propto A\wedge B$$ .

Circles

The interpolant of circles has a range of properties that are useful and clearly intrinsically tied to the geometry of spheres and point pairs. Initially we will consider the case of two circles in space that both lie on the surface of a common sphere. In past work it has been shown that circles with two common points interpolate directly without requiring re-projection and the interpolant lies on their common sphere Doran [2003] Dorst et al. [2007]. Here, as with the point pairs, we can show that this is true for a broader class of circles:

Theorem 3 If circles and together define the caps of a common sphere then $% latex2html id marker 16855 $ \langle S\rangle _4 \propto \langle C_1C_2 \rangle _4 = 0$$ where is of the form shown in equation (3.3) and thus any interpolant object $C_3 = \alpha C_1 + (1-\alpha)C_2$ is a blade without requiring re-projection to the blade manifold.

This can be proved by considering each circle $C_i$

as the intersection of a plane $P_i$

and a sphere $\Sigma_i$ . Forming this intersection via the dual (where $X^* = XI_5$

and

is the 5D space pseudoscalar as before), we have:

$\displaystyle C_1$	$\displaystyle = (\Sigma_1^\wedge P_1^)I_5$
$\displaystyle C_2$	$\displaystyle = (\Sigma_2^\wedge P_2^)I_5$
$\displaystyle \langle C_1C_2\rangle _4$	$\displaystyle = -\langle (\Sigma_1^\wedge P_1^)(\Sigma_2^\wedge P_2^)\rangle _4$

Since $(\Sigma_1^*\wedge P_1^*)$ and $(\Sigma_2^*\wedge P_2^*)$ are both bivectors:

$\displaystyle \langle C_1C_2\rangle _4 = -\Sigma_1^*\wedge P_1^*\wedge \Sigma_2^*\wedge P_2^*$

and so if $% latex2html id marker 16882 $ \Sigma_2 \propto \Sigma_1$$ :

$% latex2html id marker 16883 $\displaystyle \langle C_1C_2\rangle _4 \propto -\Sigma_1^*\wedge P_1^*\wedge \Sigma_1^*\wedge P_2^* = 0$$

We can additionally find the unique common sphere by finding the join of the circles or by reverting to linear algebra techniques:

Conjecture 1 If circles and together define the caps of a common sphere $\Sigma$ then $\langle C_1 \Sigma \rangle _3, \langle C_2 \Sigma \rangle _3 = 0$ . $\Sigma$ can be found by the following process:

First we define:

$\displaystyle \underaccent{\bar}{A}= \begin{bmatrix} \underaccent{\bar}{M}_3... ...nt{\bar}{C}_1\\ \underaccent{\bar}{M}_3\underaccent{\bar}{C}_2 \end{bmatrix}$

where is the truncated identity matrix that performs selection of grade 3 elements from a vector of coefficients and $\underaccent{\bar}{C}_1$ and $\underaccent{\bar}{C}_2$ are the matrices that perform the left geometric product of and respectively with a vector of coefficients. We can then find the $\Sigma$ for $\underaccent{\bar}{A}\underaccent{\bar}{\Sigma} = 0$ where $\underaccent{\bar}{\Sigma}$ is a vector of canonical blade coefficients limited to only the 4-vector blades. In the case that and are the same radius then $% latex2html id marker 16922 $ \Sigma \propto (C_1 + C_2)((C_1 + C_2)\wedge n_{\infty})I_5$$ .

The case for circles of the same radius is visualised in Figure 3.4.

It is also the case that the interpolant lies on the surface of the common sphere:

Theorem 4 If circles and together define the caps of a common sphere then all interpolant circles $C_3 = \alpha C_1 + (1-\alpha)C_2$ (which we have shown to be blades) also lie on the surface of the sphere $\Sigma$ common to both.

We can prove this by considering the outer product of the interpolant circle with $D$

, an arbitrary point on the common sphere $\Sigma$ :

$% latex2html id marker 16940 $\displaystyle \Sigma \propto D\wedge C_3 = \alpha D\wedge C_1 + (1-\alpha)D\wedge C_2$$

Figure 3.5 shows an example of this interpolation.

**Fig. 3.4:** The half way circle shown in red in this figure is the equator of the sphere through both and if they define a common sphere and have the same radius

**Fig. 3.5:** The interpolation of circles and is a blade and lies on the surface of a sphere if and define a common sphere.

Thus far we have dealt exclusively with circles on a common sphere. In the case in which $C_1$ and $C_2$ do not lie on the same sphere we can again look at how the interpolants behave by considering the form of the (scalar + 4-vector) that we use to project the interpolant back to the blade manifold. In the case of the geometric product between grade 4 and 3 objects we see from Table 3.1 that we produce both 1 and 3-vector grades, however the 1-vector part of the geometric product comes only from the inner product between the 4-vector and the circles. To maintain grade after the multiplication the 4-vector must therefore be the object that has an inner product of zero with both circles. This object is the sphere into whose surface both circles plunge orthogonally Dorst et al. [2007]:

Theorem 5 If circles and together do not lie on a common sphere then the 4-vector from our blade projection equation $% latex2html id marker 16989 $ \langle S\rangle _4 \propto \langle C_1C_2 \rangle _4$$ is itself a blade and geometrically represents the sphere through which both circles plunge orthogonally. ie. $C_1\cdot \langle C_1C_2 \rangle _4 = 0$ . This property means all interpolant circles after projection to the blade manifold also plunge through $\langle C_1C_2 \rangle _4$ orthogonally. ie. $C_3\cdot \langle C_1C_2 \rangle _4 = 0$ .

The intersections of the interpolant circles with the sphere $\langle S\rangle _4$ produce a set of point pairs. Intuition would suggest that these point pairs have properties tied to the interpolation of the point pairs generated by the original two circles $C_1$ and $C_2$ and indeed we can numerically verify that this is the case:

Conjecture 2 If circles and together do not lie on a common sphere then the intersection point pair $P_\alpha$ formed by the meet of the circle interpolant for a given value of $\alpha$ with the orthogonal sphere $\langle C_1C_2 \rangle _4$ ie. $% latex2html id marker 17018 $ P_\alpha \propto C_\alpha\vee\langle C_1C_2 \rangle _4$$ is the same as the re-projected interpolant $\Pi_\alpha$ of the point pairs formed from the meet of and with $\langle C_1C_2 \rangle _4$ .

Figure 3.6 shows the interpolation of two non co-spherical circles as well as the sphere these circles define and the intersection point pairs they generate.

**Fig. 3.6:** The interpolation of the two black circles which are not spherical caps, intersect orthogonally with a single sphere (shown in blue). The point pairs formed from the two intersection points are shown in yellow.

Lines

When looking at lines we can attempt to use some of the same techniques that we used for circles. First consider the form of $% latex2html id marker 17032 $ \langle S\rangle _4 \propto \langle L_1L_2\rangle _4$$ . For lines $% latex2html id marker 17034 $ \langle L_1L_2\rangle _4 \propto I_5n_{\infty}$$ , giving the form of the projection of $(X_1 + X_2)$ as:

$\displaystyle X_3 = (\mu + \nu I_5n_{\infty})(X_1 + X_2)$

(22)

where $\mu$ and $\nu$ are scalars. While neat, this form $I_5n_{\infty}$ does not on its own provide information on the properties of the interpolated line. Instead we consider the interpolation of the dual of the lines, and to understand this interpolation we must take a short detour via screw theory.

Screw Theory

Screw theory was developed by Sir Robert Stawell Ball in 1900 in his seminal work `A treatise on the theory of screws' Ball [1900]. His original applications were kinematics and one of the most important theorems in the area, Chasles' theorem, states that the most general rigid body displacement can be described by a screw transformation. More recently screw theory, and the highly related study of dual quaternions, has been applied to robotics, computational geometry and multibody dynamics K. Davidson et al. [2004]; Kavan et al. [2008]; Müller [2018].

Screw transformations consist of a translation along an axis and a rotation around that axis. To parameterise a screw we define the direction of the screw axis via a unit vector $\hat{m}$ , a point on the screw axis $p$ and a screw pitch $h$ . The pitch represents how far to move in the direction of the screw axis for each complete revolution about the axis.

Bivector representation of a line

A line in CGA is represented as a 3-vector, or dually as a bivector:

$\displaystyle L^* = \hat{m} I_3 + (p\wedge \hat{m})I_3 n_{\infty}$

(23)

This bivector formulation is equivalent to the Plücker coordinates of the line.

In Dorst and Valkenburg [2011] the authors describe the orbit of simple bivectors that describe motion. We can visualise the orbit of the dual line bivector by exponentiating the bivector to a rotor and applying it to a test point. Figure 3.7 shows the orbit of the point at the origin about a line. The motion is a circle about the line.

**Fig. 3.7:** The orbit of a line

The bivector representation of a screw

To represent a screw we will couple the rotational motion of the dual of a line with a translation in the direction of that same line. The bivector $T$

that transforms along the 3d vector $t$

is:

$\displaystyle T = tn_{\infty}$

(24)

If $t = h \hat{m}$ where $h$

is a scalar, ie. the translation is in the screw axis direction, the rotors formed from the bivectors in equations (3.9) and (3.10) commute.

It then follows that the rotor formed from the addition of the bivectors in equations (3.9) and (3.10) can be split into the rotor representing translation along the axis and the rotor representing rotation about the axis – as required for a screw. We therefore have a screw, $\mathbbm{S}$ , whose action on the point at the origin is shown in Figure 3.8.

$\displaystyle \mathbbm{S} = \hat{m} I_3 + (p\wedge \hat{m})I_3 n_{\infty}+ h\hat{m}n_{\infty}$

(25)

**Fig. 3.8:** The action of a basic screw formed by the summation of commuting bivectors applied to the point at the origin is shown in black, it forms a screw motion about the screw axis shown in red.

Hestenes and Sobczyk Hestenes et al. [1985] p81 gives an expression for decomposing any bivector into two commuting blades. In the case of our screw bivector these blades represent the dual of the screw axis $L^*$

and a translational bivector $T$

in the direction of the screw axis. ie. given a screw bivector $\mathbbm{S}$ we can decompose it as:

$\displaystyle \mathbbm{S} = L^* + T$

Adding dual lines

The addition of dual lines produces a bivector. Visualising the action of this bivector allows us to see that it is in fact also a screw transformation. Consider the addition of two dual lines:

$\displaystyle L^*_+ = L_1^* + L_2^*$

we can write this elementwise as

$\displaystyle L^*_+ = m_1I_3 + \psi_1 n_{\infty}+ m_2I_3 + \psi_2 n_{\infty}$

where $\psi_i$ is $(p_i\wedge m_i)I_3$ . We then rearrange to give something proportional to the expression in equation (3.11):

$\displaystyle L^*_+$

$\displaystyle = (m_1+ m_2 )I_3 + (\psi_1 + \psi_2 )n_{\infty}= m I_3 + (p\wedge m)I_3 n_{\infty}+ hmn_{\infty}$

Where clearly $m = m_1 + m_2$

. If we divide this by $\vert m\vert$ we have the general form of a normalised screw

$\displaystyle \mathbbm{S} = \hat{m} I_3 + (p\wedge \hat{m})I_3 n_{\infty}+ h\hat{m}n_{\infty}$

Gathering like terms, specifically those without an $n_{\infty}$ component, leads us to the conclusion that our screw axis direction $\hat{m}$ must simply be proportional to the addition of the directions of the two lines. Using this fixed axis direction we can extract the coefficient $h$

(the pitch) of the translation bivector parallel to the screw axis:

$\displaystyle m$	$\displaystyle = m_1 + m_2$
$\displaystyle L_+^*$	$\displaystyle = m I_3 + (p\wedge m)I_3 n_{\infty}+ hmn_{\infty}$
$\displaystyle L_+^*\cdot n_0$	$\displaystyle = hm + (p\wedge m)I_3$
$\displaystyle h$	$\displaystyle = \frac{(L_+^*\cdot n_0) \cdot m}{\vert m\vert^2}$

With this coefficient known we now have all the pieces in place for a full decomposition of the dual line addition bivector $L_+^*$

$\displaystyle L_+^*$	$\displaystyle = L_1^* + L_2^*$
$\displaystyle m$	$\displaystyle = -\langle L_+^* \rangle _{(e_{12},e_{13},e_{23})}I_3$
$\displaystyle h$	$\displaystyle = \left(\frac{(L_+^*\cdot n_0) \cdot m}{\vert m\vert^2}\right)$
$\displaystyle \mathbbm{S}$	$\displaystyle = L^* + T$
$\displaystyle T$	$\displaystyle = h\,\hat{m}n_{\infty}$
$\displaystyle L^*$	$\displaystyle = \mathbbm{S} - T$

Figure 3.9 shows the decomposition of the addition of two lines into its component parts.

**Fig. 3.9:** The addition of two dual lines, here shown in black, produces a screw, whose action on the point at the origin is again shown in black. The screw can be decomposed into two commuting bivectors, whose actions are shown in red and blue. The bivector whose action is shown in red is the dual of the screw axis line, also shown in red.

Relationship to object manifold reprojection

We can also analyse the screw multiplied by its own reverse, comparing this formulation with our object manifold reprojection to get the projection coefficient $S$

in terms of the screw parameters:

$\displaystyle {\mathbbm{S}}\tilde{\mathbbm{S}}$	$\displaystyle = -(L^)^2 + T\tilde{T} + T\tilde{L^} + L^*\tilde{T}$
	$\displaystyle = 1 - 2TL^*$
	$\displaystyle = 1 + 2hn_{\infty}I_3 = 1 + 2hI_5n_{\infty}$
	$\displaystyle \equiv (kS^-)^2$
$\displaystyle (kS^-)$	$\displaystyle = \mu + \nu I_5 n_{\infty}$
$\displaystyle (kS^-)^2$	$\displaystyle = \mu^2 + 2\mu\nu I_5 n_{\infty}$
$\displaystyle \mu$	$\displaystyle = 1, \,\,\, \nu = h$

Using this calculated value we can see how the projection coefficient acts on the addition of lines:

$% latex2html id marker 17145 $\displaystyle L^* \propto$$	$\displaystyle (1 - h I_5 n_{\infty})(m I_3 + (p\wedge m)I_3 n_{\infty}+ hmn_{\infty})$
	$\displaystyle = m I_3 + (p\wedge m)I_3 n_{\infty}+ hmn_{\infty}- hmn_{\infty}$
	$\displaystyle = m I_3 + (p\wedge m)I_3 n_{\infty}$

This is in fact the same line as is formed from the decomposition of the screw bivector into the screw axis bivector and pitch translation bivector. In other words, the addition of lines and reprojection to the line manifold extracts the axis of the screw formed from the addition of their duals. This axis has a direction equal to a linear interpolation of the axes of the original two lines and, as it is a mirror object, passes through the point exactly half way between the lines.

Planes

All 4-vectors are blades. Thus, for planes and spheres it is impossible to construct an invalid geometric object by addition. For planes we can analyse the form of the interpolant by again looking at the dual of a plane.

The dual of the plane can be written as:

$\displaystyle P^* = \hat{m} + dn_{\infty}$

where $\hat{m}$ is the 3D vector normal to the plane and $d$

is the perpendicular distance of the plane from the origin. Thus the interpolation of duals of two planes can be written as:

$\displaystyle \alpha P_1^* + (1-\alpha)P_2^*= \alpha\hat{m_1} + \alpha d_1n_{\infty}+ (1-\alpha)\hat{m_2} + (1-\alpha)d_2n_{\infty}$

which, when we collect like terms, is already in the form of a dual plane $P_3^*$

$\displaystyle P_3^* = \alpha \hat{m_1} + (1-\alpha)\hat{m_2} + (\alpha d_1 + (1-\alpha)d_2)n_{\infty}$

(26)

this dual plane has a normal vector that is the interpolation of the normal vectors of the original two planes and has a perpendicular distance from the origin that is also simply an interpolation of the perpendicular distance from the origin of the original two planes. An important feature of this plane interpolation is that, as noted in Cameron [2007], provided the two planes intersect, the interpolant plane always passes through the line of intersection (the meet) of the two original planes. This is visualised in Figure 3.10. In the case that the planes to not intersect (or more formally are said to intersect at infinity) the interpolation will smoothly translate one plane to the other keeping the normal fixed, the parallel vs anti-parallel cases are explored in Figures 3.11 and 3.12.

**Fig. 3.10:** The interpolant of two planes (green to red) always passes through the meet line (black) of the two original planes.

**Fig. 3.11:** The interpolant of two parallel planes smoothly moves between the start and end points (green to red) while maintaining the direction of the normal.
**Fig. 3.12:** The interpolant of two anti-parallel planes (green to red) must go via infinity due to the sign change. Care must be taken with the orientation of objects when designing algorithms using these interpolations.

Spheres

The interpolant of spheres has been studied before in Cameron [2007] and Doran [2003]. As with planes, all interpolants of spheres are valid objects as $\langle \Sigma_1\Sigma_2\rangle _4 = 0$ and have the property of making contact with the meet of the spheres at all points during the interpolation. We can see the form of the interpolant sphere by considering its dual form:

$\displaystyle I_5\Sigma_3 = \alpha I_5\Sigma_1 + (1-\alpha)I_5\Sigma_2$

The dual form of a sphere can be decomposed into the sum of the conformal centre point $P$

and negative half the radius squared times $n_{\infty}$ :

$\displaystyle I_5\Sigma = P - \frac{1}{2}\rho^2n_{\infty}$

the interpolation of the dual of two spheres is therefore

$\displaystyle I_5\Sigma_3 = \alpha P_1 + (1-\alpha)P_2 - \frac{1}{2}(\alpha\rho_1^2 + (1-\alpha)\rho_2^2)n_{\infty}$

For two concentric spheres ie. $P_2 = P_1$

we can therefore see that the interpolation between them will remain centred in the same place and will simply have a radius $\rho_3$ which varies as $\rho_3^2 = \alpha\rho_1^2 + (1-\alpha)\rho_2^2$ .

As we have seen previously, the interpolation of two conformal points $P_1$ and $P_2$ is of the form

$\displaystyle \alpha P_1 + (1-\alpha)P_2 = F(\alpha p_1 + (1-\alpha)p_2) - \alpha(1-\alpha)(P_1\cdot P_2)n_{\infty}$

we can therefore also write the interpolation of two non-concentric spheres as:

$\displaystyle I_5\Sigma_3 = F(\alpha p_1 + (1-\alpha)p_2) - \alpha(1-\alpha)(P_... ...ot P_2)n_{\infty}- \frac{1}{2}(\alpha\rho_1^2 + (1-\alpha)\rho_2^2 )n_{\infty}$

Collecting like factors shows that the centre point of the interpolated sphere moves linearly along the line joining $p_1$

and

$\displaystyle I_5\Sigma_3 = F(\alpha p_1 + (1-\alpha)p_2) -\left( \alpha(1-\alp... ...cdot P_2) + \frac{1}{2}(\alpha\rho_1^2 + (1-\alpha)\rho_2^2) \right)n_{\infty}$

Furthermore, writing the dot product of points in terms of their euclidean vectors we can see that the radius of the sphere varies along its interpolation path

$\displaystyle I_5\Sigma_3 = F(\alpha p_1 + (1-\alpha)p_2) -\frac{1}{2}\left( -\... ...a(1-\alpha)(p_1-p_2)^2 + \alpha\rho_1^2 + (1-\alpha)\rho_2^2 \right)n_{\infty}$

and so the radius $\rho_3$ varies as

$\displaystyle \rho_3^2 = -\alpha(1-\alpha)(p_1-p_2)^2 + \alpha\rho_1^2 + (1-\alpha)\rho_2^2$

For fixed values of $\rho_1$ and $\rho_2$ this implies $\rho_3^2$ varies as $-(p_1 - p_2)^2$

and so the further apart the two spheres are the smaller the radius of the interpolant. To find turning points we differentiate with respect to $\alpha$

$\displaystyle \frac{d\rho_3^2}{d\alpha} = (-1+2\alpha)(p_1-p_2)^2 + \rho_1^2 -\rho_2^2$

setting this to zero yields a single turning point at

$\displaystyle \alpha = \frac{\rho_2^2 -\rho_1^2 + (p_1-p_2)^2}{2(p_1-p_2)^2}$

Considering the second derivative

$\displaystyle \frac{d^2\rho_3^2}{d\alpha^2} = 2(p_1-p_2)^2$

we see that this is always positive and so the stationary point is a minimum.

For the case that the surfaces of $\Sigma_1$ and $\Sigma_2$ are just touching we have the condition

$\displaystyle (p_1 - p_2)^2 = (\rho_1 + \rho_2)^2 = \rho_1^2 + \rho_2^2 + 2\rho_1\rho_2$

returning to the first derivative in this case

$\displaystyle \alpha = \frac{2\rho_2(\rho_1 + \rho_2)}{2(\rho_1 + \rho_2)^2} = \frac{\rho_2}{(\rho_1 + \rho_2)}$

this value of $\alpha$ is the point at which the centre of the interpolant sphere lies on the surface of both spheres. At this point the squared radius is zero:

$\displaystyle \rho_3^2 = -\frac{\rho_2}{(\rho_1 + \rho_2)}\frac{\rho_1}{(\rho_1... ...{\rho_2}{(\rho_1 + \rho_2)}\rho_1^2 + \frac{\rho_1}{(\rho_1 + \rho_2)}\rho_2^2$

$\displaystyle \rho_3^2 = \frac{-\rho_2\rho_1(\rho_1 + \rho_2) + \rho_2\rho_1^2 + \rho_1\rho_2^2}{(\rho_1 + \rho_2)} = 0$

Pulling the spheres further apart from this point so that they no longer intersect will therefore produce a sphere with negative radius, an imaginary sphere. These results are already known Cameron [2007] and are here included for completeness.

Applications

The ability to interpolate geometric objects suggests a wide variety of applications in the areas of computer vision and graphics. There are many traditional algorithms in vision that rely solely on point information from images and ignore lines and other, potentially useful, geometric primitives. Many of these algorithms have been non trivial to translate into the framework of CGA due to having to specify transformations between objects explicitly rather than implicitly via the objects themselves. The ability to average geometric objects directly suggests immediate applications in clustering of objects extracted from real data, interpolation to produce surfaces and other areas for problems where we might normally use linear algebra.

Higher order spline interpolation through objects

With the ability to construct arbitrary linear combinations of blades we naturally might wonder about the applications of this to spline generation through control objects. Figure 3.13 shows an example of interpolating through different control objects with different orders of spline. As expected, higher order interpolation produces smoother surfaces through our objects.

**Fig. 3.13:** Interpolation through control objects.
Top: Circles. Bottom: Point pairs.
Interpolation type: (a,d) Linear, (b,e) Quadratic, (c,f) Cubic
$\begin{subfigure}[b]{0.3\textwidth} \includegraphics[width=\textwidth]{circle_spline_1.png} \end{subfigure}$ $\begin{subfigure}[b]{0.3\textwidth} \includegraphics[width=\textwidth]{circle_spline_2.png} \end{subfigure}$ $\begin{subfigure}[b]{0.3\textwidth} \includegraphics[width=\textwidth]{circle_spline_3.png} \end{subfigure}$ $\begin{subfigure}[b]{0.3\textwidth} \includegraphics[width=\textwidth]{point_pair_spline_1.png} \end{subfigure}$ $\begin{subfigure}[b]{0.3\textwidth} \includegraphics[width=\textwidth]{point_pair_spline_2.png} \end{subfigure}$ $\begin{subfigure}[b]{0.3\textwidth} \includegraphics[width=\textwidth]{point_pair_spline_3.png} \end{subfigure}$

Recursive scene simplification by averaging conformal objects

**Fig. 3.14:** A 3d line model before (above) and after (below) recursive scene simplification.
$\begin{subfigure}[b]{1.0\textwidth}{ \includegraphics[width=\textwidth]{before_filter.png} } \end{subfigure}$ $\begin{subfigure}[b]{1.0\textwidth}{ \includegraphics[width=\textwidth]{after_filter.png} } \end{subfigure}$

When extracting geometric primitives from triangulated CAD models from point cloud data or from images, there are often many objects that lie close to each other in space. Line segment detectors, for example, will often extract long lines as multiple line segments that need stitching together. We would like a way of simplifying these noisy models by collapsing objects that are close together into a single object. One way to do this is via a recursive filtering algorithm as follows:

Set a minimum cost threshold for difference between objects
Compute the cost between all objects of the same grade in the scene
If all costs are above the threshold then terminate the algorithm
Average the two objects with the smallest cost
Return to step (2)

This leads to a simplified model that retains the core features of the original model. For comparison of objects $X_i$ and $X_j$ we use the cost function $C_{ij}$ for a rotor $R_{ij}$ as defined in Eide and Lasenby [2018]:

$\displaystyle C_{ij} = \langle (R_{ij} - 1)(\tilde{R}_{ij} - 1){\rangle }_0 + \langle R_{ij\parallel} \tilde{R}_{ij\parallel}{\rangle }_0$

(27)

where $R_{\parallel}=R\cdot e$ , and gives the component of $R$

having $n_{\infty}$ as a factor and $R_{ij}$ is the rotor that takes $X_i$

as described in Lasenby et al. [2018]. An example of this algorithm working on simulated lines is shown in Figure 3.14.

This algorithm is simply one way to perform scene simplification and it has a high computational complexity making it run slowly for large numbers of objects, but is included here as an example of one potential area the averaging of object methodology may be applied to.

-means clustering of conformal objects

One of the most fundamental and simple clustering algorithms is known as $k$ -means clustering Duda et al. [2001]. Consider a 3d scene composed of $k$ geometric objects of a given grade. We have multiple noisy observations for each object and so would like to fit $k$ centroids to these clusters to represent the “true” objects in the world.

**Fig. 3.15:** Three clusters of 3d lines correctly segmented by the algorithm.

The steps for implementing this clustering are given below:

Randomly assign objects to be the initial positions of the cluster centroids, leave all other objects unassigned
Assign each object in the scene to the centroid closest under our given cost metric, again we use the cost function given in equation (3.13)
If this is not our first iteration and no objects have changed assignment then terminate the algorithm
The centroid of each cluster is moved to the mean of the objects assigned to it, where mean is defined as the sum of the objects in the cluster projected back onto the blade manifold
Go to step (2)

Figures 3.15 and 3.16 show the successful application of this algorithm on simulated data – each line or circle has been associated with the cluster (indicated by colour) to which it is most likely to belong. One of the key advantages of using the averaging of objects and correction back to a blade for this algorithm is that it is computationally cheap. A typical approach in GA to this kind of problem might involve attempting to find the mean of a given cluster by optimisation of our cost function through a space parameterising our centroid objects. Here we can simply average the objects in each cluster making it feasible to cluster very large numbers of conformal objects quickly.

**Fig. 3.16:** Three clusters of 3d circles correctly segmented by the algorithm. The black circles here are the final computed cluster centroids.

Closest point to two non intersecting lines (least squares sense)

Consider two non-intersecting non-coplanar lines in 3d space, $L_1$

and

. We wish to find the point $P$

that lies closest to both in a least squares sense. First we will construct two orthogonal intermediary lines $L_+ = S_+(L_1 + L_2)$

and

where

represents the projection of a 3-vector $X$

back onto the line manifold. $L_+$

and

both lie half way between the two original skew lines but intersect at right angles . The intersection of these lines is the point $P$

that lies half way between the original lines. To extract this point of intersection we can follow the formula given in Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]:

$\displaystyle Q = (L_-n_0L_-)n_{\infty}(L_-n_0L_-) + L_+(L_-n_0L_-)n_{\infty}(L_-n_0L_-)L_+$

$\displaystyle P = \frac{-Q n_{\infty}Q}{(Q n_{\infty}Q)\cdot n_{\infty}}$

Conclusions

This chapter has shown how we are able to add multiples of conformal objects by factoring the resulting multivector into a scalar plus 4-vector term and a valid geometric object. We have then investigated the form of this multivector for each grade of conformal object. Using the ideas of interpolating and averaging objects, a range of applications are suggested with relevance in computer vision and computer graphics.

Exploring Novel Surface Representations via an Experimental Ray-Tracer in CGA

Arithmetic! Algebra! Geometry! Grandiose trinity! Luminous triangle! Whoever has not known you is without sense!Comte de Lautreamont

Abstract

Conformal Geometric Algebra (CGA) provides a unified representation of both geometric primitives and conformal transformations, and as such holds significant promise in the field of computer graphics. In this chapter we implement a simple ray tracer in CGA with a Blinn-Phong lighting model, before putting it to use to examine ray intersections with surfaces generated from the direct interpolation of geometric primitives. General surfaces formed from these interpolations are rendered using analytic normals. In addition, special cases of point-pair interpolation, which might find use in graphics applications, are described and rendered. A closed form expression is found for the derivative of the square root of a scalar plus 4-vector element with respect to a scalar parameter. This square root derivative is used to construct an expression for the derivative of a pure-grade multivector projected to the blade manifold. The blade manifold projection provides an analytical method for finding the normal line to the interpolated surfaces and its use is shown in lighting calculations for the ray tracer and in generating vertex normals for exporting the evolved surfaces as polygonal meshes.

**Fig. 4.1:** Three images rendered with varying lighting positions. These demonstrate the lighting model, multiple light source capability and recursive tracing of rays for reflections.
$\includegraphics[width=\textwidth]{FirstFig.eps}$

Introduction

Tubular and ribbon surfaces have wide interest in fields such as neuronal modelling and streamline visualisation. The need to represent vast networks of tubular data efficiently and render these surfaces in a visually pleasing way has led to a range of different parametric representations, fitting methods and rendering techniques Bauer and Polthier [2007]; Peternell and Pottmann [1997]; Petrovic et al. [2007]. Conformal Geometric Algebra (CGA) encodes circles and line-segments, as well as planes, spheres, infinite lines and the geometric transformations between them, as natural elements of an algebra Dorst et al. [2007]; Hestenes [2001]; Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]. Given its representational power for curved surfaces and simple encoding of complicated operations, CGA appears to hold great promise in the field of Computer Graphics. Indeed several ray-tracers/path-tracers/sphere-marchers using CGA have been implemented in the past Breuils et al. [2018,2019a]; Deul et al. [2009]; De Keninck [2019]; Dorst et al. [2007]; Hildenbrand [2007]; Wareham and Lasenby [2011]. More recently the design of more intricate surfaces has been investigated with rotors Colapinto [2017] and direct-interpolation of geometric primitives as described in Chapter 3 and published in Hadfield and Lasenby [2019]. In this chapter we will press some of these techniques into use to describe tubes and ribbons as well as to develop the techniques required to render them. Figure 4.1 shows an example of output from the CGA ray-tracer we describe in this chapter.

Conformal Geometric Algebra, CGA

The ray-tracer used in this chapter is constructed using CGA and all algebraic expressions given will be in terms of elements of this algebra. CGA adds two more basis vectors, $e$

and $\bar e$ , to the original basis vectors of 3D Euclidean space, giving a complete basis for the 5D space with the following signature: $e_1^2=e_2^2=e_3^2=e^2=1$

and $\bar e^2 = -1$ . These extra basis vectors are used to define two null vectors: $n_{\infty} = e+\bar e \equiv n$ and $n_0 = \frac{\bar{e}-e}{2} \equiv -\frac{\bar n}{2}$ – note that the $(n,\bar{n})$ notation was that originally used when Hestenes first introduced this model in Hestenes et al. [1985]. The mapping from a 3D vector, $x$

, to its corresponding CGA vector, $X$

, is given by:

$\displaystyle X = F(x) = \frac{1}{2}\left(x^2n+ 2x-\bar{n}\right) \equiv \frac{1}{2}x^2n_{\infty}+ x + n_0 .$

(28)

All vectors formed from such a mapping are null. CGA is chosen for the construction of the ray-tracer since we seek neat expressions for describing intersections, reflections and lighting models, made possible in CGA since rays and scene objects are both elements of the algebra. More background on CGA can be found in Dorst et al. [2007]; Hestenes [2001]; Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004] as well as in the introduction to this thesis.

Camera Model and Ray Casting

A pinhole camera model is used with the geometry shown in Figure 4.2. It is defined by a rotor $R_{MV}$ (where $MV$

indicates model view) incorporating rotation and translation that takes the camera from the origin to its pose in space, a focal length $f$

and two bounds $x_{\textrm{max}}$ and $y_{\textrm{max}}$ on the size of the image plane.

**Fig. 4.2:** The camera is defined by a focal length, a transformation from the origin, and bounds on the image plane.
width=

**Fig. 4.3:** An image from the ray-tracer containing examples of disks, spheres and planes.
$\includegraphics[width=1.0\textwidth]{226Achawalallobjectsimage.eps}$ width=

We take $(i,j)=(0,0)$ to be at the bottom left hand corner of the image. For an image of width $w$ and height $h$ , the world coordinates of the point $P_{ij}$ at the centre of pixel $(i,j)$ are given by:

$\begin{equation*}\begin{aligned} W_{i,j} &= F(fe_{2}-\frac{x_{\max}}{2}(1-(2i/w)... .../h))e_{3}),\\ P_{ij} &=R_{MV}W_{i,j}\tilde R_{MV}. \end{aligned}\end{equation*}$

We then generate the ray from the camera centre, $L_{ij}$ , that passes through $P_{ij}$ , via the expression

$\displaystyle L_{ij}=X_0 \wedge P_{ij}\wedge n_{\infty}$

where

is the origin transformed by the model-view rotor $R_{MV}$ to the position of the camera.

Ray Geometries for Basic Objects

Initially we will start with some basic objects representable as blades in CGA. The ray-tracer will thus initially concentrate on rendering planes, spheres and circles/discs, an example of which is shown in Figures 4.1 and 4.3.

Ray-Object Intersections

In order to compute intersections between blades, the meet operator ( $\vee$ ) is used. We will, for the proposes of this chapter, always take the meet with respect to the full 5D space rather than to the join of the blades. Thus, if $X$

is an

-grade blade, $Y$

is an

-grade blade, and the number of basis vectors in the algebra is $n$

, then :

$\displaystyle X\vee Y = \left<XY\right>_{2n-r-s}I_5,$

(30)

where $\left<Z\right>_m$ indicates the $m$

-grade component of the multivector $Z$

, and

represents the 5D pseudoscalar of the algebra.

Planes

A plane is a 4-blade and a ray is a 3-blade so the meet gives a 2-blade. If the meet itself is 0, the line lies in the plane. If the meet squared is 0, there is no finite intersection. Otherwise, the intersection point, $X$ , of a line $L$ with a plane $\Phi$ , satisfies the following: $L\vee\Phi= \lambda X\wedge n_{\infty}$ Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004], where $\lambda$ is a scalar. When extracting the 3D intersection point $x$ , we need to account for the sign and magnitude of the line in our extraction, we can do this via the constant of proportionality $\lambda$ :

$\displaystyle L\vee\Phi = \lambda X \wedge n_{\infty}= \lambda x\wedge n_{\infty}- \lambda n_{\infty}\wedge n_0.$

(31)

Therefore, $x$

can be extracted from the $e_{i}e$ and $e_{i}\bar{e}$ coefficients (for $i\in\{1,2,3\}$ ) by dividing by $\lambda$ , the $e\bar{e}$ coefficient.

Spheres

Spheres are also 4-blades and so once again, taking the meet with a ray gives a 2-blade, $F$ . With spheres, there can be zero, one or two points of intersection corresponding to the cases where ${F}^2<0$ , ${F}^2=0$ and ${F}^2>0$ respectively.

If $F = A \wedge B$ (with $A$ and $B$ null vectors) and ${F}^2\ge 0$ , the points can be extracted from the point pair/blade, $F$ , by the following formula Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]:

$\begin{equation*}\begin{aligned} \sigma_a A &= \left(1 - \frac{F}{\sqrt{-F\tilde... ...= \frac{-\sigma_b B}{(\sigma_b B)\cdot n_{\infty}}, \end{aligned}\end{equation*}$

where $\sigma_a, \sigma_b$ are scalar constants. If we define $F = A \wedge B$ with $F$

oriented in the same direction as our ray $L$

, then

is the point closest to the origin of the ray, $P_0$

, as long as our sphere is `in front of' the ray source. To ensure the alignment we can pre-normalise our sphere $S$

via the following expression:

$\displaystyle S \; \longrightarrow \; -\frac{S}{S^{*}\cdot n_{\infty}}.$

(33)

For any given sphere, its dual can square to a positive or negative number; however, by carrying out the normalisation in equation 4.6, we ensure that all spheres, $S$ , satisfy $S^{*}\cdot n_{\infty}= -1$ . If the meet of a ray $L$ and a normalised sphere $S$ , is then formed from $L\vee S$ , as in equation 4.3, the resulting bivector will be ordered as $A\wedge B$ where $A$ is the point that the ray hits first in its orientation.

In figure 4.4, the meet of the ray (direction as shown) from a point $P_0$ with the smaller sphere, would result in the point pair $A_2\wedge B_2$ . For the larger sphere in figure 4.4, the point pair resulting from the meet will be $A_1\wedge B_1$ . We extract the points from the point pair and form the distance between these points and $P_0$ (via taking the inner product). We then see that for the smaller sphere the distance of the first point is less than that of the second point, whereas for the larger sphere, the distance of the first point is larger than that of the second point – which will therefore lead us to label the larger sphere as being `behind' the point $P_0$ . This allows us to perform bounces only with spheres that are in front of the ray origin point $P_0$ .

**Fig. 4.4:** Image showing the positions of intersection points with spheres.
$\includegraphics[width=0.5\textwidth]{226AchawalSphereIntersection.eps}$

Circles/Discs

Circles are 3-blades and so the meet with a ray gives a 1-vector, $Y$

If itself is zero the ray (or line) lies in the plane of the circle and either does not intersect the circle or intersects the circle in one or two points.
If , the ray does not lie in the plane of the circle and passes through the circle disc but does not intersect.
If the line does not lie in the plane of the circle and passes outside the circle disc without intersecting.
If (and $Y \ne 0$ ) the ray intersects the circumference of the circle.

Figure 4.5 shows an example of each case along with a geometric interpretation of the form of the meet. If $Y \ne 0$ and there is an intersection, the plane containing the circle is formed by taking the wedge product between the circle and $n_{\infty}$ , $C \wedge n_{\infty}$ , and the intersection point is then extracted from the ray and this plane.

**Fig. 4.5:** The 5D meet of a ray, (in red) and a circle (black) results in a vector whose dual is the sphere (grey), the properties of this sphere vary with the relative positioning of and . In the top case the ray passes outside of the circle, passes orthogonally through both the circle and the ray and squares to a negative number as we would expect from a standard CGA sphere. In the middle case the ray hits the perimeter of the circle and the meet squares to zero, in this case the dual sphere is the special case of zero radius, it represents the intersection point itself. In the bottom case the ray passes inside the circle. Here the sphere squares to a positive scalar implying it is now an imaginary sphere and in fact the circle passes through the sphere's antipodal points.
width=

If $Y = 0$ , so that the ray lies in the plane of the circle, we need to work in 2D, so that our `meet' will result from taking the 2-part of the geometric product and dualising (with respect to the plane of the circle) to give a bivector. If the bivector has negative square there are two intersections at points $A$ and $B$ , so the bivector is $A\wedge B$ . If the bivector has positive square, there is no intersection. If the bivector squares to zero there is one intersection at $A$ , and the bivector is $a\wedge n_0$ , where $A=F(a)$ . In all cases, the intersection points are easily extracted.

Extracting Normals and Reflecting Rays

Extracting the normal to the surface of an object at a ray intersection point, ${X}$

, and the reflection of that ray at $X$

, are two fundamental building blocks in our ray tracer.

For a plane $\Phi$ which intersects with a ray, we can compute the reflection $L'$ of an incident ray $L$ (we assume $\Phi$ and $L$ have been normalised such that $\Phi^2=-1, \,L^2=1$ ) with the plane by simple sandwiching: $L' = \Phi L\Phi$ . The resulting line is oriented correctly, passes through the intersection point and $L'^2 = 1$ . For the case of a sphere $S$ , we use the following formula from Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]:

$% latex2html id marker 17594 $\displaystyle L' \propto -\left({X}\cdot(SLS)\right)\wedge n_{\infty},$$

(34)

where

is the first point of intersection. Here, $SLS$

is an example of an inversion, where the incoming ray/line, $L$

, is inverted in the sphere to give a circle which passes through the two points of intersection and the origin of the sphere (note we only have a meaningful reflection if there are two points of intersection). The tangent line to this circle at the first point of intersection, $X$

, is the reflected ray, $L'$

. Figure 4.6 illustrates this geometrical construction. Note that one can also form the tangent plane at $X$

and reflect $L$

in this plane; this is performed by the following formula:

$% latex2html id marker 17610 $\displaystyle L' \propto \Phi_X L\Phi_X, \quad \Phi_X = (X\cdot S)\wedge n_{\infty}.$$

(35)

**Fig. 4.6:** A ray in blue hits a sphere. The tangent plane at the point of impact is shown in green. The inversion of the ray in the sphere produces the red circle. The reflection of the ray in the tangent plane gives the red line. The reflected ray is also the tangent to the red circle at the point of impact.

For a circle/disc $C$ , we first form the plane $C \wedge n_{\infty}= \Phi$ in which it lies. If the ray intersects the disc (see section 4.4.1), the reflected ray can then be found using the same formula as for the plane reflection case: $% latex2html id marker 17619 $ L' \propto \left( C \wedge n_{\infty}\right) L \left( C \wedge n_{\infty}\right) $$ . Note that these expressions specifically give the (correctly oriented) reflected ray that passes through the point of intersection of the incident ray and the object, rather than a parallel ray at the origin.

We end this section with two very useful constructions which we will put to use later in the chapter. Firstly, consider the reflected ray, $L'$ , and the incident ray, $L$ , both normalised such that they square to 1. The normal line to the surface of the sphere, $N$ , can be simply found:

$% latex2html id marker 17627 $\displaystyle N \propto (L' - L ).$$

The tangent line, $T$

, (in the plane containing incident and reflected rays) can similarly be found from the sum

$% latex2html id marker 17631 $\displaystyle L_T \propto L' + L .$$

Figure 4.7 shows a graphical example of these constructions for the reflection of a ray in a sphere.

**Fig. 4.7:** An incident ray (black) hits a sphere from the right hand side of the figure. The reflected ray (also black) scatters to the top right corner of the figure. For and normalised such that , the normal (red) to the surface is proportional to . The tangent line (pink) in the plane containing the incident and reflected rays can be found from $% latex2html id marker 17651 $ L_T \propto L' + L$$ .
$\includegraphics[scale=0.6]{raytrace.PNG}$

Ray Tracing Evolved Circles

**Fig. 4.8:** Polynomial interpolation through circular control objects. (a) linear, (b) quadratic, (c) cubic.
$\includegraphics[scale=0.5]{226AchawalSurfaces.eps}$

We will now turn to an interesting class of surface that arises from the direct interpolation of CGA circles Hadfield and Lasenby [2019], examples of which are shown in Figure 4.8. In order to generate such a surface, a direct interpolation is first performed between two boundary circles, $C_1$

and

both of which are normalised such that $C_1^2 = C_2^2 =1 $

. Our interpolation is of the form:

$\displaystyle C_\alpha^{'} = \alpha C_1 + (1-\alpha) C_2,$

(36)

where we take $\alpha$ moving between 0 and 1, which moves us from $C_2$

. The result of this interpolation is not itself a valid circle and needs to be `projected' onto a blade via multiplication by a projector, which we shall call $\mathcal{S}$ . This projector has only scalar and 4-vector parts and its construction is detailed in Hadfield and Lasenby [2019] and outlined in the following.

Consider a quantity $\Sigma = \langle{\Sigma}\rangle_0 + \langle{\Sigma}\rangle_4$ . We then define the quantity $[[\Sigma]] = \sqrt{\langle{\Sigma}\rangle_0^2 - \langle{\Sigma}\rangle_4^2}$ , and with this the principal square root Dorst and Valkenburg [2011] of the scalar + 4-vector, $\Sigma$ , can be found as:

$\displaystyle \sqrt{\Sigma} = \frac{\Sigma + [[\Sigma]]}{\sqrt{2}\sqrt{\langle{... ...langle{\Sigma}\rangle_4}{\sqrt{2}\sqrt{\langle{\Sigma}\rangle_0 + [[\Sigma]]}}.$

(37)

With this square root we can then form:

$\displaystyle k{\mathcal{S}}^- = \sqrt{-C_\alpha^{'} \tilde C_\alpha^{'} },$

(38)

where ${\mathcal{S}}^-$ is ${\mathcal{S}}$ with the sign of the 4-vector part reversed and $\frac{1}{k}=\mathcal{S}^{-}\mathcal{S}$ . We then construct $k{\mathcal{S}}$ by reversing the sign of the 4-vector part, ( $k{\mathcal{S}} = \left<k\mathcal{S}^-\right>_0 - \left<k\mathcal{S}^-\right>_4$ ), and use this to produce the following expression for the projector $\mathcal{S}$ and interpolated circle $C_\alpha$ :

$\displaystyle C_{\alpha} = \frac{k\mathcal{S}}{(k\mathcal{S})(k\mathcal{S}^{-})}C_\alpha^{'} \equiv \mathcal{S}C_\alpha^{'}, \quad \alpha \in [0,1].$

(39)

Given that these surfaces may find genuine applications in computer graphics and CAD, it is desirable to explore their properties with respect to the ray tracing framework. Specifically, for a given ray and scene object, the geometric constructions of interest for lighting models are the point of intersection between a ray and a surface, and the surface normal at that specific intersection point.

In order to render this surface, we first show how to extract the intersection point with a given ray and then how to construct the surface normal at this point.

Intersection Point of Ray and Interpolated Surface

We saw earlier that the intersection of a ray with a circle produces the 1-vector $Y$ . If $Y = 0$ the ray lies in the plane of the circle and if $Y \ne 0$ and $Y^2 = 0$ there is one intersection. Therefore (in the case where the meet is not zero) to find the intersection point between our interpolated surface and a ray $L$ , we need to find a value of $\alpha$ for which:

$\displaystyle (C_\alpha \vee L)^2 = 0 \;\;\implies \;\; \left< C_\alpha L \right>_4^2 = 0.$

The system must also be tested for the case of $Y = 0$

; if an $\alpha$ exists such that $(C_\alpha \vee L)=0$ , the ray may intersect $C_\alpha$ once, twice or not at all.

**Fig. 4.9:** Left: An image showing an example interpolated surface and a ray passing through it, the circles in blue show the circles which have a meet squared of 0 with the incident ray, the red circle shows where the meet squared is minimised. Right: A plot showing the value of the meet squared as a function of $\alpha$ for this case.
$\includegraphics[width=0.8\textwidth]{226AchawalskewIntFinal.eps}$

Figure 4.9 provides a simple visual illustration of one example of the shape of this curve as a function of $\alpha$ . While this example shown in the figure is particularly smooth, experiments indicate that in the general case this function is not well approximated by low order polynomials.

Non-linear Intersection Point Finder

As it is in general difficult to extract a closed form expression for the solution to $(C_\alpha \vee L)^2 = 0$ , it is necessary to design an iterative algorithm to find the roots of the equation. Our implemented algorithm works as follows:

Check for intersection with a sphere enclosing the entire surface
Calculate the value of $(C_\alpha \vee L)^2$ at N intermediate values of $\alpha$
Record where $(C_\alpha \vee L)^2$ changes sign between successive evaluated values of $\alpha$
Locally approximate $(C_\alpha \vee L)^2$ as a quadratic equation in the region of the sign change and solve to get the value of $\alpha$ at the intersection point

Computing the intermediate objects in the surface can be done once per scene and reused for all rays calculated in that scene. To generate the enclosing sphere again we reuse the intermediary objects, in the following way:

Given and , form intermediate circles, $C_{\alpha}$ , and then calculate the bounding sphere which is given by $S_\alpha = I_5C_\alpha(C_\alpha \wedge n_{\infty})$
Construct a sphere that contains all intermediate circle bounding spheres by successive application of a two sphere bounding algorithm

Again the enclosing sphere of the object can be calculated once per scene and used for all rays. Any two sphere bounding algorithm can be used, here we chose the algorithm from Hildenbrand and Hitzer [2008] which is summarised as follows:

Ensure both spheres are normalised according to equation 4.6
Construct the line joining the centres of both spheres $L = (S_1I_5)\wedge (S_2I_5)\wedge n_{\infty}$
Intersect the line with the first sphere to produce a point pair $% latex2html id marker 17782 $ L\vee S_1 = F_1 \propto A_1\wedge B_1$$ and extract using equation (4.5)
Intersect the line with the second sphere to produce a point pair $% latex2html id marker 17786 $ L\vee S_2 = F_2 \propto A_2\wedge B_2$$ and extract using equation (4.5)
Check if $B_2\cdot(S_1I_5) > 0$ , if so encloses and so is the bounding sphere
Check if $A_1\cdot(S_2I_5) > 0$ , if so encloses and so is the bounding sphere
If neither original sphere is enclosed by the other, the new bounding sphere is given by $\frac{1}{2}(A_1 + B_2)I_5$

This iterative algorithm for the most part performs perfectly satisfactorily. When compared against specially constructed test cases for which the intersection points are known it produces negligible error. The main downside to this solution is that it is not mathematically guaranteed to give correct results especially in the case of small numbers of intermediary objects. In practice we can pre-compute large numbers of intermediary objects before rendering, allowing us to get good approximations to the function of interest. Having said that, the more intermediate objects that are created, the more computationally expensive the process is, as the root finder has to evaluate our function at each one for each ray.

`Closed Form' Solution for the Intersection of a Ray and an Evolved Circle Surface

In this section we will use $C'_{\alpha}$ to be the interpolated circle; however we emphasise that the process outlined here will also hold for the intersection of rays with other evolved objects. The intersection of the ray, $L$

, and the surface, $\mathcal{S}C'_\alpha$ (see equation 4.12) occurs when:

$\displaystyle (L\vee[\mathcal{S}C'_\alpha])^2 = 0.$

Writing $\Sigma = -C'_{\alpha}\tilde{C}'_{\alpha}$ , this can be rewritten as:

$\displaystyle \left(\frac{ L\vee\left[\langle \sqrt{\Sigma}\rangle _0C'_\alpha ... ...angle _0 - \langle \sqrt{\Sigma}\rangle _4 \right)\sqrt{\Sigma}} \right)^2 = 0.$

(40)

The denominator of this expression is never infinite (other than in the uninteresting case of $\Sigma=0$ ) and so does not contribute roots. Thus we can write:

$\displaystyle \left( L\vee[\langle \sqrt{\Sigma}\rangle _0C'_\alpha] - L\vee[\langle \sqrt{\Sigma}\rangle _4C'_\alpha] \right)^2 = 0.$

Now expanding $\sqrt{\Sigma}$ as:

$\displaystyle \sqrt{\Sigma} = \frac{\Sigma + [[\Sigma]]}{\sqrt{2}\sqrt{\langle{\Sigma}\rangle_0 + [[\Sigma]]}}$

means we can write:

$\displaystyle \left( L\vee\left[\left\langle \frac{\Sigma + [[\Sigma]]}{\sqrt{2... ...\Sigma}\rangle_0 + [[\Sigma]]}}\right\rangle _4C'_\alpha\right] \right)^2 = 0.$

Again the denominator of the square root function is simply a scalar which is never infinite, thus it contributes no roots and we can write:

$\displaystyle \left( L\vee\left[\langle \Sigma + [[\Sigma]]\rangle _0C'_\alpha ... ...vee\left[\langle \Sigma + [[\Sigma]] \rangle _4C'_\alpha\right] \right)^2 = 0.$

The quantity $[[\Sigma]]$ is a scalar and so distributing the grade selection operators gives us:

$\displaystyle \left( L\vee [\langle \Sigma \rangle _0 C'_\alpha]+ [[\Sigma]]L\v... ..._\alpha - L\vee\left[\langle \Sigma \rangle _4C'_\alpha\right] \right)^2 = 0.$

Expanding this leads to:

0	$\displaystyle =[L\vee (\langle \Sigma\rangle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha) ]^2 + [[\Sigma]]^2(L\vee C'_\alpha)^2$
	$\displaystyle \phantom{=}+[[\Sigma]]\left\{ (L\vee (\langle \Sigma\rangle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha))(L\vee C'_\alpha)\right.$
	$\displaystyle \phantom{=}\left. + (L\vee C'_\alpha)(L\vee (\langle \Sigma\rangle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha)) \right\}.$	(41)

To make progress on solving this we recall that $\Sigma = -C'_\alpha\tilde{C'_\alpha}$ and that for our linear interpolation of circles we have defined $C'_\alpha$ as:

$\displaystyle C'_\alpha = \alpha C_1 + (1 - \alpha)C_2 = \alpha(C_1 - C_2) + C_2.$

As $C'_\alpha$ has terms in $\alpha$ of order 1 we would expect $\Sigma$ to have terms of order $2$

. Continuing on this train of thought one might suspect that it is possible to re-write equation (4.14) as a simple polynomial in $\alpha$ . However, it is easy to see that this is not possible due to $[[\Sigma]]$ , which is a scalar polynomial in $\alpha$ enclosed entirely in a square root:

$\displaystyle [[\Sigma]] = \sqrt{\langle \Sigma\rangle _0^2 - \langle \Sigma\rangle _4^2} .$

Thus in order to solve equation (4.14) we need to rearrange:

	$\displaystyle [L\vee (\langle \Sigma\rangle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha) ]^2 + [[\Sigma]]^2(L\vee C'_\alpha)^2$
	$\displaystyle \hspace*{10ex}= -[[\Sigma]] \left\{L\vee (\langle \Sigma\rangle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha)\right\}(L\vee C'_\alpha)$
	$\displaystyle \phantom{=}\hspace*{10.75ex} + (L\vee C'_\alpha)\left\{L\vee (\la... ...Sigma\rangle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha)\right\} .$	(42)

We can then square both sides of the equation, eliminating the square root in the process:

	$\displaystyle \left([L\vee (\langle \Sigma\rangle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha) ]^2 + [[\Sigma]]^2(L\vee C'_\alpha)^2 \right)^2$
	$\displaystyle \hspace*{10ex} = [[\Sigma]]^2 \left[ \left\{L\vee (\langle \Sigma... ...a) - L\vee (\langle \Sigma\rangle _4C'_\alpha)\right\}(L\vee C'_\alpha) \right.$
	$\displaystyle \phantom{=}\hspace*{10ex} \left. + (L\vee C'_\alpha)\left\{L\vee ... ...gle _0C'_\alpha) - L\vee (\langle \Sigma\rangle _4C'_\alpha)\right\} \right]^2.$	(43)

Expanding this out will give a polynomial in $\alpha$ – it turns out that this is a scalar polynomial due to the fact that $\langle \Sigma\rangle _4 C'_\alpha$ has only trivector components. This polynomial can then be solved with any numerical polynomial solver such as finding the eigenvalues of the companion matrix Horn and Johnson [2012].

For this case of linear evolution of circles we will get a polynomial of order 12, implying 12 potential roots. In reality 6 of these roots are extraneous, generated by the process of squaring to handle the square root term in $[[\Sigma]]$ . Some of the 6 remaining roots may be imaginary, some may be outside of the range $0 \leq \alpha \leq 1$ and some will be spurious roots corresponding to $\mathcal{S} = 0$ . To filter out the valid roots we simply take all roots between 0 and 1 and evaluate $(L\vee[\mathcal{S}C'_\alpha])^2$ at these positions, selecting the roots for which $\vert(L\vee[\mathcal{S}C'_\alpha])^2\vert < \epsilon$ for some small $\epsilon$ threshold where $\epsilon > 0$ (in our experiments $\epsilon = 10^{-6}$ works satisfactorily).

An interesting point to note here is that we could extend this intersection finding method to $C'$ being higher order functions of $\alpha$ , so long as $(L\vee[\mathcal{S}C'_\alpha])^2 = 0$ . Generating such higher order splines through geometric primitives is described in section 4.8.

A Comment on Rendering Speed

For this chapter our raytracer was implemented in Python with the Clifford Library Hadfield et al. [sent]. It is simply an investigative tool used as a framework in which to conduct basic research into the shapes and properties of surfaces as well as the algorithms used to render them. Of course in a production computer graphics environment, a higher performance language such as C/C++/GLSL would be required and the trade off of accuracy for speed with regard to the number of intermediary objects would need to be closely analysed. Such an analysis would require very careful benchmarking and comparison across multiple modern computer architectures and as such is beyond the scope of this chapter.

Analytic Form for Normals

Given the $\alpha$ for which the ray intersects the surface, we have both the interpolated circle, $C_\alpha$ , and the point of intersection $X$

. Using the result from Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004], which is also used in equation 4.7, we extract a tangential line $L_C$

in the plane of the circle at $X$

$\displaystyle L_C = (X\cdot C_{\alpha}) \wedge n_{\infty}.$

(44)

We would now like an analytic form for the tangent to the surface corresponding to evolving the surface through an increment of $\alpha$ , postulating that this will be orthogonal to $L_C$

: some future work remains to understand how these two tangent vectors are related to the directions of principal curvature. Clearly $\frac{dC(\alpha)}{d\alpha} \equiv {\dot C}_{\alpha}$ will be a key quantity in deriving this additional tangent vector. A first observation is that the circle and its derivative will be orthogonal to one another, i.e. ${\dot C} \cdot C = 0$ , and that the geometric product is minus itself under reversion, i.e. ${\dot C} \, C = - C \, {\dot C}$ (note that here, and in what follows, we will drop the $\alpha$ subscript on $C$

). This follows from the fact that $C^2 = C\cdot C = 1$ (our circles are all normalised), so that:

$\displaystyle \frac{d}{d\alpha}\left(C\cdot C\right)=C\cdot {\dot C} + {\dot C}... ..., \quad \frac{d}{d\alpha}\left(C\, C\right) = C\, {\dot C} + {\dot C} \, C = 0.$

(45)

Since $C \cdot \dot{C} = - C\cdot \dot C$ and they are both scalars, this tells us $\dot C \cdot C = 0$ . Using the fact that $\tilde{C}=-C$ , we see that $C \, \dot{C} = - C \, \dot C = -(C \, \dot{C})^{\tilde{}}$ . As there are no 6-vector parts, this indicates that the product can only have bivector parts (this is a standard construct in many areas, the most obvious being rigid body dynamics Doran and Lasenby [2003]). Let us call this bivector, $\Omega_C$ :

$\displaystyle \Omega_C = C\dot{C}.$

(46)

Using the analogy with rigid body dynamics, we think of this bivector as the angular velocity bivector of the circles as they evolve under the parameter $\alpha$ . We note here that a similar construction would be possible for the other main objects that we use in CGA, since they are all normalised to 1 or 0. The null vectors representing points, $X$ , have a constant `length' due to normalisation, so as with the circles, we can differentiate wrt $\alpha$ to see that $X$ and $\dot X$ are orthogonal, ie $X\cdot \dot{X} = 0$ . If we were to define the `velocity', $\dot{X}$ to be the inner product with the angular velocity bivector given in equation 4.19:

$\displaystyle \dot X = X \cdot \Omega_C = X \cdot (C\dot C),$

(47)

the condition $X\cdot \dot{X} = 0$ is satisfied since $X\cdot(X\cdot B)=(X\wedge X)\cdot B = 0$ . Thus, given an $X$

on the surface, lying on a circle with parameter $\alpha$ , the $\dot{X}$ defined above will preserve its length and is, we claim, the tangential direction required. In order to show this, the first thing we must do is establish that if we evolve $X$

according to this rule, generating a quantity we call $X(\alpha)$ , then $X(\alpha)$ should lie on $C_{\alpha}$ , for all $\alpha$ , i.e.,

$\displaystyle X(\alpha) \wedge C_{\alpha} = 0.$

Differentiating this (and again dropping the subscript $\alpha$ for clarity) and using $\dot{X}=X\cdot (C\dot{C})$ , gives

$\displaystyle \dot{X}\wedge C + X\wedge \dot{C}= 0$
$\displaystyle \implies \dot{X}\wedge C \equiv \left(X\cdot (C\dot{C})\right) \wedge C = - X\wedge \dot{C}.$

We now expand this expression using standard expansion results ( $a\cdot (A_r\wedge B_s) = (a\cdot A_r)\wedge B_s + (-1)^r A_r\wedge (a\cdot B_s)$ ):

$\displaystyle \left(X\cdot (C\dot{C})\right) \wedge C$	$\displaystyle = X\cdot \left((C\dot{C})\wedge C\right) - (C\dot{C}) \wedge (X\cdot C)$
	$\displaystyle = - \frac{1}{2}\langle C\dot{C}(XC+CX)\rangle _4$
	$\displaystyle = \frac{1}{2}\langle \dot{C}C(XC+CX)\rangle _4.$	(48)

The first term on the right hand side of the first line of this expansion is zero as $(C\dot{C})\wedge C = \langle C\dot{C}C\rangle _5=\langle -C^2\dot{C}\rangle _5 = \langle -\dot{C}\rangle _5 = 0$ . Since $X$

lies on

and so $X\wedge C=0$ , we see that $XC=CX$

which means that

$\displaystyle \frac{1}{2}\langle \dot{C}C(XC+CX)\rangle _4 = \langle \dot{C}C^2X\rangle _4 = \langle \dot{C}X\rangle _4 = \dot{C}\wedge X = -X\wedge \dot{C}$

giving $\dot{X}\wedge C = -X\wedge \dot{C}$ as required, so the proposed evolution is compatible with the constraint.

If we therefore assume that $\dot{X}$ is the direction we want, we can calculate the tangent line in this direction via:

$\displaystyle L_T = \dot X \wedge X \wedge n_{\infty}.$

(49)

The fact that lines $L_C$ and $L_T$ are perpendicular can be verified by showing that the quantity $L_T L_C$ has only a bivector part (see Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004] for a discussion of when intersecting lines are orthogonal – if two lines $L_1$ and $L_2$ intersect at a point, then $\langle L_1L_2\rangle _4 =0$ . In addition, if they are orthogonal, $\langle L_1L_2\rangle _0=0$ ). If this is the case, $L_T L_C$ will reverse to minus itself.

To show this, we need to consider the reverse of $(\dot X \wedge X \wedge n_{\infty})((X\cdot C)\wedge n_{\infty})$ . We will need the facts that that $XC=CX$ , $\dot{X}C=-C\dot{X}$ , $C\dot{C}=-\dot{C}C$ , $X\dot{X} = -\dot{X}X$ and $\tilde{C}=-C$ . We have shown all of these identities earlier in this section. We also need an additional fact, which is that $\dot{X}$ anticommutes with $C$ . To see this we use another standard result ( $a\wedge (A_r\cdot B_s) = (a\cdot A_r)\cdot B_s + (-1)^{r} A_r \cdot(a\wedge B_s)$ ):

$\displaystyle \dot{X}\cdot C = \left( X\cdot (C\dot{C})\right) \cdot C = X\wedge \left( (C\dot{C})\cdot C\right) - (C\dot{C})\cdot(X\wedge C).$

(50)

The first term on the RHS of this equation is zero as $(C\dot{C})\cdot C=\langle C\dot{C}C\rangle _1 = \langle -C^2 \dot{C}\rangle _1 = 0$ , and the second term on the RHS is also zero as $X$

lies on

so $X\wedge C=0$ . Thus $\dot{X}\cdot C=0$ and $\dot{X}$ therefore anticommutes with $C$

as required.

We are now in a position to expand out $(\dot X \wedge X \wedge n_{\infty})((X\cdot C)\wedge n_{\infty})$ :

	$\displaystyle (\dot X \wedge X \wedge n_{\infty})((X\cdot C)\wedge n_{\infty})$
	$\displaystyle \hspace*{10ex}=\left( (\dot{X}X)\wedge n_{\infty}\right) \left( (XC)\wedge n_{\infty}\right)$
	$\displaystyle \hspace*{10ex}=\frac{1}{4}\left[ \left( \dot{X}Xn_{\infty}+ n_{\infty}\dot{X}C\right) \left( XCn_{\infty}+ n_{\infty}XC\right) \right]$
	$\displaystyle \hspace*{10ex}=\frac{1}{4}\left[ \dot{X}Xn_{\infty}XCn_{\infty}+ n_{\infty}\dot{X}Xn_{\infty}XC\right].$	(51)

In the above we have used the facts that $\dot{X}\wedge X = \dot{X}X$ , $XC = X\cdot C$ and $X^2=0$ . Note that the term $Xn_{\infty}X$ in the final line of equation 4.24 can be written as $2(X\cdot n_{\infty})X$ (from the standard reflection formula and the fact that $C^2=0$ ). Reversing the final line of equation 4.24 and using the commutation and anticommutation relations discussed, it is easy to show that the reverse of $L_T L_C$ is indeed minus itself, implying it has only a bivector part, as required, meaning the lines are orthogonal. Note that this result relies crucially on the fact that $\dot{X}$ and $C$ anticommute, which is a good indication that $\dot{X}$ lies in the right direction.

Given these two orthogonal tangent lines $L_C$ and $L_T$ , we can construct the plane tangent to the surface at $X$ by computing the join of the two lines. Or, we can bypass the plane entirely and compute the surface normal line directly as:

$\displaystyle N = \left<L_TL_C\right>_2I_5.$

(52)

Calculating the Derivative of the Object Manifold Projection

To calculate $\dot{C}$ we must differentiate the projection onto the blade manifold of our interpolated object with respect to our evolution parameter $\alpha$ . We will continue to work with circles but note that the process works with the general case where $C'_\alpha$ is any pure-grade multivector which is a function of a scalar parameter $\alpha$ . Let the projection of $C'_\alpha$ onto the blade manifold be given by:

$\displaystyle C_\alpha = \mathcal{S}C'_\alpha$

where $\mathcal{S}$ is our blade projector. The differential of this with respect to $\alpha$ is given by:

$\displaystyle \frac{\partial C_\alpha}{\partial \alpha} = \frac{\partial \mathc... ...ial \alpha} C'_\alpha + \mathcal{S}\frac{\partial C'_\alpha}{\partial \alpha} .$

(53)

Thus, any closed form expression for the derivative on the manifold will first require a closed form for the derivative of the projector $\frac{\partial \mathcal{S}}{\partial \alpha}$ . Recall from equation 4.13 that we can write the projector $\mathcal{S}$ as a function of $\sqrt{\Sigma}$ , where $\Sigma = -C'_\alpha\tilde{C}_{\alpha}$ , and so

$\displaystyle C_\alpha = \frac{ \langle \sqrt{\Sigma}\rangle _0 - \langle \sqrt... ...ngle \sqrt{\Sigma}\rangle _4 \right)\sqrt{\Sigma}} C'_\alpha \equiv SC'_\alpha.$

(54)

Thus to find an expression for $\frac{\partial \mathcal{S}}{\partial \alpha}$ we will first need one for $\frac{\partial \sqrt{\Sigma}}{\partial \alpha}$ .

Closed Form Derivative of the Square Root Operation

The closed form for the derivative of the principal square root function can be found by repeated application of the chain and product rules:

$\displaystyle \frac{\partial [[\Sigma]]}{\partial \alpha} = \frac{ \langle{\fra... ...4\langle{\frac{\partial \Sigma}{\partial \alpha} }\rangle_4 } { 2 [[\Sigma]] },$

(55)

where we are using the fact that $\left\langle \frac{\partial \Sigma}{\partial \alpha} \right\rangle _g= \frac{\partial \langle \Sigma \rangle _g}{\partial \alpha}$ .

$\displaystyle \frac{\partial }{\partial \alpha} \left( \frac{1}{\sqrt{2}\sqrt{\langle{\Sigma}\rangle_0 + [[\Sigma]]}} \right)$	$\displaystyle = \frac{-1}{2\sqrt{2}}(\langle{\Sigma}\rangle_0 + [[\Sigma]])^{-\... ...rtial \alpha} }\rangle_0 + \frac{\partial [[\Sigma]]}{\partial \alpha} \right),$
$\displaystyle \frac{\partial \sqrt{\Sigma}}{\partial \alpha}$	$\displaystyle = \left( \frac{\partial \Sigma}{\partial \alpha} + \frac{\partial... ...t)\left( \frac{1}{\sqrt{2}\sqrt{\langle{\Sigma}\rangle_0 + [[\Sigma]]}} \right)$
	$\displaystyle \phantom{=}+\left( \Sigma + [[\Sigma]] \right)\frac{\partial }{\p... ...ac{1}{\sqrt{2}\sqrt{\left\langle{\Sigma}\right\rangle_0 + [[\Sigma]]}} \right).$	(56)

Thus the derivative of $\sqrt{\Sigma}$ is a function only of $\Sigma$ and $\frac{\partial \Sigma}{\partial \alpha}$ which in turn can be written in terms of $C'_\alpha$ and $\frac{\partial C'_\alpha}{\partial \alpha}$ :

$\displaystyle \Sigma = -C'_\alpha\tilde{C'_\alpha}, \qquad \frac{\partial \Sig... ...lde{C'_\alpha} - C'_\alpha\tilde{\frac{\partial C'_\alpha}{\partial \alpha} }.$

Closed Form Derivative of the Projector

With our square root derivative in place we can proceed to finding the derivative of the projector $\mathcal{S}$ . Recall, $\mathcal{S}$ is given by:

$\displaystyle \mathcal{S} = \frac{ \langle \sqrt{\Sigma}\rangle _0 - \langle \... ...t{\Sigma}\rangle _0 - \langle \sqrt{\Sigma}\rangle _4 \right)\sqrt{\Sigma}} .$

We can again differentiate this with repeated applications of the chain and product rule:

$\displaystyle \frac{\partial \left(\langle \sqrt{\Sigma}\rangle _0 - \langle \s... ...\left\langle \frac{\partial \sqrt{\Sigma}}{\partial \alpha} \right\rangle _4,$

	$\displaystyle \frac{\partial \left(\left[\left(\langle \sqrt{\Sigma}\rangle _0 ... ...t{\Sigma}\rangle _4 \right)\sqrt{\Sigma}\right]^{-1}\right)}{\partial \alpha} =$
	$\displaystyle \hspace*{5ex}\left[\left(\left\langle \frac{\partial \sqrt{\Sigma... ...a}\rangle _4 \right)\frac{\partial \sqrt{\Sigma}}{\partial \alpha} \right] \ast$
	$\displaystyle \hspace*{45ex}\left[-\left(\langle \sqrt{\Sigma}\rangle _0 - \langle \sqrt{\Sigma}\rangle _4 \right)\sqrt{\Sigma}\right]^{-2}$

and so finally we have our closed form expression for the projector derivative:

$\displaystyle \frac{\partial \mathcal{S}}{\partial \alpha}$	$\displaystyle =\left(\left\langle \frac{\partial \sqrt{\Sigma}}{\partial \alpha... ...a}\rangle _0 - \langle \sqrt{\Sigma}\rangle _4 \right)\sqrt{\Sigma}\right]^{-1}$
	$\displaystyle \phantom{=}\hspace*{5ex}+\left(\langle \sqrt{\Sigma}\rangle _0 - ... ...t{\Sigma}\rangle _4 \right)\sqrt{\Sigma}\right]^{-1}\right)}{\partial \alpha} .$	(57)

Now consider $C'_\alpha$ to be an interpolated circle of the form $C'_\alpha = C'_\alpha = \alpha C_1 + (1-\alpha)C_2$ . The derivative of this with respect to $\alpha$ is a constant:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = C_1 - C_2.$

This derivative is the final piece required for equation (4.26), giving us a completely closed form for $\frac{\partial C_\alpha}{\partial \alpha} \equiv \dot{C}$ .

An important point to note here is that this blade projection derivative is grade-agnostic and so can be used for objects other than just evolved circles.

Ray Tracing Evolved Point Pairs

**Fig. 4.10:** Polynomial interpolation through point-pair control objects. From left to right: linear, quadratic, cubic.

We will return to the actual ray tracing of circles later (see Figures 4.16, 4.17), but first we turn our attention to point pairs. Due to the mathematical similarities between circles and point-pairs in CGA Dorst et al. [2007], as well as the practical desire to represent ribbon-like surfaces, we can apply similar ray-tracing methods to surfaces formed from the interpolation of point-pair bivectors representing line segments. If $P_1$ and $P_2$ are point-pairs which represent a line segment, we form a surface via:

$\displaystyle P_{\alpha} = \mathcal{S}P'_{\alpha} = \mathcal{S}(\alpha P_1 + (1-\alpha)P_2),$

where again, $\mathcal{S}$ is a scalar plus 4-vector which maps the interpolated bivector onto a 2-blade. Figure 4.10 gives examples of such surfaces.

Closed Form Solution for the Intersection of a Ray and an Evolved Point-Pair Surface

To find the intersection point of a ray and these surfaces we again form the meet of a ray $L$

and the form of the evolved point-pair $P_\alpha$ . The result is a scalar quantity that can be written as:

$\displaystyle L\vee P_\alpha \equiv (L^{*}\wedge P_\alpha^{*}) I_5.$

In the case that the ray and the line that passes through both of the points in the point-pair in the surface (also known as the carrier line) intersect, this will give an answer of zero:

$\displaystyle L\vee P_\alpha = 0.$

As with the evolved circle surfaces we will attempt to construct this as a simple polynomial in $\alpha$ . We start with:

$\displaystyle L\vee [\mathcal{S}P'_\alpha] = 0.$

Expressing $\mathcal{S}$ in terms of $\Sigma$ where as before, $\Sigma = -P'_{\alpha}\tilde{P}'_{\alpha}$ , gives

$\displaystyle L\vee\left[\frac{ \langle \sqrt{\Sigma}\rangle _0 - \langle \sqrt... ...- \langle \sqrt{\Sigma}\rangle _4 \right)\sqrt{\Sigma}} P'_\alpha\right] = 0.$

As before, the denominator cannot usefully be zero, giving:

$\displaystyle L\vee[\langle \sqrt{\Sigma}\rangle _0P'_\alpha] - L\vee[\langle \sqrt{\Sigma}\rangle _4P'_\alpha] = 0,$
$\displaystyle L\vee\left[\left\langle \frac{\Sigma + [[\Sigma]]}{\sqrt{2}\sqrt{... ...rt{\langle{\Sigma}\rangle_0 + [[\Sigma]]}}\right\rangle _4P'_\alpha\right] = 0,$
$\displaystyle \left( L\vee [\langle \Sigma \rangle _0 P'_\alpha]+ [[\Sigma]]L\vee P'_\alpha - L\vee\left[\langle \Sigma \rangle _4P'_\alpha\right] \right)^2 = 0.$

as again, the denominator is not zero.We now take the term containing $[[\Sigma]]$ (a scalar), to the RHS;

$\displaystyle L\vee [\langle \Sigma \rangle _0 P'_\alpha] - L\vee\left[\langle \Sigma \rangle _4P'_\alpha\right] = -[[\Sigma]]L\vee P'_\alpha.$

Squaring then gives us:

$\displaystyle (L\vee [\langle \Sigma \rangle _0 P'_\alpha] - L\vee\left[\langle \Sigma \rangle _4P'_\alpha\right])^2 = [[\Sigma]]^2(L\vee P'_\alpha)^2.$

allowing us to form a simple scalar polynomial in $\alpha$ (we can see that this produces a scalar equation since $\langle \Sigma \rangle _4P'_\alpha$ has only bivector parts):

$\displaystyle (L\vee [\langle \Sigma \rangle _0 P'_\alpha] - L\vee\left[\langle \Sigma \rangle _4P'_\alpha\right])^2 - [[\Sigma]]^2(L\vee P'_\alpha)^2 = 0,$

(58)

which can again be solved with a fast numerical polynomial solver.

This intersection equation for linearly interpolated point pairs is of order 6, implying there are up to 6 potential hitting points. Again the same process can be used to filter the roots as was done for the roots of the circle intersection equation.

Bounding Sphere and Normal Calculation

We saw earlier that the meet will be zero if the ray hits anywhere along the carrier line of the point-pair $L_C=P\wedge n_{\infty}$ . Assuming the carrier line and ray do meet, the point of intersection can be extracted via the method outlined in the last section of Chapter 3 and published in Hadfield and Lasenby [2019]. Given that the carrier line of $P_{\alpha}$ (for some $\alpha$ ) and the ray intersect at a point $X$ , we can then check if the intersection point is within the bounding sphere $S = P_1\wedge P_2$ of the surface by ensuring:

$\displaystyle S^*\cdot X = [(P_1\wedge P_2)I_5]\cdot X > 0.$

Since the endpoints of all interpolated point-pairs will lie on the surface of $S$

(see Hadfield and Lasenby [2019]), the above condition ensures there is an intersection with the line segment and not just the carrier line of the point-pair. To find the normal to the point-pair surface we can simply use exactly the same argument, and in fact the same code, as we did before for the evolved circles case but this time extracting $L_C$

as:

$\displaystyle L_C = P\wedge n_{\infty}.$

Figure 4.11 shows an example of rendering a surface composed of interpolated point-pairs.

**Fig. 4.11:** Ray tracing evolved point-pairs. Left: The scene to be rendered, in blue is a representation of the point-pair surface to be rendered, the camera frustum is shown in black, the camera axis is shown in red. Right: The resultant rendered surface of interpolated point-pairs.
$\includegraphics[width=0.4\textwidth]{pntpairsurf.PNG}$

Special Cases of Evolved Point-Pairs

A special case of the evolved point-pairs occurs when they are co-planar and form chords of a circle, Figure 4.12 shows two examples of this. As proved in Hadfield and Lasenby [2019] this special case results in the 4-vector part of the projector becoming zero implying the interpolation requires no re-projection back to the object manifold, i.e.,

$\displaystyle P_\alpha = \alpha P_1 + (1-\alpha)P_2 .$

In this case the intersection of the carrier line with the ray can be found by looking for a point at which:

$\displaystyle P_\alpha\vee L = \alpha P_1\vee L + (1-\alpha)P_2\vee L = 0.$

Re-arranging gives an expression for $\alpha$ :

$\displaystyle \alpha = \frac{P_2\vee L}{ P_2\vee L - P_1\vee L }.$

If the $\alpha$ derived from the above expression is between 0 and 1 then there is an intersection of the ray with the carrier line.

The disappearing 4-vector part of the projector, which is proportional to $P_1\wedge P_2$ , allows the ray-tracer to detect these cases and thus reduces the computational expense of a ray-surface intersection considerably.

**Fig. 4.12:** Two examples of point-pair interpolations for which all intermediary objects are blades without requiring projection to the object manifold.

Triangular Facets from Evolved Point-Pairs

The intersection of co-circular point-pairs also allows us to examine the intersection of rays with triangular facets. Consider a ray $L$

and a set of three points $A$

which together form a triangular facet. First we will form a set of normalised point pairs:

$\displaystyle P_1 = \frac{A\wedge C}{\vert A\wedge C\vert}, \quad P_2 = \frac{A... ... B}{\vert A\wedge B\vert}, \quad P_3 = \frac{C\wedge B}{\vert C\wedge B\vert}.$

We can then check if the ray intersects the facet by computing two scalar quantities

$\displaystyle \alpha = \frac{P_2\vee L}{ P_2\vee L - P_1\vee L }, \quad \quad \beta = \frac{P_3\vee L}{ P_3\vee L - P_2\vee L }.$

If both $\alpha$ and $\beta$ are between 0 and 1 then the ray hits the facet. Figure 4.13 shows an example of rendering a triangular facet using this technique. It is of course possible to combine together multiple triangular facets and thus make meshes. Note that the line-facet intersection problem is not new in CGA, an alternative solution is already known via a reciprocal frame construction equivalent to barycentric coordinates and is well demonstrated in the raytracers of Dorst et al. [2007]; Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004].

**Fig. 4.13:** Ray tracing a triangular facet.
$\includegraphics[width=0.45\textwidth]{triangular_scene.PNG}$

Bézier Curves and Hermite Splines through Geometric Primitives

So far we have restricted our mathematics to linear interpolation of objects but have hinted that higher order interpolations are possible. A commonly used family of higher order interpolating curves are the Bézier curves Bezier [1986], which in the cubic case and with specific first order endpoint conditions are known as Hermite curves.

Linear Interpolation as a Linear Bézier Curve

The simplest form of Bézier curve is simply a linear interpolation between two vectors. If we replace the vectors with $k$

-blades and couple with the projection to the blade manifold we have the exact same linear interpolation, although this time with $\alpha$ going in the other direction. Adopting a notation of $C_0$

as the first object and $C_1$

as the second:

$\displaystyle C'_\alpha = (1-\alpha)C_0 + \alpha C_1, \quad C_\alpha = \mathcal{S}C'_\alpha.$

In sections 4.5.2 and 4.6, our analysis to extract surface normals was based on having an expression for the derivative of the pure grade multivector as a function of $\alpha$ . For the case of the linear interpolation the solution is constant:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = C_1 - C_0 .$

Quadratic Bézier Curve

With three multivectors we can specify a quadratic function of $\alpha$ :

$\displaystyle C'_\alpha = (1 - \alpha)^{2}C_0 + 2(1 - \alpha)\alpha C_1 + \alpha^{2}C_2.$

This is known as a quadratic Bézier curve. Again we can take derivatives:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = 2 (1 - \alpha) (C_1 - C_0) + 2 \alpha (C_2 - C_1).$

Cubic Bézier Curves

With four control multivectors we get the the most commonly used form of Bézier curve, the cubic Bézier curve:

$\displaystyle C'_\alpha = (1-\alpha)^3 C_0 + 3(1-\alpha)^2 \alpha C_1 + 3(1-\alpha)\alpha^2 C_2 + \alpha^3 C_3.$

Again we can take derivatives allowing us to extract surface normals:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = 3(1-\alpha)^2(C_1 - C_0) + 6(1-\alpha)\alpha(C_2 - C_1) + 3\alpha^2(C_3 - C_2).$

Figure 4.14 shows examples of orders 1,2 and 3 Bézier interpolation through circles, along with the control objects used to generate the surface.

**Fig. 4.14:** Projected multivector Bézier curves of progressively higher order. The surfaces are shown in grey while the control objects are shown in red.
[Order 1] [Order 2] [Order 3]

$N^{\text {th}}$ Order Bézier Curve

More generally we can say that an $N^{\text {th}}$ order multivector Bézier curve is of the form

$\displaystyle C'_\alpha = \sum_{i=0}^N b_{i,N} C_i,$

where

$\displaystyle b_{i,N} = \begin{cases} {N\choose i}\alpha^i(1 - \alpha)^{N - i},\;\; & 0\leq i \leq N,\\ 0, \;\; & \text{otherwise}, \end{cases}$

are known as the Bernstein polynomials. The derivative of our $N^{\text {th}}$ order Bézier curve is:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = \sum_{i=0}^N \frac{\partial b_{i,N}}{\partial \alpha} C_i$ where $\displaystyle \quad \frac{\partial b_{i,N}}{\partial \alpha} = N(b_{i-1,N-1} - b_{i,N-1}).$

If we re-arrange our coefficients the Bézier curve derivative can also be written in the form:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = N \sum_{i=0}^{N-1} b_{i,N-1} (C_{i+1} - C_i).$

Rational Bézier Curves

A rational Bézier curve adds weights $w_i$

to the polynomials allowing them to represent a broader class of curves:

$\displaystyle C'_\alpha = \frac{ \sum_{i=0}^N b_{i,N} C_{i}w_i } { \sum_{i=0}^N b_{i,N} w_i }.$

Again, a closed form for their derivatives with respect to $\alpha$ can be calculated:

$\begin{multline} \frac{\partial C'_\alpha}{\partial \alpha} = \frac{1}{ [\sum_{i... ...^N \frac{\partial b_{i,N}}{\partial \alpha} w_i \right] \right). \end{multline}$

Thus we can additionally represent projected multivector rational Bézier curves and calculate analytic normals to the evolved surfaces formed.

Hermite Cubic Curves and Splines

Hermite cubic curves are another common form of interpolating curve. They are defined by control points, $C_i$

, (where we use the notation $C$

for points, but will see shortly that these can be replaced by objects) and associated tangent vectors, $V_i$

, at each end of the curve, for $\alpha \in [0,1]$ :

$\displaystyle C'_\alpha = (2\alpha^3-3\alpha^2+1)C_0 + (\alpha^3-2\alpha^2+\alpha)V_0 + (-2\alpha^3+3\alpha^2)C_1 +(\alpha^3-\alpha^2)V_1.$

The derivative of the curve is:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = (6\alpha^2-6\alpha)C_0 + (3\alpha^2-4\alpha+1)V_0 + (-6\alpha^2+6\alpha)C_1 +(3\alpha^2-2\alpha)V_1.$

Cubic Hermite curves can be converted to cubic Bezier curves and vice-versa. As with Bezier curves, putting multivectors and multivector derivatives instead of the control points and tangents will give us a multivector valued curve.

A very common use of Hermite curves is in the construction of Hermite splines; these are piece-wise constructions in which multiple Hermite curves are placed end to end, sharing tangent vectors and control points at each endpoint. By constructing a curve in this way, a C1 continuous piece-wise curve is designed that passes through the control points exactly.

When moving the spline generation process to the multivector domain we must check whether the blade projection introduces problems with C1 continuity on the manifold. To check C1 continuity we need to evaluate the curve derivative either side of a junction between curves in the spline. Consider the form of the derivative:

$\displaystyle \frac{\partial C_\alpha}{\partial \alpha} = \frac{\partial \mathc... ...ial \alpha} C'_\alpha + \mathcal{S}\frac{\partial C'_\alpha}{\partial \alpha} .$

(60)

Let us now evaluate this at the endpoint of the $n^{\text{th}}$ curve in a piece-wise spline where $\alpha=1$ and of curve ( $n+1$

) where $\alpha=0$ . First we note that on both curves at these points, $\mathcal{S} = 1$ because the curve passes through a blade control object which requires no projection. Additionally we see that by definition of the Hermite spline at that point, the derivative in pure-grade space is shared across both curves as is the control point:

$\displaystyle \frac{\partial C'_\alpha}{\partial \alpha} = V_{n,n+1}, \quad C'_\alpha= C_{n,n+1},$

where $V_{n,n+1}$ and $C_{n,n+1}$ are the derivative (or `tangent') and control objects respectively of the curve that are shared between segments $n$

and

Thus for the derivative to evaluate to the same on either side of the boundary, we only need to check that $\frac{\partial \mathcal{S}}{\partial \alpha}$ is the same either side of the boundary. Considering the equations in section 4.6.1 we can see that $\frac{\partial \mathcal{S}}{\partial \alpha}$ is a function only of $C'$ and $\frac{\partial C'}{\partial \alpha}$ which are both constant across the junction. Thus the curve is C1 continuous on the manifold as required. With assurances that the spline is continuous across the boundaries we are free to chose any means of generating tangents in the pure-grade space that we like.

One such mechanism for generating tangents for Hermite splines comes from Kochanek and Bartels Kochanek and Bartels [1984]. The Kochanek–Bartels (KB) spline is an interpolating spline with three scalar design parameters $t,b,c$ known as tension, bias and continuity respectively. For given control objects $C_i, C_{i+1}$ the corresponding tangents $V_i, V_{i+1}$ can be calculated using the control objects in the spline $C_{i-1}$ and $C_{i+2}$ which lie previous to, and after, the curve in the order of the spline:

$\displaystyle V_i$	$\displaystyle = \frac{(1-t)(1+b)(1+c)}{2}(C_i-C_{i-1}) + \frac{(1-t)(1-b)(1-c)}{2}(C_{i+1}-C_i),$
$\displaystyle V_{i+1}$	$\displaystyle = \frac{(1-t)(1+b)(1-c)}{2}(C_{i+1}-C_{i}) + \frac{(1-t)(1-b)(1+c)}{2}(C_{i+2}-C_{i+1}).$

Setting all three scalar parameters to a value of 0 produces the commonly used Catmull–Rom spline Catmull and Rom [1974]. Figure 4.15 shows an example of a KB spline of multivector geometric primitives.

**Fig. 4.15:** A Kochanek–Bartels spline of evolved circles meshed, textured, and rendered with smooth shading.

**Fig. 4.16:** Left: A scene composed of only an evolved surface in blue and a camera. Right: The rendering of the scene from the camera.

Examples of Ray Tracing Simple Objects and Evolved Surfaces

Putting together the material from previous sections we can now raytrace both simple objects and evolved surfaces. Figure 4.3 shows an example of simple objects, spheres, planes and disks being rendered. Figure 4.16 shows an example of an evolved surface being rendered on its own. The class of surfaces that are able to be generated with the interpolation of circles is large and Figure 4.17 shows a more unusual surface being rendered in a scene with a sphere and a plane.

**Fig. 4.17:** Left: A scene composed of a ground plane in cyan, an evolved surface in blue and a sphere in red. Right: The rendering of the scene from the camera.

Meshing Evolved Surfaces

Most graphics pipelines in modern computers use triangular meshes with some form of interpolation of vertex normals for approximating the look of curved surfaces. In light of this it is clearly desirable to be able to convert from an explicitly parameterised evolved surface to a mesh approximation of that surface.

To produce a mesh approximation we first need to generate a set of points that are in some sense evenly spaced and lie on the surface itself. To do this we will begin by producing a set of evenly spaced points on the first object $C_1$ and then transform these points along a small step in $\alpha$ to give a second set of points. Continuing in this way we can cover the surface entirely. An appropriate transformation for this task needs to preserve the relative spacing of the points on the objects in order to produce a good quality mesh. TRS (Translation Rotation Scaling) rotors have this property and can map circles to circles, spheres to spheres and point pairs to point pairs (these quantities are sometimes known as rounds). A TRS rotor that takes one object $C_1$ , to another, $C_2$ , can be calculated with the following process:

Calculate and the translation rotors that bring and respectively to the origin
Apply and to and respectively, bringing them to the origin and producing and
Calculate the rotation rotor $R_{12}$ between the blades $C'_1\wedge n_{\infty}$ and $C'_2 \wedge n_{\infty}$ (if and are spheres then we do not need a rotation rotor so set $R_{12} = 1$ )
Calculate the difference in scale between the objects by extracting their relative sizes
Use the scale to generate a dilation rotor $D_{12}$ that scales to the same size as
Compose the final TRS rotor $Z_{12}$ that takes to as:

$\displaystyle Z_{12} = \tilde{T_2}D_{12}R_{12}T_1.$

Armed with our transformation, we now simply need to generate a set of starting points on the first object. First consider the case of evolved circles. We can produce a set of $N$ evenly spaced points on the unit circle in the $e_1$ , $e_2$ plane by $N$ successive rotations about the origin of a point $X$ lying initially at $X_0 = F(e_1)$ yielding $X_n$ for $n \in 0,\ldots, N$ i.e., for a fixed rotor $R_\theta$ :

$\displaystyle X_n = (R_\theta)^nX_0(\tilde{R}_\theta)^n$

where $\theta$ is chosen so that $N+1$

uniformly spaced points cover the whole circle. With the TRS rotor $Z_{01}$ that maps from the unit circle at the origin to the first object $C_1$

it is possible to transform our points to the first object:

$\displaystyle U_{n1} = Z_{01}X_{n}\tilde{Z}_{01}.$

Our process then becomes one of stepping sequentially through $\alpha$ from 1 to 0 in small increments of $\delta\alpha$ and transforming our points along the way. We will use $Z_\alpha$ to refer to the TRS rotor that maps $C_{\alpha}$ to $C_{\alpha-\delta\alpha}$ and so can write the $n^{\text{th}}$ point at $\alpha$ as:

$\displaystyle U_{n(\alpha-\delta\alpha)} = Z_{\alpha}U_{n\alpha}\tilde{Z}_{\alpha}.$

By doing this we have effectively constructed a mapping from a coordinate system in the 2D plane of $\alpha$ and $n$

to the surface manifold. This is useful as it lets us generate the mesh in the 2D plane of $\alpha$ and $n$

and map the vertex positions directly to 3D.

**Fig. 4.18:** A linear interpolation surface of evolved circles meshed and rendered with flat shading.
$\includegraphics[width=1.0\textwidth]{mesh_circles.PNG}$

Modern graphics engines allow users to write shaders that interpolate vertex normals in smart ways, giving the illusion of curved surfaces over flat facets. In our ray tracing experiments we have already identified how to calculate the normal to the evolved surface at any point on the surface provided that $\alpha$ is known at the point. The vertex normals are calculated using the formulae in the above sections. While Figure 4.18 shows a surface of evolved circles meshed and rendered using flat shading with ganja.js De Keninck [2020], Figure 4.15 shows a tubular KB spline surface, meshed, textured, and shaded with a smooth vertex normal interpolation scheme.

Summary and Conclusions

In this chapter we have outlined the basic workings of a CGA ray tracer that can render geometric primitives as well as more advanced interpolated surfaces defined by two circles or point-pairs and an evolution parameter, $\alpha$ . Integral to ray-tracing these evolved surfaces is the derivation of analytic intersection points and normal vectors.

REFORM: Rotor Estimation From Object Resampling and Matching

We are stuck with technology when what we really want is just stuff that works.Douglas Adams, The Salmon of Doubt

Abstract

In this chapter we tackle the problem of correspondence and rotor estimation between models composed of geometric primitives of different types. We frame this problem as searching for the rotor that takes a query model to a reference model. The situations that we consider are those in which our query model: contains additional primitives not present in the reference; is missing primitives that are present in the reference. We will also look at cases in which there are a large number of primitives per model. These are all common issues facing any SLAM-type (Simultaneous Localisation And Mapping) systems. To overcome these problems we introduce an inter-object rotor magnitude-based matching function and a subsampled iterative rotor estimation and matching algorithm. We title the finished algorithm: Rotor Estimation From Object Resampling and Matching - REFORM. REFORM builds on ideas from the RANSAC (RAndom SAmple Consensus) Fischler and Bolles [1981] and ICP (Iterative Closest Point) Besl and McKay [1992]Segal et al. [2009] algorithms and extends these to multivector correspondence. It is easily parallelisable and designed for good convergence performance with models of real objects.

Introduction

A fundamental problem in computer vision is the correspondence problem. How do we match features from one image to another? This correspondence problem also appears when dealing with 3D data; given a reference model of an object and a query model of the same object how do we match objects, identify discrepancies and extract the transformation between the models? Our reference might be, for example, a CAD model, and our query model might represent the output of fitting primitives to LIDAR data or structure-from-motion point clouds. Many authors have tackled the problem of rotor estimation between groups of pre-matched geometric objects Eide and Lasenby [2018] Valkenburg and Dorst [2011]Tingelstad and Egeland [2017]De Keninck and Dorst [2019] and others have applied conformal geometric algebra to 3D registration of point and sphere clouds Kleppe and Egeland [2018]Bayro-Corrochano and Rivera-Rovelo [2009]. In this chapter we tackle the problem of registration and rotor estimation for primitives of any grade.

The objects we work with here will be CGA objects unless explicitly stated otherwise. We will use the standard extension of the 3D geometric algebra, where our 5D CGA space is made up of the standard spatial basis vectors $\{e_i\}$ $i=1,2,3$ , plus two additional basis vectors, $e$ and $\bar{e}$ with signatures, $e^2=1$ , $\bar{e}^2=-1$ . Two null vectors can therefore be defined as: $n_{\infty}=e + \bar{e}$ and $n_0 = \frac{e - \bar{e}}{2}$ . The mapping of a 3D vector $x$ to its conformal representation $X$ is given by $X = F(x) = \frac{1}{2}(x^2 n_{\infty} + 2x - 2n_0)$ .

Proximity-based matching

**Fig. 5.1:** Black: 22 lines extracted from a CAD model of a table. Red: A transformation of the original model.

**Fig. 5.2:** Using a direct proximity match between objects in the example scene, the green lines are correctly matched, the black lines are incorrectly matched. In this case the method produces 11 out of 22 correct matches.

Our first attempt at matching models made from a collection of geometric objects comes simply from considering their locality in space. For cases in which our query model is a small displacement (where displacement here will refer to rotation and translation) from the reference model, we would expect that simply assigning each object in the query model to its closest object in the reference model would give us a good number of correct matches.

Several authors have proposed cost functions between objects Valkenburg and Dorst [2011] Tingelstad and Egeland [2017], and while many of these are extremely effective for extracting motors between circles and other round elements, they tend to fail to extract the transformation between parallel lines and planes. To counteract this problem we choose the cost function described in Eide and Lasenby [2018] (the properties of this cost function are further explored in Eide and Lasenby [2018]).

Consider first two arbitrary objects in 3D space represented as $O_i$ and $O_j$ in our conformal model. As in Lasenby et al. [2018] we will extract the rotor $R_{ij}$ that takes one object $O_i$ to another $O_j$ . Note that the objects will have an orientation (sign), and the rotor extraction will be orientation dependent. Once we have our rotor $R$ between our conformal objects, the next step is to use this rotor to define a cost $C$ as a function of this rotor:

$\displaystyle C(R) = \langle (R-1)(\tilde{R} - 1)\rangle _0 + \langle (R\cdot e)(R\cdot e)\,{\tilde{}} \,\,\rangle _0$

(61)

where $\langle X \rangle _r$ indicates the $r$

-grade part of $X$

. Equipped with this idea of closeness of objects, for a given $i$

, a query object $O_i$

is assigned to each of the reference objects $O_j$

(i.e. this is done for all $j$

), assuming the model and query sets are spatially close. For each object pair we form the rotor, $R_{ij}$ that takes the query object to the reference object. The minimum cost assignment is then taken as the correct match, $M_i$

, for that query object

$\displaystyle M_i = \argmin_j{[C(R_{ij})]}$

Repeating this for all $i$

, we define the total cost of this specific matching by summing the costs of each object-to-object match

$\displaystyle C_{\text{total}} = \sum_i{C(R_{iM_i})}$

The lower this cost, the better the models are matched. Figure 5.1 shows an example scene constructed of two line-based models extracted from a CAD drawing, one model is in black and the other in red, the vertices of the models are also shown but are not used for matching. Figure 5.2 shows the result of performing proximity matching on the models, the lines in green are correctly matched and those in black are incorrectly matched. In this scene 11 of 22 lines are correctly matched by proximity matching.

Finding the rotor between two sets of matched objects

Given a set of matches for all object-pairs (under the assumption that the matching is correct) we need a method for finding the rotor between the two sets of objects. One technique for doing this is to optimise over our possible rotors, via minimisation of a cost function. Typically in CGA we parameterise and optimise over rotors in bivector space. Using the above cost metric it is shown in Eide and Lasenby [2018] that given correct matching we are able to perform non-linear convex optimisation and produce the correct rotation and displacement rotor. The downside of estimated gradient non-linear convex optimisation methods is that they typically require many cost function evaluations to reach the minimum, and when we have large numbers of objects in each model the optimisation can be very slow.

Here we propose an alternative algorithm, 1, based on directly using the rotors that we calculate between matched objects as part of the proximity matching procedure:

$\begin{algorithm} % latex2html id marker 5557 [H] \SetAlgoLined \KwResult{$R_e$}... ...k on itself\; } } \caption{Direct rotor estimation algorithm} \end{algorithm}$

This algorithm does not require the computing of an explicit cost function, it is heuristic driven and has not been proven to converge. In practice however we have found it to perform well. In the case of a fully correct matching, the rotor that is found, for both the non-linear optimisation algorithm and this direct algorithm, is indeed the rotor that takes our query model to our reference model. In the case of a partially incorrect initial matching, the rotor that is produced typically takes the query model closer to the reference model but does not produce the true rotor as shown in Figures 5.3 and 5.4.

**Fig. 5.3:** The blue model is the estimated transformation from the set of red lines to the set of black lines (see Figure 5.1) given the starting proximity match using the non linear optimisation method. Given an initial proximity matching (see Figure 5.2), the rotor found by non-linear optimisation still puts the models in close proximity even if the initial matching is not perfect.

**Fig. 5.4:** As in Figure 5.3, the blue model is the estimated transformation from the red to the black given the starting proximity match, in this case using the direct rotor estimation algorithm. The direct rotor estimation algorithm in practice produces rotors of similar quality to the non-linear optimiser.

Iterative matching and rotor estimation

Armed with rotor estimation techniques for correctly matched reference and query models we will move to more difficult situations. Consider the general case where the query and reference models are not in close proximity. In this situation we first make an initial guess at the object matches and estimate the rotor between the query and reference models using the methods described in the previous section. If our initial matching was not completely correct we will not estimate the correct rotor between the objects, the resultant rotor will have some error but will likely be relatively close to the true rotor. If we transform our query model by the estimated rotor we can use proximity matching between the transformed query model and the reference to get a better set of object matches. The process is then repeated so that the number of incorrect matches decreases with each iteration and the process converges. The iterative algorithm is summarised in the following:

Each object in the query model is given a match in the reference model (there are a number of ways of making this initial guess)
Calculate the rotor between the models assuming the current matches are correct, this can be done by running an optimisation algorithm to completion or by using the direct method mentioned in the previous section.
Transform the query model by applying the rotor calculated in the previous step
Each object in the transformed query model is compared to each object in the reference model, the match with the minimum cost according to our chosen cost function is accepted
If there is no change in the matches terminate the algorithm otherwise go back to step (2)

This algorithm correctly handles partially incorrect initial matching between models, and iterates towards the answer in relatively few steps. It is also deterministic, each step is a function only of the current state and it has fixed termination criteria that clearly indicate when it has completed. In its current state this algorithm is an extension of the well known Iterative Closest Point (ICP) algorithm Besl and McKay [1992]Segal et al. [2009] routinely used for point cloud registration. As with the ICP algorithm, a significant problem arises when we consider cases in which large fractions of the initial matches are incorrect, resulting in convergence to an incorrect set of correspondences. With our geometric algorithm we additionally see local minima arise when models contain many parallel lines or planes and computationally we run into trouble when models contain a very large number of geometric objects. In these cases the algorithm may fail to converge to the true rotor and instead become stuck in a local minimum even though some matches are correct. Real manufactured objects or buildings typically contain many parallel faces and lines and as such we need a way to overcome these limitations. Figure 5.5 shows an example of the previously studied scene stuck in a local minima, in this case there are 17 of 22 lines correctly matched but the algorithm will not progress further.

**Fig. 5.5:** As our model contains a lot of symmetry the iterative algorithm is prone to getting stuck in local minima. As in Figure 5.1 the black model is our reference and the blue is our estimate, the red lines are not shown for clarity. Here the blue model is at the final output of the iterative matching algorithm. 17 of 22 matches are correct but the algorithm is stuck in a local minima.

Incorporating sampling

To counteract the local minima issue, we modify our procedure to incorporate sampling in a RANSAC-like Fischler and Bolles [1981] algorithm. This particular approach is chosen as it is readily adapted to parallel processing and is well suited to handling large numbers of incorrect matches. After each matching stage in the previous algorithm we randomly and uniformly sample $m$ lots of $k$ matches. Each of these $m$ match sets then propagates through the rotor estimation algorithm and each produces a candidate rotor for the model matching and a cost associated with that rotor for these $k$ matches. The rotor produced by the sample with the minimum cost is then chosen and used to transform the entire query model. This repeats for a fixed number of iterations or until some cost threshold is reached.

The full REFORM algorithm is now summarised as follows:

Each object in the query model is given a match in the reference model (there are a number of ways of making this initial guess)
Given our matches, randomly select multiple sample subsets
For a given sampled subset calculate the rotor that leads to minimum total cost between the subset objects as in equation (5.1)
Accept the rotor from the sample that gives the minimum total cost between the subset objects
Update our query model position by applying the estimated rotor
Each object in the query model is compared to each object in the reference model, the match with the minimum cost is accepted
Check termination criteria, go back to step 2.

The disadvantage of moving to a sampling-based model is that we no longer have fixed termination criteria – just because the matches have not changed over multiple sampling and optimisation steps, does not mean they will not change as a result of the next one. On the other hand, the rotor estimation and cost calculation for each sample is independent of every other sample allowing for easy parallelisation. The subsampling also allows the algorithm to jump out of local minima by sampling correct matches whose effect would normally be swamped by the mass of incorrect matches. A CUDA implementation of the algorithm has been written, leveraging the massive parallelisation capability afforded by modern graphics cards and is incorporated in the clifford python package Hadfield et al. [sent].

**Fig. 5.6:** Sets of synthetic random lines and planes in red along with their transformation in black to be matched. REFORM handles both in the same framework and correctly extracts the rotor between them.

Matching scenes of mixed geometric primitives

3D models of objects are typically constructed from collections of geometric objects, planes, lines and points. While traditional matching techniques typically use points from meshes Gelfand et al. [2005] or points derived from the intersection of planes/lines Bosche [2011], REFORM allows us to incorporate multiple types of 3D object together into the same matching and rotation/translation estimation framework, Figure 5.6 shows an example of two matching synthetic models composed of both lines and planes, in this example REFORM handles both types of object transparently.

Conclusions

In this chapter we have presented an algorithm for registering models composed of geometric primitives. This algorithm extends the range of traditional matching and registration algorithms from point cloud only techniques to incorporate higher grade geometric objects. The solution is available in the clifford Hadfield et al. [sent] python package with both CPU and GPU implementations.

Kinematics, Dynamics and Robotics

In this Part of the Thesis we develop techniques that relate to robotics. Geometric Algebra has found relatively widespread use in the analysis of robots. GA simplifies the intersection of geometric shapes, often used in forward and inverse kinematics, and GA provides a neat framework for the manipulation of Lie algebras and groups, again often useful in forward and inverse kinematics as well as in control problems. Here we specifically look at the embedding of screw theory into Geometric Algebra and how this embedding can be used for dynamics and multi-body kinematics.

Screw Theory in Geometric Algebra for Constrained Rigid Body Dynamics

We may always depend on it that algebra, which cannot be translated into good English and sound common sense, is bad algebra.William Kingdon Clifford

Abstract

Screw Theory and Geometric Algebra (GA) are mathematical frameworks that have found wide use in the analysis of robotic mechanisms. Here we consider an embedding of screw theoretic wrenches and twists into the motor bivectors of two common GAs, the Plane-based (also known as Projective) and Conformal Geometric Algebras. We start with statics, considering the representations of forces and moments and how the products of GA map to the products of Screw Theory. Moving on from statics we construct an inertia tensor equivalent based on the concept of the principal screws of inertia and show how to transform this inertia tensor between frames of reference. We then look at the problem of geometrically constrained dynamics in two different ways, first via the familiar concept of virtual work, and secondly via a novel idea of multivector pinning between frames. Finally we consider the problem of integrating screw motions directly in the motor bivector space, describing kinematic equations for several alternative se(3) Lie algebra to Lie group mappings. Our goal in this chapter is to kill two birds with one stone: explicitly work through how Screw Theory embeds into commonly used GAs and use Screw Theoretic ideas to show the similarity between the various approaches to statics and dynamics in GA. Along the way we will focus heavily on the geometry driving the problems we encounter, with the hope that this will shed light both on both GA and Screw Theory.

Introduction

Screw Theory

Screw theory was originally developed by Sir Robert Stawell Ball (1840-1913) and described in his seminal text `A Treatise on the Theory of Screws' Ball [1900]. Screw Theory has found wide adoption in the study of 3D mechanisms as it allows a simple, unified treatment of both the rotational and translational motion of rigid bodies. The Screw Theory literature draws on multiple sources of mathematics but Lie theory and Projective Geometry feature heavily. Modern Screw Theory appears in several different forms and is referred to by various different names by different authors. From the `spatial vector algebra' of Roy Featherstone Featherstone [2008,2010] to the `Linear Line Complexes' of Helmutt Pottmann Pottmann and Wallner [2001] the diversity of terminology reflects the diversity of fields in which Screw Theory has found applications and mathematical foundations. This chapter is an attempt to articulate a particularly elegant embedding of Screw Theory into the coordinate free language of Geometric Algebra (GA) and how this novel embedding might allow us to push the boundaries of Screw Theory in new directions.

CGA

Conformal Geometric Algebra (CGA) is a real Clifford algebra with signature $Cl(4,1,0)$

. It is very popular for its embedding of conformal transformations (and so also the Euclidean transformations) as well as its blade representations of many geometric primitives. CGA is especially well suited for robotics as the direct circle and sphere intersections appear in many forward and inverse kinematic problems while the motions involved are typically Euclidean Hadfield and Wieser [2020]; Hadfield et al. [2020]. CGA was first suggested by Hestenes Hestenes et al. [1985] and indeed the same author made the first steps in uniting it with Screw Theory Hestenes and Fasse [2002]; Hestenes [2010]. CGA has since been applied to a huge variety of problems, for more information about this algebra the reader is directed to the excellent book of Dorst, Fontijne and Mann Dorst et al. [2007].

PGA

While CGA has been a very successful idea it is by no means the most computationally efficient way to embed rigid body dynamics into GA. In many applications a practitioner is willing to lose some of the expressiveness and mathematical niceties of CGA to gain an algebra with fewer elements that nonetheless does what they desire, and crucially, does it computationally faster. The Plane-based (or Projective) Geometric Algebra (PGA) aims to strike a balance between efficiency and expressiveness. With a signature of $Cl(3, 0, 1)$

it contains the so called `flat' elements of CGA, flat points, direction bivectors, lines, planes and the rigid body rotors, but loses point-pairs, circles, spheres, dilation and inversion rotors as first class citizens of the algebra. The use of the degenerate metric also necessitates the use of an alternate dual map rather than a simple multiplication by the pseudoscalar, in practice this dualisation operation can be implemented by the right complement extended over the canonical basis vectors by linearity. For more information on PGA the reader is directed towards the 2019 SIGGRAPH course by Gunn and De Keninck Gunn and De Keninck [2019a,b] and Leo Dorst's `A Guided Tour to the Plane-Based Geometric Algebra PGA' Dorst [2020].

Forces, moments and static equilibrium

While this chapter is ostensibly about dynamics we will first cover force representations and static equilibrium to make sure we are developing sensible mathematics before we put things in motion.

For a rigid body to be in static equilibrium certain conditions have to be met. Specifically, there can be no net moment (torque) and no net force acting on it.

What is a force?

We will now do a quick survey of the different force representations, and in doing so, think about what an effective force representation requires.

3D vector representations

In many approaches to mechanics, forces are formalised as a simple 3D vector aligned with the direction of the force and a magnitude equal to the intensity of the force. In order to calculate the resultant force on a rigid body acted upon by many incident forces we take the vector sum of the force vectors acting on the body. To represent a moment in the same formalism we use a 3D vector parallel to the axis of the moment with magnitude equal to its intensity and with the convention that a positive turning force acts counter-clockwise about this axis. To calculate the moment of a force about a specific point on the body we need to combine the incident force with its incident position. If we label the force vector $f$

, its incident point $a$

and the point about which we aim to calculate the moment as $p$

then we can calculate the moment $b$

with the cross product:

$\displaystyle b = (a-p) \times f .$

6D representations

The 3D vector representation is efficient and easy to understand, however it does not encode all the information that exists in our idea of a force. Specifically, it does not inherently include any localisation information about the force. When we want to calculate the moment due to the force about a specific point we need to include the additional information of a point on the line of action of the force. Therefore to be able to use this representation effectively we actually require two 3D vectors, one the 3D force vector $f$

and one determining a point $a$

through which the force acts. We can stack these two 3D vectors on top of each other to make a 6D vector which we will call $\mathcal{F}'$ :

$\displaystyle \mathcal{F}' = \begin{bmatrix} f \\ a \end{bmatrix}{}.$

This object $\mathcal{F}'$ now contains all the information required to use the force in applications. It is however not unique. We could choose multiple values of $a$

all of which lie on the line and they would all give different 6D representations of the force while still fundamentally behaving the same way physically. To solve this problem we need to instead use a representation of the force that is unique. If we change our representation to:

$\displaystyle \mathcal{F} = \begin{bmatrix} f \\ a\times f \end{bmatrix}{}$

then we get such a formulation, as all points on the line of action of the force will have the same result for $a\times f$ . This 6D representation is known as the Plücker coordinates of the line Pottmann and Wallner [2001].

In contrast to a force, a moment is not localised to a specific point or line in the world, in the terminology of physics and engineering this type of non-localised vector is known as a free vector. While respecting their non-localised nature we would like to be able to put the 3D moments in the same 6D box as forces to be able to keep the representations consistent and maybe provide some meaningful operations on them. When putting our 3D moment vector into the 6D box we are left with several different options. One option would be to put the moment vector in the top 3 positions of the 6D vector, such that they align with the $f$ component of the force, or maybe we should put them in the lower 3 to make them align with the $a\times f$ part of the force representation.

Initially, which choice to make seems non-obvious. To help us gain some insight we consider the problem of two anti-parallel co-planar forces of equal magnitude acting on a rigid body. First, label the forces themselves:

$\displaystyle \mathcal{F}_1 = \begin{bmatrix} f \\ a_1\times f \end{bmatrix}{}, \,\, \mathcal{F}_2 = \begin{bmatrix} -f \\ a_2\times -f \end{bmatrix}{}.$

In this situation, as the forces are in opposite directions and of equal magnitude, the total resultant force on the body is $f -f = 0$

. But what about the total moment $b_t$

of the forces about a point $p$

on the body?

$\displaystyle b_t = b_1 + b_2 = (a_1-p) \times f + (a_2-p) \times -f.$

As the cross product is distributive we can expand this out and rearrange it giving:

$\displaystyle b_t = a_1\times f - a_2\times f - p\times f + p\times f ,$

$\displaystyle b_t = (a_1 - a_2)\times f .$

This is an interesting result. The moment is independent of the position $p$

on the rigid body that we have taken moments about and is equal simply to the sum of the second half of our 6D force representation. This concept of two anti-parallel forces giving rise to a pure moment is the basis of the term couple or force couple from mechanical engineering. If we were to write an external moment in the form:

$\displaystyle \mathcal{B}_e = \begin{bmatrix} 0 \\ b_e \end{bmatrix}{}.$

We could exploit the linear independence of the top and bottom of the 6D vector to state our static equilibrium condition as:

$\displaystyle \sum_i \mathcal{F}_i + \sum_i \mathcal{B}_i = 0 .$

This is a neat result and gives us hope that continuing down the route of 6D representations will lead us to other interesting things.

In fact continuing down this route leads us to the concept of constructing an algebra over these 6D objects, known as the algebra of screws and more generally into the territory of the field known as Screw Theory. The basic element of such a `screw algebra' is a 6D vector made up of the 6D representation of a force line and the 6D representation of a moment with axis parallel to that force line. Or, more formally,

$\displaystyle \mathcal{W} = \begin{bmatrix} v \\ u\times v \end{bmatrix}{} + \... ...\ hv \end{bmatrix}{} = \begin{bmatrix} v \\ u\times v + hv \end{bmatrix}{}$

where

are 3D vectors and $h$

is a scalar. Chasles' theorem states that we can, in fact, decompose any 6D vector into this form. In the language of Screw Theory pure forces and pure moments are both special cases of the more general wrenches, which is the name for the 6d representation of a screw representing a combination of force and moment. Screw Theory goes on to define a couple of useful products between screws. The first product is the `reciprocal scalar product' of screws which, for two screws $S$

and

of the form

$\displaystyle \mathcal{S} = \begin{bmatrix} A \\ B \end{bmatrix}{}, \,\, \mathcal{T} = \begin{bmatrix} C \\ D \end{bmatrix}{}$

gives a scalar:

$\displaystyle \mathcal{S} \circ \mathcal{T} = A\cdot D + B\cdot C$

where here the $\cdot$ operator is the standard inner product between vectors. The second product we will consider is the screw cross product, for screws $S$

and

the screw cross product gives:

$\displaystyle \mathcal{S} \times \mathcal{T} = \begin{bmatrix} A \times C \\ A \times D + B \times C \end{bmatrix}{}$

where $A \times C$ and $A \times D + B \times C$ are computed with the standard 3D cross product. The result of the cross product of two screws is itself another screw.

Representations of wrenches in CGA and PGA

So far we have dealt entirely with the realm of forces and moments and their representation in 6D screw theory where they are known as wrenches. We have not yet touched on the main topic of the thesis, Geometric Algebra (GA). Our goal in this section will be to describe a direct mapping from our 6D screw theory representation to two of the most commonly used GAs, namely CGA and PGA.

Forces as dual lines in CGA

In the previous section on the screw representation we saw how we could describe a force as a directed line with a magnitude using a 6D vector. Let us now consider how we might go about representing a force line in CGA. As before consider a force line $\mathcal{F}$ expressed as 6D Plücker coordinates:

$\displaystyle \mathcal{F} = \begin{bmatrix} f \\ a\times f \end{bmatrix}{}.$

We could represent this same force in CGA as the outer product of two conformal points and infinity:

$\displaystyle F =$ up $\displaystyle \left(a\right)\wedge$ up $\displaystyle \left(a+f\right)\wedge n_{\infty}$

(62)

where we use the notation up $\left(\right)$ to represent the mapping of a 3DGA point to a conformal point, ie. $X =$

up $\left(x\right) = \frac{1}{2}x^2n_{\infty}+ x + n_0$ . The object $F$

has a magnitude equal to the intensity of the force:

$\displaystyle F^2 = \vert f\vert^2 .$

Now consider the dual form of this CGA line:

$\displaystyle FI_5 = fI_3 - (a\wedge f)I_3n_{\infty}.$

(63)

This looks very similar to the 6D Plücker representation, in fact if we consider it closely we can see that it is the same. First, look at the term in front, $fI_3$

. This is the 3D dual to a 3D vector, giving a bivector, specifically the Euclidean bivector orthogonal to $f$

. Now consider the second half of the formula: $-(a\wedge f)I_3n_{\infty}$ . $-(a\wedge f)I_3$ is the 3D dual to a Euclidean bivector, ie. a vector equal to $a\times f$ , just like the lower 3 slots of the 6D screw representation. As $(a\wedge f)I_3$ is a Euclidean vector this makes $-(a\wedge f)I_3n_{\infty}$ the form of a CGA `direction bivector'. An important thing to note here is that the front and back part of this formula are linearly independent, just like in the 6D screw representation, ie.

$\displaystyle \lambda_1 F_1I_5 + \lambda_2 F_2I_5 = (\lambda_1 f_1+\lambda_2 f_2)I_3 - (\lambda_1 a_1\wedge f_1 + \lambda_2 a_2\wedge f_2)I_3n_{\infty}.$

Again as with the 6D representation we can consider two anti-parallel forces and extract the representation of a force couple or moment. Setting $\lambda_1 = \lambda_2$ and $f_1 = f = -f_2$

in the above leaves us with:

$\displaystyle F_1I_5 + F_2I_5 = -(a_1\wedge f - a_2\wedge f)I_3n_{\infty}$

which we could re-write as:

$\displaystyle B = bn_{\infty}$

where

is a 3D vector. As previously mentioned this is in the form of a CGA `direction bivector'. These bivectors have the interesting property of being invariant to the action of translation rotors, effectively making them free vectors in the physics sense. Physically, two anti-parallel forces create a force couple, a pure moment, and so we will take objects of the form $bn_{\infty}$ to be representations of moments in our scheme. This seems apt as, physically, a pure moment is often thought of as a free vector.

Let us now look at the products of GA. Consider two objects of the form:

$\displaystyle S = aI_3 - bn_{\infty}, \,\,\,\, T = cI_3 - dn_{\infty}.$

The geometric product between them gives:

$\displaystyle ST = aI_3cI_3 - bn_{\infty}cI_3 - aI_3dn_{\infty}$

$\displaystyle = -ac - (bc + ad)I_3n_{\infty}.$

In general this is a mixed grade object, if we expand out the geometric product between vectors it will give us some insight what these grades signify geometrically:

$\displaystyle -ac - (bc + ad)I_3n_{\infty}= [-a\cdot c -a\wedge c] - [(b\cdot c + a \cdot d + b\wedge c + a\wedge d)I_3n_{\infty}].$

Now collect the terms by grade:

$\displaystyle ST = [-a\cdot c]+ [-a\wedge c -(b\wedge c + a\wedge d)I_3n_{\infty}] + [-(b\cdot c + a \cdot d )I_3n_{\infty}],$

$\displaystyle ST = \langle ST\rangle _0 + \langle ST\rangle _2 + \langle ST\rangle _4$

where

$\displaystyle \langle ST\rangle _0 = -a\cdot c ,$

$\displaystyle \langle ST\rangle _2 = -a\wedge c - (b\wedge c + a\wedge d)I_3n_{\infty},$

$\displaystyle \langle ST\rangle _4 = -(b\cdot c + a \cdot d )I_3n_{\infty}.$

From this grade-based breakdown it is quite easy to see how the different parts of the result relate to the various products of Screw Theory. First, the reciprocal product of screws, which in our 6D representation produced a scalar. In our CGA formulation this scalar maps to the coefficient of $I_3n_{\infty}$ in the result of our geometric product, ie.

$\displaystyle \langle ST\rangle _4 = -(\mathcal{S}\circ \mathcal{T})I_3n_{\infty}.$

and

are in fact bivectors we can also write this in terms of the outer product:

$\displaystyle S\wedge T = -(\mathcal{S}\circ \mathcal{T})I_3n_{\infty}.$

Next we consider the bivector part of the geometric product result, alongside the cross product of screws:

$\displaystyle \langle ST\rangle _2 = -a\wedge c - (b\wedge c + a\wedge d)I_3n_{\infty},$

$\displaystyle \mathcal{S} \times \mathcal{T} = \begin{bmatrix} A \times C \\ A \times D + B \times C \end{bmatrix}{} .$

If we rewrite the GA formula here:

$\displaystyle \langle ST\rangle _2 = -a\wedge c - (b\wedge c + a\wedge d)I_3n_{\infty}= mI_3 - gn_{\infty}$

then equating terms gives:

$\displaystyle m = (a\wedge c)I_3$

and

$\displaystyle g = (b\wedge c)I_3 + (a\wedge d)I_3 .$

The form of the standard vector cross product in 3D GA is, for vectors $a$

and

, given by:

$\displaystyle a\times b = -(a\wedge b)I_3 .$

This means we can represent the cross product of screws via the negative of the bivector part of the geometric product result. In fact, we can write the bivector part of the result in terms of the `commutator product' of geometric algebra which is also, unsurprisingly perhaps, written using the ` $\times$ ' notation.

$\displaystyle S\times T = \frac{1}{2}(ST - TS) = \langle ST\rangle _2 .$

For readers well versed in Screw Theory, Clifford/Geometric Algebra and Lie Theory this is a well known result. The commutator product equips the motor bivectors with a Lie bracket just as the cross product equips the screws with one. In fact we can go further into the idea of the motor bivector as an element of $se(3)$

and note that the grade 0 element of the geometric product, equal to the GA dot product between the bivectors, is proportional to the Killing form Minguzzi [2013], $K$

, of the Lie algebra ie. $\langle ST \rangle _0 = S\cdot T = -\frac{1}{4}K(\mathcal{S},\mathcal{T})$ .

Forces as lines in PGA

In PGA a line is computed as the intersection of two planes. As by default in PGA we are in a so called `inner product null space' or IPNS we perform the intersection of these two planes by the outer product. We first break this down component-wise for the intersection of two planes:

$\displaystyle P_1 = m_1 + d_1e_0 ,$

$\displaystyle P_2 = m_2 + d_2e_0 ,$

$\displaystyle L = P_1\wedge P_2 = (m_1 + d_1e_0)\wedge (m_2 + d_2e_0)$

$\displaystyle = m_1\wedge m_2 + (d_2m_1 - d_1m_2)\wedge e_0$

which we can now re-write this in terms that look more familiar:

$\displaystyle L = m_lI_3 - (a\wedge m_l)I_3e_0 .$

In this form the line looks very similar to the CGA representation of the line. The line squares to a scalar, there is a section that is a Euclidean bivector and a section that is a null bivector. In fact the only thing that we have changed is the form of the null element. In CGA we typically use the null element $n_{\infty}$ constructed from the sum of two orthogonal basis vectors, one squaring to $+1$

and one to $-1$

. Here in PGA the null element $e_0$

is itself a basis vector. For our 6D force line representation we therefore have exactly the same mapping as we did in CGA:

$\displaystyle F = fI_3 - (a\wedge f)I_3e_0$

and so moments appear as:

$\displaystyle B = be_0 .$

We would also expect the various products of Screw Theory to still appear as the various grades of the result of a geometric product between two elements of this algebra. Indeed they are the same as the CGA results, but of course replacing $I_3n_{\infty}$ with $I_3e_0$

which is now the pseudo-scalar of the algebra.

Force and moment representations in the GA literature

Different authors have taken different approaches when considering the form of forces in conformal and projective geometric algebra. In Lasenby, Lasenby and Doran's `Rigid Body Dynamics and Conformal Geometric Algebra' Lasenby et al. [2011] the conformal force representation is taken to be of the form

$\displaystyle F = f + \alpha n_{\infty}$

where

is the normal Euclidean 3D force vector and $\alpha$ is a scalar multiplier. This force formulation is then wedged with $n_{\infty}$ and combined with their equation (1.42) to define the equations of motion. When wedged with $n_{\infty}$ this force takes the form:

$\displaystyle F\wedge n_{\infty}= f\wedge n_{\infty}$

this is in the form of a direction bivector.

This paper further goes on to specify in equation (1.50) that the form of the moment bivector that a force of this type generates about a point is:

$\displaystyle M = X\wedge F.$

We can break this up into its constituent parts as follows:

$\displaystyle M = \left(\frac{1}{2}x^2n_{\infty}+ x + n_0\right)\wedge \left( f + \alpha n_{\infty}\right)$

$\displaystyle = \frac{1}{2}x^2(n_{\infty}\wedge f) + (x\wedge f) + (n_0\wedge f) + \alpha (x\wedge n_{\infty}) + \alpha(n_0\wedge n_{\infty}).$

This moment formulation is then wedged with $n_{\infty}$ and combined with their equation (1.42) to define the equations of motion. We now consider, as they do in their equation (1.55), the form of the moment wedged with $n_{\infty}$ :

$\displaystyle M\wedge n_{\infty}= \frac{1}{2}x^2(n_{\infty}\wedge f)\wedge n_{\... ...ge n_{\infty})\wedge n_{\infty}+ \alpha(n_0\wedge n_{\infty})\wedge n_{\infty}$

$\displaystyle = x\wedge f\wedge n_{\infty}+ n_0\wedge f\wedge n_{\infty}.$

Now take the dual of this quantity with respect to 5D pseudoscalar:

$\displaystyle (M\wedge n_{\infty})I_5 = fI_3 - (x\wedge f)I_3 n_{\infty}.$

This is the bivector form of a CGA line in the direction of the 3D vector $f$

and passing through the point $x$

`Rigid Body Dynamics and Conformal Geometric Algebra' Lasenby et al. [2011] therefore uses a mixture of 1-vectors and bivectors to represent forces but in the end their static equilibrium conditions are in the form of bivectors and trivectors, indeed they are left with something very similar to our formulation in section 6.2.2. The main difference in fact is related to which of the two orthogonal elements of the wrench formulation we take to be a moment and which to be a line. Lasenby et al. effectively choose the $f\wedge n_{\infty}$ part to be a force and here we choose the dual line section to be the force.

In fact conceptually our approach here of having lines as forces is the same as Charles Gunn's approach in his paper `On the Homogeneous Model of Euclidean Geometry' Gunn [2011a] and his PhD Thesis Gunn [2011b]. To make this connection more explicit consider the form of a PGA `ideal line', used to represent moments in Gunn's formulation:

$\displaystyle L = e_0(\alpha_1e_1 + \alpha_2e_2 + \alpha_3e_3) = e_0a$

where $a = \alpha_1e_1 + \alpha_2e_2 + \alpha_3e_3$ . The properties of this object are that it squares to 0, ie. it is null but that the magnitude and direction of the 3D vector $a$

are the magnitude and axis of the moment. This is identical to the $a\wedge n_{\infty}$ formulation for bivector moments that we described in section 6.2.2.

As discussed in Anthony Lasenby's `Rigid Body Dynamics in a Constant Curvature Space and the ‘1D-up’ Approach to Conformal Geometric Algebra' Lasenby [2011], both `On the Homogeneous Model of Euclidean Geometry' Gunn [2011a] and `Rigid Body Dynamics and Conformal Geometric Algebra' Lasenby et al. [2011] use the same form for rotors and generalised instantaneous velocities, known in the Screw Theory literature as twists.

Screw transformations, instantaneous twists, and the motor manifold

Time derivatives of frame transformations

Before looking at dynamics in detail, we will define the notation used and state various definitions.

We will consider a world that contains a single rigid body. A frame is rigidly attached to the rigid body and the body moves through space such that a time varying rotor $R$ will transform an arbitrary fixed point $X'$ in the body frame into the corresponding point $X$ in the world frame:

$\displaystyle X = RX'\tilde{R}.$

We will take time derivatives. As $X'$

is fixed:

$\displaystyle \dot{X} = \dot{R}X'\tilde{R} + RX'\dot{\tilde{R}}.$

Now substitute in $X'\tilde{R} = \tilde{R}X$ and $RX' = XR$

$\displaystyle \dot{X} = \dot{R}\tilde{R}X + XR\dot{\tilde{R}}.$

A rotor by definition has the property that $R\tilde{R} = 1$ . If we differentiate this constraint with respect to time

$\displaystyle \dot{R}\tilde{R} + R\dot{\tilde{R}} = 0$

which means we can write:

$\displaystyle \dot{X} = (\dot{R}\tilde{R})X - X(\dot{R}\tilde{R})$

which is twice the anti-commutator of $X$

and $\dot{R}\tilde{R}$ . If $X$

is a 1-vector and as $\dot{R}\tilde{R}$ is a bivector we can write:

$\displaystyle \dot{X} = -2X\cdot (\dot{R}\tilde{R}).$

If we now label our bivector quantity:

$\displaystyle \Omega_w = -2\dot{R}\tilde{R}$

(64)

it allows us to write:

$\displaystyle \dot{X} = -\frac{1}{2}\left(\Omega_w X - X \Omega_w\right) = X \times \Omega_w$

(65)

where $\times$ represents the commutator product. Note that in this form with the commutator product no assumptions are made of the grade or other properties of $X$

. If we do restrict $X$

to 1-vectors however, this allows us to write:

$\displaystyle \dot{X} = X \cdot \Omega_w .$

We can further re-arrange equation (6.3) to get an equation for the relationship between the rotor $R$ and its time derivative in terms of this quantity $\Omega_w$ :

$\displaystyle \dot{R} = -\frac{1}{2}\Omega_w R .$

(66)

This quantity $\Omega_w$ is actually our generalised instantaneous screw velocity, expressed in the world frame. Geometrically it is a screw and we can transform it just like any other screw between frames. We can therefore write $\Omega_w = R\Omega\tilde{R}$ and change equation (6.5) to:

$\displaystyle \dot{R} = -\frac{1}{2}R\Omega\tilde{R} R = -\frac{1}{2}R\Omega$

(67)

where $\Omega$ is the velocity bivector expressed in the body frame. For reference the reverse of this quantity is:

$\displaystyle \dot{\tilde{R}} = -\frac{1}{2}\tilde{\Omega}\tilde{R} .$

In the screw theory literature $\Omega$ is known as a `twist' or `velocity screw'. To take further time derivatives we can just use the chain rule:

$\displaystyle \ddot{R} = -\frac{1}{2}R\dot{\Omega} - \frac{1}{2}\dot{R}\Omega$

(68)

and we will also calculate the reverse for reference:

$\displaystyle \ddot{\tilde{R}} = -\frac{1}{2}\dot{\tilde{\Omega}}\tilde{R} - \frac{1}{2}\tilde{\Omega}\dot{\tilde{R}} .$

Momentum and inertia

Screw momentum

In traditional 3D dynamics formulations we specify that the resultant force is the rate of change of linear momentum and the resultant moment is the rate of change of angular momentum. In a screw formulation we can specify that, for a body under the influence of multiple external forces $W_i$

, the resultant wrench $W_r$

is the rate of change of screw momentum $\Psi$ with time:

$\displaystyle W_r = \sum W_i = \frac{\partial \Psi}{\partial t} .$

We can, of course, write this whether we are working in the 6D vector space, CGA or PGA.

Mapping from screw velocity to screw momentum

In 3D dynamics we are used to the idea of converting between linear velocity and linear momentum via multiplication or division by the mass of the rigid body. For a body of mass $m$

and linear velocity $v_l$

the linear momentum is simply $\rho_l$ :

$\displaystyle \rho_l = mv_l .$

When it comes to angular velocity, $v_a$

, and angular momentum, $\rho_a$ , however we have a more complicated relationship. In fact, for a body centred and axis aligned reference frame, the two are related by a diagonal matrix known as the inertia tensor that we label here as $\mathcal{M}_a$ :

$\displaystyle \rho_a = \mathcal{M}_av_a$

where $\mathcal{M}_a$ can be expressed in the form of $\gamma_i$ , the second moments of volume:

$\displaystyle \mathcal{M}_a = m \begin{bmatrix} \gamma_1 & 0 & 0 \\ 0 & \gamma_2 & 0 \\ 0 & 0 & \gamma_3 \\ \end{bmatrix} .$

While it is clear that our screw equivalent of the inertia tensor should also be a linear function that somehow combines the linear and rotational aspects of the above and maps between screw velocity and screw momentum it is not immediately obvious how we should go about constructing such a function. Let us write the matrix version of this linear function for now as $\mathcal{Q}$ and the GA version as $Q$ :

$\displaystyle \Psi = Q[\Omega_w] .$

To solve the problem of determining this function we will construct a little thought experiment. Imagine a rigid body that is initially at rest at the origin but which is then acted on by a wrench $W$ . The velocity of objects attached to the frame of the body are given by:

$\displaystyle \dot{X} = X \times \Omega_w = X \times Q^{-1}[\Psi] .$

By differentiating this equation we can calculate the acceleration of objects in the frame:

$\displaystyle \ddot{X} = \dot{X} \times Q^{-1}[\Psi] + X \times Q^{-1}[\dot{\Psi}] .$

(69)

In this case the body is initially at rest implying $\Omega_w = 0$ allowing us to eliminate the first term and leaving us with:

$\displaystyle \ddot{X} = X \times Q^{-1}[\dot{\Psi}] = X \times Q^{-1}[W] .$

(70)

We therefore have a direct mapping from the wrench acting on the body to the initial acceleration of objects in the frame and we can now use our physical intuition about the world to guide us to a compatible form of $Q^{-1}$ and hence $Q$

To calculate the form of this linear function $Q$ in GA we will first have to define what is known in the GA literature as a reciprocal frame. An important thing to note here is that this is different to the concept in screw theory of a screw and a twist being reciprocal Gallardo-Alvarado [2016] which we will come to later when considering virtual work and power. What we mean by a reciprocal frame here is a set of reciprocal bases $S^i$ that are matched to the motor bivectors $S_i$ such that when the GA inner product is taken between the two matched elements the result is 1 but otherwise 0, ie.:

$\begin{displaymath} S^i\cdot S_j = \begin{cases} 1 \quad i = j \\ 0 \quad i\ne j \end{cases}\end{displaymath}$

These reciprocal frame constructs are fairly common in the GA literature and are especially useful for this type of problem.

As $Q^{-1}$ is a linear function operating on a motor bivector $W$ in CGA it can be written as:

$\displaystyle Q^{-1}[W] = \sum_i \left[ (S^i\cdot W) \sum_k c_{ik} S_k \right]$

(71)

where $c_{ki}$ are scalar coefficients, $S^i$

is the reciprocal frame for the motor bivectors and $S_k$

are the motor bivectors themselves. We will take $S_k$

and

as follows:

$\displaystyle S_k = [e_1I_3,\,\, e_2I_3,\,\, e_3I_3,\,\, e_1\wedge n_{\infty},\,\, e_2\wedge n_{\infty},\,\, e_3\wedge n_{\infty}] ,$

$\displaystyle S^i = [-e_1I_3,\,\, -e_2I_3,\,\, -e_3I_3,\,\, e_1\wedge n_0,\,\, e_2\wedge n_0,\,\, e_3\wedge n_0]$

and will seek to determine $c_{ik}$ . This formulation is equivalent to a square matrix formulation of $Q^{-1}$ where the coefficients $c_{ik}$ are the elements in the matrix. What we need to do here is to find out exactly what these $c_{ik}$ coefficients have to be.

Firstly we will probe the response of the function to the action of a force passing through the centre of mass of the object, which lies at the origin. We can write such a force as a dual line in CGA as:

$\displaystyle W = FI_5 = fI_3 - (a\wedge f)I_3n_{\infty}$

where in this case $a$

is 0 as we are at the origin, giving:

$\displaystyle W = FI_5 = fI_3 .$

The mapping that embeds a 3D point $x$

in CGA is given by:

$\displaystyle X = \frac{1}{2}x^2n_{\infty}+ x + n_0$

which, after differentiation twice leads to:

$\displaystyle \ddot{X} = \frac{1}{2}(\ddot{x}x + 2\dot{x}\dot{x} + x\ddot{x})n_{\infty}+ \ddot{x} .$

As the body is initially stationary, $\dot{x}$ is zero and so we have for time $t=0$

$\displaystyle \ddot{X} = \frac{1}{2}(\ddot{x}x + x\ddot{x})n_{\infty}+ \ddot{x} = (\ddot{x} \cdot x) n_{\infty}+ \ddot{x} .$

For the case of a force of magnitude $\vert f\vert$ applied to the centre of mass of a rigid body we would expect all points on the body to accelerate linearly at $\frac{\vert f\vert}{m}$ in the direction of $f$

, or more formally $\ddot{x} = \frac{f}{m}$ for all $x$

. We can encode this expectation by substituting $\ddot{X}$ and $X$

into equation (6.9):

$\displaystyle (\ddot{x} \cdot x)n_{\infty}+ \ddot{x} = \left(\frac{1}{2}x^2n_{\infty}+ x + n_0\right) \times Q^{-1}[W]$

(72)

and we can then apply our specific case:

$\displaystyle \left(\frac{f}{m} \cdot x\right)n_{\infty}+ \frac{f}{m} = \left(\... ...+ x + n_0\right) \times \sum_i \left[ (S^i\cdot W) \sum_k c_{ik} S_k \right].$

Looking closely at

$\displaystyle \sum_i \left[ (S^i\cdot W) \sum_k c_{ik} S_k \right]$

we can see that $W=fI_3$

is a Euclidean bivector and, when dotted with each of the reciprocal frame elements in turn it is non zero only for $[-e_1I_3,\, -e_2I_3,\, -e_3I_3]$ in which cases the dot product gives only $f_i$

, the component of the force in the $e_i$

direction. In other words we could write:

$\displaystyle \sum_i \left[ (S^i\cdot W) \sum_k c_{ik} S_k \right] = \sum_j f_j \left[\sum_k c_{jk} S_k \right] = \sum_j \sum_k f_j c_{jk} S_k$

where $j\in(1,2,3)$ and $f_j = e^j\cdot f$ . We can combine this with $\ddot{X} = X \cdot \dot{\Omega}_w$ to give:

$\displaystyle \left(\frac{f}{m} \cdot x\right)n_{\infty}+ \frac{f}{m} = \left(\frac{1}{2}x^2n_{\infty}+ x + n_0\right) \cdot \sum_j \sum_k f_j c_{jk} S_k$

$\displaystyle \left(\frac{f}{m} \cdot x\right)n_{\infty}+ \frac{f}{m} = \sum_j \sum_k f_j c_{jk} \left(\frac{1}{2}x^2n_{\infty}+ x + n_0\right) \cdot S_k$

Let us now examine the terms on the right hand side of the above equation:

$\displaystyle \left(\frac{1}{2}x^2n_{\infty}+ x + n_0\right) \cdot S_k = \frac{1}{2}x^2n_{\infty}\cdot S_k + x\cdot S_k + n_0\cdot S_k$

For $S_k \in e_n\wedge n_{\infty}, n \in [1,2,3]$ we have the results that $n_0\cdot S_k = e_n$ , $n_{\infty}\cdot S_k = 0$ , $x\cdot (e_n\wedge n_{\infty}) = x_nn_{\infty}$ . For $S_k \in [e_1I_3,\,\, e_2I_3,\,\, e_3I_3]$ we have the results that $n_0\cdot S_k = 0$ , $n_{\infty}\cdot S_k = 0$ .

$\begin{multline} \left(\frac{f}{m} \cdot x\right)n_{\infty}+ \frac{f}{m} = \sum... ...j4} x_1 n_{\infty}+ c_{j5} x_2 n_{\infty}+ c_{j6} x_3 n_{\infty}) \end{multline}$

Now we can consider the above equation component-wise and extract the required transformation coefficients. Firstly, consider the coefficients of $n_{\infty}$ :

$\displaystyle \left(\frac{f}{m} \cdot x\right)n_{\infty}= \sum_j f_j(c_{j4} x_1 n_{\infty}+ c_{j5} x_2 n_{\infty}+ c_{j6} x_3 n_{\infty})$

Expand the left side into a sum of elements of the force vector:

$\displaystyle \left(\frac{f}{m} \cdot x\right)n_{\infty}= \sum_j \frac{f_jx_j}{m}n_{\infty}.$

Comparing terms we see we have arrived at the following solution:

$\displaystyle c_{14}, c_{25}, c_{36} = \frac{1}{m}$

$\displaystyle c_{24}, c_{34}, c_{15}, c_{35}, c_{16}, c_{26} = 0$

Now that we have dealt with the $n_{\infty}$ terms we can return to the Euclidean bivectors:

$\displaystyle \frac{f}{m} = \sum_j f_j (c_{j1} x \cdot (e_1I_3) + c_{j2}x \cdot (e_2I_3) + c_{j3} x \cdot (e_3I_3) + c_{j4} e_1 + c_{j5} e_2 + c_{j6} e_3) .$

We have just calculated the $c_{j4},c_{j5},c_{j6}$ terms in the right hand side of this equation and so by substituting these in we get something of the form:

$\displaystyle \frac{f}{m} = \frac{f}{m} + \sum_j f_j (c_{j1} x \cdot (e_1I_3) + c_{j2}x \cdot (e_2I_3) + c_{j3} x \cdot (e_3I_3))$

which implies that $c_{11}, c_{12}, c_{13}, c_{21}, c_{22}, c_{23}, c_{31}, c_{32}, c_{33} = 0$ . A visual representation of the argument made here can be seen in Figure 6.1.

**Fig. 6.1:** The multivector field generated by the commutator product of a twist and a field of multivector points can be visualised and provides a visual verification of the action of specific types of twist. Here we are visualising a twist in the form of a `direction bivector' which produces a uniform linear velocity field. This implies that our inertia tensor, or indeed some preprocessing step to the inertia tensor, must map from lines through the origin to direction bivectors, a fact proved in the mathematical content of this section.

**Fig. 6.2:** As in Figure 6.1 here we are visualising the multivector field generated by the commutator product of a twist and a field of points. Here the twist that is visualised is a line through the origin. The field that is generated is rotational about the line. Again as with the translational field of Figure 6.1 this implies specific constraints on our inertia tensor and shows the inherent reciprocal nature of momentum screws/wrenches and velocity twists in the Screw Theory formulation of dynamics.

So far we have dealt with all the translational elements of the transformation, however we have only prescribed 18 of the total 36 (6 $\times$ 6) degrees of freedom of the problem. To continue calculating the required form of the transformation we will now analyse the effect of a torque applied to the rigid body. Firstly we will again need the general form of the acceleration of a point due to a wrench:

$\displaystyle (\ddot{x} \cdot x) n_{\infty}+ \ddot{x} = \left(\frac{1}{2}x^2n_{\infty}+ x + n_0\right) \cdot \sum_i \sum_k (S^i\cdot W) c_{ik} S_k .$

(74)

We will now apply a moment to the rigid body, again at rest at the origin. From standard 3D kinematics we know that if we have a body rotating about its centre of mass and about an axis $\hat{v}_a$ with angular speed $\vert v_a\vert$ , ie. $v_a = \vert v_a \vert\hat{v}_a$ then, provided that the centre of mass has no linear velocity we can calculate the linear velocity of a point on the body as:

$\displaystyle \dot{x} = v_a \times x = (\mathcal{M}_a^{-1} \rho_a) \times x$

and the acceleration of a point as:

$\displaystyle \ddot{x} = \dot{v}_a \times x + v_a \times \dot{x}$

where $\times$ is the traditional cross product operation. In this thought experiment the body is initially at rest allowing us to remove the $\dot{x}$ term and leaving us with:

$\displaystyle \ddot{x} = \dot{v}_a \times x = (\mathcal{M}_a^{-1} \dot{\rho_a}) \times x = (\mathcal{M}_a^{-1} t) \times x$

where

is the 3d torque vector. The GA equivalent of the traditional 3d cross product for vectors $a$

and

is $-I_3(a\wedge b)$ . Which means we can write:

$\displaystyle (\mathcal{M}_a^{-1} \rho_a) \times x = -I_3((\mathcal{M}_a^{-1}t) \wedge x) .$

We can represent this same moment in our CGA formulation as the bivector wrench $W = tn_{\infty}$ . Again we look at:

$\displaystyle \sum_i \sum_k (S^i\cdot W) c_{ik} = \sum_i \sum_k (S^i\cdot (tn_{\infty})) c_{ik}$

and from this we note that $(S^i\cdot (tn_{\infty})) = 0$ for $S^i \in [-e_nI_3]$ and that $(S^i\cdot (tn_{\infty})) = t_{i-3}$ for $i \in [4,5,6]$ . For this case we can therefore write equation (6.13) as:

$\displaystyle (\left(-I_3((\mathcal{M}_a^{-1}t) \wedge x)\right) \cdot x) n_{\i... ..._{\infty}+ x + n_0\right) \cdot \sum_{i \in [4,5,6]} \sum_k t_{i-3} c_{ik} S_k$

The left hand side of this equation collapses simply to $\ddot{x}$ and we bring $X$

inside the summation again on the right:

$\displaystyle -I_3((\mathcal{M}_a^{-1}t) \wedge x) = \sum_{i \in [4,5,6]} \sum_k t_{i-3} c_{ik} \left(\frac{1}{2}x^2n_{\infty}+ x + n_0\right) \cdot S_k.$

We can now use the same results as noted previously to simplify the right hand side of this equation, leading to:

$\displaystyle -I_3((\mathcal{M}_a^{-1}t) \wedge x) = \sum_{i \in [4,5,6]}$	$\displaystyle t_{i-3} (c_{i1} x \cdot (e_1I_3) + c_{i2}x \cdot (e_2I_3) + c_{i3} x \cdot (e_3I_3)$
	$\displaystyle + c_{i4} e_1 + c_{i5} e_2 + c_{i6} e_3 + c_{i4} x_1 n_{\infty}+ c_{i5} x_2 n_{\infty}+ c_{i6} x_3 n_{\infty})$

Equating like terms eliminates all the $n_{\infty}$ terms, ie. $c_{i4}, c_{i5}, c_{i6} = 0$ for $i \in [4,5,6]$ . This leaves us with:

$\displaystyle -I_3((\mathcal{M}_a^{-1}t) \wedge x) = \sum_{i \in [4,5,6]} t_{i-3} (c_{i1} x \cdot (e_1I_3) + c_{i2}x \cdot (e_2I_3) + c_{i3} x \cdot (e_3I_3)) .$

We can break up the left side of this component-wise:

$\displaystyle \sum_{i \in [4,5,6]} -I_3((\mathcal{M}_a^{-1}t_{i-3}e_{i-3}) \wed... ...} (c_{i1} x \cdot (e_1I_3) + c_{i2}x \cdot (e_2I_3) + c_{i3} x \cdot (e_3I_3))$

Now assume that our torque vector is aligned with one of the principal axes of the rigid body, say $e_1$

, this implies $t_{i-3}$ is only non-zero for $i=4$

$\displaystyle -I_3((t_{1}\mathcal{M}_a^{-1}e_{1}) \wedge x) = t_{1} (c_{41} x \cdot (e_1I_3) + c_{42}x \cdot (e_2I_3) + c_{43} x \cdot (e_3I_3))$

$\displaystyle \frac{t_{1}}{m\gamma_1}\left(x \wedge e_{1} \right)I_3 = t_{1} (c_{41} x \cdot (e_1I_3) + c_{42}x \cdot (e_2I_3) + c_{43} x \cdot (e_3I_3))$

Which is true because $M_a^{-1}e_1 = \frac{1}{m\gamma_1}e_1$ . Noting that $\left(x \wedge e_{1} \right)I_3 = x \cdot (e_1I_3)$ allows us to identify the $c$

parameters: $c_{41} = \frac{1}{m\gamma_1}$ , $c_{42} = 0$ , $c_{43} = 0$ .

Of course our choice of $e_1$ as the principal axis to which our torque vector aligns was arbitrary, we could equally have chosen one of the other two principal axes, $e_2$ and $e_3$ with their respective $\gamma_2$ and $\gamma_3$ . As a result of this symmetry we can directly identify our final $c$ parameters: $c_{51} = 0$ , $c_{52} = \frac{1}{m\gamma_2}$ , $c_{53} = 0$ , $c_{61} = 0$ , $c_{62} = 0$ , $c_{63} = \frac{1}{m\gamma_3}$ .

Figure 6.2 gives a graphical illustration of why this form of mapping is required for rotations and torques.

The Screw Inertia Tensor

For our 6D vector representation had we not just done the maths of the previous section we might have been tempted to stack the $v_l$ and $v_a$ on top of each other to produce a combined inertia tensor like:

$\displaystyle \Psi = \begin{bmatrix} \rho_l \\ \rho_a \end{bmatrix} = \begin... ...... & \mathcal{M}_a \end{bmatrix} \begin{bmatrix} v_l \\ v_a \end{bmatrix} .$

This however would ignore the fundamental conceptual differences in the action of the top 3 and bottom 3 elements of the 6D representation as a wrench vs as a twisting motion generator. Instead, we need to use the coefficients $c_{ik}$ of the last section to produce a matrix that performs a `flip' in the relative positions of the parts of the vectors while also applying the required scaling along each principal axis. Putting the calculated coefficients in place we can see that the matrix that comes out for $Q$

looks like this:

$\displaystyle \Psi = \begin{bmatrix} \rho_l \\ \rho_a \end{bmatrix} = \begin... ...a & 0..& 0.. & 0.. \end{bmatrix} \begin{bmatrix} v_a \\ v_l \end{bmatrix} .$

We could also achieve this `flip' effect via a diagonal matrix and a permutation matrix:

$\displaystyle \Psi = \begin{bmatrix} \rho_l \\ \rho_a \end{bmatrix}= \begin{bm... ... \\ I(3\times 3) & 0... \end{bmatrix}\begin{bmatrix} v_a \\ v_l \end{bmatrix}$

(75)

where $I(3\times 3)$ is the $3\times 3$ identity matrix.

When writing the inertia tensor in CGA it is convenient to do a little relabelling for ease of reading. The reciprocal frame of the motor bivectors in CGA is as follows:

$\displaystyle -e_1I_3,\,\, -e_2I_3,\,\, -e_3I_3,\,\, e_1\wedge n_0,\,\, e_2\wedge n_0,\,\, e_3\wedge n_0.$

Which we can break into the two groups:

$\displaystyle l_i = e_iI_3, \,\,\,\, l^i = -e_iI_3$

$\displaystyle t_i = e_i\wedge n_{\infty}, \,\,\,\, t^i = e_i\wedge n_0$

and so with these we can write the inertia tensor:

$\displaystyle \Psi = Q(\Omega) = m\sum_{i=1}^{i=3} \left[ (\Omega\cdot t^i)l_i + \gamma_i(\Omega\cdot l^i)t_i \right]$

and the inverse inertia tensor:

$\displaystyle Q^{-1}(\Psi) = \Omega = \frac{1}{m}\sum_{i=1}^{i=3}\left[ \frac{1}{\gamma_i}(\Psi\cdot t^i)l_i + (\Psi\cdot l^i)t_i \right] .$

This inertia tensor performs the same kind of `flip' that we saw for the 6D screw representation. Effectively our inertia tensor does the following mapping:

$\displaystyle e_1I_3$	$\displaystyle \,\,\rightarrow\,\, e_1\wedge n_{\infty},$	$\displaystyle e_1\wedge n_{\infty}$	$\displaystyle \,\,\rightarrow\,\, e_1I_3$
$\displaystyle e_2I_3$	$\displaystyle \,\,\rightarrow\,\, e_2\wedge n_{\infty},$	$\displaystyle e_2 \wedge n_{\infty}$	$\displaystyle \,\,\rightarrow\,\, e_2I_3$
$\displaystyle e_3I_3$	$\displaystyle \,\,\rightarrow\,\, e_3\wedge n_{\infty},$	$\displaystyle e_3\wedge n_{\infty}$	$\displaystyle \,\,\rightarrow\,\, e_3I_3$

The reciprocal frame construction of the inertia tensor that we have discussed so far works for many algebras but for degenerate metric algebras such as PGA the fact that we have an element squaring to zero means this setup does not work. Instead we need to do something a little different. The degenerate metric approach to reciprocal frames Gunn [2020] is to consider some blade which we will label $x^i$

that wedges with a given blade of magnitude $b_i$

to produce the pseudoscalar with magnitude $b_i$

ie. $b_ix_i\wedge x^i = b_iI$ . Clearly, despite us labelling it $x^i$

, this object is not quite the same as the reciprocal frame, although it allows us to perform the same function of coordinate free coefficient selection producing the magnitude $b_i$

as the scalar coefficient of the pseudoscalar. Let us now identify this pseudo-reciprocal frame for the PGA bivectors:

$\displaystyle e_1\wedge e_2$	$\displaystyle \,\,\rightarrow\,\, -e_3\wedge e_0,$	$\displaystyle e_1\wedge e_0$	$\displaystyle \,\,\rightarrow\,\, -e_2\wedge e_3$
$\displaystyle e_1\wedge e_3$	$\displaystyle \,\,\rightarrow\,\, e_2 \wedge e_0,$	$\displaystyle e_2 \wedge e_0$	$\displaystyle \,\,\rightarrow\,\, e_1\wedge e_3$
$\displaystyle e_2\wedge e_3$	$\displaystyle \,\,\rightarrow\,\, -e_1\wedge e_0,$	$\displaystyle e_3\wedge e_0$	$\displaystyle \,\,\rightarrow\,\, -e_1\wedge e_2$

Comparing the PGA pseudo-reciprocal frame mapping with that of our CGA-inertia tensor mapping it is immediately clear that they are equivalent up to a minus sign. We will now define a function to perform this PGA mapping and will call it $J$

has the following action:

$\displaystyle J(b_ix_i) = \langle b_ix_i\wedge x^i\rangle _{e_{1230}}x^i = b_ix^i$

where the syntax $\langle A \rangle _{e_{1230}}$ returns the scalar coefficient of $e_1\wedge e_2\wedge e_3\wedge e_0$ in $A$

. We can extend this operation to combinations of basis elements by linearity so that for $X = \sum_i b_ix_i$ :

$\displaystyle X^J = J(X) = J\left(\sum_i b_ix_i\right) = \sum_i J(b_ix_i).$

As with our CGA reciprocal frame let's now write our PGA pseudo-reciprocal frame in two groups:

$\displaystyle l_i = e_iI_3, \,\,\,\, l^i = e_i\wedge e_0$

$\displaystyle t_i = e_i\wedge e_0, \,\,\,\, t^i = e_iI_3$

This means we can write our PGA inertia tensor as:

$\displaystyle \Psi = Q(\Omega) = -m\sum_{i=1}^{i=3} \left[ \langle \Omega\wedge... ...angle _{e1230}l^i + \gamma_i\langle \Omega\wedge t^i\rangle _{e1230}t^i\right]$

and its inverse inertia tensor:

$\displaystyle Q^{-1}(\Psi) = \Omega = -\frac{1}{m}\sum_{i=1}^{i=3}\left[ \fra... ...e l^i\rangle _{e1230}l^i + \langle \Psi\wedge t^i\rangle _{e1230}t^i \right] .$

We could also apply the $J$

map first to first `flip' the input and apply a component-wise scaling $A$

to the result:

$\displaystyle \Psi = Q(\Omega) = A[J(\Omega)] .$

This

map first formulation is conceptually the same as the screw formulation with a flip permutation matrix as in equation (6.14).

Motor Bivectors as the Principal Screws of Inertia

We can visualise the motor bivectors as a frame of screws attached to the origin. These motor bivectors are a version of the Principal Screws of Inertia, specifically they are the principal screws in a Plücker and Hunt sense as opposed to Ball's original formulation of the principal screws. Effectively what we are doing in the inertia tensor is considering these principal screws as wrenches and analysing the impact of them on the motion of the body.

By considering the bivectors as localised screws we can begin to build intuition about them and their properties. First of all, we will consider how they, and their reciprocal frame, transform under the action of rigid body rotors. The Euclidean bivectors are in the form of a dual line through the origin and are affected by rotors exactly as lines are. The direction type motor bivectors ( $e_i\wedge n_{\infty}$ ) are, as we mentioned previously, invariant to translation rotors but are affected by rotation rotors. For the general case of the rigid body rotor the direction bivectors therefore are affected only by the rotational aspect of the rotor. This is in keeping with the view of these bivectors as dual lines at infinity, known in the projective geometry world as `ideal' lines.

The transformation properties of the motor bivectors suggests an inroad on the problem of defining non-axis aligned inertia tensors which often comes up in analysis problems when we wish to transform the frame of a body to be about a known axis of rotation. We can phrase this specific problem as follows. Consider a body with a known, axis aligned, inertia tensor $Q_B$ and a world frame momentum $\Psi$ . We can get the velocity of the body by transforming the world frame momentum back to the body frame, applying the inverse inertia tensor, and transforming back:

$\displaystyle \Omega = RQ_B^{-1}(\tilde{R}\Psi R)\tilde{R} .$

We will wrap this $RQ_B^{-1}(\tilde{R}\Psi R)\tilde{R}$ operation up as an inertia tensor in its own right and label it $Q_W^{-1}$ , we can therefore write:

$\displaystyle \Omega = Q_W^{-1}(\Psi) = RQ_B^{-1}(\tilde{R}\Psi R)\tilde{R} .$

(76)

The problem of defining non-axis aligned inertia tensors is, given we know $Q_B^{-1}$ and $R$

, exactly what form does $Q_W^{-1}$ take?

The form of $Q_B^{-1}$ is:

$\displaystyle Q_B^{-1}(\Psi_B) = \Omega_B = \frac{1}{m}\sum_{i=1}^{i=3}\left[ \frac{1}{\gamma_i}(\Psi_B \cdot t^i)l_i + (\Psi_B \cdot l^i)t_i \right] .$

Based on the fact we can transform screws and their reciprocal frames with rotors as usual, we will guess the form of $Q_W^{-1}$ to be as follows:

$\displaystyle Q_W^{-1}(\Psi) = \Omega = \frac{1}{m}\sum_{i=1}^{i=3}\left[ \frac... ...lde{R}))(Rl_i\tilde{R}) + (\Psi \cdot (Rl^i\tilde{R}))(Rt_i\tilde{R}) \right] .$

(77)

Note this is just the same as $M_B^{-1}$ except that we have transformed all of the motor bivectors and reciprocals by $R$

, effectively transforming the frame of principal screws to be centred and aligned with the rigid body but expressed in the world frame. We can check whether our guess is correct as follows. Substitute $\Psi = (R\Psi_B\tilde{R})$ :

$\displaystyle Q_W^{-1}(\Psi) = \frac{1}{m}\sum_{i=1}^{i=3}\left[ \frac{1}{\gamm... ...ilde{R}) + ((R\Psi_B\tilde{R}) \cdot (Rl^i\tilde{R}))(Rt_i\tilde{R}) \right] .$

Noting that $(R\Psi_B\tilde{R})\cdot (Rt^i\tilde{R}) \equiv \Psi_B\cdot t^i$ and $(R\Psi_B\tilde{R})\cdot (Rl^i\tilde{R}) \equiv \Psi_B\cdot l^i$ we can write:

$\displaystyle Q_W^{-1}(\Psi) = \frac{1}{m}\sum_{i=1}^{i=3}\left[ \frac{1}{\gamm... ...(\Psi_B\cdot t^i)(Rl_i\tilde{R}) + (\Psi_B \cdot l^i)(Rt_i\tilde{R}) \right] .$

We can then take the rotor application outside of the summation by factorisation and we have arrived at the required relation of equation (6.15):

$\displaystyle Q_W^{-1}(\Psi)$	$\displaystyle = R\left(\frac{1}{m}\sum_{i=1}^{i=3}\left[ \frac{1}{\gamma_i}(\Psi_B\cdot t^i)l_i + (\Psi_B \cdot l^i)t_i) \right]\right)\tilde{R}$
	$\displaystyle = RQ_B^{-1}(\Psi_B)\tilde{R} = RQ_B^{-1}(\tilde{R}\Psi R)\tilde{R} .$

Of course by allowing the movement of the frame of the inertia tensor with the body it will no longer be constant and we therefore might also want the form of the time derivative for potential applications:

$\displaystyle \dot{Q_W^{-1}}(\Psi) + Q_W^{-1}(\dot{\Psi}) = \dot{R}Q_B^{-1}(\ti... ...tilde{R}} + RQ_B^{-1}(\dot{\tilde{R}}\Psi R + \tilde{R}\Psi \dot{R})\tilde{R} .$

(78)

In this subsection we have analysed the principal screws as CGA objects and our analysis of the translation of the inertia tensor was phrased using the transformation of the reciprocal frame. In reality we could equally have done our analysis from the PGA perspective as well, with our pseudo-reciprocal frame taking the place of the reciprocal frame and so our equivalent transformed inverse inertia tensor would appear as:

$\displaystyle Q_W^{-1}(\Psi) = -\frac{1}{m}\sum_{i=1}^{i=3}\left[ \frac{1}{\g... ... + \langle \Psi\wedge (Rt^i\tilde{R})\rangle _{e1230}(Rt^i\tilde{R}) \right] .$

Unconstrained rigid body dynamics

**Fig. 6.3:** A cuboid is simulated spinning about its 2nd principal axis of inertia while translating linearly. Due to the intermediate axis theorem small instabilities in the rotation build quickly causing rapid flips in orientation. Despite these rapid flips the linear motion of the centre of mass is unaffected. Blue: the path of the centre of mass, Green, Red, Orange: the path of several vertices on the cuboid as it undergoes a flip in orientation.

Equipped with forces, moments, momentum, velocities and inertia tensors we are now at a position where we can formulate the equations of motion and simulate them. We will start by considering the dynamics of an unconstrained rigid body moving under the influence of external forces and moments. We can write the state of our rigid body at a time $t$ as:

$\displaystyle Y_t = \begin{bmatrix} R_t\\ \Psi_t \end{bmatrix}$

and its first time derivative is:

$\displaystyle \dot{Y_t} = \begin{bmatrix} \dot{R_t}\\ \dot{\Psi}_t \end{bmatrix} = \begin{bmatrix} -\frac{1}{2}R_t\Omega_t\\ R_tW_{bt}\tilde{R}_t \end{bmatrix}$

where $\Psi_t$ is the momentum bivector at time $t$

expressed in the world frame and $W_{bt}$ is the resultant external wrench acting on the body expressed in the body frame. From this point on we will drop the $t$

subscript and simply state that all variables are functions of time. From our discussion in the previous section we know that we can further expand $\Omega$ using the inverse inertia tensor $Q^{-1}$ :

$\displaystyle \Omega = Q^{-1}[\tilde{R}\Psi R] .$

Re-writing the time derivative of the state with this equation for $\Omega$ gives:

$\displaystyle \dot{Y} = \begin{bmatrix} \dot{R}\\ \dot{\Psi} \end{bmatrix} = \... ...matrix} -\frac{1}{2}RQ^{-1}[\tilde{R}\Psi R]\\ RW_{b}\tilde{R} \end{bmatrix}.$

Constrained dynamics via virtual power

Unconstrained dynamics, while important, do not allow us to represent all the types of motion that we see in the real world around us. In many practical situations we are faced with the problem of constrained motion. Consider modelling a rigid body that can move dynamically under external forces but is constrained so that one or more points lie on a surface or a situation where a rigid body is constrained such that it can translate but not rotate. These are the types of problem we will attack here.

To impose a constraint on our dynamics model we will use the concept of a reaction wrench. The reaction wrench provides a combined external force and moment that acts on the rigid body in addition to the other external wrenches and, in doing so, forces the body to move in a way that respects the constraints. We will write $W_b$ as the sum of external wrenches, $S$ , plus some reaction wrench, $F$ , caused by the constraints. As we already know $S$ , all we need to calculate $F$ is the value of $W_b$ required to keep the constraints valid.

In traditional constrained dynamics work the concepts of virtual work and virtual power are widespread. In the virtual work/virtual power literature constraints are enforced by imagining several independent virtual reaction forces and moments at the constraint position and ensuring that any velocity of the body produces zero power against these forces/moments. In the screw framework that we have developed, the virtual power, $P$ , produced by a virtual world frame wrench, $T$ , when the body moves with a body frame screw velocity $\Omega$ is given by:

$\displaystyle P = \Omega\wedge(\tilde{R}TR)$

and is of the form

$\displaystyle P = pI_3n_{\infty}$

where

is a virtual scalar power. Differentiating this gives:

$\displaystyle \dot{P} = \dot{\Omega}\wedge(\tilde{R}TR) + \Omega\wedge\left(\dot{\tilde{R}}TR + \tilde{R}(\dot{T}R + T\dot{R})\right) .$

We can now substitute in our dynamics equation for $\dot{\Omega}$ :

$\displaystyle \dot{\Omega} = Q^{-1}[\dot{\tilde{R}}\Psi R + W_b + \tilde{R}\Psi\dot{R}] ,$

$\displaystyle \dot{P} = Q^{-1}[\dot{\tilde{R}}\Psi R + W_b + \tilde{R}\Psi\dot{... ...+ \Omega\wedge\left(\dot{\tilde{R}}TR + \tilde{R}(\dot{T}R + T\dot{R})\right).$

Setting the virtual power and rate of change of virtual power to 0 gives us the virtual power condition for our constraint. First:

$\displaystyle 0 = \Omega\wedge(\tilde{R}TR)$

tells us that the virtual wrench must be parallel to the screw velocity. Setting the rate of change of virtual power to be zero allows us to write:

$\displaystyle Q^{-1}[\dot{\tilde{R}}\Psi R + W_b + \tilde{R}\Psi\dot{R}]\wedge(... ... - \Omega\wedge\left(\dot{\tilde{R}}TR + \tilde{R}(\dot{T}R + T\dot{R})\right)$

$\displaystyle Q^{-1}[W_b]\wedge(\tilde{R}TR) + Q^{-1}[\dot{\tilde{R}}\Psi R + \... ... - \Omega\wedge\left(\dot{\tilde{R}}TR + \tilde{R}(\dot{T}R + T\dot{R})\right)$

$\displaystyle Q^{-1}[W_b]\wedge(\tilde{R}TR) = - Q^{-1}[\dot{\tilde{R}}\Psi R +... ... - \Omega\wedge\left(\dot{\tilde{R}}TR + \tilde{R}(\dot{T}R + T\dot{R})\right)$

is a static constraint we can specify that $\dot{T} = 0$ leaving us with:

$\displaystyle Q^{-1}[W_b]\wedge(\tilde{R}TR) = - Q^{-1}[\dot{\tilde{R}}\Psi R +... ...tilde{R}TR) - \Omega\wedge\left(\dot{\tilde{R}}TR + \tilde{R}T\dot{R}\right) .$

Which can again be solved for $W_b$

and hence

If we specify the way that $T$ varies with time we can add curved surface constraints. Consider a situation in which a rigid body is constrained such that one point $A$ always touches a sphere centred at point $V$ . Given the point is always touching the sphere we know that $T$ must always be parallel to the line joining $A$ and $V$ , we would therefore write:

$\displaystyle T = A\wedge V \wedge n_{\infty}.$

Taking a time derivative of this we see:

$\displaystyle \dot{T} = \dot{A}\wedge V \wedge n_{\infty}.$

is driven by the rotor $R$

, ie:

$\displaystyle A = RA_0\tilde{R}$

we get:

$\displaystyle \dot{A} = \dot{R}A_0\tilde{R} + RA_0\dot{\tilde{R}}$

and so:

$\displaystyle \dot{T} = (\dot{R}A_0\tilde{R} + RA_0\dot{\tilde{R}})\wedge V \wedge n_{\infty}.$

We can then directly substitute this into:

$\displaystyle Q^{-1}[W_b]\wedge(\tilde{R}TR) = - Q^{-1}[\dot{\tilde{R}}\Psi R +... ... - \Omega\wedge\left(\dot{\tilde{R}}TR + \tilde{R}(\dot{T}R + T\dot{R})\right)$

and so calculate $W_b$

. To constrain this same point to a circle we would add an additional planar constraint, ie. the point must lie on the plane in which the circle lies and on the sphere of which the circle is the equator.

Constrained dynamics by pinned multivectors

Consider a geometric primitive represented by multivector $U$ in the body frame and the same geometric primitive represented by multivector $V$ when expressed in the world frame. These two multivectors can be related by the rotor $R$ :

$\displaystyle V = RU\tilde{R}$

or equivalently:

$\displaystyle U = \tilde{R}VR .$

Taking first and second derivatives gives us the expressions:

$\displaystyle \dot{U} = \dot{\tilde{R}}VR + \tilde{R}(\dot{V}R + V\dot{R}) ,$

(79)

$\displaystyle \ddot{U} = \ddot{\tilde{R}}VR + \dot{\tilde{R}}(\dot{V}R + V\dot{... ...}}(\dot{V}R + V\dot{R}) + \tilde{R}(\ddot{V}R + 2\dot{V}\dot{R} + V\ddot{R})$

$\displaystyle = \ddot{\tilde{R}}VR + 2\dot{\tilde{R}}(\dot{V}R + V\dot{R}) + \tilde{R}\ddot{V}R + 2\tilde{R}\dot{V}\dot{R} + \tilde{R}V\ddot{R} .$

(80)

The next step is to think about what these expressions mean physically. Essentially we have two `views' of the same object, one in body space and one in world space. For example we can imagine the $U$

is a point attached to our rigid body and $V$

is a point in the world that that point is also attached to. In a sense we are `pinning' the rigid body to $V$

by its extremity $U$

Lets consider first the case that both of these `views' of the object are fixed, ie. the position and orientation of $U$ cannot change with respect to the coordinate system of the body and the position and orientation of $V$ cannot change with respect to the origin. Mathematically we are stating that $\dot{U}, \ddot{U}, \dot{V}, \ddot{V} = 0$ . If we substitute these values into (6.19) for the time derivatives we end up with the following equation:

$\displaystyle 0 = \ddot{\tilde{R}}VR + \tilde{R}V\ddot{R} + 2\dot{\tilde{R}}V\dot{R} .$

(81)

This equation is a constraint on the second time derivative of $R$

that will ensure that $U$

and

do not vary with time. We can go a step further here and substitute our expression for $\ddot{R}$ from equation (6.7), leading to:

$\displaystyle 0 = -\frac{1}{2}\dot{\tilde{\Omega}}\tilde{R}VR - \frac{1}{2}\til... ...{\Omega} - \frac{1}{2}\tilde{R}V \dot{R}\Omega + 2\dot{\tilde{R}}V\dot{R} .$

Now if we substitute in $\dot{R} = -\frac{1}{2}R\Omega$ and $\dot{\tilde{R}} = -\frac{1}{2}\tilde{\Omega}\tilde{R}$ :

$\displaystyle 0 = -\frac{1}{2}\dot{\tilde{\Omega}}\tilde{R}VR + \frac{1}{4}\til... ...c{1}{4}\tilde{R}VR\Omega\Omega + \frac{1}{2}\tilde{\Omega}\tilde{R}VR\Omega .$

Simplify and gather, substituting $U = \tilde{R}VR$ :

$\displaystyle 0 = -\frac{1}{2}\dot{\tilde{\Omega}}U - \frac{1}{2}U\dot{\Omega}... ...ilde{\Omega}^2U + \frac{1}{4}U\Omega^2 + \frac{1}{2}\tilde{\Omega}U\Omega .$

Now separate the terms with $\dot{\Omega}$ :

$\displaystyle \frac{1}{2}\dot{\tilde{\Omega}}U + \frac{1}{2}U\dot{\Omega} = ... ...\tilde{\Omega}^2U + \frac{1}{4}U\Omega^2 + \frac{1}{2}\tilde{\Omega}U\Omega .$

As $\Omega, \dot{\Omega}$ are bivectors their reverse is just a negation:

$\displaystyle -\frac{1}{2}\dot{\Omega}U + \frac{1}{2}U\dot{\Omega} = \frac{1}{4}\Omega^2U + \frac{1}{4}U\Omega^2 + \frac{1}{2}\Omega U\Omega .$

(82)

So far, we have done a lot of algebra but so far appear to be no closer to calculating our reaction wrench. If we calculate an expression for $\dot{\Omega}$ however, we start to make headway towards a solution:

$\displaystyle \dot{\Omega} = Q^{-1}[\dot{\tilde{R}}\Psi R + \tilde{R}\dot{\Psi}R + \tilde{R}\Psi\dot{R}]$

using $W_b = \tilde{R}\dot{\Psi}R$ we can also write:

$\displaystyle = Q^{-1}[\dot{\tilde{R}}\Psi R + \tilde{R}\Psi\dot{R}] + Q^{-1}[W_b]$

(83)

and so we can now substitute in on the left hand side of equation (6.21) for $\dot{\Omega}$ :

$\displaystyle LHS =-\frac{1}{2}\left( Q^{-1}[\dot{\tilde{R}}\Psi R + \tilde{R}\... ...( Q^{-1}[\dot{\tilde{R}}\Psi R + \tilde{R}\Psi\dot{R}] + Q^{-1}[W_b] \right) .$

Now we separate out the terms with $W_b$

$\displaystyle = -\frac{1}{2}Q^{-1}[W_b]U + \frac{1}{2}UQ^{-1}[W_b] + \left( -\f... ...]U + \frac{1}{2}U Q^{-1}[\dot{\tilde{R}}\Psi R + \tilde{R}\Psi\dot{R}] \right)$

and take all terms not containing $W_b$

onto the right side of the equation. We now have something of the form:

$\displaystyle -\frac{1}{2}Q^{-1}[W_b]U + \frac{1}{2}UQ^{-1}[W_b] =$ Some function of $\displaystyle R, \Psi, U.$

We can rewrite this to use the commutator product:

$\displaystyle (U\times Q^{-1}[W_b]) =$ Some function of $\displaystyle R, \Psi, U.$

(84)

If we now decide to write our total bivector wrench, $W_b$

, as the sum of external wrenches, $S$

, plus some reaction wrench, $F$

, caused by the constraints:

$\displaystyle (U\times Q^{-1}[F]) = -(U\times Q^{-1}[S]) +$ Some function of $\displaystyle R, \Psi, U$

we now have a constraint expression that fixes the reaction wrench $F$

as a function of the state of the system and the forces applied to it.

For a given $R, \Psi, U$ this constraint is linear in $F$ and can be solved for $F$ so long as we provide a correct basis for the constraint wrench. An important point to make here is that this discussion has been entirely algebra agnostic. This framework works equally well for CGA, PGA or indeed many other geometric algebras, a topic that we will return to later on.

Geometric objects as constraints

Now that we have identified a means of enforcing constraints via pinned geometric primitives let us have a look at exactly what constraints are imposed by specific choices of this pinned multivector for CGA and PGA.

Point constraint

If we want to pin a specific point in our rigid body to a point in the world we can set an invariant point constraint. In this case $U$

and

are both points and, due to the rotational symmetry of a point our reaction wrench $F$

can only support translational forces and no moments. In this case $F$

has 3 degrees of freedom, corresponding to each of the translational principal screws of inertia. This point constraint can be set in both CGA and PGA.

Point-pair constraint

Consider a bivector $G$

of the form:

$\displaystyle G = A\wedge B$

where

and

are CGA points. This object is invariant only under the action of rotation rotors about an axis parallel to the line joining $A$

and

. In this case $F$

has 5 degrees of freedom corresponding to 3 translational forces and 2 moments. This point-pair constraint is specific to CGA.

Direction constraint

Consider a so called `direction' bivector in CGA of the form:

$\displaystyle D = d\wedge n_{\infty}$

where

is a 3D vector. This object is invariant under the action of all translation rotors and is invariant to rotation rotors with axis of rotation parallel to $d$

, ie. rotors of the form:

$\displaystyle R = e^{-\frac{\theta}{2} dI_3}$

where

is the pseudoscalar of 3D space. In this case $F$

only has 2 degrees of freedom corresponding to two moments with axes perpendicular to d. For a PGA equivalent of this constraint a bivector of the form:

$\displaystyle D = d\wedge e_0$

can be used to achieve the same thing.

Flat point constraint

Consider a so called `flat point' bivector in CGA of the form:

$\displaystyle D = d\wedge n_{\infty}- n_{\infty}\wedge n_0$

where

is a 3D vector. This object is invariant under the action of all rotation rotors about the point $d$

, but is not invariant to translation. Thus, under the action of the rigid body rotors it behaves like a CGA 1-vector point. For PGA this constraint can again be implemented by a standard PGA point.

Line constraint

A line is invariant to translation along the line and rotation about the line axis. Thus we would expect to be able to support reaction forces orthogonal to the line, and moments with axes orthogonal to the line axis, ie. $F$

has 4 degrees of freedom. In PGA a line is a bivector and is formed by the intersection of two planes, in CGA a line can be represented directly as the wedge of two points and $n_{\infty}$ or dually as given in equation (6.2). Both the dual and direct CGA form work fine for pinning.

**Fig. 6.4:** Left: A physical pendulum moves under the effect of gravity and with a starting linear momentum. It is constrained such that a line, coincident with one end of the pendulum shown in blue, is pinned between the body and world reference frames. The symmetry of the line leads to constrained motion along and about the line.
Right: A spinning cone is affected by gravity but is constrained such that its end point, shown in blue, does not move. Precession and nutation are observable in the movement of the centre of mass, shown in green, and a point on the rim of the cone, shown in red.

Circle constraint

A circle is invariant to rotation about the axis of the circle only. Thus we would expect to be able to support reaction forces in all directions and all moments other than the one with axis parallel to the circle axis, in other words $F$

has 5 degrees of freedom. In retrospect this should be unsurprising as the dual to a 3-vector CGA circle is an imaginary point-pair bivector, and we have already seen that point-pairs have this same form of invariance.

**Fig. 6.5:** A physical pendulum moves under gravity and is constrained such that one end of it, shown in blue, is always in contact with the surface of an object. The green trace shows the midpoint of the pendulum and the red shows the free end. Left: a sphere. Right: a circle.

Plane constraint

A plane is invariant to translation in plane and rotation about axes parallel to its normal direction. Thus it can support one direction of force and two directions of moments giving 3 degrees of freedom to $F$

. Both CGA and PGA can use planes for pinning.

Sphere constraint

A sphere is invariant to all rotation about its centre, but not to translation. Thus it acts like a point at the sphere centre and $F$

can support translation reaction forces only giving it 3 degrees of freedom. CGA can represent spheres directly as the outer product of 4 points. In PGA this type of constraint would have to be enforced by pinning a point at the centre of the sphere.

Pinning parametric multivectors paths

So far in our construction of multivector pinning constraints we have assumed that the objects we are pinning are static in both the world and body frame. When working with constrained dynamics in the real world we often want to pin parts of our rigid body to moving things in the real world, such as a manipulator attached to the moving end-point of a robot, or a flywheel fixed in a moving vehicle. Consider once again equation (6.19):

$\displaystyle \ddot{U} = \ddot{\tilde{R}}VR + 2\dot{\tilde{R}}(\dot{V}R + V\dot{R}) + \tilde{R}\ddot{V}R + 2\tilde{R}\dot{V}\dot{R} + \tilde{R}V\ddot{R}.$

In the previous section we enforced static multivector constraints by setting $\dot{U},\ddot{U},\dot{V},\ddot{V}$ to zero, rearranging to isolate the $\ddot{R}$ terms and solving the resultant linear equation for $W_b$

. Now we will relax the static constraint and consider the cases when $U, V$

are known time varying multivector functions, ie when $\dot{U},\ddot{U},\dot{V},\ddot{V} \ne 0$ .

First note we can still rearrange to separate terms in $\ddot{R}$ :

$\displaystyle \ddot{\tilde{R}}VR + \tilde{R}V\ddot{R} = - 2\dot{\tilde{R}}(\dot{V}R + V\dot{R}) - \tilde{R}\ddot{V}R - 2\tilde{R}\dot{V}\dot{R} + \ddot{U}.$

In fact, if we continue as in our previous analysis by breaking up $\ddot{R}$ as a function of $\dot{\Omega}$ and extracting $W_b$

$\displaystyle \ddot{\tilde{R}}VR + \tilde{R}V\ddot{R} = (-\frac{1}{2}\dot{\tild... ...ilde{R}})VR + \tilde{R}V(-\frac{1}{2}R\dot{\Omega} - \frac{1}{2}\dot{R}\Omega)$

$\displaystyle = -\frac{1}{2}\dot{\tilde{\Omega}}\tilde{R}VR -\frac{1}{2}\tilde... ...frac{1}{2}\tilde{\Omega}\dot{\tilde{R}}VR - \frac{1}{2}\tilde{R}V\dot{R}\Omega$

$\displaystyle = \dot{\Omega} \times (\tilde{R}VR ) - \frac{1}{2}\tilde{\Omega}\dot{\tilde{R}}VR - \frac{1}{2}\tilde{R}V\dot{R}\Omega$

$\displaystyle = Q^{-1}[W_b]\times (\tilde{R}VR) + Q^{-1}[\dot{\tilde{R}}\Psi R ... ...frac{1}{2}\tilde{\Omega}\dot{\tilde{R}}VR - \frac{1}{2}\tilde{R}V\dot{R}\Omega$

$\displaystyle = Q^{-1}[F]\times (\tilde{R}VR) + Q^{-1}[S]\times (\tilde{R}VR) +... ...frac{1}{2}\tilde{\Omega}\dot{\tilde{R}}VR - \frac{1}{2}\tilde{R}V\dot{R}\Omega$

and so:

$\displaystyle (\tilde{R}VR) \times Q^{-1}[F] = - \ddot{U} + Q^{-1}[S]\times (... ...e{R}}(\dot{V}R + V\dot{R}) + \tilde{R}\ddot{V}R + 2\tilde{R}\dot{V}\dot{R}$

$\displaystyle + Q^{-1}[\dot{\tilde{R}}\Psi R + \tilde{R}\Psi\dot{R}] \times (\t... ...ac{1}{2}\tilde{\Omega}\dot{\tilde{R}}VR - \frac{1}{2}\tilde{R}V\dot{R}\Omega.$

If we substitute $U = \tilde{R}VR$ we have eventually got to a position where:

$\displaystyle (U\times Q^{-1}[F]) =$ Some function of $\displaystyle R, \Psi, U, \dot{U}, \ddot{U}, V, \dot{V}, \ddot{V}, S .$

Again this is a linear function in $F$

and so solvable as long as it is of sufficient rank.

What this means practically is that we can set $U$ and $V$ to follow any desired path we like in their respective spaces and extract the reaction forces and moments acting on the body that are required to keep them pinned to each other.

Pinning linear functions of parametric multivector paths

In the previous two sections we dealt directly with transformations that pin static multivectors or time varying multivector paths directly to each other in space. In many practical situations what we would really like to pin is a linear function of one multivector to another. For example we could pin the outer product of a point in the body frame and a plane in the world frame to 0, effectively forcing them to be coincident without specifying anything about their relative orientation (unlike in the transformed plane invariant case). Mathematically we can express our linear function constraint as $A[\,]$ and time derivatives as:

$\displaystyle U = A\left[\tilde{R}VR\right], \,\,\,\,\,\,\,\, \dot{U} = A\left[\dot{\tilde{R}}VR + \tilde{R}(\dot{V}R + V\dot{R})\right] ,$
$\displaystyle \ddot{U} = A\left[ \ddot{\tilde{R}}VR + 2\dot{\tilde{R}}(\dot{V}R... ... + \tilde{R}\ddot{V}R + 2\tilde{R}\dot{V}\dot{R} + \tilde{R}V\ddot{R} \right] .$

Once again we can rearrange:

$\displaystyle A\left[ \ddot{\tilde{R}}VR + \tilde{R}V\ddot{R}\right] = A\left[ ... ...V\dot{R}) - \tilde{R}\ddot{V}R - 2\tilde{R}\dot{V}\dot{R} \right] + \ddot{U}$

leading to an equation of the form:

$\displaystyle A\left[(\tilde{R}VR)\times Q^{-1}[F]\right] = A\left[-((\tilde{R}... ...n of }R, \Psi, U, \dot{U}, \ddot{U}, V, \dot{V}, \ddot{V} \right] - \ddot{U} .$

Again this is linear and solvable as before. Figure 6.5 shows the simulation with the Clifford Python library Hadfield et al. [sent] of two cases in which the linear function is the outer product with one end of a physical pendulum.

Mapping Screw Velocity to Lie Algebra Velocity

Throughout this chapter so far we have represented the derivative of the state of the body on the motor manifold as $\dot{R}$ . In practice numerical integration schemes which integrate $\dot{R}$ will undoubtedly accumulate errors and so wander off the motor manifold. Depending on the application this may or may not be a problem Boyle [2017]. We can, however, nicely side step the problem by directly mapping the velocity screw to a velocity in a suitable se(3) Lie algebra that generates the rotor $R$

. If we label the generator for the current position and orientation in the Lie algebra as $\Phi$ which maps to the current rotor with a function $Q$

then we are interested in finding a function $K$

that does the following:

$\displaystyle R = Q(\Phi), \,\,\,\,\,\, \dot{\Phi} = K(\Omega, \Phi) .$

We will refer to this function $K$

as the `kinematic equation' for the given Lie algebra to Lie group mapping and will have a look at its form for a few choices of $Q$

Exponential Mapping and the Bortz Equation

One commonly used mapping for SE(3) is the exponential mapping:

$\displaystyle R = e^{\Phi}.$

Conveniently, the kinematic equation for the exponential mapping of the motor bivectors has already been derived in Candy [2012] and a screw Lie algebra version in Selig [2004] and appears again (with a corrected typo) at the end of Section 4.5 in Selig [2005]. The result in Selig [2004,2005] is derived via idempotents and nilpotents of the adjoint matrix representation of the se(3) Lie algebra but it is readily translatable into our own notation:

$\displaystyle \dot{\Phi} = \Omega_w + \frac{1}{2}\Phi\times\Omega_w + \left(\fr... ...heta\vert(\cos{\vert\theta\vert} - 1)} \right)\Phi\times [\Phi \times \Omega_w]$
$\displaystyle + \left(\frac{1}{\vert\theta\vert^4} + \frac{\vert\theta\vert + \... ...t} - 1)} \right) \Phi \times [\Phi \times [\Phi \times [\Phi \times \Omega_w]]]$	(85)

where $\vert\theta\vert = \sqrt{\Omega_w\cdot \tilde{\Omega}_w}$ .

The exponential se(3) kinematic equation is also the subject of Section 5.3 of Liam Candy's PhD thesis Candy [2012], and while the set up of the problem is certainly correct we were unable to make their equation 5.41 work in our implementations. Suspecting simply a typo somewhere in their derivations of the derivatives we can calculate the derivatives ourselves and check the final result. We start with the following setup, following mostly along the lines of Candy [2012]. Our objective is to calculate $\dot{\Phi}$ as a function of $\Phi$ and $\Omega_b$ or $\Omega_w$ . First, we note that it is possible to write a motor bivector in the form:

$\displaystyle \Phi = \alpha B + tn_{\infty}$

(86)

where

is a rotation bivector and $t$

is a 3DGA vector. We then follow Candy [2012] choosing to define two quantaties:

$\displaystyle d = 2\langle R\tilde{R}_\alpha \rangle _2\cdot n_0, \,\,\,\,\,\, \Omega_b = \Omega_\alpha + vn_{\infty}$

where:

$\displaystyle R = \exp{\frac{-\Phi}{2}}, \quad R_\alpha = \cos{\frac{\alpha}{2}} - B\sin{\frac{\alpha}{2}}$

$\displaystyle d_\parallel =\frac{1}{2}(d + BdB), \quad d_\bot =\frac{1}{2}(d - BdB)$

$\displaystyle t = t_\bot + t_\parallel, \quad t_\bot = d_\bot, \quad t_\paralle... ... \frac{\alpha}{2}} = \frac{\alpha d_\parallel R_\alpha}{2\sin\frac{\alpha}{2}}$

We can then write an expression for $\dot{\Phi}$ as the time derivative of (6.25):

$\displaystyle \dot{\Phi} = \dot{\alpha}B + \alpha\dot{B} + \dot{t}n_{\infty}.$

This is made up of the following components:

$\displaystyle \dot{\alpha} = -B\cdot\Omega_\alpha, \quad \dot{B} = \frac{1}{2}\... ...alpha + (\Omega_\alpha\cdot B)B) - \frac{1}{2}\langle B\Omega_\alpha\rangle _2$

$\displaystyle \dot{t} = \dot{t}_\bot + \dot{t}_\parallel$

where:

$\displaystyle \dot{t}_\parallel = \frac{\dot{\alpha} d_\parallel R_\alpha}{2\s... ..._\parallel R_\alpha)\frac{\dot{\alpha}\cos\frac{\alpha}{2}}{2(\cos\alpha - 1)}$

$\displaystyle \dot{t}_\bot = \dot{d}_\bot = \frac{1}{2}(\dot{d} - \dot{B}dB - B\dot{d}B - Bd\dot{B})$

$\displaystyle \dot{d} = 2\langle \dot{R}\tilde{R}_\alpha + R\tilde{\dot{R}}_\al... ...d \dot{d}_\parallel = \frac{1}{2}(\dot{d} + \dot{B}dB - B\dot{d}B - Bd\dot{B})$

$\displaystyle \dot{R}_\alpha = -\frac{1}{2}R_\alpha \Omega_\alpha, \quad \dot{R} = - \frac{1}{2}R\Omega_b$

If we restrict ourselves to so(3) these equations will produce the same answer as that of the Bortz equation Bortz [1971] familiar to practitioners from the field of strapdown inertial navigation.

Cayley Mapping

An alternative mapping to the exponential that is simple and potentially useful is the Cayley mapping Hestenes and Fasse [2002]; Tingelstad and Egeland [2018]. For small rotations and translations the Cayley mapping approximates the exponential however it diverges somewhat as we move further from the origin. Matrix versions of the Cayley map are well known and the kinematic equations for the matrix versions of this map have been studied before in the aeronautics literature Sinclair [2005]. We have failed to find any previous attempts at the GA version of the kinematic equation however we can reuse much of the logic of the matrix derivation, simply substituting transposes for tildes. Start with the expression for the mapping:

$\displaystyle R = (1 - \Phi)(1 + \Phi)^{-1} ,$

(87)

$\displaystyle R(1 + \Phi) = 1 - \Phi .$

Take the time derivative:

$\displaystyle \dot{R}(1 + \Phi) + R\dot{\Phi} = - \dot{\Phi} .$

Rearrange:

$\displaystyle (1 + R)\dot{\Phi} = -\dot{R}(1 + \Phi) = \frac{1}{2}\Omega_w R(1 + \Phi) = \frac{1}{2}\Omega_w(1 - \Phi) .$

Now we will seek a reformulation of $1 + R$

$\displaystyle 1 + R = (1 + \Phi)(1 + \Phi)^{-1} + (1 - \Phi)(1 + \Phi)^{-1}$

$\displaystyle = (1 + \Phi + 1 - \Phi)(1 + \Phi)^{-1} = 2(1 + \Phi)^{-1} .$

This allows us to write:

$\displaystyle (1 + R)\dot{\Phi} = 2(1 + \Phi)^{-1}\dot{\Phi} = \frac{1}{2}\Omega_w(1 - \Phi)$

and so we are left with:

$\displaystyle \dot{\Phi} = \frac{1}{4}(1 + \Phi)\Omega_w(1 - \Phi) .$

(88)

Outer Exponential Mapping

The final mapping that we will consider is the so called `outer exponential' mapping as presented in Tingelstad and Egeland [2018]. This mapping is defined as taking a Taylor series of the exponential but replacing geometric products with wedge products, for an algebra with maximum grade 5 or below this can be written as:

$\displaystyle R = \exp_{\wedge}(\Phi) = \frac{1 + \Phi + \frac{1}{2}\langle \Phi ^2\rangle _4}{\sqrt{1 - \langle \Phi ^2\rangle }} .$

(89)

Again we were unable to find a GA kinematic equation for this mapping in the existing literature and so present our own as follows. First we take a time derivative of the outer exponential:

$\displaystyle R = \frac{1 + \Phi + \frac{1}{2}\langle \Phi ^2\rangle _4}{\sqrt{... ...Phi ^2\rangle _4\right)\left(1 - \langle \Phi ^2\rangle \right)^{-\frac{1}{2}}$

$\displaystyle \dot{R} = \left(\dot{\Phi } + \frac{1}{2}\langle \dot{\Phi }\Phi ... ...t{\Phi }\rangle \right)\left(1 - \langle \Phi ^2\rangle \right)^{-\frac{3}{2}}$

$\displaystyle = \frac{\dot{\Phi } + \frac{1}{2}\langle \dot{\Phi }\Phi + \Phi \... ...gle \dot{\Phi }\Phi + \Phi \dot{\Phi }\rangle }{1 - \langle \Phi ^2\rangle } .$

Then, given that in (6.3) we have defined:

$\displaystyle \Omega_w = -2\dot{R}\tilde{R}= 2R\dot{\tilde{R}}$

we can therefore write:

$\displaystyle \Omega_w = \frac{-2\left(\dot{\Phi } + \frac{1}{2}\langle \dot{\P... ...gle \dot{\Phi }\Phi + \Phi \dot{\Phi }\rangle }{1 - \langle \Phi ^2\rangle } .$

Noticing that:

$\displaystyle \langle \dot{\Phi }\Phi + \Phi \dot{\Phi }\rangle \equiv 2\langle \dot{\Phi }\Phi \rangle$

$\displaystyle \langle \dot{\Phi }\Phi + \Phi \dot{\Phi }\rangle _4 \equiv 2\langle \dot{\Phi }\Phi \rangle _4$

and rewriting our equation for $\Omega_w$ as:

$\displaystyle \left(1 - \langle \Phi ^2\rangle \right)\Omega_w R = -2\left(1 - ... ...\Phi }\rangle _4 \right) - \langle \dot{\Phi }\Phi + \Phi \dot{\Phi }\rangle R$

gives us:

$\displaystyle \left(1 - \langle \Phi ^2\rangle \right)\Omega_w R = -2\left(1 - ... ...ngle \dot{\Phi }\Phi \rangle _4 \right) - 2\langle \dot{\Phi }\Phi \rangle R .$

We can now equate grade 0 elements:

$\displaystyle \langle \left(1 - \langle \Phi ^2\rangle \right)\Omega_w R \rangle = - 2\langle \dot{\Phi }{\Phi }\rangle \langle R\rangle$

and so:

$\displaystyle \langle \dot{\Phi }{\Phi }\rangle = \frac{\left(1 - \langle \Phi ^2\rangle \right)\langle \Omega_w R \rangle } {- 2\langle R\rangle } .$

We can also equate grade 2 elements:

$\displaystyle \langle \left(1 - \langle \Phi ^2\rangle \right)\Omega_w R \rangl... ...c{1}{2}}\dot{\Phi } - 2\langle \dot{\Phi }{\Phi }\rangle \langle R\rangle _2 .$

This, after some rearrangement, leaves us with the form of the kinematic equation for the outer exponential:

$\displaystyle \dot{\Phi } = \frac{ \left(1 - \langle \Phi ^2\rangle \right)\lan... ...Phi }\rangle \langle R\rangle _2}{-2(1-\langle \Phi ^2\rangle )^{\frac{1}{2}}}$

$\displaystyle = \frac{ \left(1 - \langle \Phi ^2\rangle \right)\langle \Omega_w... ...angle \frac{ \langle R\rangle _2 }{-(1-\langle \Phi ^2\rangle )^{\frac{1}{2}}}$

$\displaystyle = \frac{ \left(1 - \langle \Phi ^2\rangle \right)\langle \Omega_w... ...ngle }\frac{ \langle R\rangle _2 }{-(1-\langle \Phi ^2\rangle )^{\frac{1}{2}}}$

$\displaystyle = -\frac{1}{2}(1-\langle \Phi ^2\rangle )^{\frac{1}{2}}\langle \O... ...{1}{2}}\langle \Omega_w R \rangle \langle R\rangle _2} { 2\langle R\rangle } .$

We can write this neatly as:

$\displaystyle \dot{\Phi } = \frac{1}{2}\sqrt{1-\langle \Phi ^2\rangle } \left[-... ... \langle \Omega_w R \rangle \langle R\rangle _2} { \langle R\rangle } \right] .$

(90)

Conclusions

In this chapter we have looked at forces, moments, free and constrained dynamics in both CGA and PGA. As well as considering how to apply virtual power as a constraint mechanism in our GA formulations we have constructed a novel technique for constrained dynamics in GA via the concept of multivector pinning. While in this chapter we have only considered two algebras, CGA and PGA, the techniques are expected to work across the board for algebras with easily representable line elements and motor bivectors. Using other higher dimensional algebras such as Cl(4,4) Du et al. [2017], Cl(8,2) Easter and Hitzer [2017] or even Cl(9,6) Breuils et al. [2018] with this technique in the future should allow for easy configurations of exotic constraints such as pinning dynamic objects to the surface of quadrics.

The Kinematics of Multi-body Systems in Geometric Algebra

Wealth, beauty and fame are transient. When those are gone, little is left except the need to be useful.Gene Tierney

Abstract

Screw Theory is a framework for analysing articulated mechanisms and performing statics and dynamics calculations that has found much success in the kinematic analysis of mechanisms. In this chapter we consider the embedding of Screw Theory into another extremely powerful framework for robotics, namely Geometric Algebra (GA). We start by rederiving well known results for the accumulation of twists along kinematic chains within our GA framework before turning our attention to the analysis of kinematic pairs. We derive an elegant representation of kinematic pairs via bilinear functions of basic geometric primitives and use this to describe the most common types of robotic joints. We then address multi-body systems and, using the Delta robot as a case study, we compare the screw theoretic approach to a direct differentiation method for extracting the Jacobians of the system.

Introduction

Modern manufacturing is increasingly mechanised and automated. The drive for automation has led to a need for advanced modelling capabilities and many modern analysis frameworks have been developed as a result. Perhaps the most successful of these frameworks for the analysis for 3D mechanisms is known as Screw Theory. Screw Theory is, unsurprisingly perhaps, concerned with the study of `screws'. In its modern form Screw Theory was first described by Sir Robert Stawell Ball in his `Treatise on the Theory of Screws' Ball [1900], however the mathematical roots that underlie this remarkable field come from Projective Geometry Pottmann and Wallner [2001] and the study of Lie groups and Lie algebras Selig [2004]. More recent proponents of Screw Theory include Hunt Hunt [1991], Selig Selig [2005], Davidson Davidson and Hunt [2004], Martins Selig and Martins [2014]; Tischler et al. [2000], Featherstone Featherstone [2008], Gallardo Alvarado Gallardo-Alvarado [2016], Pottmann Pottmann and Wallner [2001], Minguzzi Minguzzi [2013], Lipkin Lipkin [2005] and indeed many others.

Alongside the development of Screw Theory, and indeed building on some of the same fundamental mathematics, we have seen the rise of the Clifford/Geometric Algebra (GA) as a robotics modelling framework Aristidou [2010]; Fu et al. [2013]; Hildenbrand et al. [2008,2019]; Kim et al. [2015]; Kleppe and Egeland [2016]; Tichý [2020]; Zamora and Bayro-Corrochano [2004]. With modern computing capabilities, the high level description of geometry that these algebras afford allows researchers elegant, concise and coordinate-free descriptions of physically intricate mechanisms and constraints. Many of the modern applications of GA are in the analysis of conformal Dorst and Valkenburg [2011] and Euclidean motions Gunn [2011b], however there have been only a few attempts to properly embed the tools of modern Screw Theory into GA Hestenes and Fasse [2002]; Tingelstad and Egeland [2018]. This chapter, along with the previous one, is an attempt to lay out the overlap between the two fields with language and ideas familiar to both Screw Theory and GA practitioners.

Geometric Algebra

In this chapter we will, by default, work with Conformal Geometric Alegbra (CGA). CGA adds two more basis vectors, $e$

and $\bar e$ , to the original basis vectors of 3D Euclidean space, giving a complete basis for the 5D space with the following signature: $e_1^2=e_2^2=e_3^2=e^2=1$

, to its corresponding CGA vector, $X$

, is given by:

$\displaystyle X = F(x) = \frac{1}{2}\left(x^2n+ 2x-\bar{n}\right) \equiv \frac{1}{2}x^2n_{\infty}+ x + n_0 .$

(91)

which is often referred to as: up $\left(x\right) = X$ . We can invert this mapping quite easily so long as we remember to normalise the CGA point such that $X\cdot n_{\infty}= -1$ :

$\displaystyle x = F^{-1}(X) = \gadown {(X)} = \frac{-(X\wedge E_0) E_0}{X\cdot n_{\infty}}$

where $E_0 = n_{\infty}\wedge n_0$ . There are many excellent expositions on CGA in the literature and so we will refrain from a lengthy introduction of all the features of the algebra in this chapter, instead simply preferring to remind the reader of immediately relevant facts about the framework as we go along. If the reader is looking for a more thorough coverage of CGA we would recommend they turn to the excellent book Geometric Algebra for Computer Science by Dorst, Fontijne and Mann Dorst et al. [2007].

One of the important sections of the CGA framework that this chapter will deal with is the set of bivectors that form the generators of rotors that perform Euclidean motion, we will refer to these bivectors as the motor bivectors. The motor bivector basis set contains 6 elements, reflecting the 6 degrees of freedom present in rigid body motion. A common choice for the motor bivector basis in CGA is the following ordered set:

$\displaystyle m_i \in \{e_1I_3, e_2I_3, e_3I_3, e_1n_{\infty}, e_2n_{\infty}, e_3n_{\infty}\}$

where

represents the 3D pseudo-scalar $e_1\wedge e_2 \wedge e_3$ . This set would have a corresponding reciprocal frame as follows:

$\displaystyle m^i \in \{-e_1I_3, -e_2I_3, -e_3I_3, e_1n_0, e_2n_0, e_3n_0\}$

such that $m_i\cdot m^j = 1$ if $i=j$

and $m_i\cdot m^j = 0$ if $i\ne j$ . The full set of motor bivectors is then a linear combination of this motor bivector basis. Introducing 6 scalar parameters $\gamma_i$ we can write a general motor bivector $T$

as:

$\displaystyle T = \sum^{6}_{i=1}\gamma_i m_i.$

An alternative and increasingly popular GA framework to work in is the Plane-based or Projective Geometric Algebra (PGA) Dorst [2020]; Gunn [2011b]. This algebra has signature Cl(3,0,1), meaning it has 3 basis vectors that square to $+1$

, 0 basis vectors that square to $-1$

and one null basis vector that squares to 0. It is a subalgebra of CGA that contains only the `flat' elements and the Euclidean rotors Hrdina et al. [2021]; Lasenby [2011]. Restricting the algebra in this way is useful for certain applications that do not require the round elements of CGA and especially as the reduced dimensionality can produce significant speedups for some numerical packages. PGA also contains the motor bivectors with its null basis vector, $e_0$

taking the place of $n_{\infty}$ . Due to its degenerate signature PGA cannot use reciprocal frames in the same way as CGA but this restriction can be neatly sidestepped via pseudo-reciprocal frames as described in Gunn [2020] and explicitly worked through for the motor bivectors in Chapter 6. While PGA is not a direct focus of this specific chapter a sharp eyed reader will note that almost all of the formulae in the Screw Theory framework developed here will work with no, or only minor, modifications in PGA.

Whichever GA you choose, the motor bivectors, when exponentiated, produce rotors that can implement rigid body motions. In general it is possible to split $T$ into two commuting bivectors, one a generator of rotational motion about a line in the world and one a generator of translational motion along that line. This combination of rotational and translational motion leads to the identification of the motor bivectors as `screws' and if we were to take the 6 scalar parameters $\gamma_i$ and arrange them in a 6x1 vector we would get the vector space that is the core subject matter of the field of Screw Theory. This chapter is one of a pair which focus on the embedding of Screw Theory into GA and is the less theoretical of the two, here we will remind readers of any information from the previous chapter when needed but will focus more on the practical realities of using Screw Theory ideas in GA.

Twists in Kinematic Chains

**Fig. 7.1:** Rigid body is shown in blue. Rotor transforms from a fixed principal axis aligned frame to an arbitrary frame fixed to the body. Rotor transforms from the arbitrary fixed frame to the world frame and rotor transforms directly from the principal axis aligned body frame to the world frame. Likewise for body shown in green.
$\includegraphics[width=0.8\textwidth]{frames.pdf}$

Consider a rigid body $i$ . The rotor $R_i$ transforms between the world frame and an arbitrary frame fixed to, but not principal axis aligned with, the body. Specifically it transforms objects from the fixed frame to the world frame. The rotor $V_i$ then transforms from a fixed principal axis aligned frame to that arbitrary fixed frame. The combined rotor that transforms directly from the principal axis aligned frame to the world frame is written $S_i$ . Consider another body, $j$ , again with a rotor $R_j$ from fixed frame to world frame, $V_j$ from principal axis aligned fixed frame to arbitrary fixed frame and $S_j$ from principal axis aligned fixed to world frame . Figure 7.1 shows a visual depiction of this. The rotor that transforms between fixed frames of limb $i$ to limb $j$ is labelled $R_{ij}$ :

$\displaystyle R_{ij}R_i = R_j.$

(92)

Using $\sim$ to represent the standard geometric algebra reversion operator we can write:

$\displaystyle R_{ij} = R_j\tilde{R_i}.$

(93)

Now imagine that these bodies form part of a jointed kinematic chain with a sequence of limbs. The position and orientation of each limb relative to the world origin is defined by the rotor $S_i$

, which can in turn be written $S_i = R_iV_i$

Each of the fixed frames attached to the limbs has a combined rotational and translational velocity in the world frame that we can represent as a motor bivector which we will label $\Omega_i$ . These motor bivectors representing combined rotational and translational velocity are known in the Screw Theory literature as a `velocity screw' or a `twist'. From standard results we know we can write the time derivative of these rotors as:

$\displaystyle \dot{S}_i = -\frac{1}{2}\Omega_iS_i.$

(94)

We can decompose this rotor into $R_i$

and

$\displaystyle \dot{S}_i = -\frac{1}{2}\Omega_iR_iV_i.$

(95)

We can also explicitly take the time derivative of $S_i$

, which simplifies as $\dot{V_i} = 0$ due to both body frames being fixed relative to one another:

$\displaystyle \dot{S}_i = \dot{R_i}V_i + R_i\dot{V_i} = \dot{R_i}V_i.$

(96)

This implies:

$\displaystyle \dot{R_i}V_i = -\frac{1}{2}\Omega_iR_iV_i.$

(97)

If we right multiply by $\tilde{V_i}$ we can see that $\Omega_i$ is the bivector velocity of both $S_i$

and

$\displaystyle \dot{R}_i = -\frac{1}{2}\Omega_iR_i.$

(98)

This should not surprise us as both body frames are fixed to one another.

Taking the time derivative of equation (7.2) leads to a formula relating the relative velocity screws of the body in the chain:

$\displaystyle \dot{R}_{ij}R_i + R_{ij}\dot{R}_i = \dot{R}_j$

(99)

$\displaystyle -\frac{1}{2}\Omega_{ij}R_{ij}R_i -\frac{1}{2}R_{ij}\Omega_iR_i = -\frac{1}{2}\Omega_jR_j$

(100)

$\displaystyle -\frac{1}{2}\Omega_{ij}R_j -\frac{1}{2}R_{ij}\Omega_i\tilde{R}_{ij}R_j = -\frac{1}{2}\Omega_jR_j$

(101)

$\displaystyle \Omega_{ij} + R_{ij}\Omega_i\tilde{R}_{ij} = \Omega_j .$

(102)

This is a potentially convenient result but there is a particular case that proves to be of special interest. Our results so far are not a function of the rotors $V_i$ at all, just of $R_i$ . We could choose any decomposition of $S_i$ into $R_i$ and $V_i$ at a given time point and our equations will still be valid. The special case we will concentrate on is to choose to instantaneously align the fixed frames with the world frame, in other words we choose at this instant: $R_i = 1$ , $V_i = S_i$ and therefore $R_{ij}=1$ . Under these conditions equation (7.12) simplifies to the following:

$\displaystyle \Omega_{ij} + \Omega_i = \Omega_j.$

(103)

This equation is particularly convenient if we are able to measure the relative velocity screw of one limb with respect to another as we can simply accumulate the relative limb velocities along the chain to arrive at the final global limb velocity screw with respect to the world frame. This is well known in traditional Screw Theory and forms the basis of many practical techniques for analysing robots.

A direct result of the definition of $\Omega$ is that the time derivative of a geometric object $X$ is given by the commutator product of the object with $\Omega$ :

$\displaystyle \dot{X} = \frac{1}{2}(X\Omega - \Omega X) = X \times \Omega.$

(104)

This is a particularly helpful result as it allows us to define geometry associated with our kinematic chains and calculate how they evolve through time as the kinematic chain moves.

Geometrically Constrained Kinematic Pairs

Consider a pair of adjacent limbs, $i$

and

, in a kinematic chain. These two limbs are connected by some form of joint and together they are known as a kinematic pair Reuleaux [1876]. Exactly how they are constrained to move relative to one another is defined by the type of joint. Let us first consider a specific type of joint, one defined by shared geometry.

In a shared geometry joint there exists a piece of geometry, which we will label $X$ , that is fixed relative to both frames in the kinematic pair. Practically what this means is that when acted on by the velocity of $i$ or $j$ this object $X$ must have the same $\dot{X}$ , or more formally:

$\displaystyle \dot{X} = X \times \Omega_i = X \times \Omega_j.$

(105)

Due to the linearity of the commutator product we can write this as:

$\displaystyle 0 = X \times (\Omega_i - \Omega_j) = -X\times \Omega_{ij}.$

(106)

This is an extremely useful result. The quantity $\Omega_{ij}$ is the relative velocity screw of one limb relative to the other and the restriction of its commutator with X being zero means that the shared geometry of the joint must be invariant to the effect of the relative velocity.

We can represent many types of joint as shared geometry joints, or combinations of shared geometry joints. In practice however we often want to represent slightly more complex joints compactly. In many commonly used GA frameworks we can represent more advanced types of joint with some form of bilinear mapping between two relevant objects, one in each frame, that we know remains constant for the joint geometry. More formally, for an object $X$ in frame $i$ and an object $Y$ in frame $j$ :

$\displaystyle B(X, Y) = H$

(107)

where

is a multivector that is constant for the joint, and in many cases is simply 0. Taking time derivatives once again gives us a linear constraint on $\Omega_i$ and $\Omega_j$ :

$\displaystyle B(\dot{X}, Y)$	$\displaystyle + B(X, \dot{Y}) = 0,$	(108)
$\displaystyle B((X\times \Omega_i),Y)$	$\displaystyle + B(X, (Y\times \Omega_j)) = 0.$	(109)

So far at no point have we specified what objects $X$

and

are, nor the form of the bilinear mapping $B$

(7.17). In practice, common forms of $B$

include the outer product:

$\displaystyle X\wedge Y = 0, \quad (X\times \Omega_i)\wedge Y + X\wedge (Y\times \Omega_j) = 0$

the inner product:

$\displaystyle X\cdot Y = 0, \quad (X\times \Omega_i)\cdot Y + X\cdot (Y\times \Omega_j) = 0$

and the meet taken with respect to a specific fixed subspace:

$\displaystyle X\vee Y = 0, \quad (X\times \Omega_i)\vee Y + X\vee (Y\times \Omega_j) = 0.$

The Geometry of Real Joints

It is all well and good to state some neat theoretical results but it is more useful to robotics practitioners to outline how to describe actual joints in our framework. A good place to start is with the most common types of kinematic pair.

Spherical Joint

A spherical joint, often known as a ball joint, is one in which two bodies are free to move about a single common fixed point. Often in practice these joints are implemented with one body designed to have a spherical cavity (known as the socket) and another with a ball on the end which is held captive within the socket. To implement a spherical joint within our shared geometry constraint we can use a sphere shared between both bodies in the kinematic pair. The sphere can be of any radius and in many practical applications (and in some GAs which do not naturally embed non-zero radius spheres as objects, such as PGA) it makes sense to use a sphere of zero radius, or, in other words, a point.

Cylindrical Joint

A cylindrical joint is a mechanism that allows two bodies to translate parallel to a shared axis and rotate about the same axis simultaneously. In practice this joint appears often as a cylindrical cuff free to slip over and around a rod. To implement the cylindrical joint within our framework we can use the geometry of a shared line. A line is invariant to rotation about its axis and translation along its axis. In CGA we can use a line either in its trivector form as the direct wedge product of two conformal points and $n_{\infty}$ , $A\wedge B\wedge n_{\infty}$ or in its bivector form $\hat{m}I_3 - (p\wedge \hat{m})I_3n_{\infty}$ which is equivalent to the PGA line formulation $\hat{m}I_3 - (p\wedge \hat{m})I_3e_0$ where $\hat{m}$ is the direction of the line and $p$

is a point on the line.

Planar Joint

A planar joint is one that allows two bodies to slide over one another in a specific plane of motion and rotate about the normal to that plane. An example of this in the real world would be a wheeled trolley that is free to move in any direction on a flat surface. To implement the planar joint within our framework we can use a shared plane between the kinematic pairs. Again in CGA this can be done with either a 4-vector plane, $A\wedge B\wedge C\wedge n_{\infty}$ , or with a dual 1-vector plane $m + dn_{\infty}$ which is again the same as the PGA plane $m + de_0$

where

is the normal to the plane and $d$

is the distance of the plane from the origin.

Revolute Joint

A revolute joint, also known as a pin joint or hinge joint, is one that allows rotation about an axis but, unlike the cylindrical joint, does not allow translation parallel to the axis. The revolute joint is one of the most common mechanical joints in use in the wild and so having a convenient mechanism to handle it is important. To implement a revolute joint within our framework with CGA we can use a shared circle. A circle is invariant only to rotation about the axis passing through its centre and normal to the plane in which it lies. The radius of this circle is again arbitrary and choices of zero radius, or radius appropriate to the mechanical proportions of the robot at hand, are likely reasonable. In CGA this circle comes in the form of a trivector $A\wedge B\wedge C$ or in the form of the bivector dual to that. In CGA a point-pair $P_1\wedge P_2$ that is aligned with the axis can also be used.

Prismatic Joint

A prismatic joint is one in which two bodies are constrained to move relative to one another parallel to a fixed line without any rotation. The prismatic joint is the first of the classic robotic joints that we cannot represent with a single geometric primitive in 3D conformal GA. Instead we are forced to use some form of compound constraint. There are several options here, two parallel line constraints, a line and a plane constraint, a line and a direction bivector constraint and two planar constraints are all constructs that would work.

Universal Joint

A universal joint is a mechanism which consists of two orthogonal revolute joints with axes that are incident at a fixed point. It is often used to transmit rotational motion between two shafts whose axes are not parallel but are incident.

We can represent a universal joint in our framework as a joint consisting of two equal radius circles. The planes of these circles must lie orthogonal to each other and both circles must have the same centre point. If we choose two circles, $X$ and $Y$ , the orthogonality constraint can be written as $X\cdot Y = \langle XY\rangle _0 = 0$ and the same centre point constraint for circles of the same radius can be written $\langle X Y \rangle _4 = 0$ Hadfield and Lasenby [2019]. We can combine these constraints neatly with the anti-commutator product

$\displaystyle X\bar{\times}Y = \frac{1}{2}(XY + YX).$

When operating on two circles, the anti-commutator product produces grade 0 and grade 4 elements only. Using this notation and taking derivatives our kinematic constraint can therefore be described as:

$\displaystyle (X\times \Omega_i)\bar{\times}Y + X\bar{\times}(Y\times \Omega_j) = 0.$

(110)

Of course we could have chosen to model this joint with two revolute joints and an additional body to represent the internals of the joint but this would introduce additional and unnecessary variables into our problem definition.

The Kinematic Constraint Matrix and the Jacobian Matrix

So far we have only looked at individual kinematic pairs, we will now address full multi-body systems. For an articulated robot with $N$

limbs we write the velocity state as:

$\displaystyle T = \begin{bmatrix} \Omega_1\\ \Omega_2\\ ... \\ \Omega_N \end{bmatrix}$

(111)

For a system with $M$

constraints, there exists an $M\times N$ matrix $C$

with linear functions as elements that embodies the combination of all the linear constraints on the velocity state. $C$

acts on the velocity state giving a result of M zero valued multivectors:

$\displaystyle 0 = CT.$

(112)

We will name $C$

the kinematic constraint matrix and readers familiar with the Screw Theory literature will recognise it as a cousin of Davies' method Davies [1981]. We can choose some ordered $L$

dimensional basis in which to represent the multivectors $\Omega_i$ in which case $T$

becomes an $NL\times 1$ dimensional vector, $C$

becomes an $ML\times NL$ dimensional matrix and the zero vector is $ML\times 1$ dimensional. To find the various coefficients of $C$

we can simply set each of the coefficients of $T$

to 1 in turn and all the rest to zero and calculate the output of each of our constraints for that input, calculating its output representation in the chosen multivector basis.

Consider a robot with velocity state $T$ , now specify some of the limbs as inputs, with known $\Omega$ and some of the limbs as outputs with unknown $\Omega$ . We can write this as follows:

$\displaystyle T = T_{k} + T_{u}$

(113)

where the subscripts $k$

and

refer to known and unknown $\Omega$ respectively. By linearity, applying the constraint matrix produces:

$\displaystyle 0 = CT = CT_{k} + CT_{u},$

$\displaystyle CT_{u} = -CT_{k} .$

(114)

This is in the standard form $Ax = b$

and can be solved with the pseudo-inverse in our chosen basis to produce a valid set of unknown limb velocities given the set of known velocities.

Designating a linear map $K$ from the input degrees of freedom, $x$ , to the known velocities of the system, and a linear map $U$ from the unknown velocities to the output degrees of freedom, $y$ , allows us to calculate a matrix that is equivalent to the Jacobian matrix of the system:

$\displaystyle T_k = K(x), \quad U(T_u) = y,$

$\displaystyle y = (-UC^+CK)x$

(115)

where

is the pseudo-inverse of $C$

The matrix $C^+C$ describes the relationship between the known and unknown twists and, so long as we can form the various components of equation (7.25) we should be able to use this Screw Theory inspired framework to calculate the Jacobian matrix for any robot whose joints can be modelled using the bilinear mappings of Section 7.5.

Case Study: The Delta Robot

The Delta robot Clavel [1990,1991] was invented in 1985 by Raymond Clavel at EPFL after being inspired by a visit to a chocolate packing factory Pessina [2012]. It has since become a particularly popular robot in industrial settings due to its good precision coupled with high speed and acceleration.

The Delta robot is a specific type of robot known as a parallel manipulator. Parallel manipulators, also known as parallel robots, are a class of robots that feature end-effectors driven by multiple underactuated parallel kinematic chains Gallardo-Alvarado [2016]; Merlet [2006]. Typically a parallel robot is designed such that all actuators remain fixed to the support structure of the robot thereby minimising the mass of the moving parts of the robot and enabling very fast accelerations. Indeed this goal of high speed/fast acceleration has been the primary driving force in the development of parallel robots for industry and today architectures such as the Delta robot are widespread in many high precision, high throughput manufacturing applications. Parallel robots, while practically very useful, are often significantly more difficult to analyse than their serial cousins due to the end-point position being a function of the configuration of multiple kinematic chains.

Here we will do a case study of the Delta robot, analysing it with Geometric Algebra, calculating Jacobians with our Screw Theory based framework, and then finally comparing our screw theory setup with a direct differentiation approach to get the Jacobians.

Geometry of a Delta Robot

**Fig. 7.2:** Left: the 3D geometry of the delta robot. Right: The geometry of a single arm in plane.

Since its inception, there have been many variants of the Delta robot Pierrot et al. [1991]. In this chapter we will assume the simple robot described in this section and illustrated in Figure 7.2. The static part of the robot is a base plate to which three motors are rigidly attached, we will assume a space in which the origin is at the centre of this plate. Each motor shaft is rigidly attached to an `upper arm' of length $l$ ; we will number each upper arm $i \in [1,2,3]$ . The connection point of the motor and upper arm will be labelled $B_i$ . The arm can only rotate in plane about the motor axis as the motor shaft and upper arm are rigidly connected. We will refer to the other end of this upper arm as the `elbow point' and will label it $W_i$ . At the elbow point each arm is rigidly attached to a central point of a horizontal rod we will refer to as the `elbow rod'. At each end of the elbow rod a ball joint connects to a `forearm' piece. The two forearm pieces for each arm are the same length and, at the other end from the elbow rod, are connected to a rigid plate that we will refer to as the end-effector plate. The point half-way between where the two forearm rods connect to the end-effector plate is labelled $X_i$ . We will label the point at the centre of the end-effector plate $Y$ . Assuming the robot is infinitely stiff, the end plate is constrained, due to this specific arrangement of the forearms, to always remain parallel to the base plate and to have its in-plane orientation fixed as well. The Delta robot is therefore a purely translational mechanism.

Calculating the Robot Pose

The first task in analysing the Delta robot is to calculate its pose for a given set of inputs. Typically we are attempting to map from the end-effector to the actuators (inverse kinematics) or from the actuators to the end-effector (forward kinematics). Along the way we are interested in the positions of the limbs and joints that define the physical structure of the robot. To solve these pose calculation problems for the Delta robot initially we will use a simplified geometry for each leg.

Inverse Kinematics

The inverse kinematic problem for the Delta robot is summarised as follows: To what angle relative to the base should we move the upper arms given we want the centre of the end-effector plate to be in a specific position in 3D space?

**Fig. 7.3:** The geometry of the inverse kinematic problem. There are three spheres, one for each arm of the robot, centred at the conformal points . Each sphere intersects with a circle centred at allowing the extraction of the conformal elbow point .

To solve this problem we need to work backwards from the 3D end-effector plate position $y$ to the motor angles $\theta_i$ considering the geometry of the robot as we go. Starting at the end-effector plate the 3D points $x_i$ are translationally offset in plane in the direction $s_i$ giving $x_i = y + r_es_i$ where $r_e$ is the radius of the end-effector plate. Due to the geometry of the robot the elbow point $W_i$ is constrained to lie on a sphere with radius equal to the length of the forearms $\rho$ centred at this point $x_i$ . We can write the dual form of this sphere as:

$\displaystyle \Sigma_i^* = X_i - \frac{1}{2}\rho^2n_{\infty}.$

The elbow point is also simultaneously constrained to lie on a circle of radius $l$

centred at the motor shaft to upper-arm joint, $B_i$

. We can represent this circle in its dual form $C_i^*$

in CGA, where $C_i$

is the intersection of a sphere of radius $l$

centred at the position $B_i$

, with dual form $\left( B_i - \frac{1}{2}l^2n_\infty \right)$ , and the plane through the origin, $B_i$

and

which has dual form $I_3(s_i\wedge e_3)$ . Here $e_3$

is the vertical unit vector and for the most common orientation of a delta robot, points vertically downward from the base plate. In CGA we calculate the intersection of objects via the `meet' operator, as both operands are in their dual form however, here we simply need an outer product:

$\displaystyle C_i^* = \left( B_i - \frac{1}{2}l^2n_\infty \right)\wedge(I_3(s_i\wedge e_3)).$

So long as $y$

is within the reachable volume of the robot there are two possible solutions for this pair of constraints. These two solutions lie at the intersection points of the sphere and circle and the `meet' operation of CGA provides us with a direct means to calculate these intersection points. The sphere and circle are in the dual form (1 and 2-vectors respectively), and so the point-pair bivector resulting from their meet is calculated as simply their outer product followed by multiplication with the 5D pseudo-scalar, $I_5$

$\displaystyle T_i = (C_i^*\wedge\Sigma_i^*)I_5.$

The desired individual solution can be extracted from this point-pair object by relying on the oriented nature of the algebra and relying on a projection operator Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]:

$\displaystyle P_i = \frac{1}{2}\left(1 + \frac{T_i}{\sqrt{T_i^2}}\right),$
$\displaystyle W_i = -\tilde{P_i}(T_i\cdot n_\infty).$

We can then convert from the CGA to the 3D vector point:

$\displaystyle w_i = \gadown {(W_i)}$

and so, with a little trigonometry we can extract the motor angles:

$\displaystyle \theta_i = \operatorname{atan2}({z_i\cdot e_3, z_i\cdot s_i}) , \hspace{0.5cm} z_i = w_i - r_bs_i$

where

is the radius of the base-plate. Figure 7.3 illustrates the geometry of the inverse kinematic problem graphically. As mentioned we have relied on the oriented nature of the algebra to extract the solution of interest from the point-pair. This solution has the elbow position being as far from the $e_3$

axis as possible and is normally the only feasible position of the elbow in a real Delta robot as the other typically causes self intersection. If the other solution is desired the same projection formulae can be used, simply substituting $\tilde{T}_i$ for $T_i$

Forward Kinematics

The forward kinematic problem is, in some sense, the opposite of the inverse kinematic one. Our goal here is to calculate the 3D vector position of the end-effector plate $y$ given the motor angles $\theta_i, i \in [1,2,3]$ .

**Fig. 7.4:** The geometry of the forward kinematic problem. Each motor connects to an upper arm at position . The upper arms end in the elbow point . Each elbow point has an associated pseudo-elbow point and forearm constraint sphere. All three constraint spheres meet at the centre of the end-point plate, .

To solve the forward kinematic problem we will consider the robot one arm at a time. For a given arm motor angle $\theta_i$ the 3D position of the elbow point $w_i$ can be calculated as:

$\displaystyle w_i = (r_b + l\cos(\theta_i))s_i + l\sin(\theta_i)e_3.$

For each arm we will now define a pseudo-elbow point, $a_i$

which is offset horizontally from the true elbow point by the radius of the end effector plate and in the direction of the origin.

$\displaystyle a_i = (r_b - r_e + l\cos(\theta_i))s_i + l\sin(\theta_i)e_3.$

The equivalent CGA point is then:

$\displaystyle A_i = \frac{1}{2}a_i^2n_\infty + a_i + n_0.$

Given the geometry of the robot, these pseudo-elbow points all lie a distance equal to the length of the robot's forearms, $\rho$ , from the centre of the end-point plate $Y$

. Geometrically these fixed distance constraints manifest themselves as spheres, which we will label $\Sigma_i$ , on which the centre of the end-point plate can lie:

$\displaystyle \Sigma_i^* = A_i - \frac{1}{2}\rho^2n_\infty.$

Each arm contributes one constraint sphere and the intersection of the three spheres produces a point-pair, $T$

, that represents the two possible configurations of the end-plate:

$\displaystyle T = I_5\bigwedge_{i=1}^{i=3}\Sigma_i^*$

where the $\bigwedge$ notation implies an outer product of all elements following it.

Practically only one of these possible solutions is feasible, the solution which places $Y$ at a greater position along the $e_3$ axis. As we have been careful throughout our equations to ensure our signs are correct we can exploit the oriented nature of CGA to extract the 3D position of this point, $y$ with the use of a single projector:

$\displaystyle P = \frac{1}{2}\left(1 + \frac{T}{\sqrt{T^2}}\right),$
$\displaystyle y = \gadown {Y} , \hspace{0.5cm} Y = -\tilde{P}(T\cdot n_\infty) . \hspace{0.5cm}$

Figure 7.4 illustrates the geometry of the forward kinematic problem.

Full Geometry and Kinematic Constraint Matrix of the Delta Robot

Knowing static kinematic solutions is useful but to do more advanced analysis of the Delta robot mechanism we need to look at velocities. Given we have, up to this point, simplified the pose of the robot in both the forward and inverse cases, we will turn our attention to analysis of the full geometry and possible movements of the limbs. First we will do our analysis with our Screw Theory based framework and then we will compare it with a direct differentiation approach.

Looking at the Delta robot there are two types of joint present. The motor shaft to motor body connection is a revolute joint, the arm to arm connections are spherical joints, and the arm to end-effector platform are also spherical joints. The revolute joints of each motor connection can neatly be represented by a piece of shared geometry, namely a circle with normal parallel to the motor axis. Any radius of circle could be chosen, however as the limbs attached to the motor have a length of $l$ it makes sense for us to use a radius of $l$ for our circle. We will label these circles $P_1, P_4, P_7$ and we can construct them identically to those in section 7.7.2:

	$\displaystyle P_1 = C_1^* = \left( B_1 - \frac{1}{2}l^2n_\infty \right)\wedge(I_3(s_1\wedge e_3)),$
	$\displaystyle P_4 = C_2^* = \left( B_2 - \frac{1}{2}l^2n_\infty \right)\wedge(I_3(s_2\wedge e_3)),$
	$\displaystyle P_7 = C_3^* = \left( B_3 - \frac{1}{2}l^2n_\infty \right)\wedge(I_3(s_3\wedge e_3)).$

The remaining joints are all spherical/ball joints and can therefore be conveniently represented by shared spheres. Again, the radius of these spheres is irrelevant. In this case we will choose a radius of zero, which makes the shared spheres into shared points. We will label each of the limbs of the Delta robot according to the diagram in Figure 7.6. According to this labelling scheme we then have the following mapping from simplified to full robot geometry:

	$\displaystyle P_2 =$ up $\displaystyle \left(\gadown {W_1} + r_ta\right),$
	$\displaystyle P_3 =$ up $\displaystyle \left(\gadown {W_1} - r_ta\right),$
	$\displaystyle P_5 =$ up $\displaystyle \left(\gadown {W_2} + r_tb\right),$
	$\displaystyle P_6 =$ up $\displaystyle \left(\gadown {W_2} - r_tb\right),$
	$\displaystyle P_8 =$ up $\displaystyle \left(\gadown {W_3} + r_tc\right),$
	$\displaystyle P_9 =$ up $\displaystyle \left(\gadown {W_3} - r_tc\right),$
	$\displaystyle P_{10} =$ up $\displaystyle \left(\gadown {X_1} + r_ta\right),$
	$\displaystyle P_{11} =$ up $\displaystyle \left(\gadown {X_1} - r_ta\right),$
	$\displaystyle P_{12} =$ up $\displaystyle \left(\gadown {X_2} - r_tb\right)$

where:

$\displaystyle r_t$	$\displaystyle = \frac{r_e}{\tan(\pi/6)},$
$\displaystyle a = -(s_1\wedge e_3)I_3, \quad b$	$\displaystyle = -(s_2\wedge e_3)I_3, \quad c = -(s_3\wedge e_3)I_3.$

due to the geometry of the end plate as shown in Figure 7.5.

**Fig. 7.5:** The geometry of the end plate allows us to convert from the simplified end plate geometry to the full end plate geometry $P_{10}, P_{11}, P_{12}$ .
$\begin{figure}\begin{center} \par \begin{tikzpicture} \draw[gray, thick] (0,0) -... ...[anchor=east] {$\frac{\pi}{6}$}; \par \end{tikzpicture}\end{center} \end{figure}$

**Fig. 7.6:** A diagrammatic representation of the Delta robot and the relevant objects and screw quantities.
$\includegraphics[width=0.8\textwidth]{DeltaScrews.pdf}$

For the motor shaft connections the circles are shared between the fixed world and the arm pieces. As it is fixed, the velocity screw of the world $\Omega_W$ is 0. Considering equation 7.16 we can substitute in our circles $X_i$ , for $X$ and 0 for $\Omega_W$ :

$\displaystyle 0 = X_i \times (\Omega_W - \Omega_i) = -X_i\times \Omega_{i}.$

(116)

Putting this together with the shared point ball joints the full set of kinematic constraints for the Delta robot are as follows:

	$\displaystyle 0 = P_1 \times\Omega_1,$
	$\displaystyle 0 = P_4 \times\Omega_4,$
	$\displaystyle 0 = P_7 \times\Omega_7,$
	$\displaystyle 0 = P_2 \times(\Omega_2 - \Omega_1),$
	$\displaystyle 0 = P_3 \times(\Omega_3 - \Omega_1),$
	$\displaystyle 0 = P_5 \times(\Omega_5 - \Omega_4),$
	$\displaystyle 0 = P_6 \times(\Omega_6 - \Omega_4),$
	$\displaystyle 0 = P_8 \times(\Omega_8 - \Omega_7),$
	$\displaystyle 0 = P_9 \times(\Omega_9 - \Omega_7),$
	$\displaystyle 0 = P_{10} \times(\Omega_{10} - \Omega_2),$
	$\displaystyle 0 = P_{11} \times(\Omega_{10} - \Omega_3),$
	$\displaystyle 0 = P_{11} \times(\Omega_{10} - \Omega_5),$
	$\displaystyle 0 = P_{12} \times(\Omega_{10} - \Omega_6),$
	$\displaystyle 0 = P_{12} \times(\Omega_{10} - \Omega_8),$
	$\displaystyle 0 = P_{10} \times(\Omega_{10} - \Omega_9).$

The forward kinematic problem has known input velocities $\Omega_1, \Omega_4, \Omega_7$ and output velocity $\Omega_{10}$ . The inverse kinematic problem has known input velocity $\Omega_{10}$ and output velocities $\Omega_1, \Omega_4, \Omega_7$ .

From Constraint Matrix to Jacobian Matrices

In order to convert from a constraint matrix to a Jacobian matrix we will need to define the linear maps mentioned in Section 7.6, $K$ that maps from the input degrees of freedom to the known velocities, and $U$ that maps from the unknown velocities to the output degrees of freedom. We will start by defining three normalised lines $\hat{L}_1, \hat{L}_4, \hat{L}_7$ whose dual forms are directly proportional to the twists $\Omega_1, \Omega_4, \Omega_7$ . These lines are proportional to the twists as they are the axes about which pivot the limbs attached to the motors. We can form these lines from the dual circles $P_1, P_4, P_7$ :

$\displaystyle L_1 = P_1\wedge n_{\infty}, \quad L_4 = P_4\wedge n_{\infty}, \quad L_7 = P_7\wedge n_{\infty},$

$\displaystyle \hat{L}_1^* = \frac{L_1I_5}{\sqrt{-L_1\tilde{L_1}}}, \quad \hat{L... ...{-L_4\tilde{L_4}}}, \quad \hat{L}_7^* = \frac{L_7I_5}{\sqrt{-L_7\tilde{L_7}}}.$

With these dual lines we can map between the motor angular speeds $\dot{\theta}_1, \dot{\theta}_2, \dot{\theta}_3$ and the twists $\Omega_1, \Omega_4, \Omega_7$ .

$\displaystyle \Omega_1 = \dot{\theta}_1\hat{L}_1^*, \quad \Omega_4 = \dot{\theta}_2\hat{L}_4^*, \quad \Omega_7 = \dot{\theta}_3\hat{L}_7^*.$

For the forward kinematic problem we therefore can construct a map $K_F$

as follows:

$\displaystyle \begin{bmatrix} \Omega_1 \\ ... \\ \Omega_4 \\ ... \\ \Omega_... ...bmatrix} \dot{\theta}_1 \\ \dot{\theta}_2 \\ \dot{\theta}_3 \\ \end{bmatrix}$

and as

for the inverse kinematic case we can construct $U_I$

from negative the transpose of $K_F$

At the end-effector in the forward kinematic problem we need to get the translational effect of $\Omega_{10}$ on the central point of the end-plate, we therefore need to calculate $\dot{y}$ . Given we have already calculated $Y$ , getting $\dot{y}$ is a relatively straightforward task. First, we form the line $Y\wedge \dot{Y}\wedge n_{\infty}$ which has the orientation and magnitude of $\dot{y}$ , then we extract the 3D direction and magnitude from this line. To extract the 3D element we can use the dual form of the line and extract the euclidean bivector components which when multiplied with the negative 3DGA pseudoscalar gives the direction of the line. Or, symbolically:

$\displaystyle L^* = (Y\wedge(\Omega_{10}\times Y)\wedge n_{\infty})I_5,$

$\displaystyle \dot{y} = -(L^*\wedge E_0)E_0 I_3.$

This is a general method to extract the linear velocity of a point from the CGA point and its velocity screw. In the case of the Delta robot of course we know that the end-plate is constrained to move only translationally and so $\Omega_{10}$ will be a purely translational bivector of the form $-\dot{y}\wedge n_{\infty}$ . We can therefore extract $\dot{y}$ directly from $\Omega_{10}$ :

$\displaystyle \dot{y} = n_0\cdot \Omega_{10} .$

To turn $\dot{y}$ into individual components we can simply dot it with each basis vector in turn $\dot{y}_{e_1} = \dot{y}\cdot e_1$ etc.

With this in mind we can now construct the known input map for the inverse problem $K_I$ :

$\displaystyle \begin{bmatrix} ... \\ \Omega_{10} \\ ... \end{bmatrix}= \begin... ...in{bmatrix} \dot{y}_{e_1} \\ \dot{y}_{e_2} \\ \dot{y}_{e_3} \\ \end{bmatrix}$

and the unknown output map for the forward problem $U_f$

$\displaystyle \begin{bmatrix} \dot{y}_{e_1} \\ \dot{y}_{e_2} \\ \dot{y}_{e_3}... ...]) & ... \end{bmatrix}\begin{bmatrix} ... \\ \Omega_{10} \\ ... \end{bmatrix}$

Calculating the Jacobian with Direct Differentiation

The Inverse Jacobian

For our direct differentiation method we will start with our simplified inverse kinematic solution and simply differentiate the expressions directly.

First we will write the 3D end-point plate position as a linear combination of basis vectors with coefficients denoted $\alpha_j ,\,\,\, j \in {1,2,3}$ :

$\displaystyle y = \alpha_1e_1 + \alpha_2e_2 + \alpha_3e_3.$

Our goal is to calculate the partial derivative of each motor angle with respect to each of these $\alpha$ coefficients. Taking partial derivatives of $y$

with respect to one of the $\alpha_j$ coefficients trivially gives:

$\displaystyle \frac{\partial y}{\partial \alpha_j} = e_j.$

For now we do not need to worry about which $\alpha$ parameter we are taking derivatives with respect to, so we will leave the derivative of the end-point written as $\frac{\partial y}{\partial \alpha}$ . Using this notation, our ultimate goal in this section is to find an equation for the partial derivative of a given motor angle $\theta_i$ with respect to $\alpha$ , $\frac{\partial \theta_i}{\partial \alpha}$ . To find $\frac{\partial \theta_i}{\partial \alpha}$ we select a specific robot arm $i$

and work back through its joints from the end-point.

The first joint position of interest is $x_i$ , we saw in section 7.7.2 that:

$\displaystyle x_i = y + r_es_i.$

Taking partial derivatives gives:

$\displaystyle \frac{\partial x_i}{\partial \alpha} = \frac{\partial y}{\partial \alpha}.$

The 3D point $x_i$

can then be represented as the CGA point $X_i$

$\displaystyle X_i = \frac{1}{2}x_i^2n_\infty + x_i + n_0.$

The derivative of this CGA point is then easily found:

$\displaystyle \frac{\partial X_i}{\partial \alpha} = \left( \frac{\partial x_i}... ...rtial \alpha}\cdot x_i \right)n_\infty + \frac{\partial x_i}{\partial \alpha}.$

We then form the dual constraint sphere:

$\displaystyle \Sigma_i^* = X_i - \frac{1}{2}\rho^2n_{\infty}$

which, as the radius is fixed, has partial derivative:

$\displaystyle \frac{\partial \Sigma_i^*}{\partial \alpha} = \frac{\partial X_i}{\partial \alpha} .$

As we saw in the previous section, the intersection of this dual constraint sphere $\Sigma_i^*$ and the dual circle $C_i^*$

centred on the motor shaft produces a point-pair $T_i$

that represents the two possible elbow positions for that arm:

$\displaystyle C_i^*$	$\displaystyle = \left( B_i - \frac{1}{2}l^2n_\infty \right)\wedge(I_3(s_i\wedge e_3)),$
$\displaystyle T_i$	$\displaystyle = (\Sigma_i^\wedge C_i^)^*.$

The outer product and dual operations are both linear, which means that taking derivatives is particularly easy here:

$\displaystyle \frac{\partial T_i}{\partial \alpha} = \left(\frac{\partial \Sigma_i^*}{\partial \alpha} \wedge C_i^*\right)^*.$

Of course the elbow can only actually be in one position which we can extract via a projection operation:

$\displaystyle P_i = \frac{1}{2}\left(1 + \frac{T_i}{\sqrt{T_i^2}}\right) , \hsp... ... T_i\frac{\frac{\partial T_i}{\partial \alpha} \cdot T_i}{\sqrt{T_i^2}}\right),$
$\displaystyle W_i = -\tilde{P_i}(T_i\cdot n_\infty),$
$\displaystyle \frac{\partial W_i}{\partial \alpha} = -\frac{\partial \tilde{P_i... ...) -\tilde{P_i}\left(\frac{\partial T_i}{\partial \alpha} \cdot n_\infty\right).$

We can then convert from the CGA to the 3D vector point:

$\displaystyle w_i = \frac{-(W_i\wedge E_0) E_0}{W_i\cdot n_{\infty}},$
$\displaystyle \frac{\partial w_i}{\partial \alpha} = \frac{-(\frac{\partial W_i... ...frac{\partial W_i}{\partial \alpha} \cdot n_{\infty})}{(W_i\cdot n_{\infty})^2}$

and use this to form the derivative of the motor angles with respect to $\alpha$ :

$\displaystyle z_i = w_i - r_bs_i , \hspace{0.5cm} \frac{\partial z_i}{\partial \alpha} = \frac{\partial w_i}{\partial \alpha} ,$
$\displaystyle \theta_i = \operatorname{atan2}({z_i\cdot e_3, z_i\cdot s_i}),$
$\displaystyle \frac{\partial \theta_i}{\partial \alpha} = \frac{z_i\cdot s_i}{\... ...i\cdot e_3) \left(\frac{\partial z_i}{\partial \alpha}\cdot s_i\right)}{z_i^2}.$

This finally gives us an expression for the derivative of the motor angle with respect to the $\alpha$ of the endpoint. Typically in engineering scenarios we would construct a matrix of the partial derivatives with respect to $\alpha_j, \,\, j\in {1,2,3}$ , known as the Jacobian matrix:

$\displaystyle J^* = \begin{bmatrix} \frac{\partial \theta_1}{\partial \alpha_1}... ...\partial \alpha_2} & \frac{\partial \theta_3}{\partial \alpha_3} \end{bmatrix}$

This matrix can then be used to convert an end-point velocity vector to a set of motor velocities:

$\displaystyle \begin{bmatrix} \frac{\partial \theta_1}{\partial t} \\ \frac{\p... ...l \alpha_2}{\partial t} \\ \frac{\partial \alpha_3}{\partial t} \end{bmatrix}$

As it is the Jacobian matrix for the inverse kinematic problem, this matrix is specifically labelled the inverse Jacobian matrix.

The Forward Jacobian

Many problems in robotics require us to take derivatives of the forward kinematic equations. Specifically, we need to know the end-point plate velocity as a function of the motor speeds.

Our forward kinematic solution begins with calculating the position of the elbow point for a given arm $i$ :

$\displaystyle w_i = (r_b + l\cos(\theta_i))s_i + l\sin(\theta_i)e_3 , \hspace{... ...c{\partial w_i}{\partial \theta_i} = -l\sin(\theta_i)s_i + l\cos(\theta_i)e_3.$

With the elbow point we can then calculate the pseudo-elbow point:

$\displaystyle a_i = w_i - r_es_i , \hspace{0.5cm} \frac{\partial a_i}{\partial \theta_i} = \frac{\partial w_i}{\partial \theta_i} .$

We then convert the pseudo-elbow to a CGA point:

$\displaystyle A_i = \frac{1}{2}a_i^2n_{\infty}+ a_i + n_0 , \hspace{0.5cm} \fr... ...l \theta_i}\cdot a_i \right)n_\infty + \frac{\partial a_i}{\partial \theta_i}.$

The forearm length dual constraint sphere can then be constructed about the pseudo-elbow point

$\displaystyle \Sigma_i^* = A_i - \frac{1}{2}\rho^2n_\infty , \hspace{0.5cm} \f... ...tial \Sigma_i^*}{\partial \theta_i} = \frac{\partial A_i}{\partial \theta_i} .$

The intersection of all three constraint spheres, one from each arm, produces the point pair on which the solution lies:

$\displaystyle T = (\Sigma_1\vee \Sigma_2 \vee \Sigma_3) \equiv I_5(\Sigma_1^*\wedge \Sigma_2^*\wedge \Sigma_3^*).$

We can take derivatives of this point-pair with respect to each of the motor angles:

$\displaystyle \frac{\partial T}{\partial \theta_1} = I_5\left(\frac{\partial \S... ...\wedge \frac{\partial \Sigma_2^}{\partial \theta_2} \wedge \Sigma_3^*\right),$
$\displaystyle \frac{\partial T}{\partial \theta_3} = I_5\left(\Sigma_1^* \wedge \Sigma_2^\wedge \frac{\partial \Sigma_3^}{\partial \theta_3} \right).$

We can re-write these derivatives as follows:

$\displaystyle \frac{\partial T}{\partial \theta_i} = (-1)^{i-1}I_5\left(\frac{\partial \Sigma_i^*}{\partial \theta_i} \wedge C^*\right) ,$ where $\displaystyle \,\,\, C^* = \bigwedge_{j \in {1,2,3} \,\,\, j\ne i}\Sigma_j^*.$

(117)

Practically, when we take partial derivatives with respect to one $\theta$ at a time we are effectively freezing two of the motors in position and moving the third. Geometrically, this process forces the end-point plate to move along a circle formed by the intersection of the two constraint spheres centred at the pseudo-elbow points of the frozen motors.

**Fig. 7.7:** With two limbs frozen the end-point plate is constrained to move such that its centre always lies on the circle (shown in red) formed from the intersection of the other two limbs' constraint spheres. The numbers displayed on the figure are the motor angles in radians.

Figure 7.7 shows the geometric significance of Equation 7.27. To get the end-point plate position we again extract one end of the point-pair $T$ :

$\displaystyle P = \frac{1}{2}\left(1 + \frac{T}{\sqrt{T^2}}\right) , \hspace{1c... ..._i} - T\frac{\frac{\partial T}{\partial \theta_i} \cdot T}{\sqrt{T^2}}\right),$

$\displaystyle Y = -\tilde{P}(T\cdot n_\infty) , \hspace{0.5cm} \frac{\partial Y... ...y) -\tilde{P}\left(\frac{\partial T}{\partial \theta_i} \cdot n_\infty\right).$

Finally we convert our end-point back to a 3D point:

$\displaystyle y = \gadown {Y},$
$\displaystyle \frac{\partial y}{\partial \theta_i} = \frac{-(\frac{\partial Y}{... ...\frac{\partial Y}{\partial \theta_i} \cdot n_{\infty})}{(Y\cdot n_{\infty})^2}.$

We can write the end-point plate position as:

$\displaystyle y = \alpha_1e_1 + \alpha_2e_2 + \alpha_3e_3 , \hspace{0.5cm} \fra... ...a_2}{\partial \theta_i} e_2 + \frac{\partial \alpha_3}{\partial \theta_i} e_3.$

With $\frac{\partial y}{\partial \theta_i}$ we are therefore in a position to build the forward Jacobian matrix:

$\displaystyle J = \begin{bmatrix} \frac{\partial y}{\partial \theta_1} \cdot e_... ...ta_2} \cdot e_3 & \frac{\partial y}{\partial \theta_3} \cdot e_3 \end{bmatrix}$

The inverse Jacobian matrix and the forward Jacobian matrix are, as the names suggest, inverse to each other.

$\displaystyle JJ^* = I.$

Comparing Direct Differentiation to Screw Theory

Both the direct differentiation method and our Screw Theory inspired method can be shown to give us numerically identical Jacobian matrices, however the screw theoretic approach is significantly easier for a practitioner to use as it does not require explicitly taking derivatives. In this example of the Delta robot we have the unusual luxury of a simple closed form solution for the pose of the robot in both the forward and inverse case and so the direct differentiation method is easy to compute either manually as we have done here or using automatic differentiation. With other robot architectures we may not have this blessing and would have to rely on non-linear optimisation or algebraic geometry methods to find a pose that satisfies the constraints. Table 7.1 shows a brief qualitative comparison of the two methods, as ever in practice the choice of which is best to use will be down to the problem itself and what tools and computational resources are available.

Table 7.1: A qualitative comparison of the direct differentiation method vs the Screw Theory based technique for analysing the Jacobians of the system.

	Direct differentiation	Screw Theory based
Does not require explicit derivatives of pose calculation	❌	✅
Can be implemented directly with automatic differentiation	✅	❌
Encapsulates information about the entire system	❌	✅

Conclusions and Future Work

In this chapter we have embedded screw theoretic concepts within Geometric Algebra and used this embedding to analyse kinematic pairs and full multi-body systems. The combination of Screw Theory and GA is a particularly potent mix for robotics, allowing clean expressions of geometric constraints and neat representations of kinematic limitations. There are many potential avenues for future work in this area, one promising route would be in the use of higher dimensional Clifford Algebras Breuils et al. [2019b] to represent complex contact surfaces in joints, another might be to expand the allowable motions to include expansion and shear, allowing us to form an extended screw theory for soft robotic modelling. More immediately there is the issue of characterising the computational cost of these Screw Theory based methods. A comparison of the compute speed of these techniques with the direct differentiation method is beyond the scope of this chapter but as it is of practical importance it will likely be a focus of follow up work on this topic.

$\displaystyle ab$	$\displaystyle = (\alpha_1e_1 + \alpha_2e_2)(\beta_1e_1 + \beta_2e_2)$
	$\displaystyle = \alpha_1\beta_1e_1e_1 + \alpha_1\beta_2e_1e_2 + \alpha_2\beta_1e_2e_1 + \alpha_2\beta_2e_2e_2$
	$\displaystyle = \alpha_1\beta_1 - \alpha_2\beta_2 + (\alpha_1\beta_2 - \alpha_2\beta_1)e_1e_2$

Applications of Geometric Algebra in Mathematical Engineering

Abstract:

Direct linear interpolation of geometric objects in conformal geometric algebra

Ray-Object Intersections

Circles/Discs

Extracting Normals and Reflecting Rays

Analytic Form for Normals

Calculating the Derivative of the Object Manifold Projection

Closed Form Derivative of the Square Root Operation

Bézier Curves and Hermite Splines through Geometric Primitives

Representations of wrenches in CGA and PGA

The Screw Inertia Tensor