# Applications of Geometric Algebra in Mathematical Engineering

Key words and phrases: LaTeX PhD Thesis Engineering University of Cambridge

#### Abstract:

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

# Introduction and Background

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

## Introduction

There are many reasons that are given in the literature for studying Geometric Algebra (GA). It is claimed to produce many elegant formulae, it is described as imbuing linear algebra with intuitive geometric meaning and it is claimed to be both computationally efficient and easy to program up. During the course of my PhD I have found many of these claims to be true but from my perspective the real benefit of the GA framework is in bootstrapping' the learning of multiple fields. GA is a very broadly applicable framework for expressing problems from various topics and as a result it allows a practitioner to leverage their knowledge from previous domains when coming to a new problem or field. This bootstrapping advantage is one of the main reasons I have enjoyed this PhD, with GA in my toolbox it feels like large swathes of previously inaccessible research are now within reach.

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:

• Chapter 2 has been published in Advances in Applied Clifford Algebras' as Calculating the Rotor Between Conformal Objects' Lasenby et al. [2018]
• Chapter 3 has been published in Advances in Applied Clifford Algebras' as Direct Linear Interpolation of Geometric Objects in Conformal Geometric Algebra' Hadfield and Lasenby [2019]
• Chapter 4 has been published in Advances in Applied Clifford Algebras' as Exploring Novel Surface Representations via an Experimental Ray-Tracer in CGA' Hadfield et al. [2021]
• Chapter 5 has been published in Advances in Applied Clifford Algebras' as REFORM: Rotor Estimation From Object Resampling and Matching' Hadfield et al. [2019]
• Chapter 6 has been submitted to Advances in Applied Clifford Algebras' as Screw Theory in Geometric Algebra for Con- strained Rigid Body Dynamics'
• Chapter 7 has been submitted to Advances in Applied Clifford Algebras' as The Kinematics of Multi-body Systems in Geometric Algebra'.

## Basics of GA

Geometric Algebras, also known as Clifford Algebras, are a class of algebras that allow operations between objects of mixed grade. Here we will give a quick introduction to the topic and state some of the required background and results that will be put to use throughout the thesis. As the literature is fairly extensive we would suggest a reader not familiar with the subject, or looking for proofs of some of the more fundamental results, to consult one of the excellent books on the topic Doran and Lasenby [2003]Dorst et al. [2007]. 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.

### Defining a Specific Subalgebra

Typically when working in GA we choose to work in a specific closed subalgebra. Often we will define this algebra with 3 integers 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.

### The Products of GA

The fundamental product of GA is the 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 geometric product of these two vectors is:

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:

A linear combination of blades of different grades is known as a multivector. A linear combination of blades of the same grade is known as a pure grade multivector. For an algebra with orthogonal basis vectors the highest grade blade that can be constructed is also of grade and is known as the pseudoscalar of the algebra and often denoted .

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, :

All other products are defined in terms of the geometric product and grade selection. In this thesis we will refer to the inner product which we will represent by , the outer product that we will represent with , and the commutator product that we will represent with . The inner product that we will use throughout this thesis is defined as follows:

Note that the choice of the inner product with scalars being 0 or not varies among researchers and software systems. We have chosen this inner product to match the one used in the clifford' python package Hadfield et al. [sent]. The outer product has the following standard definition:

In the same way we use the shorthand to imply the geometric product of several elements so we use the notation to imply the wedge product of several multivectors. The commutator product is typically denoted with a cross, , and is defined in terms of the geometric product:

### Reciprocal and Pseudo-reciprocal Frames

A particularly useful tool for constructing proofs in GA is the concept of reciprocal frames. Given an ordered set of multivectors, , 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 [2020]. 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.

### Linear Multivector Mappings

One of the most commonly used types of mapping of a multivector is a linear mapping (also known as a linear function). Linear mappings are exceptionally useful in applied mathematics and are typically represented as matrices with a fixed basis. GA by default does not define a fixed basis for its multivectors, it does however define products that act in a coordinate-free way. If we have a linear mapping of a multivector 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:

For a degenerate metric algebra we can do the same with a pseudo-reciprocal frame, followed by the selection of the coefficient of the pseudoscalar:

### Multivector Reverse and Multivector Inverses

The reverse is one of the fundamental operations of GA. The reverse of a multivector is defined as:

Like many operations in GA the reverse is linear and thus distributes over addition:

A topic that has been the subject of much study in recent years is multivector inverses. The goal of multivector inversion is this: for a given multivector, , to find an multivector, , such that , we would then describe as the inverse of . Not all multivectors are invertible, however for those that are we can invert them with the matrix based techniques of Perwass Perwass [2009], the specific closed form solutions of Hitzer Hitzer and Sangwine [2019] or the general closed form solutions of Shirokov Shirokov [2021]. For invertible multivectors in low dimensional algebras these techniques should all give the same results but for the sake of clarity we will state that unless explictly mentioned we use the inverse as described by Shirokov in Shirokov [2021].

### Rotors

An aspect of the GA framework that has found widespread use is the inclusion of geometrical transformations as elements of the algebra. The primary mechanism by which we represent geometrical transformations are with a class of multivectors known as rotors.

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.

Rotor application is grade preserving and can be extended over the outer product producing the result Doran and Lasenby [2003]:

The signature of the algebra determines what types of geometric transformations are available in the rotor group of that algebra.

## 3D CGA

### Homogeneous Point Embedding

The Conformal Geometric Algebra (CGA) was first described by Hestenes in Hestenes et al. [1985] 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:

which is sometimes referred to as    up. An important consequence of this specific point embedding is for a given point :

and more generally for a point and point we get the result:

The embedding is a homogeneous representation of the point, ie. and for scalar both represent the same 3D point. To invert the mapping we therefore first have to normalise the point such that:

we can then extract the components of via the dot product:

or via:

where . The mapping from a conformal point back to the 3D point it represents is sometimes referred to as .

### Geometric Primitives

One of the primary motivations for shifting to a higher dimensional algebra to represent 3D space is that it allows us to represent a richer set of possible entities. Some of the most immediately useful objects that CGA allows us to represent are the geometric primitives, these primitives are represented by the blades of the algebra and can be constructed with the outer product of several points.

These geometric primitives can also be parameterised in their dual form with the classic parameters that we would associate with a geometric primitive of that type. We will introduce these primitives for the specific 3D CGA in more detail as we work our way through the thesis but here we will simply mention that they come in two types, primitives formed from the outer product with , sometimes known as flats', and those without sometimes known as rounds' Dorst et al. [2007].

### Intersections

In many contexts where we are working with geometric primitives we are interested in computing their intersections. In CGA we can use the meet to compute intersections; the meet is represented by the 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.

### Transformations

The rotor driven transformations of CGA include all conformal transformations. The full group of conformal transformations has 8 degrees of freedom, but many applications stick to using only the 6 degree of freedom rigid body motion subgroup. Additionally the primitives of CGA support reflection/inversion via the double sided formula:

where is a given blade. The reflection operations are practically useful in many situations and particularly find application in computer graphics.

### Duality

The duality operation in CGA is performed by multiplication with the pseudoscalar. In non-degenerate metric algebras such as CGA primitives behave reasonably under rotor transformations whether represented dually or directly, and as such the idea of privileging one representation over another makes very little sense. In practice we tend to use one representation over another when it fits the way in which we choose to parameterise a problem. In this thesis we will switch between dual and direct primitive representations fairly fluidly throughout.

## 3D PGA

### Homogeneous Point Embedding

The 3D Plane-based (or Projective) Geometric Algebra (PGA), like the 3D CGA uses a homogeneous embedding, representing points as trivectors. PGA has a signature of . The specific embedding used in 3D PGA is:

where and . The inverse mapping is:

where represents the (linear) duality operation of PGA that we will describe shortly.

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.

### Geometric Primitives

PGA is a subalgebra of CGA, it is in fact the subalgebra of the flat elements of CGA. This means that it contains points, lines and planes as primitives which have grades 3, 2 and 1 respectively.

### Intersections

Intersections in PGA can again be computed with the same collection of operators that we use in CGA. PGA however is a more opinionated' algebra, it chooses a default space for its geometric primitives and the intersection between these primitives are always computed with the outer product.

### Transformations

The bivectors of PGA are exactly the Euclidean motion generators and as such the rotor driven transformations in PGA are the rigid body motions. This is generally a useful thing, one of the advantages of PGA is its simplicity; if all you need is flat objects and the Euclidean group there is not much point creating the whole conformal group and then just playing with part of it.

### Duality

Duality in degenerate metric algebras is a tricky business. As the pseudoscalar now squares to 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:

where we use the notation to mean taking the coefficient of the pseudoscalar. This linear operator acts as a dual for PGA and, especially when dealing with lines, is a convenient and useful construction.

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.

## Other GA Frameworks

Several other GA frameworks have been studied over the last few years, and some of them will be mentioned in this thesis. Many of these constructions have been designed to extend CGA in some way to include more geometric primitives and transformations, most notably there has been a push to include higher dimensional primitives such as conics Hrdina et al. [2018] and quadrics Breuils et al. [2018]; Easter and Hitzer [2017] as well as a goal of representing richer groups of transformations Du et al. [2017].

# Computer Graphics, Computer Vision and Visualisation

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

# Calculating the rotor between conformal objects

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

## Introduction

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

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

## Related Work

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

## Conformal Geometric Algebra

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 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.

## A Rotor between Objects

Suppose we wish to find the rotor (rotation, translation, dilation) which takes an object to an object (where and are conformal -blades representing the lines/circles/planes/ spheres/point pairs). If we firstly take lines as an example, conventionally we would translate along the common perpendicular and then rotate about the intersection point – which requires a series of 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:

Note, that . We motivate our approach by considering the quantity which is in some sense the average' object; ie, if we reflect in , we should get some function of (we assume for convenience that , ie ):

 (1)

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

we propose to use the spinor quantity to form . As above (but now with included) gives ;

 (2)

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

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

so that is our required rotor and . To find such an we can use the square root formula given in Dorst and Valkenburg [2011] or simply equate scalar and 4-vector parts of the equation . We do the latter first in order to see how the particular form of our scalar plus 4-vector behaves and then confirm that it agrees with the formula in Dorst and Valkenburg [2011]:

Since and , we have:

From equating 4-vector parts we see that so that, provided ;

If we simply have if is positive, which it is for lines, planes, circles and point pairs. can take negative values for some sphere cases. If we then find from the equation which equates scalar parts:

where , since the 4-vectors always square to give zero or a negative scalar. This is a quadratic in :

 (3)

with solutions given by:

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:

1. If :

 (4)

 (5)

2. If and

 (6)

3. If and ,

 (7)

where and .

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

where and are as given previously, ie and takes the form in equations 2.4,2.6,2.7 depending on the nature of . We will see in Chapter 3 that the quantity projects the -vector obtained from the addition of the two blades and onto an -blade and therefore an object – the object being that in which we reflect in to get .

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

where . Here, our , so that (taking the solution corresponding to the sign):

and

giving , as required (taking ).

### Lines

Conformal lines take the form , 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:

 (8)

### Planes

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:

 (9)

where .

### Circles

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

Let us start with two conformal circles, 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 .

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

 (10)

 (11)

with .

### Spheres

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:

 (12)

where , if and if . is the same sphere as , so in a sense it does not matter whether we take to or to – this additional complexity occurs with spheres as they lack any intrinsic orientation, which is not the case for lines, planes, circles and point pairs.

### Point Pairs

In the conformal setting, point pairs take the form 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 .

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

 (13)

 (14)

with .

### Lines to Circles: Planes to Spheres

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.

## The Non-Uniqueness of the Recovered Rotors

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

To investigate this further we extract the bivector, , 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.

## Conclusion

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

# Direct linear interpolation of geometric objects in conformal geometric algebra

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

## Abstract

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

## Introduction

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

The objects we work with here will be CGA objects unless explicitly stated otherwise. We will use the standard extension of the 3d geometric algebra, where our 5D CGA space is made up of the standard spatial basis vectors , 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 .

## Motivation

In our conformal representation of space the blades of a given grade 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?

## Linearly interpolating conformal points

The result of linear combinations of conformal points is well known Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]. Consider two arbitrary points in 3d space 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 :

 (15)

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 ):

## Linearly interpolating higher grade conformal objects

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

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

 (16)

where 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 [2003] and point pairs Dorst et al. [2007] 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:

 (17)

where 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 .

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

Let , where is a valid object (squaring to 1). Now define , ie. it only contains 0 and 4 grade coefficients and is an element of the conformal algebra. Then, defining , the square root can be found as:

 (18)

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

 (19)

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

 (20)

It now simply remains to isolate via multiplication by where . Since is a scalar, we have

 (21)

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.

## Techniques for understanding interpolant properties

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

From this we immediately see that for the multiplication to be grade preserving we require and to give only objects of grade where is the grade of and . Table 3.1 shows the resultant grades from the geometric product of pure grade objects and Table 3.2 shows the resultant grades from the inner product. These tables are presented here for reference and will be returned to when dealing with individual grade blades.

Table 3.1: Resulting grades from the geometric product of pure grade objects. The grade 4 row is highlighted here as it is of specific interest for the form of the in equation (3.3).
 0 1 2 3 4 5 0 0 1 2 3 4 5 1 1 0,2 1,3 2,4 3,5 4 2 2 1,3 0,2,4 1,3,5 2,4 3 3 3 2,4 1,3,5 0,2,4 1,3 2 Gray 4 4 3,5 2,4 1,3 0,2 1 5 5 4 3 2 1 0

Table 3.2: Resulting grades from the inner product of pure grade objects. In the case of the inner product of a multivector and a scalar by definition the result is always 0 rather than some other scalar valued function of the scalar and multivector. Here again the grade 4 row is highlighted here as it is of specific interest for the form of the in equation (3.3).
 0 1 2 3 4 5 0 0 0 0 0 0 0 1 0 0 1 2 3 4 2 0 1 0 1 2 3 3 0 2 1 0 1 2 Gray 4 0 3 2 1 0 1 5 0 4 3 2 1 0

## Point pairs

We start with point pairs. Previous work Dorst et al. [2007] has shown that when an end point is shared between point pairs 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:

Theorem 1   If point pairs and are both chords of a common circle the interpolant point pairs are blades and also have end points lying on as and .

Note, since , the projector is a scalar. The common circle itself is the join' of the two original point pairs and can be computed with the algorithms supplied in Chapter 21 of Dorst, Fontijne and Mann Dorst et al. [2007]. Figure 3.2 shows two cases of the interpolation of 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:

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

as we see that

Now we just need to prove that is a scalar multiple of . From equation (3.3) we know that

it is therefore sufficient to prove that:

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

As we can write this as:

and this can further simplified using the fact that:

As is a scalar:

As is a scalar we see that the proof is completed

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

Theorem 2   For non-coplanar point pairs and , all interpolant point pairs lie on the surface of the sphere .

## Circles

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

Theorem 3   If circles and together define the caps of a common sphere then where is of the form shown in equation (3.3) and thus any interpolant object is a blade without requiring re-projection to the blade manifold.

This can be proved by considering each circle as the intersection of a plane and a sphere . Forming this intersection via the dual (where and is the 5D space pseudoscalar as before), we have:

Since and are both bivectors:

and so if :

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

Conjecture 1   If circles and together define the caps of a common sphere then . can be found by the following process:

First we define:

where is the truncated identity matrix that performs selection of grade 3 elements from a vector of coefficients and and are the matrices that perform the left geometric product of and respectively with a vector of coefficients. We can then find the for where is a vector of canonical blade coefficients limited to only the 4-vector blades. In the case that and are the same radius then .

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

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

Theorem 4   If circles and together define the caps of a common sphere then all interpolant circles (which we have shown to be blades) also lie on the surface of the sphere common to both.

We can prove this by considering the outer product of the interpolant circle with , an arbitrary point on the common sphere :

Figure 3.5 shows an example of this interpolation.

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. [2007]:

Theorem 5   If circles and together do not lie on a common sphere then the 4-vector from our blade projection equation is itself a blade and geometrically represents the sphere through which both circles plunge orthogonally. ie. . This property means all interpolant circles after projection to the blade manifold also plunge through orthogonally. ie. .

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:

Conjecture 2   If circles and together do not lie on a common sphere then the intersection point pair formed by the meet of the circle interpolant for a given value of with the orthogonal sphere ie. is the same as the re-projected interpolant of the point pairs formed from the meet of and with .

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

## Lines

When looking at lines we can attempt to use some of the same techniques that we used for circles. First consider the form of . For lines , giving the form of the projection of as:

 (22)

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

### Screw Theory

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

Screw transformations consist of a translation along an axis and a rotation around that axis. To parameterise a screw we define the direction of the screw axis via a unit vector , 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.

### Bivector representation of a line

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

 (23)

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

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

### The bivector representation of a screw

To represent a screw we will couple the rotational motion of the dual of a line with a translation in the direction of that same line. The bivector that transforms along the 3d vector is:

 (24)

If 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.

 (25)

Hestenes and Sobczyk Hestenes et al. [1985] p81 gives an expression for decomposing any bivector into two commuting blades. In the case of our screw bivector these blades represent the dual of the screw axis and a translational bivector in the direction of the screw axis. ie. given a screw bivector we can decompose it as:

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

we can write this elementwise as

where is . We then rearrange to give something proportional to the expression in equation (3.11):

Where clearly . If we divide this by we have the general form of a normalised screw

Gathering like terms, specifically those without an component, leads us to the conclusion that our screw axis direction must simply be proportional to the addition of the directions of the two lines. Using this fixed axis direction we can extract the coefficient (the pitch) of the translation bivector parallel to the screw axis:

With this coefficient known we now have all the pieces in place for a full decomposition of the dual line addition bivector :

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

### Relationship to object manifold reprojection

We can also analyse the screw multiplied by its own reverse, comparing this formulation with our object manifold reprojection to get the projection coefficient in terms of the screw parameters:

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

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

## Planes

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

The dual of the plane can be written as:

where is the 3D vector normal to the plane and is the perpendicular distance of the plane from the origin. Thus the interpolation of duals of two planes can be written as:

which, when we collect like terms, is already in the form of a dual plane :

 (26)

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

## Spheres

The interpolant of spheres has been studied before in Cameron [2007] and Doran [2003]. As with planes, all interpolants of spheres are valid objects as 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:

The dual form of a sphere can be decomposed into the sum of the conformal centre point and negative half the radius squared times :

the interpolation of the dual of two spheres is therefore

For two concentric spheres ie. we can therefore see that the interpolation between them will remain centred in the same place and will simply have a radius which varies as .

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

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

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

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

and so the radius varies as

For fixed values of and this implies varies as and so the further apart the two spheres are the smaller the radius of the interpolant. To find turning points we differentiate with respect to

setting this to zero yields a single turning point at

Considering the second derivative

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

For the case that the surfaces of and are just touching we have the condition

returning to the first derivative in this case

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

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

## Applications

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

### Higher order spline interpolation through objects

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

### Recursive scene simplification by averaging conformal 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:

1. Set a minimum cost threshold for difference between objects
2. Compute the cost between all objects of the same grade in the scene
3. If all costs are above the threshold then terminate the algorithm
4. Average the two objects with the smallest cost

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

 (27)

where , and gives the component of having as a factor and is the rotor that takes to as described in Lasenby et al. [2018]. An example of this algorithm working on simulated lines is shown in Figure 3.14.

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

### -means clustering of conformal objects

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

The steps for implementing this clustering are given below:

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

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

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

Consider two non-intersecting non-coplanar lines in 3d space, 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. [2004]:

## Conclusions

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

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

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

## Abstract

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

## Introduction

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

## Conformal Geometric Algebra, CGA

The ray-tracer used in this chapter is constructed using CGA and all algebraic expressions given will be in terms of elements of this algebra. CGA adds two more basis vectors, 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. [1985]. The mapping from a 3D vector, , to its corresponding CGA vector, , is given by:

 (28)

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

## Camera Model and Ray Casting

A pinhole camera model is used with the geometry shown in Figure 4.2. It is defined by a rotor (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.

 width=

 width=

We take to be at the bottom left hand corner of the image. For an image of width and height , the world coordinates of the point at the centre of pixel are given by:

We then generate the ray from the camera centre, , that passes through , via the expression

where is the origin transformed by the model-view rotor to the position of the camera.

## Ray Geometries for Basic Objects

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

### Ray-Object Intersections

In order to compute intersections between blades, the meet operator () 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 :

 (30)

where indicates the -grade component of the multivector , and represents the 5D pseudoscalar of the algebra.

#### Planes

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

 (31)

Therefore, can be extracted from the and coefficients (for ) by dividing by , the coefficient.

#### Spheres

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. [2004]:

where are scalar constants. If we define with oriented in the same direction as our ray , then is the point closest to the origin of the ray, , as long as our sphere is in front of' the ray source. To ensure the alignment we can pre-normalise our sphere via the following expression:

 (33)

For any given sphere, its dual can square to a positive or negative number; however, by carrying out the normalisation in equation 4.6, we ensure that all spheres, , 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 .

#### Circles/Discs

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

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

Figure 4.5 shows an example of each case along with a geometric interpretation of the form of the meet. If and there is an intersection, the plane containing the circle is formed by taking the wedge product between the circle and , , and the intersection point is then extracted from the ray and this plane.

 width=

If , so that the ray lies in the plane of the circle, we need to work in 2D, so that our meet' will result from taking the 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.

### Extracting Normals and Reflecting Rays

Extracting the normal to the surface of an object at a ray intersection point, , 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. [2004]:

where is the first point of intersection. Here, is an example of an inversion, where the incoming ray/line, , is inverted in the sphere to give a circle which passes through the two points of intersection and the origin of the sphere (note we only have a meaningful reflection if there are two points of intersection). The tangent line to this circle at the first point of intersection, , is the reflected ray, . Figure 4.6 illustrates this geometrical construction. Note that one can also form the tangent plane at and reflect  in this plane; this is performed by the following formula:

 (35)

For a circle/disc , we first form the plane in which it lies. If the ray intersects the disc (see section 4.4.1), the reflected ray can then be found using the same formula as for the plane reflection case: . Note that these expressions specifically give the (correctly oriented) reflected ray that passes through the point of intersection of the incident ray and the object, rather than a parallel ray at the origin.

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

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

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

## Ray Tracing Evolved Circles

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

 (36)

where we take moving between 0 and 1, which moves us from to . The result of this interpolation is not itself a valid circle and needs to be projected' onto a blade via multiplication by a projector, which we shall call  . This projector has only scalar and 4-vector parts and its construction is detailed in Hadfield and Lasenby [2019] and outlined in the following.

Consider a quantity . We then define the quantity , and with this the principal square root Dorst and Valkenburg [2011] of the scalar + 4-vector, , can be found as:

 (37)

With this square root we can then form:

 (38)

where 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 :

 (39)

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

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

### Intersection Point of Ray and Interpolated Surface

We saw earlier that the intersection of a ray with a circle produces the 1-vector . 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:

The system must also be tested for the case of ; if an exists such that , the ray may intersect once, twice or not at all.

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.

#### Non-linear Intersection Point Finder

As it is in general difficult to extract a closed form expression for the solution to , it is necessary to design an iterative algorithm to find the roots of the equation. Our implemented algorithm works as follows:
1. Check for intersection with a sphere enclosing the entire surface
2. Calculate the value of at N intermediate values of
3. Record where changes sign between successive evaluated values of
4. Locally approximate as a quadratic equation in the region of the sign change and solve to get the value of at the intersection point
Computing the intermediate objects in the surface can be done once per scene and reused for all rays calculated in that scene. To generate the enclosing sphere again we reuse the intermediary objects, in the following way:
1. Given and , form intermediate circles, , and then calculate the bounding sphere which is given by
2. Construct a sphere that contains all intermediate circle bounding spheres by successive application of a two sphere bounding algorithm
Again the enclosing sphere of the object can be calculated once per scene and used for all rays. Any two sphere bounding algorithm can be used, here we chose the algorithm from Hildenbrand and Hitzer [2008] which is summarised as follows:
1. Ensure both spheres are normalised according to equation 4.6
2. Construct the line joining the centres of both spheres
3. Intersect the line with the first sphere to produce a point pair and extract using equation (4.5)
4. Intersect the line with the second sphere to produce a point pair and extract using equation (4.5)
5. Check if , if so encloses and so is the bounding sphere
6. Check if , if so encloses and so is the bounding sphere
7. If neither original sphere is enclosed by the other, the new bounding sphere is given by

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

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

In this section we will use 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:

Writing , this can be rewritten as:

 (40)

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

Now expanding as:

means we can write:

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

The quantity is a scalar and so distributing the grade selection operators gives us:

 0 (41)

To make progress on solving this we recall that and that for our linear interpolation of circles we have defined as:

As has terms in of order 1 we would expect to have terms of order . Continuing on this train of thought one might suspect that it is possible to re-write equation (4.14) as a simple polynomial in . However, it is easy to see that this is not possible due to , which is a scalar polynomial in enclosed entirely in a square root:

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

 (42)

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

 (43)

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

For this case of linear evolution of circles we will get a polynomial of order 12, implying 12 potential roots. In reality 6 of these roots are extraneous, generated by the process of squaring to handle the square root term in . 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.

#### A Comment on Rendering Speed

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

### Analytic Form for Normals

Given the 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. [2004], which is also used in equation 4.7, we extract a tangential line in the plane of the circle at :

 (44)

We would now like an analytic form for the tangent to the surface corresponding to evolving the surface through an increment of , 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:

 (45)

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

 (46)

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

 (47)

the condition is satisfied since . Thus, given an on the surface, lying on a circle with parameter , the defined above will preserve its length and is, we claim, the tangential direction required. In order to show this, the first thing we must do is establish that if we evolve according to this rule, generating a quantity we call , then should lie on , for all , i.e.,

Differentiating this (and again dropping the subscript for clarity) and using , gives

We now expand this expression using standard expansion results ( ):

The first term on the right hand side of the first line of this expansion is zero as . Since lies on and so , we see that which means that

giving as required, so the proposed evolution is compatible with the constraint.

If we therefore assume that is the direction we want, we can calculate the tangent line in this direction via:

 (49)

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

To show this, we need to consider the reverse of . We will need the facts that that , , , and . We have shown all of these identities earlier in this section. We also need an additional fact, which is that anticommutes with . To see this we use another standard result ( ):

 (50)

The first term on the RHS of this equation is zero as , and the second term on the RHS is also zero as lies on so . Thus and therefore anticommutes with as required.

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:

 (52)

## Calculating the Derivative of the Object Manifold Projection

To calculate 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:

where is our blade projector. The differential of this with respect to is given by:

 (53)

Thus, any closed form expression for the derivative on the manifold will first require a closed form for the derivative of the projector . Recall from equation 4.13 that we can write the projector as a function of , where , and so

 (54)

Thus to find an expression for we will first need one for .

### Closed Form Derivative of the Square Root Operation

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

 (55)

where we are using the fact that .

 (56)

Thus the derivative of is a function only of and which in turn can be written in terms of and :

### Closed Form Derivative of the Projector

With our square root derivative in place we can proceed to finding the derivative of the projector . Recall, is given by:

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

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

 (57)

Now consider to be an interpolated circle of the form . The derivative of this with respect to is a constant:

This derivative is the final piece required for equation (4.26), giving us a completely closed form for .

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

## Ray Tracing Evolved Point Pairs

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

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

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

To find the intersection point of a ray and these surfaces we again form the meet of a ray and the form of the evolved point-pair . The result is a scalar quantity that can be written as:

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

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

Expressing in terms of where as before, , gives

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

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

Squaring then gives us:

allowing us to form a simple scalar polynomial in (we can see that this produces a scalar equation since has only bivector parts):

 (58)

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

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

### Bounding Sphere and Normal Calculation

We saw earlier that the meet will be zero if the ray hits anywhere along the carrier line of the point-pair . Assuming the carrier line and ray do meet, the point of intersection can be extracted via the method outlined in the last section of Chapter 3 and published in Hadfield and Lasenby [2019]. Given that the carrier line of (for some ) and the ray intersect at a point , we can then check if the intersection point is within the bounding sphere of the surface by ensuring:

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

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

### Special Cases of Evolved Point-Pairs

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

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

Re-arranging gives an expression for :

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

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

### Triangular Facets from Evolved Point-Pairs

The intersection of co-circular point-pairs also allows us to examine the intersection of rays with triangular facets. Consider a ray and a set of three points , , which together form a triangular facet. First we will form a set of normalised point pairs:

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

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

## Bézier Curves and Hermite Splines through Geometric Primitives

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

### Linear Interpolation as a Linear Bézier Curve

The simplest form of Bézier curve is simply a linear interpolation between two vectors. If we replace the vectors with -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:

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

With three multivectors we can specify a quadratic function of :

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

### Cubic Bézier Curves

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

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

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

 [Order 1] [Order 2] [Order 3]

### Order Bézier Curve

More generally we can say that an order multivector Bézier curve is of the form

where

are known as the Bernstein polynomials. The derivative of our order Bézier curve is:

where

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

### Rational Bézier Curves

A rational Bézier curve adds weights to the polynomials allowing them to represent a broader class of curves:

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

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

### Hermite Cubic Curves and Splines

Hermite cubic curves are another common form of interpolating curve. They are defined by control points, , (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 :

The derivative of the curve is:

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

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

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

 (60)

Let us now evaluate this at the endpoint of the curve in a piece-wise spline where and of curve () where . First we note that on both curves at these points, because the curve passes through a blade control object which requires no projection. Additionally we see that by definition of the Hermite spline at that point, the derivative in pure-grade space is shared across both curves as is the control point:

where and are the derivative (or tangent') and control objects respectively of the curve that are shared between segments and .

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 [1984]. The Kochanek–Bartels (KB) spline is an interpolating spline with three scalar design parameters known as tension, bias and continuity respectively. For given control objects the corresponding tangents can be calculated using the control objects in the spline and which lie previous to, and after, the curve in the order of the spline:

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

## Examples of Ray Tracing Simple Objects and Evolved Surfaces

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

## Meshing Evolved Surfaces

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

To produce a mesh approximation we first need to generate a set of points that are in some sense evenly spaced and lie on the surface itself. To do this we will begin by producing a set of evenly spaced points on the first object 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:

• Calculate and the translation rotors that bring and respectively to the origin
• Apply and to and respectively, bringing them to the origin and producing and
• Calculate the rotation rotor between the blades and (if and are spheres then we do not need a rotation rotor so set )
• Calculate the difference in scale between the objects by extracting their relative sizes
• Use the scale to generate a dilation rotor that scales to the same size as
• Compose the final TRS rotor that takes to as:

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 :

where is chosen so that uniformly spaced points cover the whole circle. With the TRS rotor that maps from the unit circle at the origin to the first object it is possible to transform our points to the first object:

Our process then becomes one of stepping sequentially through from 1 to 0 in small increments of and transforming our points along the way. We will use to refer to the TRS rotor that maps to and so can write the point at as:

By doing this we have effectively constructed a mapping from a coordinate system in the 2D plane of and to the surface manifold. This is useful as it lets us generate the mesh in the 2D plane of and and map the vertex positions directly to 3D.

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

## Summary and Conclusions

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

# REFORM: Rotor Estimation From Object Resampling and Matching

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

## Abstract

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

## Introduction

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

The objects we work with here will be CGA objects unless explicitly stated otherwise. We will use the standard extension of the 3D geometric algebra, where our 5D CGA space is made up of the standard spatial basis vectors , 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 .

## Proximity-based matching

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

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

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

 (61)

where indicates the -grade part of . Equipped with this idea of closeness of objects, for a given , a query object is assigned to each of the reference objects (i.e. this is done for all ), assuming the model and query sets are spatially close. For each object pair we form the rotor, that takes the query object to the reference object. The minimum cost assignment is then taken as the correct match, , for that query object

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

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

## Finding the rotor between two sets of matched objects

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

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

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.

## Iterative matching and rotor estimation

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

1. Each object in the query model is given a match in the reference model (there are a number of ways of making this initial guess)

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

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

## Incorporating sampling

To counteract the local minima issue, we modify our procedure to incorporate sampling in a RANSAC-like Fischler and Bolles [1981] algorithm. This particular approach is chosen as it is readily adapted to parallel processing and is well suited to handling large numbers of incorrect matches. After each matching stage in the previous algorithm we randomly and uniformly sample 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:

1. Each object in the query model is given a match in the reference model (there are a number of ways of making this initial guess)

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

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

## Matching scenes of mixed geometric primitives

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

## Conclusions

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

# Kinematics, Dynamics and Robotics

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

# Screw Theory in Geometric Algebra for Constrained Rigid Body Dynamics

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

## Abstract

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

## Introduction

### Screw Theory

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

### CGA

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

### PGA

While CGA has been a very successful idea it is by no means the most computationally efficient way to embed rigid body dynamics into GA. In many applications a practitioner is willing to lose some of the expressiveness and mathematical niceties of CGA to gain an algebra with fewer elements that nonetheless does what they desire, and crucially, does it computationally faster. The Plane-based (or Projective) Geometric Algebra (PGA) aims to strike a balance between efficiency and expressiveness. With a signature of it contains the so called flat' elements of CGA, flat points, direction bivectors, lines, planes and the rigid body rotors, but loses point-pairs, circles, spheres, dilation and inversion rotors as first class citizens of the algebra. The use of the degenerate metric also necessitates the use of an alternate dual map rather than a simple multiplication by the pseudoscalar, in practice this dualisation operation can be implemented by the right complement extended over the canonical basis vectors by linearity. For more information on PGA the reader is directed towards the 2019 SIGGRAPH course by Gunn and De Keninck Gunn and De Keninck [2019a,b] and Leo Dorst's A Guided Tour to the Plane-Based Geometric Algebra PGA' Dorst [2020].

## Forces, moments and static equilibrium

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

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

### What is a force?

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

#### 3D vector representations

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

#### 6D representations

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

This object now contains all the information required to use the force in applications. It is however not unique. We could choose multiple values of all of which lie on the line and they would all give different 6D representations of the force while still fundamentally behaving the same way physically. To solve this problem we need to instead use a representation of the force that is unique. If we change our representation to:

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

In contrast to a force, a moment is not localised to a specific point or line in the world, in the terminology of physics and engineering this type of non-localised vector is known as a free vector. While respecting their non-localised nature we would like to be able to put the 3D moments in the same 6D box as forces to be able to keep the representations consistent and maybe provide some meaningful operations on them. When putting our 3D moment vector into the 6D box we are left with several different options. One option would be to put the moment vector in the top 3 positions of the 6D vector, such that they align with the 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 this situation, as the forces are in opposite directions and of equal magnitude, the total resultant force on the body is . But what about the total moment of the forces about a point on the body?

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

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

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

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

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

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

gives a scalar:

where here the operator is the standard inner product between vectors. The second product we will consider is the screw cross product, for screws and the screw cross product gives:

where and are computed with the standard 3D cross product. The result of the cross product of two screws is itself another screw.

### Representations of wrenches in CGA and PGA

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

### Forces as dual lines in CGA

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

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

 up   up (62)

where we use the notation up to represent the mapping of a 3DGA point to a conformal point, ie.    up. The object has a magnitude equal to the intensity of the force:

Now consider the dual form of this CGA line:

 (63)

This looks very similar to the 6D Plücker representation, in fact if we consider it closely we can see that it is the same. First, look at the term in front, . This is the 3D dual to a 3D vector, giving a bivector, specifically the Euclidean bivector orthogonal to . Now consider the second half of the formula: . is the 3D dual to a Euclidean bivector, ie. a vector equal to , just like the lower 3 slots of the 6D screw representation. As is a Euclidean vector this makes the form of a CGA direction bivector'. An important thing to note here is that the front and back part of this formula are linearly independent, just like in the 6D screw representation, ie.

Again as with the 6D representation we can consider two anti-parallel forces and extract the representation of a force couple or moment. Setting and in the above leaves us with:

which we could re-write as:

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

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

The geometric product between them gives:

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

Now collect the terms by grade:

where

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

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

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

If we rewrite the GA formula here:

then equating terms gives:

and

The form of the standard vector cross product in 3D GA is, for vectors and , given by:

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

For readers well versed in Screw Theory, Clifford/Geometric Algebra and Lie Theory this is a well known result. The commutator product equips the motor bivectors with a Lie bracket just as the cross product equips the screws with one. In fact we can go further into the idea of the motor bivector as an element of and note that the grade 0 element of the geometric product, equal to the GA dot product between the bivectors, is proportional to the Killing form Minguzzi [2013], , of the Lie algebra ie. .

### Forces as lines in PGA

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

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

In this form the line looks very similar to the CGA representation of the line. The line squares to a scalar, there is a section that is a Euclidean bivector and a section that is a null bivector. In fact the only thing that we have changed is the form of the null element. In CGA we typically use the null element constructed from the sum of two orthogonal basis vectors, one squaring to and one to . Here in PGA the null element is itself a basis vector. For our 6D force line representation we therefore have exactly the same mapping as we did in CGA:

and so moments appear as:

We would also expect the various products of Screw Theory to still appear as the various grades of the result of a geometric product between two elements of this algebra. Indeed they are the same as the CGA results, but of course replacing with which is now the pseudo-scalar of the algebra.

### Force and moment representations in the GA literature

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

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

this is in the form of a direction bivector.

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

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

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

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

This is the bivector form of a CGA line in the direction of the 3D vector and passing through the point .

Rigid Body Dynamics and Conformal Geometric Algebra' Lasenby et al. [2011] therefore uses a mixture of 1-vectors and bivectors to represent forces but in the end their static equilibrium conditions are in the form of bivectors and trivectors, indeed they are left with something very similar to our formulation in section 6.2.2. The main difference in fact is related to which of the two orthogonal elements of the wrench formulation we take to be a moment and which to be a line. Lasenby et al. effectively choose the 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:

where . The properties of this object are that it squares to 0, ie. it is null but that the magnitude and direction of the 3D vector are the magnitude and axis of the moment. This is identical to the formulation for bivector moments that we described in section 6.2.2.

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

## Screw transformations, instantaneous twists, and the motor manifold

### Time derivatives of frame transformations

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

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

We will take time derivatives. As is fixed:

Now substitute in and :

A rotor by definition has the property that . If we differentiate this constraint with respect to time

which means we can write:

which is twice the anti-commutator of and . If is a 1-vector and as is a bivector we can write:

If we now label our bivector quantity:

 (64)

it allows us to write:

 (65)

where represents the commutator product. Note that in this form with the commutator product no assumptions are made of the grade or other properties of . If we do restrict to 1-vectors however, this allows us to write:

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 :

 (66)

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

 (67)

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

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

 (68)

and we will also calculate the reverse for reference:

## Momentum and inertia

### Screw momentum

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

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

### Mapping from screw velocity to screw momentum

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

When it comes to angular velocity, , and angular momentum, , however we have a more complicated relationship. In fact, for a body centred and axis aligned reference frame, the two are related by a diagonal matrix known as the inertia tensor that we label here as :

where can be expressed in the form of , the second moments of volume:

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:

 (69)

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

 (70)

We therefore have a direct mapping from the wrench acting on the body to the initial acceleration of objects in the frame and we can now use our physical intuition about the world to guide us to a compatible form of 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 [2016] which we will come to later when considering virtual work and power. What we mean by a reciprocal frame here is a set of reciprocal bases that are matched to the motor bivectors such that when the GA inner product is taken between the two matched elements the result is 1 but otherwise 0, ie.:

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

As is a linear function operating on a motor bivector in CGA it can be written as:

 (71)

where are scalar coefficients, is the reciprocal frame for the motor bivectors and are the motor bivectors themselves. We will take and as follows:

and will seek to determine . This formulation is equivalent to a square matrix formulation of where the coefficients are the elements in the matrix. What we need to do here is to find out exactly what these coefficients have to be.

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

where in this case is 0 as we are at the origin, giving:

The mapping that embeds a 3D point in CGA is given by:

which, after differentiation twice leads to:

As the body is initially stationary, is zero and so we have for time :

For the case of a force of magnitude applied to the centre of mass of a rigid body we would expect all points on the body to accelerate linearly at in the direction of , or more formally for all . We can encode this expectation by substituting and into equation (6.9):

 (72)

and we can then apply our specific case:

Looking closely at

we can see that is a Euclidean bivector and, when dotted with each of the reciprocal frame elements in turn it is non zero only for in which cases the dot product gives only , the component of the force in the direction. In other words we could write:

where and . We can combine this with to give:

Let us now examine the terms on the right hand side of the above equation:

For we have the results that , , . For we have the results that , .

Now we can consider the above equation component-wise and extract the required transformation coefficients. Firstly, consider the coefficients of :

Expand the left side into a sum of elements of the force vector:

Comparing terms we see we have arrived at the following solution:

Now that we have dealt with the terms we can return to the Euclidean bivectors:

We have just calculated the terms in the right hand side of this equation and so by substituting these in we get something of the form:

which implies that . A visual representation of the argument made here can be seen in Figure 6.1.

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:

 (74)

We will now apply a moment to the rigid body, again at rest at the origin. From standard 3D kinematics we know that if we have a body rotating about its centre of mass and about an axis 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:

and the acceleration of a point as:

where is the traditional cross product operation. In this thought experiment the body is initially at rest allowing us to remove the term and leaving us with:

where is the 3d torque vector. The GA equivalent of the traditional 3d cross product for vectors and is . Which means we can write:

We can represent this same moment in our CGA formulation as the bivector wrench . Again we look at:

and from this we note that for and that for . For this case we can therefore write equation (6.13) as:

The left hand side of this equation collapses simply to and we bring inside the summation again on the right:

We can now use the same results as noted previously to simplify the right hand side of this equation, leading to:

Equating like terms eliminates all the terms, ie. for . This leaves us with:

We can break up the left side of this component-wise:

Now assume that our torque vector is aligned with one of the principal axes of the rigid body, say , this implies is only non-zero for :

Which is true because . Noting that allows us to identify the parameters: , , .

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.

### The Screw Inertia Tensor

For our 6D vector representation had we not just done the maths of the previous section we might have been tempted to stack the and on top of each other to produce a combined inertia tensor like:

This however would ignore the fundamental conceptual differences in the action of the top 3 and bottom 3 elements of the 6D representation as a wrench vs as a twisting motion generator. Instead, we need to use the coefficients of the last section to produce a matrix that performs a flip' in the relative positions of the parts of the vectors while also applying the required scaling along each principal axis. Putting the calculated coefficients in place we can see that the matrix that comes out for looks like this:

We could also achieve this flip' effect via a diagonal matrix and a permutation matrix:

 (75)

where is the identity matrix.

When writing the inertia tensor in CGA it is convenient to do a little relabelling for ease of reading. The reciprocal frame of the motor bivectors in CGA is as follows:

Which we can break into the two groups:

and so with these we can write the inertia tensor:

and the inverse inertia tensor:

This inertia tensor performs the same kind of flip' that we saw for the 6D screw representation. Effectively our inertia tensor does the following mapping:

The reciprocal frame construction of the inertia tensor that we have discussed so far works for many algebras but for degenerate metric algebras such as PGA the fact that we have an element squaring to zero means this setup does not work. Instead we need to do something a little different. The degenerate metric approach to reciprocal frames Gunn [2020] is to consider some blade which we will label that wedges with a given blade of magnitude ie to produce the pseudoscalar with magnitude ie. . Clearly, despite us labelling it , this object is not quite the same as the reciprocal frame, although it allows us to perform the same function of coordinate free coefficient selection producing the magnitude as the scalar coefficient of the pseudoscalar. Let us now identify this pseudo-reciprocal frame for the PGA bivectors:

Comparing the PGA pseudo-reciprocal frame mapping with that of our CGA-inertia tensor mapping it is immediately clear that they are equivalent up to a minus sign. We will now define a function to perform this PGA mapping and will call it . has the following action:

where the syntax returns the scalar coefficient of in . We can extend this operation to combinations of basis elements by linearity so that for :

As with our CGA reciprocal frame let's now write our PGA pseudo-reciprocal frame in two groups:

This means we can write our PGA inertia tensor as:

and its inverse inertia tensor:

We could also apply the map first to first flip' the input and apply a component-wise scaling to the result:

This map first formulation is conceptually the same as the screw formulation with a flip permutation matrix as in equation (6.14).

### Motor Bivectors as the Principal Screws of Inertia

We can visualise the motor bivectors as a frame of screws attached to the origin. These motor bivectors are a version of the Principal Screws of Inertia, specifically they are the principal screws in a Plücker and Hunt sense as opposed to Ball's original formulation of the principal screws. Effectively what we are doing in the inertia tensor is considering these principal screws as wrenches and analysing the impact of them on the motion of the body.

By considering the bivectors as localised screws we can begin to build intuition about them and their properties. First of all, we will consider how they, and their reciprocal frame, transform under the action of rigid body rotors. The Euclidean bivectors are in the form of a dual line through the origin and are affected by rotors exactly as lines are. The direction type motor bivectors ( ) 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:

We will wrap this operation up as an inertia tensor in its own right and label it , we can therefore write:

 (76)

The problem of defining non-axis aligned inertia tensors is, given we know and , exactly what form does take?

The form of is:

Based on the fact we can transform screws and their reciprocal frames with rotors as usual, we will guess the form of to be as follows:

 (77)

Note this is just the same as except that we have transformed all of the motor bivectors and reciprocals by , effectively transforming the frame of principal screws to be centred and aligned with the rigid body but expressed in the world frame. We can check whether our guess is correct as follows. Substitute :

Noting that and we can write:

We can then take the rotor application outside of the summation by factorisation and we have arrived at the required relation of equation (6.15):

Of course by allowing the movement of the frame of the inertia tensor with the body it will no longer be constant and we therefore might also want the form of the time derivative for potential applications:

 (78)

In this subsection we have analysed the principal screws as CGA objects and our analysis of the translation of the inertia tensor was phrased using the transformation of the reciprocal frame. In reality we could equally have done our analysis from the PGA perspective as well, with our pseudo-reciprocal frame taking the place of the reciprocal frame and so our equivalent transformed inverse inertia tensor would appear as:

## Unconstrained rigid body dynamics

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:

and its first time derivative is:

where is the momentum bivector at time expressed in the world frame and is the resultant external wrench acting on the body expressed in the body frame. From this point on we will drop the subscript and simply state that all variables are functions of time. From our discussion in the previous section we know that we can further expand using the inverse inertia tensor :

Re-writing the time derivative of the state with this equation for gives:

## Constrained dynamics via virtual power

Unconstrained dynamics, while important, do not allow us to represent all the types of motion that we see in the real world around us. In many practical situations we are faced with the problem of constrained motion. Consider modelling a rigid body that can move dynamically under external forces but is constrained so that one or more points lie on a surface or a situation where a rigid body is constrained such that it can translate but not rotate. These are the types of problem we will attack here.

To impose a constraint on our dynamics model we will use the concept of a reaction wrench. The reaction wrench provides a combined external force and moment that acts on the rigid body in addition to the other external wrenches and, in doing so, forces the body to move in a way that respects the constraints. We will write 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:

and is of the form

where is a virtual scalar power. Differentiating this gives:

We can now substitute in our dynamics equation for :

Setting the virtual power and rate of change of virtual power to 0 gives us the virtual power condition for our constraint. First:

tells us that the virtual wrench must be parallel to the screw velocity. Setting the rate of change of virtual power to be zero allows us to write:

If is a static constraint we can specify that leaving us with:

Which can again be solved for and hence .

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:

Taking a time derivative of this we see:

As is driven by the rotor , ie:

we get:

and so:

We can then directly substitute this into:

and so calculate . To constrain this same point to a circle we would add an additional planar constraint, ie. the point must lie on the plane in which the circle lies and on the sphere of which the circle is the equator.

## Constrained dynamics by pinned multivectors

Consider a geometric primitive represented by multivector 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 :

or equivalently:

Taking first and second derivatives gives us the expressions:

 (79)

 (80)

The next step is to think about what these expressions mean physically. Essentially we have two views' of the same object, one in body space and one in world space. For example we can imagine the is a point attached to our rigid body and is a point in the world that that point is also attached to. In a sense we are pinning' the rigid body to by its extremity .

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:

 (81)

This equation is a constraint on the second time derivative of that will ensure that and do not vary with time. We can go a step further here and substitute our expression for from equation (6.7), leading to:

Now if we substitute in and :

Simplify and gather, substituting :

Now separate the terms with :

As are bivectors their reverse is just a negation:

 (82)

So far, we have done a lot of algebra but so far appear to be no closer to calculating our reaction wrench. If we calculate an expression for however, we start to make headway towards a solution:

using we can also write:

 (83)

and so we can now substitute in on the left hand side of equation (6.21) for :

Now we separate out the terms with

and take all terms not containing onto the right side of the equation. We now have something of the form:

Some function of

We can rewrite this to use the commutator product:

 Some function of (84)

If we now decide to write our total bivector wrench, , as the sum of external wrenches, , plus some reaction wrench, , caused by the constraints:

Some function of

we now have a constraint expression that fixes the reaction wrench as a function of the state of the system and the forces applied to it.

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.

## Geometric objects as constraints

Now that we have identified a means of enforcing constraints via pinned geometric primitives let us have a look at exactly what constraints are imposed by specific choices of this pinned multivector for CGA and PGA.

### Point constraint

If we want to pin a specific point in our rigid body to a point in the world we can set an invariant point constraint. In this case 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.

### Point-pair constraint

Consider a bivector of the form:

where and are CGA points. This object is invariant only under the action of rotation rotors about an axis parallel to the line joining and . In this case has 5 degrees of freedom corresponding to 3 translational forces and 2 moments. This point-pair constraint is specific to CGA.

### Direction constraint

Consider a so called direction' bivector in CGA of the form:

where is a 3D vector. This object is invariant under the action of all translation rotors and is invariant to rotation rotors with axis of rotation parallel to , ie. rotors of the form:

where is the pseudoscalar of 3D space. In this case only has 2 degrees of freedom corresponding to two moments with axes perpendicular to d. For a PGA equivalent of this constraint a bivector of the form:

can be used to achieve the same thing.

### Flat point constraint

Consider a so called flat point' bivector in CGA of the form:

where is a 3D vector. This object is invariant under the action of all rotation rotors about the point , but is not invariant to translation. Thus, under the action of the rigid body rotors it behaves like a CGA 1-vector point. For PGA this constraint can again be implemented by a standard PGA point.

### Line constraint

A line is invariant to translation along the line and rotation about the line axis. Thus we would expect to be able to support reaction forces orthogonal to the line, and moments with axes orthogonal to the line axis, ie. 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.

### Circle constraint

A circle is invariant to rotation about the axis of the circle only. Thus we would expect to be able to support reaction forces in all directions and all moments other than the one with axis parallel to the circle axis, in other words 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.

### Plane constraint

A plane is invariant to translation in plane and rotation about axes parallel to its normal direction. Thus it can support one direction of force and two directions of moments giving 3 degrees of freedom to . Both CGA and PGA can use planes for pinning.

### Sphere constraint

A sphere is invariant to all rotation about its centre, but not to translation. Thus it acts like a point at the sphere centre and can support translation reaction forces only giving it 3 degrees of freedom. CGA can represent spheres directly as the outer product of 4 points. In PGA this type of constraint would have to be enforced by pinning a point at the centre of the sphere.

## Pinning parametric multivectors paths

So far in our construction of multivector pinning constraints we have assumed that the objects we are pinning are static in both the world and body frame. When working with constrained dynamics in the real world we often want to pin parts of our rigid body to moving things in the real world, such as a manipulator attached to the moving end-point of a robot, or a flywheel fixed in a moving vehicle. Consider once again equation (6.19):

In the previous section we enforced static multivector constraints by setting to zero, rearranging to isolate the terms and solving the resultant linear equation for . Now we will relax the static constraint and consider the cases when are known time varying multivector functions, ie when .

First note we can still rearrange to separate terms in :

In fact, if we continue as in our previous analysis by breaking up as a function of and extracting :

and so:

If we substitute we have eventually got to a position where:

Some function of

Again this is a linear function in and so solvable as long as it is of sufficient rank.

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.

## Pinning linear functions of parametric multivector paths

In the previous two sections we dealt directly with transformations that pin static multivectors or time varying multivector paths directly to each other in space. In many practical situations what we would really like to pin is a linear function of one multivector to another. For example we could pin the outer product of a point in the body frame and a plane in the world frame to 0, effectively forcing them to be coincident without specifying anything about their relative orientation (unlike in the transformed plane invariant case). Mathematically we can express our linear function constraint as and time derivatives as:

Once again we can rearrange:

leading to an equation of the form:

Again this is linear and solvable as before. Figure 6.5 shows the simulation with the Clifford Python library Hadfield et al. [sent] of two cases in which the linear function is the outer product with one end of a physical pendulum.

## Mapping Screw Velocity to Lie Algebra Velocity

Throughout this chapter so far we have represented the derivative of the state of the body on the motor manifold as . 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 [2017]. We can, however, nicely side step the problem by directly mapping the velocity screw to a velocity in a suitable se(3) Lie algebra that generates the rotor . 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:

We will refer to this function as the kinematic equation' for the given Lie algebra to Lie group mapping and will have a look at its form for a few choices of .

### Exponential Mapping and the Bortz Equation

One commonly used mapping for SE(3) is the exponential mapping:

Conveniently, the kinematic equation for the exponential mapping of the motor bivectors has already been derived in Candy [2012] and a screw Lie algebra version in Selig [2004] and appears again (with a corrected typo) at the end of Section 4.5 in Selig [2005]. The result in Selig [2004,2005] is derived via idempotents and nilpotents of the adjoint matrix representation of the se(3) Lie algebra but it is readily translatable into our own notation:

where .

The exponential se(3) kinematic equation is also the subject of Section 5.3 of Liam Candy's PhD thesis Candy [2012], and while the set up of the problem is certainly correct we were unable to make their equation 5.41 work in our implementations. Suspecting simply a typo somewhere in their derivations of the derivatives we can calculate the derivatives ourselves and check the final result. We start with the following setup, following mostly along the lines of Candy [2012]. Our objective is to calculate as a function of and or . First, we note that it is possible to write a motor bivector in the form:

 (86)

where is a rotation bivector and is a 3DGA vector. We then follow Candy [2012] choosing to define two quantaties:

where:

We can then write an expression for as the time derivative of (6.25):

This is made up of the following components:

where:

If we restrict ourselves to so(3) these equations will produce the same answer as that of the Bortz equation Bortz [1971] familiar to practitioners from the field of strapdown inertial navigation.

### Cayley Mapping

An alternative mapping to the exponential that is simple and potentially useful is the Cayley mapping Hestenes and Fasse [2002]; Tingelstad and Egeland [2018]. For small rotations and translations the Cayley mapping approximates the exponential however it diverges somewhat as we move further from the origin. Matrix versions of the Cayley map are well known and the kinematic equations for the matrix versions of this map have been studied before in the aeronautics literature Sinclair [2005]. We have failed to find any previous attempts at the GA version of the kinematic equation however we can reuse much of the logic of the matrix derivation, simply substituting transposes for tildes. Start with the expression for the mapping:

 (87)

Take the time derivative:

Rearrange:

Now we will seek a reformulation of :

This allows us to write:

and so we are left with:

 (88)

### Outer Exponential Mapping

The final mapping that we will consider is the so called outer exponential' mapping as presented in Tingelstad and Egeland [2018]. This mapping is defined as taking a Taylor series of the exponential but replacing geometric products with wedge products, for an algebra with maximum grade 5 or below this can be written as:

 (89)

Again we were unable to find a GA kinematic equation for this mapping in the existing literature and so present our own as follows. First we take a time derivative of the outer exponential:

Then, given that in (6.3) we have defined:

we can therefore write:

Noticing that:

and rewriting our equation for as:

gives us:

We can now equate grade 0 elements:

and so:

We can also equate grade 2 elements:

This, after some rearrangement, leaves us with the form of the kinematic equation for the outer exponential:

We can write this neatly as:

 (90)

## Conclusions

In this chapter we have looked at forces, moments, free and constrained dynamics in both CGA and PGA. As well as considering how to apply virtual power as a constraint mechanism in our GA formulations we have constructed a novel technique for constrained dynamics in GA via the concept of multivector pinning. While in this chapter we have only considered two algebras, CGA and PGA, the techniques are expected to work across the board for algebras with easily representable line elements and motor bivectors. Using other higher dimensional algebras such as Cl(4,4) Du et al. [2017], Cl(8,2) Easter and Hitzer [2017] or even Cl(9,6) Breuils et al. [2018] with this technique in the future should allow for easy configurations of exotic constraints such as pinning dynamic objects to the surface of quadrics.

# The Kinematics of Multi-body Systems in Geometric Algebra

Wealth, beauty and fame are transient. When those are gone, little is left except the need to be useful.Gene Tierney

## Abstract

Screw Theory is a framework for analysing articulated mechanisms and performing statics and dynamics calculations that has found much success in the kinematic analysis of mechanisms. In this chapter we consider the embedding of Screw Theory into another extremely powerful framework for robotics, namely Geometric Algebra (GA). We start by rederiving well known results for the accumulation of twists along kinematic chains within our GA framework before turning our attention to the analysis of kinematic pairs. We derive an elegant representation of kinematic pairs via bilinear functions of basic geometric primitives and use this to describe the most common types of robotic joints. We then address multi-body systems and, using the Delta robot as a case study, we compare the screw theoretic approach to a direct differentiation method for extracting the Jacobians of the system.

## Introduction

Modern manufacturing is increasingly mechanised and automated. The drive for automation has led to a need for advanced modelling capabilities and many modern analysis frameworks have been developed as a result. Perhaps the most successful of these frameworks for the analysis for 3D mechanisms is known as Screw Theory. Screw Theory is, unsurprisingly perhaps, concerned with the study of screws'. In its modern form Screw Theory was first described by Sir Robert Stawell Ball in his Treatise on the Theory of Screws' Ball [1900], however the mathematical roots that underlie this remarkable field come from Projective Geometry Pottmann and Wallner [2001] and the study of Lie groups and Lie algebras Selig [2004]. More recent proponents of Screw Theory include Hunt Hunt [1991], Selig Selig [2005], Davidson Davidson and Hunt [2004], Martins Selig and Martins [2014]; Tischler et al. [2000], Featherstone Featherstone [2008], Gallardo Alvarado Gallardo-Alvarado [2016], Pottmann Pottmann and Wallner [2001], Minguzzi Minguzzi [2013], Lipkin Lipkin [2005] and indeed many others.

Alongside the development of Screw Theory, and indeed building on some of the same fundamental mathematics, we have seen the rise of the Clifford/Geometric Algebra (GA) as a robotics modelling framework Aristidou [2010]; Fu et al. [2013]; Hildenbrand et al. [2008,2019]; Kim et al. [2015]; Kleppe and Egeland [2016]; Tichý [2020]; Zamora and Bayro-Corrochano [2004]. With modern computing capabilities, the high level description of geometry that these algebras afford allows researchers elegant, concise and coordinate-free descriptions of physically intricate mechanisms and constraints. Many of the modern applications of GA are in the analysis of conformal Dorst and Valkenburg [2011] and Euclidean motions Gunn [2011b], however there have been only a few attempts to properly embed the tools of modern Screw Theory into GA Hestenes and Fasse [2002]; Tingelstad and Egeland [2018]. This chapter, along with the previous one, is an attempt to lay out the overlap between the two fields with language and ideas familiar to both Screw Theory and GA practitioners.

## Geometric Algebra

In this chapter we will, by default, work with Conformal Geometric Alegbra (CGA). CGA adds two more basis vectors, 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. [1985]. The mapping from a 3D vector, , to its corresponding CGA vector, , is given by:

 (91)

which is often referred to as: up. We can invert this mapping quite easily so long as we remember to normalise the CGA point such that :

where . There are many excellent expositions on CGA in the literature and so we will refrain from a lengthy introduction of all the features of the algebra in this chapter, instead simply preferring to remind the reader of immediately relevant facts about the framework as we go along. If the reader is looking for a more thorough coverage of CGA we would recommend they turn to the excellent book Geometric Algebra for Computer Science by Dorst, Fontijne and Mann Dorst et al. [2007].

One of the important sections of the CGA framework that this chapter will deal with is the set of bivectors that form the generators of rotors that perform Euclidean motion, we will refer to these bivectors as the motor bivectors. The motor bivector basis set contains 6 elements, reflecting the 6 degrees of freedom present in rigid body motion. A common choice for the motor bivector basis in CGA is the following ordered set:

where represents the 3D pseudo-scalar . This set would have a corresponding reciprocal frame as follows:

such that if and if . The full set of motor bivectors is then a linear combination of this motor bivector basis. Introducing 6 scalar parameters we can write a general motor bivector as:

An alternative and increasingly popular GA framework to work in is the Plane-based or Projective Geometric Algebra (PGA) Dorst [2020]; Gunn [2011b]. This algebra has signature Cl(3,0,1), meaning it has 3 basis vectors that square to , 0 basis vectors that square to and one null basis vector that squares to 0. It is a subalgebra of CGA that contains only the flat' elements and the Euclidean rotors Hrdina et al. [2021]; Lasenby [2011]. Restricting the algebra in this way is useful for certain applications that do not require the round elements of CGA and especially as the reduced dimensionality can produce significant speedups for some numerical packages. PGA also contains the motor bivectors with its null basis vector, taking the place of . Due to its degenerate signature PGA cannot use reciprocal frames in the same way as CGA but this restriction can be neatly sidestepped via pseudo-reciprocal frames as described in Gunn [2020] and explicitly worked through for the motor bivectors in Chapter 6. While PGA is not a direct focus of this specific chapter a sharp eyed reader will note that almost all of the formulae in the Screw Theory framework developed here will work with no, or only minor, modifications in PGA.

Whichever GA you choose, the motor bivectors, when exponentiated, produce rotors that can implement rigid body motions. In general it is possible to split 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.

## Twists in Kinematic Chains

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 :

 (92)

Using to represent the standard geometric algebra reversion operator we can write:

 (93)

Now imagine that these bodies form part of a jointed kinematic chain with a sequence of limbs. The position and orientation of each limb relative to the world origin is defined by the rotor , which can in turn be written .

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:

 (94)

We can decompose this rotor into and :

 (95)

We can also explicitly take the time derivative of , which simplifies as due to both body frames being fixed relative to one another:

 (96)

This implies:

 (97)

If we right multiply by we can see that is the bivector velocity of both and :

 (98)

This should not surprise us as both body frames are fixed to one another.

Taking the time derivative of equation (7.2) leads to a formula relating the relative velocity screws of the body in the chain:

 (99)

 (100)

 (101)

 (102)

This is a potentially convenient result but there is a particular case that proves to be of special interest. Our results so far are not a function of the rotors at all, just of . We could choose any decomposition of into and at a given time point and our equations will still be valid. The special case we will concentrate on is to choose to instantaneously align the fixed frames with the world frame, in other words we choose at this instant: , and therefore . Under these conditions equation (7.12) simplifies to the following:

 (103)

This equation is particularly convenient if we are able to measure the relative velocity screw of one limb with respect to another as we can simply accumulate the relative limb velocities along the chain to arrive at the final global limb velocity screw with respect to the world frame. This is well known in traditional Screw Theory and forms the basis of many practical techniques for analysing robots.

A direct result of the definition of is that the time derivative of a geometric object is given by the commutator product of the object with :

 (104)

This is a particularly helpful result as it allows us to define geometry associated with our kinematic chains and calculate how they evolve through time as the kinematic chain moves.

## Geometrically Constrained Kinematic Pairs

Consider a pair of adjacent limbs, and , in a kinematic chain. These two limbs are connected by some form of joint and together they are known as a kinematic pair Reuleaux [1876]. Exactly how they are constrained to move relative to one another is defined by the type of joint. Let us first consider a specific type of joint, one defined by shared geometry.

In a shared geometry joint there exists a piece of geometry, which we will label , that is fixed relative to both frames in the kinematic pair. Practically what this means is that when acted on by the velocity of or this object must have the same , or more formally:

 (105)

Due to the linearity of the commutator product we can write this as:

 (106)

This is an extremely useful result. The quantity is the relative velocity screw of one limb relative to the other and the restriction of its commutator with X being zero means that the shared geometry of the joint must be invariant to the effect of the relative velocity.

We can represent many types of joint as shared geometry joints, or combinations of shared geometry joints. In practice however we often want to represent slightly more complex joints compactly. In many commonly used GA frameworks we can represent more advanced types of joint with some form of bilinear mapping between two relevant objects, one in each frame, that we know remains constant for the joint geometry. More formally, for an object in frame and an object in frame :

 (107)

where is a multivector that is constant for the joint, and in many cases is simply 0. Taking time derivatives once again gives us a linear constraint on and :

 (108) (109)

So far at no point have we specified what objects and are, nor the form of the bilinear mapping (7.17). In practice, common forms of include the outer product:

the inner product:

and the meet taken with respect to a specific fixed subspace:

## The Geometry of Real Joints

It is all well and good to state some neat theoretical results but it is more useful to robotics practitioners to outline how to describe actual joints in our framework. A good place to start is with the most common types of kinematic pair.

### Spherical Joint

A spherical joint, often known as a ball joint, is one in which two bodies are free to move about a single common fixed point. Often in practice these joints are implemented with one body designed to have a spherical cavity (known as the socket) and another with a ball on the end which is held captive within the socket. To implement a spherical joint within our shared geometry constraint we can use a sphere shared between both bodies in the kinematic pair. The sphere can be of any radius and in many practical applications (and in some GAs which do not naturally embed non-zero radius spheres as objects, such as PGA) it makes sense to use a sphere of zero radius, or, in other words, a point.

### Cylindrical Joint

A cylindrical joint is a mechanism that allows two bodies to translate parallel to a shared axis and rotate about the same axis simultaneously. In practice this joint appears often as a cylindrical cuff free to slip over and around a rod. To implement the cylindrical joint within our framework we can use the geometry of a shared line. A line is invariant to rotation about its axis and translation along its axis. In CGA we can use a line either in its trivector form as the direct wedge product of two conformal points and , 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.

### Planar Joint

A planar joint is one that allows two bodies to slide over one another in a specific plane of motion and rotate about the normal to that plane. An example of this in the real world would be a wheeled trolley that is free to move in any direction on a flat surface. To implement the planar joint within our framework we can use a shared plane between the kinematic pairs. Again in CGA this can be done with either a 4-vector plane, , 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.

### Revolute Joint

A revolute joint, also known as a pin joint or hinge joint, is one that allows rotation about an axis but, unlike the cylindrical joint, does not allow translation parallel to the axis. The revolute joint is one of the most common mechanical joints in use in the wild and so having a convenient mechanism to handle it is important. To implement a revolute joint within our framework with CGA we can use a shared circle. A circle is invariant only to rotation about the axis passing through its centre and normal to the plane in which it lies. The radius of this circle is again arbitrary and choices of zero radius, or radius appropriate to the mechanical proportions of the robot at hand, are likely reasonable. In CGA this circle comes in the form of a trivector 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.

### Prismatic Joint

A prismatic joint is one in which two bodies are constrained to move relative to one another parallel to a fixed line without any rotation. The prismatic joint is the first of the classic robotic joints that we cannot represent with a single geometric primitive in 3D conformal GA. Instead we are forced to use some form of compound constraint. There are several options here, two parallel line constraints, a line and a plane constraint, a line and a direction bivector constraint and two planar constraints are all constructs that would work.

### Universal Joint

A universal joint is a mechanism which consists of two orthogonal revolute joints with axes that are incident at a fixed point. It is often used to transmit rotational motion between two shafts whose axes are not parallel but are incident.

We can represent a universal joint in our framework as a joint consisting of two equal radius circles. The planes of these circles must lie orthogonal to each other and both circles must have the same centre point. If we choose two circles, and , the orthogonality constraint can be written as and the same centre point constraint for circles of the same radius can be written Hadfield and Lasenby [2019]. We can combine these constraints neatly with the anti-commutator product

When operating on two circles, the anti-commutator product produces grade 0 and grade 4 elements only. Using this notation and taking derivatives our kinematic constraint can therefore be described as:

 (110)

Of course we could have chosen to model this joint with two revolute joints and an additional body to represent the internals of the joint but this would introduce additional and unnecessary variables into our problem definition.

## The Kinematic Constraint Matrix and the Jacobian Matrix

So far we have only looked at individual kinematic pairs, we will now address full multi-body systems. For an articulated robot with limbs we write the velocity state as:

 (111)

For a system with constraints, there exists an matrix with linear functions as elements that embodies the combination of all the linear constraints on the velocity state. acts on the velocity state giving a result of M zero valued multivectors:

 (112)

We will name the kinematic constraint matrix and readers familiar with the Screw Theory literature will recognise it as a cousin of Davies' method Davies [1981]. We can choose some ordered dimensional basis in which to represent the multivectors in which case becomes an dimensional vector, becomes an dimensional matrix and the zero vector is dimensional. To find the various coefficients of we can simply set each of the coefficients of to 1 in turn and all the rest to zero and calculate the output of each of our constraints for that input, calculating its output representation in the chosen multivector basis.

Consider a robot with velocity state , now specify some of the limbs as inputs, with known and some of the limbs as outputs with unknown . We can write this as follows:

 (113)

where the subscripts and refer to known and unknown respectively. By linearity, applying the constraint matrix produces:

 (114)

This is in the standard form and can be solved with the pseudo-inverse in our chosen basis to produce a valid set of unknown limb velocities given the set of known velocities.

Designating a linear map 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:

 (115)

where is the pseudo-inverse of .

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.

## Case Study: The Delta Robot

The Delta robot Clavel [1990,1991] was invented in 1985 by Raymond Clavel at EPFL after being inspired by a visit to a chocolate packing factory Pessina [2012]. It has since become a particularly popular robot in industrial settings due to its good precision coupled with high speed and acceleration.

The Delta robot is a specific type of robot known as a parallel manipulator. Parallel manipulators, also known as parallel robots, are a class of robots that feature end-effectors driven by multiple underactuated parallel kinematic chains Gallardo-Alvarado [2016]; Merlet [2006]. Typically a parallel robot is designed such that all actuators remain fixed to the support structure of the robot thereby minimising the mass of the moving parts of the robot and enabling very fast accelerations. Indeed this goal of high speed/fast acceleration has been the primary driving force in the development of parallel robots for industry and today architectures such as the Delta robot are widespread in many high precision, high throughput manufacturing applications. Parallel robots, while practically very useful, are often significantly more difficult to analyse than their serial cousins due to the end-point position being a function of the configuration of multiple kinematic chains.

Here we will do a case study of the Delta robot, analysing it with Geometric Algebra, calculating Jacobians with our Screw Theory based framework, and then finally comparing our screw theory setup with a direct differentiation approach to get the Jacobians.

### Geometry of a Delta Robot

Since its inception, there have been many variants of the Delta robot Pierrot et al. [1991]. In this chapter we will assume the simple robot described in this section and illustrated in Figure 7.2. The static part of the robot is a base plate to which three motors are rigidly attached, we will assume a space in which the origin is at the centre of this plate. Each motor shaft is rigidly attached to an upper arm' of length ; we will number each upper arm . The connection point of the motor and upper arm will be labelled . The arm can only rotate in plane about the motor axis as the motor shaft and upper arm are rigidly connected. We will refer to the other end of this upper arm as the elbow point' and will label it . At the elbow point each arm is rigidly attached to a central point of a horizontal rod we will refer to as the elbow rod'. At each end of the elbow rod a ball joint connects to a forearm' piece. The two forearm pieces for each arm are the same length and, at the other end from the elbow rod, are connected to a rigid plate that we will refer to as the 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.

### Calculating the Robot Pose

The first task in analysing the Delta robot is to calculate its pose for a given set of inputs. Typically we are attempting to map from the end-effector to the actuators (inverse kinematics) or from the actuators to the end-effector (forward kinematics). Along the way we are interested in the positions of the limbs and joints that define the physical structure of the robot. To solve these pose calculation problems for the Delta robot initially we will use a simplified geometry for each leg.

#### Inverse Kinematics

The inverse kinematic problem for the Delta robot is summarised as follows: To what angle relative to the base should we move the upper arms given we want the centre of the end-effector plate to be in a specific position in 3D space?

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 elbow point is also simultaneously constrained to lie on a circle of radius centred at the motor shaft to upper-arm joint, . We can represent this circle in its dual form in CGA, where is the intersection of a sphere of radius centred at the position , with dual form , and the plane through the origin, and which has dual form . Here is the vertical unit vector and for the most common orientation of a delta robot, points vertically downward from the base plate. In CGA we calculate the intersection of objects via the meet' operator, as both operands are in their dual form however, here we simply need an outer product:

So long as is within the reachable volume of the robot there are two possible solutions for this pair of constraints. These two solutions lie at the intersection points of the sphere and circle and the meet' operation of CGA provides us with a direct means to calculate these intersection points. The sphere and circle are in the dual form (1 and 2-vectors respectively), and so the point-pair bivector resulting from their meet is calculated as simply their outer product followed by multiplication with the 5D pseudo-scalar, :

The desired individual solution can be extracted from this point-pair object by relying on the oriented nature of the algebra and relying on a projection operator Lasenby, A.N., Lasenby, J. Wareham, R.J. [2004]:

We can then convert from the CGA to the 3D vector point:

and so, with a little trigonometry we can extract the motor angles:

where is the radius of the base-plate. Figure 7.3 illustrates the geometry of the inverse kinematic problem graphically. As mentioned we have relied on the oriented nature of the algebra to extract the solution of interest from the point-pair. This solution has the elbow position being as far from the axis as possible and is normally the only feasible position of the elbow in a real Delta robot as the other typically causes self intersection. If the other solution is desired the same projection formulae can be used, simply substituting for .

#### Forward Kinematics

The forward kinematic problem is, in some sense, the opposite of the inverse kinematic one. Our goal here is to calculate the 3D vector position of the end-effector plate 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:

For each arm we will now define a pseudo-elbow point, which is offset horizontally from the true elbow point by the radius of the end effector plate and in the direction of the origin.

The equivalent CGA point is then:

Given the geometry of the robot, these pseudo-elbow points all lie a distance equal to the length of the robot's forearms, , from the centre of the end-point plate . Geometrically these fixed distance constraints manifest themselves as spheres, which we will label , on which the centre of the end-point plate can lie:

Each arm contributes one constraint sphere and the intersection of the three spheres produces a point-pair, , that represents the two possible configurations of the end-plate:

where the notation implies an outer product of all elements following it.

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:

Figure 7.4 illustrates the geometry of the forward kinematic problem.

### Full Geometry and Kinematic Constraint Matrix of the Delta Robot

Knowing static kinematic solutions is useful but to do more advanced analysis of the Delta robot mechanism we need to look at velocities. Given we have, up to this point, simplified the pose of the robot in both the forward and inverse cases, we will turn our attention to analysis of the full geometry and possible movements of the limbs. First we will do our analysis with our Screw Theory based framework and then we will compare it with a direct differentiation approach.

Looking at the Delta robot there are two types of joint present. The motor shaft to motor body connection is a revolute joint, the arm to arm connections are spherical joints, and the arm to end-effector platform are also spherical joints. The revolute joints of each motor connection can neatly be represented by a piece of shared geometry, namely a circle with normal parallel to the motor axis. Any radius of circle could be chosen, however as the limbs attached to the motor have a length of 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:

The remaining joints are all spherical/ball joints and can therefore be conveniently represented by shared spheres. Again, the radius of these spheres is irrelevant. In this case we will choose a radius of zero, which makes the shared spheres into shared points. We will label each of the limbs of the Delta robot according to the diagram in Figure 7.6. According to this labelling scheme we then have the following mapping from simplified to full robot geometry:

 up up up up up up up up up

where:

due to the geometry of the end plate as shown in Figure 7.5.

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 :

 (116)

Putting this together with the shared point ball joints the full set of kinematic constraints for the Delta robot are as follows:

The forward kinematic problem has known input velocities and output velocity . The inverse kinematic problem has known input velocity and output velocities .

### From Constraint Matrix to Jacobian Matrices

In order to convert from a constraint matrix to a Jacobian matrix we will need to define the linear maps mentioned in Section 7.6, 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 :

With these dual lines we can map between the motor angular speeds and the twists .

For the forward kinematic problem we therefore can construct a map as follows:

and as for the inverse kinematic case we can construct from negative the transpose of , .

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:

This is a general method to extract the linear velocity of a point from the CGA point and its velocity screw. In the case of the Delta robot of course we know that the end-plate is constrained to move only translationally and so will be a purely translational bivector of the form . We can therefore extract directly from :

To turn into individual components we can simply dot it with each basis vector in turn etc.

With this in mind we can now construct the known input map for the inverse problem :

and the unknown output map for the forward problem :

### Calculating the Jacobian with Direct Differentiation

#### The Inverse Jacobian

For our direct differentiation method we will start with our simplified inverse kinematic solution and simply differentiate the expressions directly.

First we will write the 3D end-point plate position as a linear combination of basis vectors with coefficients denoted :

Our goal is to calculate the partial derivative of each motor angle with respect to each of these coefficients. Taking partial derivatives of with respect to one of the coefficients trivially gives:

For now we do not need to worry about which parameter we are taking derivatives with respect to, so we will leave the derivative of the end-point written as . Using this notation, our ultimate goal in this section is to find an equation for the partial derivative of a given motor angle with respect to , . To find we select a specific robot arm and work back through its joints from the end-point.

The first joint position of interest is , we saw in section 7.7.2 that:

Taking partial derivatives gives:

The 3D point can then be represented as the CGA point :

The derivative of this CGA point is then easily found:

We then form the dual constraint sphere:

which, as the radius is fixed, has partial derivative:

As we saw in the previous section, the intersection of this dual constraint sphere and the dual circle centred on the motor shaft produces a point-pair that represents the two possible elbow positions for that arm:

The outer product and dual operations are both linear, which means that taking derivatives is particularly easy here:

Of course the elbow can only actually be in one position which we can extract via a projection operation:

We can then convert from the CGA to the 3D vector point:

and use this to form the derivative of the motor angles with respect to :

This finally gives us an expression for the derivative of the motor angle with respect to the of the endpoint. Typically in engineering scenarios we would construct a matrix of the partial derivatives with respect to , known as the Jacobian matrix:

This matrix can then be used to convert an end-point velocity vector to a set of motor velocities:

As it is the Jacobian matrix for the inverse kinematic problem, this matrix is specifically labelled the inverse Jacobian matrix.

#### The Forward Jacobian

Many problems in robotics require us to take derivatives of the forward kinematic equations. Specifically, we need to know the end-point plate velocity as a function of the motor speeds.

Our forward kinematic solution begins with calculating the position of the elbow point for a given arm :

With the elbow point we can then calculate the pseudo-elbow point:

We then convert the pseudo-elbow to a CGA point:

The forearm length dual constraint sphere can then be constructed about the pseudo-elbow point

The intersection of all three constraint spheres, one from each arm, produces the point pair on which the solution lies:

We can take derivatives of this point-pair with respect to each of the motor angles:

We can re-write these derivatives as follows:

 where (117)

Practically, when we take partial derivatives with respect to one at a time we are effectively freezing two of the motors in position and moving the third. Geometrically, this process forces the end-point plate to move along a circle formed by the intersection of the two constraint spheres centred at the pseudo-elbow points of the frozen motors.

Figure 7.7 shows the geometric significance of Equation 7.27. To get the end-point plate position we again extract one end of the point-pair :

Finally we convert our end-point back to a 3D point:

We can write the end-point plate position as:

With we are therefore in a position to build the forward Jacobian matrix:

The inverse Jacobian matrix and the forward Jacobian matrix are, as the names suggest, inverse to each other.

### Comparing Direct Differentiation to Screw Theory

Both the direct differentiation method and our Screw Theory inspired method can be shown to give us numerically identical Jacobian matrices, however the screw theoretic approach is significantly easier for a practitioner to use as it does not require explicitly taking derivatives. In this example of the Delta robot we have the unusual luxury of a simple closed form solution for the pose of the robot in both the forward and inverse case and so the direct differentiation method is easy to compute either manually as we have done here or using automatic differentiation. With other robot architectures we may not have this blessing and would have to rely on non-linear optimisation or algebraic geometry methods to find a pose that satisfies the constraints. Table 7.1 shows a brief qualitative comparison of the two methods, as ever in practice the choice of which is best to use will be down to the problem itself and what tools and computational resources are available.

Table 7.1: A qualitative comparison of the direct differentiation method vs the Screw Theory based technique for analysing the Jacobians of the system.
 Direct differentiation Screw Theory based Does not require explicit derivatives of pose calculation ❌ ✅ Can be implemented directly with automatic differentiation ✅ ❌ Encapsulates information about the entire system ❌ ✅

## Conclusions and Future Work

In this chapter we have embedded screw theoretic concepts within Geometric Algebra and used this embedding to analyse kinematic pairs and full multi-body systems. The combination of Screw Theory and GA is a particularly potent mix for robotics, allowing clean expressions of geometric constraints and neat representations of kinematic limitations. There are many potential avenues for future work in this area, one promising route would be in the use of higher dimensional Clifford Algebras Breuils et al. [2019b] to represent complex contact surfaces in joints, another might be to expand the allowable motions to include expansion and shear, allowing us to form an extended screw theory for soft robotic modelling. More immediately there is the issue of characterising the computational cost of these Screw Theory based methods. A comparison of the compute speed of these techniques with the direct differentiation method is beyond the scope of this chapter but as it is of practical importance it will likely be a focus of follow up work on this topic.

# Conclusions

## Main Contributions

In this body of work we have made the following contributions:

• We link the problem of finding a valid rotor between two geometric primitives to the direct addition of these primitives
• We introduce a technique for extracting a blade from an arbitrary pure grade multivector
• We show how linear combinations of multivectors can be used to make tubular and ribbon surfaces and give formulae and techniques for intersections of these surfaces with lines and the normal to the surface at any point
• We make explicit the mapping between lines in CGA and PGA and the screws of Screw Theory
• We show how GA can provide Screw Theory with a coordinate free screw representation and show the applications of this in rigid body dynamics
• We show how the GA screw formulation allows us to extend Screw Theory multi-body kinematics to include geometric primitives directly as kinematic constraints within joints

## Future Work

Much of this thesis has been focused on the theory and applications of the 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.

### References

Aristidou, A. (2010).
Tracking and Modelling Motion for Biomechanical Analysis.
PhD thesis, University of Cambridge.

Ball, R. S. (1900).
A treatise on the theory of screws.
Cambridge University Press Cambridge.

Bauer, U. and Polthier, K. (2007).
Parametric reconstruction of bent tube surfaces.
In 2007 International Conference on Cyberworlds (CW’07), page 465–474. IEEE.

Bayro-Corrochano, E. and Rivera-Rovelo, J. (2009).
The use of geometric algebra for 3D modeling and registration of medical data.
Journal of Mathematical Imaging and Vision, 34(1):48–60.

Besl, P. J. and McKay, N. D. (1992).
A method for registration of 3D shapes.
IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(2):239–256.

Bezier, P. (1986).
The Mathematical Basis of the UNISURF CAD System.
Butterworth-Heinemann, USA.

Bortz, J. E. (1971).
A new mathematical formulation for strapdown inertial navigation.
IEEE Transactions on Aerospace and Electronic Systems, AES-7(1):61–66.

Bosche, F. (2011).
Plane-based coarse registration of 3D point clouds with 4d models.
28th International Symposium on Automation and Robotics in Construction.

Boyle, M. (2017).
The integration of angular velocity.
Advances in Applied Clifford Algebras, 27(3):2345–2374.

Breuils, S., Nozick, V., and Fuchs, L. (2019a).
Garamon: A geometric algebra library generator.
Advances in Applied Clifford Algebras, 29(3).

Breuils, S., Nozick, V., Fuchs, L., and Sugimoto, A. (2019b).
Transverse approach to geometric algebra models for manipulating quadratic surfaces.
In Gavrilova, M., Chang, J., Thalmann, N. M., Hitzer, E., and Ishikawa, H., editors, Advances in Computer Graphics, pages 523–534, Cham. Springer International Publishing.

Breuils, S., Nozick, V., Sugimoto, A., and Hitzer, E. (2018).
Quadric conformal geometric algebra of R(9,6).
Advances in Applied Clifford Algebras, 28(2).

Cameron, J. I. (2007).
Random Disconnected Applications of Geometric Algebra in Computer Graphics and Computer Vision.
PhD thesis, University of Cambridge.

Candy, L. P. (2012).
Kinematics in conformal geometric algebra with applications in strapdown inertial navigation.
PhD thesis, University of Cambridge.

Catmull, E. and Rom, R. (1974).
A class of local interpolating splines.
In Computer Aided Geometric Design, pages 317–326. Academic Press.

Clavel, R. (1990).
Device for the movement and positioning of an element in space.
United States Patent and Trademark Office, US4976582A.

Clavel, R. (1991).
Conception d’un robot parallèle rapide à 4 degrés de liberté.
PhD thesis, École polytechnique fédérale de Lausanne.

Colapinto, P. (2011).
Spatial Computing with Conformal Geometric Algebra.
PhD thesis, University of California, Santa Barbara.

Colapinto, P. (2017).
Composing surfaces with conformal rotors.
Advances in Applied Clifford Algebras, 27(1):453–474.

Davidson, J. K. and Hunt, K. H. (2004).
Robots and screw theory: Applications of kinematics and statics to robotics.
Journal of Mechanical Design, 126(4):763–764.

Davies, T. H. (1981).
Kirchhoff’s circulation law applied to multi-loop kinematic chains.
Mechanism and Machine Theory, 16(3):171–183.

De Keninck, S. (2019).
Non-parametric realtime rendering of subspace objects in arbitrary geometric algebras.
In Gavrilova, M., Chang, J., Thalmann, N. M., Hitzer, E., and Ishikawa, H., editors, Advances in Computer Graphics, Lecture Notes in Computer Science, page 549–555. Springer International Publishing.

De Keninck, S. (2020).
ganja.js.
Zenodo. https://doi.org/10.5281/ZENODO.3635774.

De Keninck, S. and Dorst, L. (2019).
Geometric algebra levenberg-marquardt.
Advances in Computer Graphics. CGI 2019. Lecture Notes in Computer Science, 11542:511–522.

Deul, C., Burger, M., Hildenbrand, D., and Koch, A. (2009).
Raytracing point clouds using geometric algebra.
GraVisMa proceedings.

Doran, C. (2003).
Circle and sphere blending with conformal geometric algebra.
arXiv:cs/0310017.

Doran, C. and Lasenby, A. (2003).
Geometric Algebra for Physicists.
Cambridge University Press.

Dorst, L. (2020).
A guided tour to the plane-based geometric algebra PGA, version 1.15.
Bivector.net.

Dorst, L., Fontijne, D., and Mann, S. (2007).
Geometric algebra for computer science: an object-oriented approach to geometry.
Morgan Kaufmann series in computer graphics. Elsevier; Morgan Kaufmann.

Dorst, L. and Valkenburg, R. (2011).
Square root and logarithm of rotors in 3D conformal geometric algebra using polar decomposition.
In Guide to Geometric Algebra in Practice, page 81–104. Springer, London.

Du, J., Goldman, R., and Mann, S. (2017).
Modeling 3D geometry in the clifford algebra R(4,4).
Advances in Applied Clifford Algebras, 27(4):3039–3062.

Duda, R., Hart, P., and Stork, D. (2001).
Pattern Classification.
Wiley Interscience, 2nd edition edition.

Easter, R. B. and Hitzer, E. (2017).
Double conformal geometric algebra.
Advances in Applied Clifford Algebras, 27(3):2175–2199.

Eggert, D., Lorusso, A., and Fisher, R. (1997).
Estimating 3D rigid body transformations: a comparison of four major algorithms.
Machine Vision and Applications, 9(5):272–290.

Eide, E. R. and Lasenby, J. (2018).
A novel way of estimating rotors between conformal objects and its applications in computer vision.
AACA: Topical Collection AGACSE 2018, IMECC – UNICAM, Campinas, Brazil.

Featherstone, R. (2008).
Rigid body dynamics algorithms.
Springer.

Featherstone, R. (2010).
A beginner’s guide to 6D vectors (part 1).
IEEE Robotics & Automation Magazine, 17(3):83–94.

Fischler, M. A. and Bolles, R. C. (1981).
Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography.
Commun. ACM, 24(6):381–395.

Fu, Z., Yang, W., and Yang, Z. (2013).
Solution of inverse kinematics for 6R robot manipulators with offset wrist based on geometric algebra.
Journal of Mechanisms and Robotics, 5(3).

Kinematic Analysis of Parallel Manipulators by Algebraic Screw Theory.
Springer International Publishing.

Gelfand, N., Mitra, N. J., Guibas, L. J., and Pottmann, H. (2005).
Robust global registration.
Eurographics Symposiumon Geometry Processing.

Gunn, C. (2011a).
On the homogeneous model of euclidean geometry.
arXiv:1101.4542 [math].
arXiv: 1101.4542.

Gunn, C. (2020).
On-metric alternatives to reciprocal frames. the bivector.net forums. https://discourse.bivector.net/t/non-metric-alternative-to-reciprocal-frame/105/4. Accessed 12th May 2020.

Gunn, C. G. (2011b).
Geometry, Kinematics, and Rigid Body Mechanics in Cayley-Klein Geometries.
PhD thesis, Technical University Berlin.

Gunn, C. G. and De Keninck, S. (2019a).
Geometric algebra and computer graphics.
ACM SIGGRAPH 2019 Courses (SIGGRAPH '19). Association for Computing Machinery, New York, NY, USA, Article 12, 1–140.

Gunn, C. G. and De Keninck, S. (2019b).
Siggraph geometric algebra course video.

Hadfield, H., Achawal, S., and Lasenby, J. (2021).
Exploring novel surface representations via an experimental ray-tracer in CGA.
Advances in Applied Clifford Algebras, 31(16).

Hadfield, H. and Lasenby, J. (2019).
Direct linear interpolation of geometric objects in conformal geometric algebra.
Advances in Applied Clifford Algebras, 29(85).

Hadfield, H., Lasenby, J., Ramage, M., and Doran, C. (2019).
Reform: Rotor estimation from object resampling and matching.
Advances in Applied Clifford Algebras: Topical Collection AGACSE 2019, IMECC – UNICAM, Campinas, Brazil.

Hadfield, H., Wei, L., and Lasenby, J. (2020).
The forward and inverse kinematics of a Delta robot.
In Magnenat-Thalmann, N., Stephanidis, C., Wu, E., Thalmann, D., Sheng, B., Kim, J., Papagiannakis, G., and Gavrilova, M., editors, Advances in Computer Graphics, Lecture Notes in Computer Science, page 447–458. Springer International Publishing.

Hadfield, H. and Wieser, E. (2020).
Robots, ganja & screw theory.

Hadfield, H., Wieser, E., Arsenovic, A., Kern, R., and The Pygae Team (2018-Present).
www.github.com/pygae/clifford.

Hestenes, D. (2001).
Old Wine in New Bottles: A New Algebraic Framework for Computational Geometry, pages 3–17.
Birkhäuser Boston, Boston, MA.

Hestenes, D. (2010).
New Tools for Computational Geometry and Rejuvenation of Screw Theory, pages 3–33.
Springer London, London.

Hestenes, D. and Fasse, E. D. (2002).
Homogeneous Rigid Body Mechanics with Elastic Coupling, pages 197–212.
Birkhäuser Boston, Boston, MA.

Hestenes, D., Sobczyk, G., and S. Marsh, J. (1985).
Clifford algebra to geometric calculus. a unified language for mathematics and physics.
American Journal of Physics, 53:510–511.

Hildenbrand, D. (2007).
Geometric Computing in Computer Graphics and Robotics using Conformal Geometric Algebra.

Hildenbrand, D. and Hitzer, E. (2008).
Analysis of point clouds - using conformal geometric algebra.
In Proceedings of the Third International Conference on Computer Graphics Theory and Applications, page 99–106. SciTePress - Science and and Technology Publications.

Hildenbrand, D., Hrdina, J., Návrat, A., and Vašík, P. (2019).
Local controllability of snake robots based on CRA, theory and practice.
Advances in Applied Clifford Algebras, 30(1):2.

Hildenbrand, D., Zamora, J., and Bayro-Corrochano, E. (2008).
Inverse kinematics computation in computer graphics and robotics using conformal geometric algebra.
Advances in Applied Clifford Algebras, 18(3–4):699–713.

Hitzer, E. and Sangwine, S. J. (2019).
Construction of multivector inverse for clifford algebras over 2m+1-dimensional vector spaces from multivector inverse for clifford algebras over 2m-dimensional vector spaces.
Advances in Applied Clifford Algebras, 29(2):29.

Hitzer, E., Tachibana, K., Buchholz, S., and Yu, I. (2009).
Carrier method for the general evaluation and control of pose, molecular conformation, tracking, and the like.
Advances in Applied Clifford Algebras, 19(2):339–364.

Horn, R. A. and Johnson, C. R. (2012).
Matrix analysis.
Cambridge University Press, 2nd edition.

Hrdina, J., Návrat, A., and Vašík, P. (2018).
Geometric algebra for conics.
Advances in Applied Clifford Algebras, 28:1–21.

Hrdina, J., Návrat, A., and Vašík, P. (2021).
Projective geometric algebra as a subalgebra of conformal geometric algebra.
Advances in Applied Clifford Algebras, 31(18).

Hunt, K. (1991).
Kinematic geometry of mechanisms.
Robotica, 9(1).

K. Davidson, J., H. Hunt, K., and Pennock, G. (2004).
Robots and screw theory: Applications of kinematics and statics to robotics.
Journal of Mechanical Design - J MECH DESIGN, 126.

Kavan, L., Collins, S., Žára, J., and O'Sullivan, C. (2008).
Geometric skinning with approximate dual quaternion blending.
ACM Trans. Graph., 27(4):105:1–105:23.

Kim, J. S., Jeong, J. H., and Park, J. H. (2015).
Inverse kinematics and geometric singularity analysis of a 3-SPS/S redundant motion mechanism using conformal geometric algebra.
Mechanism and Machine Theory, 90:23–36.

Kleppe, A. L. and Egeland, O. (2016).
Inverse kinematics for industrial robots using conformal geometric algebra.
Modeling, Identification and Control: A Norwegian Research Bulletin, 37(1):63–75.

Kleppe, A. L. and Egeland, O. (2018).
A curvature-based descriptor for point cloud alignment using conformal geometric algebra.
Advances in Applied Clifford Algebras, 28(2):50.

Kochanek, D. H. U. and Bartels, R. H. (1984).
Interpolating splines with local tension, continuity, and bias control.
ACM SIGGRAPH Computer Graphics, 18(3):33–41.

Lasenby, A. (2011).
Rigid Body Dynamics in a Constant Curvature Space and the `1D-up' Approach to Conformal Geometric Algebra, pages 371–389.
Springer London, London.

Lasenby, A., Lasenby, R., and Doran, C. (2011).
Rigid Body Dynamics and Conformal Geometric Algebra, page 3–24.
Springer.

Lasenby, J., Hadfield, H., and Lasenby, A. (2018).
Calculating the rotor between conformal objects.
AACA: Topical Collection AGACSE 2018, IMECC – UNICAM, Campinas, Brazil, 29:102.

Lasenby, A.N., Lasenby, J. Wareham, R.J. (2004).
A covariant approach to geometry using geometric algebra,.
Cambridge University Engineering Department, Technical Report, CUED/F-INFENG/TR-483.

Lipkin, H. (2005).
Time derivatives of screws with applications to dynamics and stiffness.
Mechanism and Machine Theory, 40(3):259–273.

Merlet, J.-P. (2006).
Parallel robots.
Solid mechanics and its applications. Kluwer Academic Publishers, 2nd edition.

Minguzzi, E. (2013).
A geometrical introduction to screw theory.
European Journal of Physics, 34(3):613–632.
arXiv: 1201.4497.

Müller, A. (2018).
Screw and lie group theory in multibody dynamics.
Multibody System Dynamics, 42(2):219–248.

Perwass, C. (2009).
Geometric algebra with applications in engineering.
Geometry and computing. Springer.

Pessina, L. (2012).
Reymond clavel, creator of the delta robot, reflects on his career - sti - school of engineering.
https://sti.epfl.ch/reymond-clavel-creator-of-the-delta-robot-reflects-on-his-career/.

Peternell, M. and Pottmann, H. (1997).
Computing rational parametrizations of canal surfaces.
Journal of Symbolic Computation, 23(2–3):255–266.

Petrovic, V., Fallon, J., and Kuester, F. (2007).
Visualizing whole-brain dti tractography with gpu-based tuboids and lod management.
IEEE Transactions on Visualization and Computer Graphics, 13(6):1488–1495.

Pierrot, F., Fournier, A., and Dauchex, P. (1991).
Towards a fully-parallel 6 dof robot for high-speed applications.
In 1991 IEEE International Conference on Robotics and Automation Proceedings, page 1288–1293 vol.2.

Pottmann, H. and Wallner, J. (2001).
Computational Line Geometry.
Mathematics and Visualization. Springer-Verlag.

Reuleaux, F. (1876).
Kinematics of Machinery. Translated and Edited by Kennedy, A. B. W.
Macmillan and Company.

Segal, A., Haehnel, D., and Thrun, S. (2009).
Generalized-icp.
In Robotics: science and systems, volume 2, page 435. Seattle, WA.

Selig, J. M. (2004).
Lie groups and lie algebras in robotics.
In Byrnes, J., editor, Computational Noncommutative Algebra and Applications, pages 101–125, Dordrecht. Springer Netherlands.

Selig, J. M. (2005).
Geometric fundamentals of robotics.
Monographs in computer science. Springer, 2nd edition.

Selig, J. M. and Martins, D. (2014).
On the line geometry of rigid-body inertia.
Acta Mechanica, 225(11):3073–3101.

Shirokov, D. S. (2021).
On computing the determinant, other characteristic polynomial coefficients, and inverse in clifford algebras of arbitrary dimension.
Computational and Applied Mathematics, 40:173.

Sinclair, A. J. (2005).
Generalization of rotational mechanics and application to aerospace systems.
Texas A&M University.

Tichý, R. (2020).
Inverse kinematics for the industrial robot IRB4400 based on conformal geometric algebra.
In Mazal, J., Fagiolini, A., and Vasik, P., editors, Modelling and Simulation for Autonomous Systems, Lecture Notes in Computer Science, page 148–161. Springer International Publishing.

Tingelstad, L. and Egeland, O. (2017).
Motor estimation using heterogeneous sets of objects in conformal geometric algebra.
Advances in Applied Clifford Algebras, 27(3):2035–2049.

Tingelstad, L. and Egeland, O. (2018).
Motor parameterization.
Advances in Applied Clifford Algebras, 28(2).

Tischler, C. R., Downing, D. M., Lucas, S. R., and Martins, D. (2000).
Rigid-body inertia and screw geometry.
Proceedings of a Symposium Commemorating the Legacy, Works, and Life of Sir Robert Stawell Ball Upon the 100th Anniversary of A Treatise on the Theory of Screws, page 14.

Valkenburg, R. and Dorst, L. (2011).
Estimating motors from a variety of geometric data in 3D conformal geometric algebra.
Guide to Geometric Algebra in Practice, pages 25–45.

Wareham, R. and Lasenby, J. (2008).
Mesh vertex pose and position interpolation using geometric algebra.
In Perales, F. J. and Fisher, R. B., editors, Articulated Motion and Deformable Objects, Lecture Notes in Computer Science, page 122–131. Springer Berlin Heidelberg.

Wareham, R. J. and Lasenby, J. (2011).
Generating fractals using geometric algebra.
Advances in Applied Clifford Algebras, 21(3):647–659.

Zamora, J. and Bayro-Corrochano, E. (2004).
Inverse kinematics, fixation and grasping using conformal geometric algebra.
In 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (IEEE Cat. No.04CH37566), volume 4, page 3841–3846.

Applications of Geometric Algebra in Mathematical Engineering

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

##### Footnotes
....
There are a range of different notations in use in the literature to define the basis of CGA: is sometimes referred to in other works as , as , as and an is sometimes defined which is equal to