Direct Linear Interpolation of Geometric Objects in Conformal Geometric Algebra

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 paper 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 [6].
The objects we work with here will be CGA objects unless explicitly stated otherwise. We will use the standard extension of the 3D geometric algebra, where our 5D CGA space is made up of the standard spatial basis vectors {e i } i = 1, 2, 3, plus two additional basis vectors, e andē with signatures, e 2 = 1,ē 2 = −1. Two null vectors can therefore be defined as: n ∞ = e +ē and n 0 = e−ē 2 . The mapping of a 3D vector x to its conformal representation X is given by X = F (x) = 1 2 (x 2 n ∞ + 2x − 2n 0 ). 1

Motivation
In our conformal representation of space the blades of a given grade n represent specific classes of object and lie on a manifold within the overall subspace This article is part of the Topical Collection on Proceedings of AGACSE 2018, IMECC-UNICAMP, Campinas, Brazil, edited by Sebastià Xambó-Descamps and Carlile Lavor.

Linearly Interpolating Conformal Points
The result of linear combinations of conformal points is well known [12]. Consider two arbitrary points in 3d space a and b represented as A and B in our conformal model. Linear interpolation of these points followed by our conformal mapping produces a linear interpolation of our conformal points with an additional term in n ∞ : We can therefore get a useful interpolation of points by taking a direct linear interpolation and simply adding the final α(1 − α)(A · B)n ∞ term to the result. If Y = αA + (1 − α)B, we can recover Y = F (αa + (1 − α)b) via the following formula (assuming Y · n ∞ = −1):

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 [4], or the problem is attacked indirectly via carriers [9] or by forming the rotor between the objects, extracting the corresponding bivector, which is then interpolated [15] and applied to the first object. It is shown in [11] that we can represent the mirror object X m that reflects one object X 1 , into another X 2 , as the left multiplication of the summation of the blades by a scalar + 4-vector factor S: where β and γ are scalars and the 4-vector part of S is proportional to the anticommutator of X 1 and X 2 .
For the previously known cases in which the linear interpolation of higher grade objects gives a blade, such as with circles [4] and point pairs [5] both with common points, the factor S is a scalar and the object X m is simply 'half-way' between the objects. We can extend this notion to the cases where the addition of objects is not a blade by using our object X m , which has been corrected to being a blade, as the half-way object. We can use this idea of the half-way object to recursively subdivide the space between X 1 and X 2 allowing us to create objects that are any fraction of X 1 and X 2 . While this technique allows us to generate interpolant objects from any two objects (of the same type), it is nevertheless clumsy to represent fractional interpolant objects via repeated subdivision. This subdivision technique also provides no obvious way of performing an average of many objects. What we would really like is some way of directly dealing with the linear interpolation

Creating a Blade from a Pure Grade Multivector
Consider the general interpolant, X α = αX 1 + (1 − α)X 2 where X 1 and X 2 are blades of the same grade. We claim that we can project X α into object space in a simple and general way. First we will generalise Eq. (2) to the interpolation case: where β α and γ α are once again scalars.
Since S is of the form (scalar + 4-vector) it is self reverse. Defining S − = S 0 − S 4 we get the result that S − S is a scalar, and can therefore write X α = kS − X α where k is a scalar and k = 1 S − S . To use this decomposition we need to extract S from X α . To do this we can use the methods of [11], or as follows using the square root operator of Dorst and Valkenburg [6].
Let SX α = X α , where X α is a valid object (squaring to ± 1). Now define Σ ∈ Cl 0,4 4,1 , ie. it only contains 0 and 4 grade coefficients and is an element of the conformal algebra. Then, defining [[Σ]] = Σ 2 0 − Σ 2 4 , the square root can be found as: To use this method to find S − we multiply our non-blade object by its own reverse: this is now in a form where we can apply the above square root formula: it now simply remains to isolate X α via multiplication by kS where This result is particularly important as we have identified a way of projecting any pure grade object of the form S − X (with X a blade) back to the blade manifold. An immediate application of this is that we can now deal with arbitrary linear combinations of objects, allowing us to smoothly interpolate as well as to average and cluster geometric primitives. Additionally we can correct numerical errors that result from arithmetic operations to give true blades again. Figure 1 shows examples of interpolating various geometric objects.
As shown in [11] this method holds for all the standard normalised conformal objects of grade 2 or above (point pairs, lines, circles, planes, spheres). The direct interpolation method is potentially more computationally efficient than the bivector interpolation method, and its form indicates that it is covariant, ie, for a rotor transformation given by R, Unfortunately, we cannot apply this form of interpolation to points as we encounter a problem due to the fact that for a conformal point P , PP = 0. However, we saw in Eq. (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 S 4 X 1 and S 4 X 2 to give only objects of grade n where n is the grade of X 1 and X 2 . Table 1 shows the resultant grades from the geometric product of pure grade objects and Table 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.

Point Pairs
We start with point pairs. Previous work [5] has shown that when an end point is shared between point pairs A and B the interpolant point pairs are also blades and their end points trace out the circumference of the circle     Note, since AB 4 = A∧B = 0, the projector S is a scalar. The common circle itself is the 'join' of the two original point pairs and can be computed with the algorithms supplied in Chapter 21 of Dorst, Fontijne and Mann [5]. Figure 2 shows two cases of the interpolation of co-planar point pairs that lie on the same circle.
Turing 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 1 that we produce both 2 and 4-vector grades. The 2-vector part of the geometric product comes from the inner product between the point pairs and the 4-vector. ie. for point pairs A and B, S 4 (αA For the general case of two point pairs not lying in plane ie. A ∧ B = 0, we can show that there is only one object that behaves in this way, the sphere Σ ∝ A ∧ B, as it passes through both end points of both point pairs. This is illustrated in Fig. 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 C of point pairs A and B always has endpoints lying on the surface of the sphere Σ. We can prove this by showing that C ∧ B or C ∧ A also produces the sphere: First consider an interpolant object C and its outer product with one of the original objects, B Now we just need to prove that S 4 ((αA 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 AB 4 ≡ BA 4 we can write this as: and this can further simplified using the fact that: As AB 0 = A · B is a scalar: Vol. 29 (2019) Direct Linear Interpolation of Geometric Objects in CGA Page 9 of 25 85 is a scalar we see that the proof is completed

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 [4,5]. Here, as with the point pairs, we can show that this is true for a broader class of circles: Theorem 3. If circles C 1 and C 2 together define the caps of a common sphere then This can be proved by considering each circle C i as the intersection of a plane P i and a sphere Σ i . Forming this intersection via the dual (where X * = XI 5 and I 5 is the 5D space pseudoscalar), we have: and so if Σ 2 ∝ Σ 1 : 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 C 1 and C 2 together define the caps of a common sphere Σ then C 1 Σ 3 , C 2 Σ 3 = 0. Σ can be found by the following process: First we define: where M 3 is the truncated identity matrix that performs selection of grade 3 elements from a vector of coefficients andC 1 andC 2 are the matrices that perform the left geometric product of C 1 and C 2 respectively with a vector The case for circles of the same radius is visualised in Fig. 4. It is also the case that the interpolant lies on the surface of the common sphere: Theorem 4. If circles C 1 and C 2 together define the caps of a common sphere then all interpolant circles C 3 = αC 1 + (1 − α)C 2 (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 D, an arbitrary point on the common sphere Σ: Figure 5 shows an example of this interpolation. Thus far we have dealt exclusively with circles on a common sphere. In the case in which C 1 and C 2 do not lie on the same sphere we can again look at how the interpolants behave by considering the form of the (scalar + 4vector) that we use to project the interpolant back to the blade manifold. In the case of the geometric product between grade 4 and 3 objects we see from Table 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 [5]: The intersections of the interpolant circles with the sphere S 4 produce a set of point pairs. Intuition would suggest that these point pairs have properties tied to the interpolation of the point pairs generated by the original two circles C 1 and C 2 and indeed we can numerically verify that this is the case: Conjecture 2. If circles C 1 and C 2 together do not lie on a common sphere then the intersection point pair P α formed by the meet of the circle interpolant for a given value of α with the orthogonal sphere C 1 C 2 4 ie. P α ∝ C α ∨ C 1 C 2 4 is the same as the re-projected interpolant Π α of the point pairs formed from the meet of C 1 and C 2 with C 1 C 2 4 . Figure 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 S 4 ∝ L 1 L 2 4 . For lines L 1 L 2 4 ∝ I 5 n ∞ , giving the form of the projection of (X 1 + X 2 ) as: where μ and ν are scalars. While neat, this form I 5 n ∞ 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' [1]. 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 [3,10,13]. 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 vectorm, a point on the screw axis p and a screw pitch h. The pitch represents how far to move in the direction of the screw axis for each complete revolution about the axis.

Bivector Representation of a Line
A line in CGA is represented as a 3-vector, or dually as a bivector: This bivector formulation is equivalent to the Plücker coordinates of the line. In [6] 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 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 T that transforms along the 3D vector t is: If t = hm where h is a scalar, ie. the translation is in the screw axis direction, the rotors formed from the bivectors in Eqs. (9) and (10) commute. It then follows that the rotor formed from the addition of the bivectors in Eqs. (9) and (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, S, whose action on the point at the origin is shown in Fig. 8.
Hestenes and Sobczyk [8] p81 gives an expression for decomposing any bivector into two commuting blades. In the case of our screw bivector these blades represent the dual of the screw axis L * and a translational bivector T in the direction of the screw axis. ie. given a screw bivector S we can decompose it as:

Adding Dual Lines
The addition of dual lines produces a bivector. Visualising the action of this bivector allows us to see that it is in fact also a screw transformation. Consider the addition of two dual lines: we can write this elementwise as We then rearrange to give something proportional to the expression in Eq. (11): If we divide this by |m| we have the general form of a normalised screw S =mI 3 + (p ∧m)I 3 n ∞ + hmn ∞ Gathering like terms, specifically those without an n ∞ component, leads us to the conclusion that our screw axis directionm must simply be proportional to the addition of the directions of the two lines. Using this fixed axis direction we can extract the coefficient h (the pitch) of the translation bivector parallel to the screw axis: Vol. 29 (2019) Direct Linear Interpolation of Geometric Objects in CGA Page 15 of 25 85 Figure 9. The addition of two dual lines, here shown in black, produces a screw, whose action on the point at the origin is again shown in black. The screw can be decomposed into two commuting bivectors, whose actions are shown in red and blue. The bivector whose action is shown in red is the dual of the screw axis line, also shown in red With this coefficient known we now have all the pieces in place for a full decomposition of the dual line addition bivector L * + : Figure 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 S 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: wherem is the 3d vector normal to the plane and d is the perpendicular distance of the plane from the origin. Thus the interpolation of duals of two planes can be written as: which, when we collect like terms, is already in the form of a dual plane P * 3 : 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 [2], 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 Fig. 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 Figs. 11 and 12.
Vol. 29 (2019) Direct Linear Interpolation of Geometric Objects in CGA Page 17 of 25 85 Figure 10. The interpolant of two planes (green to red) always passes through the meet line (black) of the two original planes Figure 11. The interpolant of two parallel planes smoothly moves between the start and end points (green to red) while maintaining the direction of the normal

Spheres
The interpolant of spheres has been studied before in [2] and [4]. As with planes, all interpolants of spheres are valid objects as Σ 1 Σ 2 4 = 0 and have the property of making contact with the meet of the spheres at all points during the interpolation. We can see the form of the interpolant sphere by considering its dual form:  Figure 12. The interpolant of two anti-parallel planes (green to red) must go via infinity due to the sign change.
Care must be taken with the orientation of objects when designing algorithms using these interpolations The dual form of a sphere can be decomposed into the sum of the conformal centre point P and negative half the radius squared times n ∞ : the interpolation of the dual of two spheres is therefore For two concentric spheres ie. P 2 = P 1 we can therefore see that the interpolation between them will remain centred in the same place and will simply have a radius ρ 3 which varies as ρ 2 3 = αρ 2 1 + (1 − α)ρ 2 2 . As we have seen previously, the interpolation of two conformal points P 1 and P 2 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 p 1 and p 2 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 ρ 3 varies as For fixed values of ρ 1 and ρ 2 this implies ρ 2 3 varies as −(p 1 − p 2 ) 2 and so the further apart the two spheres are the smaller the radius of the interpolant. To find turning points we differentiate with respect to α 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 Σ 1 and Σ 2 are just touching we have the condition 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 [2] 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 Figure 13. Interpolation through control objects. Top: circles. Bottom: point pairs. Interpolation type: a, d linear, b, e quadratic, c, f cubic objects directly suggests immediate applications in clustering of objects extracted from real data, interpolation to produce surfaces and other areas for problems 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 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  Figure 14. A 3D line model before (a) and after (b) recursive scene simplification 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 This leads to a simplified model that retains the core features of the original model. For comparison of objects X i and X j we use the cost function C ij for a rotor R ij as defined in [7]: where R = R · e, and gives the component of R having n ∞ as a factor and R ij is the rotor that takes X i to X j as described in [11]. An example of this algorithm working on simulated lines is shown in Fig. 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.

k-Means Clustering of Conformal Objects
One of the most fundamental and simple clustering algorithms is known as k-means clustering [14]. Consider a 3d scene composed of k geometric objects of a given grade. We have multiple noisy observations for each object and so would like to fit k centroids to these clusters to represent the "true" objects in the world. The steps for implementing this clustering are given below: 1. Randomly assign k 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 Eq. (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 15 and 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, L 1 and L 2 . We wish to find the point P that lies closest to both in a least squares sense. First we will construct two orthogonal intermediary lines L + = S + (L 1 + L 2 ) and L − = S − (L 1 − L 2 ) where S(X) represents the projection of a 3-vector X back onto the line manifold. L + and L − both lie half way between the two original skew lines but intersect at right angles . The intersection of these lines is the point P that lies half way between the original lines. To extract this point of intersection we can follow the formula given in [12]:

Conclusions
This paper 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.