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:
Doran and Lasenby Dorst et al. . As the notation used in the field varies quite widely between research groups we will also remind the reader of our notation at the start of each chapter.
which represent the number of orthogonal basis vectors in a set which square to respectively. To define an algebra we then can write . This notation is often known as the signature of the algebra. Algebras with are known as degenerate algebras and have to be handled in specific ways in order to work correctly. We often label the orthogonal basis vectors with numbers or symbols and ascribe specific geometric or problem specific significance to the various members of the set.
geometric product. The geometric product is both associative and distributive. The orthonormal basis vectors anticommute with each other and, as we have already specified what they square to, we can compute simplified forms of the products of several vectors. Take for example the algebra , consider two vectors in this space:
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, :
, we define a reciprocal frame of multivectors, , which have the property , for all . When dealing with degenerate metric algebras (those with in the signature) however we have to use a slightly different construction, the pseudo-reciprocal frame Gunn . A pseudo-reciprocal frame has the following property , for all , where is the pseudoscalar of the algebra. The concept of reciprocal and pseudo-reciprocal frames will come up many times throughout this thesis.
we can represent this mapping as a matrix of real coefficients acting on a vector of multivector coefficients for a given basis by the following construction:
is defined as:
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 , the specific closed form solutions of Hitzer Hitzer and Sangwine  or the general closed form solutions of Shirokov Shirokov . 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 .
A rotor, often denoted , in a specific subalgebra can be constructed by the geometric product of an even number of unit 1-vectors 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.
Hestenes et al.  and has been one of the most popular GA frameworks for solving practical problems ever since. CGA has a signature of , three basis vectors and that square to and an additional pair of basis vectors and that square to and respectively. These additional basis vectors are used to form two null vectors and . These null vectors are then used to define a point embedding via the formula:
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. .
operator. The meet is taken with respect to a subspace, depending on the entities involved and the context in which this meet operation takes place, this subspace is either the full space or the smallest possible subspace that contains the operands of the meet operation. In this thesis we will state explicitly in which subspace the meet takes place.
. The specific embedding used in 3D PGA is:
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.
0 it is not possible to dualise simply by multiplication with it. Instead we use the concept of pseudo-reciprocal frames to calculate an alternative dual formulation.
For a given basis vector of PGA, , we will select a pseudo-reciprocal 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 plane-based representation of primitves where points are 3-vectors.
Hrdina et al.  and quadrics Breuils et al. ; Easter and Hitzer  as well as a goal of representing richer groups of transformations Du et al. .
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.
Eggert et al. Valkenburg and Dorst Tingelstad and Egeland De Keninck and Dorst . In Hitzer et al.  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.
, plus two additional basis vectors, and with signatures, , . Two null vectors can then be defined as: and . The mapping of a 3D vector to its conformal representation is given by . 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.
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 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 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:
. 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 4-vector, of . Since we can write the LHS of equation 2.1 as
Since all 4-vectors 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 self-reverse and commutes with and ; we can therefore write
, we have:
From equating 4-vector parts we see that
so that, provided
we simply have
is positive, which it is for lines, planes, circles and point pairs.
can take negative values for some sphere cases.
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
, with , being the conformal representations of two points lying on the line, and the point at infinity. and we normalise such that , therefore . For lines, the 4 vector part of the anticommutator takes the form , thus the square of this is always zero, which means and , which reduces equation 2.3 to and [note that it does not matter which sign we take], giving us the simpler form of the rotor as:
With planes, as with lines, there is no issue of scaling as the objects are infinite. A plane is taken to be the conformal 4-blade 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 plane-to-plane case is particularly simple as the term vanishes:
Let us start with two conformal circles, and not necessarily of the same radius. A conformal circle is a 3-blade 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 .
We start with two conformal spheres, and not necessarily of the same radius. A conformal sphere is a 4-blade 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:
where are conformal points – we can think of a point pair as a line segment. For a point pair, , clearly and gives a positive scalar. We will therefore assume that point pairs are normalised so that .
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.4-2.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  Hadfield and Lasenby  and Hadfield et al. . It is also interesting to note that the nature of the quantity was investigated first in Lasenby, A.N., Lasenby, J. Wareham, R.J. , and then in Dorst et al. , 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 . 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
Dorst and Valkenburg .
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 .
represent specific classes of object and lie on a manifold within the overall subspace of grade . 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?
The result of linear combinations of conformal points is well known Lasenby, A.N., Lasenby, J. Wareham, R.J. . 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 , or the problem is attacked indirectly via carriers Hitzer et al.  or by forming the rotor between the objects, extracting the corresponding bivector, which is then interpolated Wareham and Lasenby  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 + 4-vector factor :and are scalars and the 4-vector 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  and point pairs Dorst et al.  both with common points, the factor is a scalar and the object 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 , 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 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:and are once again scalars.
Since is of the form (scalar + 4-vector) it is self reverse. Defining we get the result that is a scalar, and can therefore write where is a scalar and .
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 non-blade object by its own reverse:
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.
(scalar + 4-vector) required to project the interpolant back to the blade manifold. As before we write our blades as:
We start with point pairs. Previous work Dorst et al.  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. . Figure 3.2 shows two cases of the interpolation of co-planar 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 + 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 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:
Doran  Dorst et al. . Here, as with the point pairs, we can show that this is true for a broader class of circles:
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 + 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. :
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: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.
Ball . 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. ; Kavan et al. ; Müller .
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  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.
that transforms along the 3d vector is: where 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, , whose action on the point at the origin is shown in Figure 3.8.
in terms of the screw parameters:
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:
The interpolant of spheres has been studied before in Cameron  and Doran . 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
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.
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 :, and gives the component of having as a factor and is the rotor that takes to as described in Lasenby et al. . 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. . 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.
and . We wish to find the point that lies closest to both in a least squares sense. First we will construct two orthogonal intermediary lines and where represents the projection of a 3-vector back onto the line manifold. and both lie half way between the two original skew lines but intersect at right angles . The intersection of these lines is the point 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. :
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 ; Peternell and Pottmann ; Petrovic et al. . 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. ; Hestenes ; Lasenby, A.N., Lasenby, J. Wareham, R.J. . 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. ; De Keninck ; Dorst et al. ; Hildenbrand ; Wareham and Lasenby . More recently the design of more intricate surfaces has been investigated with rotors Colapinto  and direct-interpolation of geometric primitives as described in Chapter 3 and published in Hadfield and Lasenby . 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.
and , to the original basis vectors of 3D Euclidean space, giving a complete basis for the 5D space with the following signature: and . These extra basis vectors are used to define two null vectors: and – note that the notation was that originally used when Hestenes first introduced this model in Hestenes et al. . The mapping from a 3D vector, , to its corresponding CGA vector, , is given by:
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. ; Hestenes ; Lasenby, A.N., Lasenby, J. Wareham, R.J.  as well as in the introduction to this thesis.
4.2. It is defined by a rotor (where indicates model view) incorporating rotation and translation that takes the camera from the origin to its pose in space, a focal length and two bounds and on the size of the image plane.
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 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.
) 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 is an -grade blade, is an -grade blade, and the number of basis vectors in the algebra is , then : indicates the -grade component of the multivector , and represents the 5D pseudoscalar of the algebra.
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, , of a line with a plane , satisfies the following: Lasenby, A.N., Lasenby, J. Wareham, R.J. , 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 :
Spheres are also 4-blades and so once again, taking the meet with a ray gives a 2-blade, . 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. :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 pre-normalise 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.
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 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 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.
, and the reflection of that ray at , are two fundamental building blocks in our ray tracer.
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. :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:
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  of the scalar + 4-vector, , can be found as:is with the sign of the 4-vector part reversed and . We then construct by reversing the sign of the 4-vector 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 1-vector . 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.
, it is necessary to design an iterative algorithm to find the roots of the equation. Our implemented algorithm works as follows:
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.
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, , and the surface, (see equation 4.12) occurs when:
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.
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.
for which the ray intersects the surface, we have both the interpolated circle, , and the point of intersection . Using the result from Lasenby, A.N., Lasenby, J. Wareham, R.J. , which is also used in equation 4.7, we extract a tangential line in the plane of the circle at : , postulating that this will be orthogonal to : some future work remains to understand how these two tangent vectors are related to the directions of principal curvature. Clearly 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. , and that the geometric product is minus itself under reversion, i.e. (note that here, and in what follows, we will drop the subscript on ). This follows from the fact that (our circles are all normalised), so that:
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:
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.  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 ( ):
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:
we must differentiate the projection onto the blade manifold of our interpolated object with respect to our evolution parameter . We will continue to work with circles but note that the process works with the general case where is any pure-grade multivector which is a function of a scalar parameter . Let the projection of onto the blade manifold be given by:
. Recall, is given by:
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.
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. , 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 and are point-pairs which represent a line segment, we form a surface via:
and the form of the evolved point-pair . The result is a scalar quantity that can be written as:
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 point-pair . 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 . 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:
4.12 shows two examples of this. As proved in Hadfield and Lasenby  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.,
The disappearing 4-vector part of the projector, which is proportional to , allows the ray-tracer to detect these cases and thus reduces the computational expense of a ray-surface intersection considerably.
and a set of three points , , which together form a triangular facet. First we will form a set of normalised point pairs:
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 , which in the cubic case and with specific first order endpoint conditions are known as Hermite curves.
-blades and couple with the projection to the blade manifold we have the exact same linear interpolation, although this time with going in the other direction. Adopting a notation of as the first object and as the second:
[Order 1] [Order 2] [Order 3]
order multivector Bézier curve is of the form
rational Bézier curve adds weights to the polynomials allowing them to represent a broader class of curves:
, (where we use the notation for points, but will see shortly that these can be replaced by objects) and associated tangent vectors, , at each end of the curve, for :
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:
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 pure-grade space that we like.
One such mechanism for generating tangents for Hermite splines comes from Kochanek and Bartels Kochanek and Bartels . 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:
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.
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 , 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 point-pairs and an evolution parameter, . Integral to ray-tracing 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
Fischler and Bolles  and ICP (Iterative Closest Point) Besl and McKay Segal et al.  algorithms and extends these to multivector correspondence. It is easily parallelisable and designed for good convergence performance with models of real objects.
Eide and Lasenby  Valkenburg and Dorst Tingelstad and Egeland De Keninck and Dorst  and others have applied conformal geometric algebra to 3D registration of point and sphere clouds Kleppe and Egeland Bayro-Corrochano and Rivera-Rovelo . 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 , 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  Tingelstad and Egeland , 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  (the properties of this cost function are further explored in Eide and Lasenby ).
Consider first two arbitrary objects in 3D space represented as and in our conformal model. As in Lasenby et al.  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: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
Eide and Lasenby  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:
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.
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 Segal et al.  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 RANSAC-like Fischler and Bolles  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 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].
Hadfield et al. [sent] python package with both CPU and GPU implementations.
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.
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
Ball . 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  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.
. 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 ; Hadfield et al. . CGA was first suggested by Hestenes Hestenes et al.  and indeed the same author made the first steps in uniting it with Screw Theory Hestenes and Fasse ; Hestenes . 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. .
faster. The Plane-based (or Projective) Geometric Algebra (PGA) aims to strike a balance between efficiency and expressiveness. With a signature of 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 .
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.
, its incident point and the point about which we aim to calculate the moment as then we can calculate the moment with the cross product:
actually require two 3D vectors, one the 3D force vector and one determining a point 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 :
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 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 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:
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,
expressed as 6D Plücker coordinates:
Let us now look at the products of GA. Consider two objects of the form:
Lasenby et al.  the conformal force representation is taken to be 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.  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 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 ‘1D-up’ Approach to Conformal Geometric Algebra' Lasenby , both `On the Homogeneous Model of Euclidean Geometry' Gunn [2011a] and `Rigid Body Dynamics and Conformal Geometric Algebra' Lasenby et al.  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 re-arrange equation (6.3) to get an equation for the relationship between the rotor and its time derivative in terms of 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:
, the resultant wrench is the rate of change of screw momentum with time:
and linear velocity the linear momentum is simply :
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: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 Gallardo-Alvarado  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: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 pseudo-reciprocal 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 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 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: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.
and are both points and, due to the rotational symmetry of a point our reaction wrench can only support translational forces and no moments. In this case 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.
of the form:
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 or dually as given in equation (6.2). Both the dual and direct CGA form work fine for pinning.
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.
. Both CGA and PGA can use planes for pinning.
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.
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.
and time derivatives as:
. In practice numerical integration schemes which integrate will undoubtedly accumulate errors and so wander off the motor manifold. Depending on the application this may or may not be a problem Boyle . 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 . If we label the generator for the current position and orientation in the Lie algebra as which maps to the current rotor with a function then we are interested in finding a function that does the following:
The exponential se(3) kinematic equation is also the subject of Section 5.3 of Liam Candy's PhD thesis Candy , 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 . 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:is a rotation bivector and is a 3DGA vector. We then follow Candy  choosing to define two quantaties:
Hestenes and Fasse ; Tingelstad and Egeland . 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 . 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:
Tingelstad and Egeland . 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:
Du et al. , Cl(8,2) Easter and Hitzer  or even Cl(9,6) Breuils et al.  with this technique in the future should allow for easy configurations of exotic constraints such as pinning dynamic objects to the surface of quadrics.
Wealth, beauty and fame are transient. When those are gone, little is left except the need to be useful.Gene Tierney
Ball , however the mathematical roots that underlie this remarkable field come from Projective Geometry Pottmann and Wallner  and the study of Lie groups and Lie algebras Selig . More recent proponents of Screw Theory include Hunt Hunt , Selig Selig , Davidson Davidson and Hunt , Martins Selig and Martins ; Tischler et al. , Featherstone Featherstone , Gallardo Alvarado Gallardo-Alvarado , Pottmann Pottmann and Wallner , Minguzzi Minguzzi , Lipkin Lipkin  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 ; Fu et al. ; Hildenbrand et al. [2008,2019]; Kim et al. ; Kleppe and Egeland ; Tichý ; Zamora and Bayro-Corrochano . 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  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 ; Tingelstad and Egeland . 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.
and , to the original basis vectors of 3D Euclidean space, giving a complete basis for the 5D space with the following signature: and . These extra basis vectors are used to define two null vectors: and – note that the notation was that originally used when Hestenes first introduced this model in Hestenes et al. . The mapping from a 3D vector, , to its corresponding CGA vector, , is given by: up. We can invert this mapping quite easily so long as we remember to normalise the CGA point such that :
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 :to represent the standard geometric algebra reversion operator we can write:
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:and :
Taking the time derivative of equation (7.2) leads to a formula relating the relative velocity screws of the body in the chain:
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:
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 :
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 . 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 , 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:
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 :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 :
, or in its bivector form which is equivalent to the PGA line formulation where is the direction of the line and is a point on the line.
, or with a dual 1-vector plane which is again the same as the PGA plane where is the normal to the plane and is the distance of the plane from the origin.
or in the form of the bivector dual to that. In CGA a point-pair that is aligned with the axis can also be used.
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 . We can combine these constraints neatly with the anti-commutator product
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.
limbs we write the velocity state as:
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:
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.
Clavel [1990,1991] was invented in 1985 by Raymond Clavel at EPFL after being inspired by a visit to a chocolate packing factory Pessina . 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 ; Merlet . 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.
Since its inception, there have been many variants of the Delta robot Pierrot et al. . 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 end-effector plate. The point half-way between where the two forearm rods connect to the end-effector plate is labelled . We will label the point at the centre of the end-effector 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 in-plane 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 end-effector plate position to the motor angles considering the geometry of the robot as we go. Starting at the end-effector plate the 3D points are translationally offset in plane in the direction giving where is the radius of the end-effector 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 end-effector 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 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 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:
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 :
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 end-effector in the forward kinematic problem we need to get the translational effect of on the central point of the end-plate, 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 end-point 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 :
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.
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.
In this body of work we have made the following contributions:
addition of GA blades representing geometric primitives. We saw in Part I how addition can be used in an intuitive way in to create intermediary objects that blend between extremal control objects and how these intermediary objects are related to transformations that take one control object to another. In the case of CGA circles and point pairs, while the properties of these interpolant objects are still relatively poorly understood. There may be promising future work that could be done in this area, investigating links between the differential geometry of these surfaces and the form of the interpolant objects.
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 multi-body 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 2023-09-05