Hugo Hadfield
Key words and phrases: LaTeX PhD Thesis Engineering University of Cambridge
Let's think the unthinkable, let's do the undoable. Let us prepare to grapple with the ineffable itself, and see if we may not eff it after all.Douglas Adams, Dirk Gently's Holistic Detective Agency
The last few years have seen somewhat of a spike in interest in GA from various disciplines but perhaps most notably from the computer graphics community. Open source GA software has advanced to a point where it is now mature and performant enough to use in real production code environments and there are growing communities of professionals and academics alike designing and implementing GA algorithms on a daily basis.
Despite enjoying the topic I think it is important to emphasise that, despite some grandiose claims in the community, GA is fundamentally one amongst several alternative tools in a mathematical toolbox. As someone new to the field, and to mathematical engineering research in general, there is a danger of myopic referencing and working siloed within the GA literature. I have tried throughout to tie this work into the broader mathematical world but the more research I have read the more my reading list has expanded and, as my supervisor has reminded me, at some point you just have to sit down and write the thesis.
This thesis has been put together from published work of the last few years. It pretty much follows the germination of one idea, namely that adding together geometric primitives in CGA can be a useful, and perhaps intuitive, thing to do. In Part I. we show how this idea came about, its relationship to rotor driven transformations between primitives, some of the properties of the addition of objects, and some of the applications in graphics, vision and visualisation. In Part II. we show how we can use the addition of lines to represent the screws of Screw Theory allowing us to construct a coordinate free basis for, and extend techniques from, that field.
Along the way specific bits of this thesis have been published as stand alone papers:
The result of the geometric product of two or more orthogonal vectors is known as a blade. We refer to the number of orthogonal vectors multiplied together to create a given blade as the grade of the blade. A blade made of the product of basis vectors is often given a subscript notation:
The second most fundamental operation on a multivector in GA after the geometric product is grade selection (also known as grade projection). Grade selection returns only the part of the multivector of a given grade. We often use the notation to denote selecting grade from multivector . With this notation we can define grade selection by considering a multivector defined as the sum of blades of several grades with scalar coefficients, :
A topic that has been the subject of much study in recent years is multivector inverses. The goal of multivector inversion is this: for a given multivector, , to find an multivector, , such that , we would then describe as the inverse of . Not all multivectors are invertible, however for those that are we can invert them with the matrix based techniques of Perwass Perwass [2009], the specific closed form solutions of Hitzer Hitzer and Sangwine [2019] or the general closed form solutions of Shirokov Shirokov [2021]. For invertible multivectors in low dimensional algebras these techniques should all give the same results but for the sake of clarity we will state that unless explictly mentioned we use the inverse as described by Shirokov in Shirokov [2021].
A rotor, often denoted , in a specific subalgebra can be constructed by the geometric product of an even number of unit 1vectors from that subalgebra. Whatever signature algebra you use rotors always have certain properties, namely and is made up of a sum of strictly even grade elements. To perform a transformation with a rotor we sandwich a vector, , with a rotor and the reverse of that rotor.
These geometric primitives can also be parameterised in their dual form with the classic parameters that we would associate with a geometric primitive of that type. We will introduce these primitives for the specific 3D CGA in more detail as we work our way through the thesis but here we will simply mention that they come in two types, primitives formed from the outer product with , sometimes known as `flats', and those without sometimes known as `rounds' Dorst et al. [2007].
PGA is an algebra designed to model only the Euclidean motions and has a point embedding that is familiar to graphics programmers. As we will see later in the thesis this point embedding is identical in fact to the embedding of 3D points as `flat points' in CGA. One of the advantages of embedding points in this way is that normalised points can be interpolated linearly without correction terms.
For a given basis vector of PGA, , we will select a pseudoreciprocal frame, , such that . We can then construct a linear operator in the following form:
This dual does however produce one of the more interesting issues with the PGA; in general the dual of the transformation of a primitive by a rotor is not the same as the transformation of the dual of the primitive, ie. . In general this leads to a selection of one specific space as the primary space for primitves in the PGA, normally the so called planebased representation of primitves where points are 3vectors.
We shape our tools, and thereafter our tools shape us.John M. Culkin
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.
Suppose we wish to find the rotor (rotation, translation, dilation) which takes an object to an object (where and are conformal 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 nontrivial 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 to be , where this comprises both rotation, translation and dilation rotors. We assume both objects are normalised such that
, where for lines, circles and point pairs, and
for planes and spheres:
Note, that
. We motivate our approach by considering the quantity
which is in some sense the `average' object; ie, if we reflect in
, we should get some function of (we assume for convenience that , ie ):
So the reflection does indeed produce a multiple, though the multiple is a scalar plus 4vector, of . Since we can write the LHS of equation 2.1 as
Since all 4vectors square to give a scalar, we can take
, such that
, is a scalar, which we call . We show later that is always positive. We now multiply both sides of equation 2.2 by to give:
We now look to split up such that , where
and and are scalars. If takes this form, it is clear that it is both selfreverse and commutes with and ; we can therefore write
Since
and
, we have:
From equating 4vector parts we see that
so that, provided
;
If
we simply have
if
is positive, which it is for lines, planes, circles and point pairs.
can take negative values for some sphere cases.
If
we then find from the equation which equates scalar parts:
As we need the solution which is guaranteed to be positive:
Recall , , , , so is always positive (as ). We can now write the explicit form of the rotor as:
Taking the positive or negative square root for 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 to via a reflection in the quantity where
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 4vector, , is given by
With planes, as with lines, there is no issue of scaling as the objects are infinite. A plane is taken to be the conformal 4blade of the form , with any 3 conformal points lying on the plane. Conformal planes square to a negative number, so we assume that planes are normalised such that , therefore . Note that .
For planes the anticommutator is a scalar and it is not hard to show that (for normalised planes) is always positive. Thus, the form for the rotor in the planetoplane case is particularly simple as the term vanishes:
(9) 
Let us start with two conformal circles, and not necessarily of the same radius. A conformal circle is a 3blade of the form , where lie on the circle. Circles square to a positive scalar, so we will assume that our circles are normalised such that and therefore . Note that .
The anticommutator, , is in general a scalar plus 4vector, so we must use the form given in equations 2.4,2.5 and little simplification is possible:
with .
We start with two conformal spheres, and not necessarily of the same radius. A conformal sphere is a 4blade of the form , where lie on the sphere. Circles square to a negative scalar, so we will assume that our spheres are normalised such that and therefore . Note that .
As for planes,
is zero, so the rotor takes a very simple form:
(12) 
Since the anticommutator will generally have both scalar and 4vector parts, we again have the general form taken from equations 2.4,2.5:
Note that in the previous rotor derivation we assumed and were blades of the same grade, but nothing further. Therefore, we should, and indeed do, find that the rotor formulae in equations 2.42.7 work for moving between lines and circles and between planes and spheres.
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, , for each recovered rotor, with e, and plot the interpolated objects for each of , , with ee, where and . Figure 2.1 shows these interpolations for each class of object.
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 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 , where is the rotor taking to . 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.
I write rhymes with addition and algebra, mental geometry.Ice T, Mind Over Matter
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 , plus two additional basis vectors, and with signatures, , . Two null vectors can therefore be defined as: and . The mapping of a 3d vector to its conformal representation is given by . ^{}
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?
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 and represented as and 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 :
We can therefore get a useful interpolation of points by taking a direct linear interpolation and simply adding the final term to the result. If , we can recover via the following formula (assuming ):
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 that reflects one object , into another , as the left multiplication of the summation of the blades by a scalar + 4vector factor :
where and are scalars and the 4vector part of is proportional to the anticommutator of 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 is a scalar and the object is simply `halfway' between the objects. We can extend this notion to the cases where the addition of objects is not a blade by using our object , which has been corrected to being a blade, as the halfway object. We can use this idea of the halfway object to recursively subdivide the space between and allowing us to create objects that are any fraction of and . 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 .
Consider the general interpolant, where and are blades of the same grade. We claim that we can project into object space in a simple and general way. First we will generalise equation (3.2) to the interpolation case:
where and are once again scalars.Since is of the form (scalar + 4vector) it is self reverse. Defining we get the result that is a scalar, and can therefore write where is a scalar and .
To use this decomposition we need to extract from . 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 , where is a valid object (squaring to 1). Now define , ie. it only contains 0 and 4 grade coefficients and is an element of the conformal algebra. Then, defining , the square root can be found as:
To use this method to find we multiply our nonblade object by its own reverse:
This is now in a form where we can apply the above square root formula: It now simply remains to isolate via multiplication by where . Since is a scalar, we have
This result is particularly important as we have identified a way of projecting any pure grade object of the form (with 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 ,
Unfortunately, we cannot apply this form of interpolation to points as we encounter a problem due to the fact that for a conformal point , . However, we saw in equation (3.1) that points can be interpolated very easily using known explicit formulae.
We start with point pairs. Previous work Dorst et al. [2007] has shown that when an end point is shared between point pairs and 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 define a circle and a fourth point lying on the circle will satisfy , this allows us to define a check to see if two point pairs are chords of the same circle. Point pairs , will satisfy if they are both chords and any additional chord, , of the same circle will satisfy and thus . This leads to:
Note, since , the projector 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 coplanar point pairs that lie on the same 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 + 4vector 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 4vector grades. The 2vector part of the geometric product comes from the inner product between the point pairs and the 4vector. ie. for point pairs and , . For the general case of two point pairs not lying in plane ie. , we can show that there is only one object that behaves in this way, the sphere , as it passes through both end points of both point pairs. This is illustrated in Figure 3.3 and suggests that the sphere is intrinsically tied to the form of the interpolant objects. Indeed we can see from the same visualisation that the interpolant of point pairs and always has endpoints lying on the surface of the sphere .
We can prove this by showing that or also produces the sphere:
We can additionally find the unique common sphere by finding the join of the circles or by reverting to linear algebra techniques:
First we define:
It is also the case that the interpolant lies on the surface of the common sphere:
Thus far we have dealt exclusively with circles on a common sphere. In the case in which and do not lie on the same sphere we can again look at how the interpolants behave by considering the form of the (scalar + 4vector) 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 3vector grades, however the 1vector part of the geometric product comes only from the inner product between the 4vector and the circles. To maintain grade after the multiplication the 4vector 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]:
The intersections of the interpolant circles with the sphere 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 and and indeed we can numerically verify that this is the case:
When looking at lines we can attempt to use some of the same techniques that we used for circles. First consider the form of . For lines , giving the form of the projection of as:
where and are scalars. While neat, this form 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 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 , a point on the screw axis and a screw pitch . The pitch represents how far to move in the direction of the screw axis for each complete revolution about the axis.
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.
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, , whose action on the point at the origin is shown in Figure 3.8.
All 4vectors 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:
(26) 
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 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:
As we have seen previously, the interpolation of two conformal points and is of the form
For the case that the surfaces of and are just touching we have the condition
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:
This leads to a simplified model that retains the core features of the original model. For comparison of objects and we use the cost function for a rotor as defined in Eide and Lasenby [2018]:
where , and gives the component of having as a factor and is the rotor that takes to 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.
One of the most fundamental and simple clustering algorithms is known as means clustering Duda et al. [2001]. Consider a 3d scene composed of geometric objects of a given grade. We have multiple noisy observations for each object and so would like to fit centroids to these clusters to represent the “true” objects in the world.
The steps for implementing this clustering are given below:
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.
Arithmetic! Algebra! Geometry! Grandiose trinity! Luminous triangle! Whoever has not known you is without sense!Comte de Lautreamont
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 linesegments, 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 raytracers/pathtracers/spheremarchers 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 directinterpolation 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 raytracer we describe in this chapter.
All vectors formed from such a mapping are null. CGA is chosen for the construction of the raytracer 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.
width= 
width= 
We take to be at the bottom left hand corner of the image. For an image of width and height , the world coordinates of the point at the centre of pixel are given by:
Initially we will start with some basic objects representable as blades in CGA. The raytracer will thus initially concentrate on rendering planes, spheres and circles/discs, an example of which is shown in Figures 4.1 and 4.3.
A plane is a 4blade and a ray is a 3blade so the meet gives a 2blade. 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, , of a line with a plane , satisfies the following: Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004], where is a scalar. When extracting the 3D intersection point , we need to account for the sign and magnitude of the line in our extraction, we can do this via the constant of proportionality :
(31) 
Spheres are also 4blades and so once again, taking the meet with a ray gives a 2blade, . With spheres, there can be zero, one or two points of intersection corresponding to the cases where , and respectively.
If (with and null vectors) and , the points can be extracted from the point pair/blade, , by the following formula Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]:
where are scalar constants. If we define with oriented in the same direction as our ray , then is the point closest to the origin of the ray, , as long as our sphere is `in front of' the ray source. To ensure the alignment we can prenormalise our sphere via the following expression: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, , satisfy . If the meet of a ray and a normalised sphere , is then formed from , as in equation 4.3, the resulting bivector will be ordered as where 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 with the smaller sphere, would result in the point pair . For the larger sphere in figure 4.4, the point pair resulting from the meet will be . We extract the points from the point pair and form the distance between these points and (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 . This allows us to perform bounces only with spheres that are in front of the ray origin point .
Figure 4.5 shows an example of each case along with a geometric interpretation of the form of the meet. If and there is an intersection, the plane containing the circle is formed by taking the wedge product between the circle and , , and the intersection point is then extracted from the ray and this plane.
width=

If , 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 2part 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 and , so the bivector is . If the bivector has positive square, there is no intersection. If the bivector squares to zero there is one intersection at , and the bivector is , where . In all cases, the intersection points are easily extracted.
For a plane which intersects with a ray, we can compute the reflection of an incident ray (we assume and have been normalised such that ) with the plane by simple sandwiching: . The resulting line is oriented correctly, passes through the intersection point and . For the case of a sphere , we use the following formula from Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]:
where is the first point of intersection. Here, is an example of an inversion, where the incoming ray/line, , 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, , is the reflected ray, . Figure 4.6 illustrates this geometrical construction. Note that one can also form the tangent plane at and reflect in this plane; this is performed by the following formula:(35) 
For a circle/disc , we first form the plane 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: . 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, , and the incident ray, , both normalised such that they square to 1. The normal line to the surface of the sphere, , can be simply found:
Consider a quantity . We then define the quantity , and with this the principal square root Dorst and Valkenburg [2011] of the scalar + 4vector, , can be found as:
With this square root we can then form: where is with the sign of the 4vector part reversed and . We then construct by reversing the sign of the 4vector part, ( ), and use this to produce the following expression for the projector and interpolated circle :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.
We saw earlier that the intersection of a ray with a circle produces the 1vector . If the ray lies in the plane of the circle and if and 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 , we need to find a value of for which:
Figure 4.9 provides a simple visual illustration of one example of the shape of this curve as a function of . 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.
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 precompute 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.
(42) 
(43) 
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 . Some of the 6 remaining roots may be imaginary, some may be outside of the range and some will be spurious roots corresponding to . To filter out the valid roots we simply take all roots between 0 and 1 and evaluate at these positions, selecting the roots for which for some small threshold where (in our experiments works satisfactorily).
An interesting point to note here is that we could extend this intersection finding method to being higher order functions of , so long as . Generating such higher order splines through geometric primitives is described in section 4.8.
(45) 
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 . 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, , have a constant `length' due to normalisation, so as with the circles, we can differentiate wrt to see that and are orthogonal, ie . If we were to define the `velocity', to be the inner product with the angular velocity bivector given in equation 4.19:
(47) 
If we therefore assume that is the direction we want, we can calculate the tangent line in this direction via:
The fact that lines and are perpendicular can be verified by showing that the quantity 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 and intersect at a point, then . In addition, if they are orthogonal, ). If this is the case, will reverse to minus itself.
To show this, we need to consider the reverse of . We will need the facts that that , , , and . We have shown all of these identities earlier in this section. We also need an additional fact, which is that anticommutes with . To see this we use another standard result ( ):
(50) 
We are now in a position to expand out :
In the above we have used the facts that , and . Note that the term in the final line of equation 4.24 can be written as (from the standard reflection formula and the fact that ). 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 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 and anticommute, which is a good indication that lies in the right direction.
Given these two orthogonal tangent lines and , we can construct the plane tangent to the surface at by computing the join of the two lines. Or, we can bypass the plane entirely and compute the surface normal line directly as:
(55) 
(56) 
(57) 
An important point to note here is that this blade projection derivative is gradeagnostic and so can be used for objects other than just evolved circles.
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 pointpairs in CGA Dorst et al. [2007], as well as the practical desire to represent ribbonlike surfaces, we can apply similar raytracing methods to surfaces formed from the interpolation of pointpair bivectors representing line segments. If and are pointpairs which represent a line segment, we form a surface via:
(58) 
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.
We saw earlier that the meet will be zero if the ray hits anywhere along the carrier line of the pointpair . 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 (for some ) and the ray intersect at a point , we can then check if the intersection point is within the bounding sphere of the surface by ensuring:
The disappearing 4vector part of the projector, which is proportional to , allows the raytracer to detect these cases and thus reduces the computational expense of a raysurface intersection considerably.
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.
[Order 1]
[Order 2]
[Order 3]

A very common use of Hermite curves is in the construction of Hermite splines; these are piecewise 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 piecewise 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:
(60) 
Thus for the derivative to evaluate to the same on either side of the boundary, we only need to check that is the same either side of the boundary. Considering the equations in section 4.6.1 we can see that is a function only of and 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 puregrade 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 known as tension, bias and continuity respectively. For given control objects the corresponding tangents can be calculated using the control objects in the spline and which lie previous to, and after, the curve in the order of the spline:
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 and then transform these points along a small step in 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 , to another, , can be calculated with the following process:
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 evenly spaced points on the unit circle in the , plane by successive rotations about the origin of a point lying initially at yielding for i.e., for a fixed rotor :
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 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.
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 pointpairs and an evolution parameter, . Integral to raytracing these evolved surfaces is the derivation of analytic intersection points and normal vectors.
We are stuck with technology when what we really want is just stuff that works.Douglas Adams, The Salmon of Doubt
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 , plus two additional basis vectors, and with signatures, , . Two null vectors can therefore be defined as: and . The mapping of a 3D vector to its conformal representation is given by .
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 and in our conformal model. As in Lasenby et al. [2018] we will extract the rotor that takes one object to another . Note that the objects will have an orientation (sign), and the rotor extraction will be orientation dependent. Once we have our rotor between our conformal objects, the next step is to use this rotor to define a cost as a function of this rotor:
where indicates the grade part of . Equipped with this idea of closeness of objects, for a given , a query object is assigned to each of the reference objects (i.e. this is done for all ), assuming the model and query sets are spatially close. For each object pair we form the rotor, that takes the query object to the reference object. The minimum cost assignment is then taken as the correct match, , for that query object
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:
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 nonlinear 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.
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:
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.
To counteract the local minima issue, we modify our procedure to incorporate sampling in a RANSAClike 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 lots of matches. Each of these 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 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:
The disadvantage of moving to a samplingbased 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].
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 multibody kinematics.
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
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.
We will now do a quick survey of the different force representations, and in doing so, think about what an effective force representation requires.
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 nonlocalised vector is known as a free vector. While respecting their nonlocalised 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 component of the force, or maybe we should put them in the lower 3 to make them align with the part of the force representation.
Initially, which choice to make seems nonobvious. To help us gain some insight we consider the problem of two antiparallel coplanar forces of equal magnitude acting on a rigid body. First, label the forces themselves:
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,
Let us now look at the products of GA. Consider two objects of the form:
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:
`Rigid Body Dynamics and Conformal Geometric Algebra' Lasenby et al. [2011] therefore uses a mixture of 1vectors 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 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:
As discussed in Anthony Lasenby's `Rigid Body Dynamics in a Constant Curvature Space and the ‘1Dup’ 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.
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 will transform an arbitrary fixed point in the body frame into the corresponding point in the world frame:
We can further rearrange equation (6.3) to get an equation for the relationship between the rotor and its time derivative in terms of this quantity :
This quantity 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 and change equation (6.5) to:(67) 
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 and the GA version as :
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 . The velocity of objects attached to the frame of the body are given by:
By differentiating this equation we can calculate the acceleration of objects in the frame:
In this case the body is initially at rest implying allowing us to eliminate the first term and leaving us with:
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 and hence .To calculate the form of this linear function 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 GallardoAlvarado [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 that are matched to the motor bivectors such that when the GA inner product is taken between the two matched elements the result is 1 but otherwise 0, ie.:
As is a linear function operating on a motor bivector in CGA it can be written as:
where are scalar coefficients, is the reciprocal frame for the motor bivectors and are the motor bivectors themselves. We will take and as follows:
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:
So far we have dealt with all the translational elements of the transformation, however we have only prescribed 18 of the total 36 (66) 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:
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 with angular speed , ie. then, provided that the centre of mass has no linear velocity we can calculate the linear velocity of a point on the body as:
Of course our choice of as the principal axis to which our torque vector aligns was arbitrary, we could equally have chosen one of the other two principal axes, and with their respective and . As a result of this symmetry we can directly identify our final parameters: , , , , , .
Figure 6.2 gives a graphical illustration of why this form of mapping is required for rotations and torques.
For our 6D vector representation had we not just done the maths of the previous section we might have been tempted to stack the and on top of each other to produce a combined inertia tensor like:
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:
As with our CGA reciprocal frame let's now write our PGA pseudoreciprocal frame in two groups:
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 ( ) 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 nonaxis 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 and a world frame momentum . 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:
The form of is:
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 as:
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 as the sum of external wrenches, , plus some reaction wrench, , caused by the constraints. As we already know , all we need to calculate is the value of 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, , produced by a virtual world frame wrench, , when the body moves with a body frame screw velocity is given by:
If we specify the way that varies with time we can add curved surface constraints. Consider a situation in which a rigid body is constrained such that one point always touches a sphere centred at point . Given the point is always touching the sphere we know that must always be parallel to the line joining and , we would therefore write:
Consider a geometric primitive represented by multivector in the body frame and the same geometric primitive represented by multivector when expressed in the world frame. These two multivectors can be related by the rotor :
Lets consider first the case that both of these `views' of the object are fixed, ie. the position and orientation of cannot change with respect to the coordinate system of the body and the position and orientation of cannot change with respect to the origin. Mathematically we are stating that . If we substitute these values into (6.19) for the time derivatives we end up with the following equation:
This equation is a constraint on the second time derivative of that will ensure that and do not vary with time. We can go a step further here and substitute our expression for from equation (6.7), leading to:
For a given this constraint is linear in and can be solved for 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.
First note we can still rearrange to separate terms in :
What this means practically is that we can set and 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.
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 as a function of and or . First, we note that it is possible to write a motor bivector in the form:
where is a rotation bivector and is a 3DGA vector. We then follow Candy [2012] choosing to define two quantaties:
Wealth, beauty and fame are transient. When those are gone, little is left except the need to be useful.Gene Tierney
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 BayroCorrochano [2004]. With modern computing capabilities, the high level description of geometry that these algebras afford allows researchers elegant, concise and coordinatefree 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.
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:
Whichever GA you choose, the motor bivectors, when exponentiated, produce rotors that can implement rigid body motions. In general it is possible to split 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 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.
Consider a rigid body . The rotor 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 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 . Consider another body, , again with a rotor from fixed frame to world frame, from principal axis aligned fixed frame to arbitrary fixed frame and 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 to limb is labelled :
Using to represent the standard geometric algebra reversion operator we can write:(93) 
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 . 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:
We can decompose this rotor into and :(95) 
(96) 
(97) 
Taking the time derivative of equation (7.2) leads to a formula relating the relative velocity screws of the body in the chain:
(99) 
(100) 
(101) 
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 at all, just of . We could choose any decomposition of into and 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: , and therefore . Under these conditions equation (7.12) simplifies to the following:
(103) 
A direct result of the definition of is that the time derivative of a geometric object is given by the commutator product of the object with :
(104) 
In a shared geometry joint there exists a piece of geometry, which we will label , that is fixed relative to both frames in the kinematic pair. Practically what this means is that when acted on by the velocity of or this object must have the same , or more formally:
(105) 
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 in frame and an object in frame :
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 and :(108)  
(109) 
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, and , the orthogonality constraint can be written as and the same centre point constraint for circles of the same radius can be written Hadfield and Lasenby [2019]. We can combine these constraints neatly with the anticommutator product
(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.
(111) 
(112) 
Consider a robot with velocity state , now specify some of the limbs as inputs, with known and some of the limbs as outputs with unknown . We can write this as follows:
(113) 
(114) 
Designating a linear map from the input degrees of freedom, , to the known velocities of the system, and a linear map from the unknown velocities to the output degrees of freedom, , allows us to calculate a matrix that is equivalent to the Jacobian matrix of the system:
The matrix 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.
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 endeffectors driven by multiple underactuated parallel kinematic chains GallardoAlvarado [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 endpoint 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.
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 ; we will number each upper arm . The connection point of the motor and upper arm will be labelled . 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 . 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 endeffector plate. The point halfway between where the two forearm rods connect to the endeffector plate is labelled . We will label the point at the centre of the endeffector plate . 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 inplane orientation fixed as well. The Delta robot is therefore a purely translational mechanism.
To solve this problem we need to work backwards from the 3D endeffector plate position to the motor angles considering the geometry of the robot as we go. Starting at the endeffector plate the 3D points are translationally offset in plane in the direction giving where is the radius of the endeffector plate. Due to the geometry of the robot the elbow point is constrained to lie on a sphere with radius equal to the length of the forearms centred at this point . We can write the dual form of this sphere as:
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 endeffector plate given the motor angles .
To solve the forward kinematic problem we will consider the robot one arm at a time. For a given arm motor angle the 3D position of the elbow point can be calculated as:
Practically only one of these possible solutions is feasible, the solution which places at a greater position along the 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, with the use of a single projector:
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 endeffector 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 it makes sense for us to use a radius of for our circle. We will label these circles and we can construct them identically to those in section 7.7.2:
up  
up  
up  
up  
up  
up  
up  
up  
up 
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 is 0. Considering equation 7.16 we can substitute in our circles , for and 0 for :
Putting this together with the shared point ball joints the full set of kinematic constraints for the Delta robot are as follows:
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, that maps from the input degrees of freedom to the known velocities, and that maps from the unknown velocities to the output degrees of freedom. We will start by defining three normalised lines whose dual forms are directly proportional to the twists . 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 :
At the endeffector in the forward kinematic problem we need to get the translational effect of on the central point of the endplate, we therefore need to calculate . Given we have already calculated , getting is a relatively straightforward task. First, we form the line which has the orientation and magnitude of , 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:
With this in mind we can now construct the known input map for the inverse problem :
First we will write the 3D endpoint plate position as a linear combination of basis vectors with coefficients denoted :
The first joint position of interest is , we saw in section 7.7.2 that:
Our forward kinematic solution begins with calculating the position of the elbow point for a given arm :
Figure 7.7 shows the geometric significance of Equation 7.27. To get the endpoint plate position we again extract one end of the pointpair :
In this body of work we have made the following contributions:
In the case of the addition of lines we have seen in Part II how linear combinations of these flat elements produces the screws of Screw Theory and the various grades of the geometric product map neatly onto its operators. GA provides a particularly elegant framework for Screw Theory and the unification of geometric primitives, screws motions and covariant products in one algebraic framework and holds great promise for concise descriptions of joint types for multibody systems. The combination of Screw Theory and Geometric Algebra is underdeveloped, given the wide interest of GA practitioners in robotics.
This document was generated using the LaTeX2HTML translator Version 2021.2 (Released July 1, 2021)
The command line arguments were:
latex2html thesis.tex split 0 html_version=4.0 no_navigation
The translation was initiated on 20230905