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

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 paper 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
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 [2,23,24]. 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 [13,16,22]. Given its representational power for curved surfaces and simple encoding of complicated operations, CGA appears to hold great promise in the field of Computer Graphics. Indeed several raytracers/path-tracers/sphere-marchers using CGA have been implemented in the past [5,6,10,11,13,19,25]. More recently the design of more intricate surfaces has been investigated with rotors [8] and direct-interpolation of geometric primitives [15]. In this paper 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 1 shows an example of output from the CGA ray-tracer we describe in this paper.

Conformal Geometric Algebra, CGA
The ray-tracer used in this paper is constructed using CGA and all algebraic expressions given will be in terms of elements of this algebra. CGA adds two more basis vectors, e andē, to the original basis vectors of 3D Euclidean space, giving a complete basis for the 5D space with the following signature: e 2 1 = e 2 2 = e 2 3 = e 2 = 1 andē 2 = −1. These extra basis vectors are used to define two null vectors: n ∞ = e +ē ≡ n and n 0 =ē −e 2 ≡ −n 2 -note that the (n,n) notation was that originally used when Hestenes first introduced this model in [17]. The mapping from a 3D vector, x, to its corresponding CGA vector, X, is given by: All vectors formed from such a mapping are null. CGA is chosen for the construction of the ray-tracer since we seek neat expressions for describing intersections, reflections and lighting models, made possible in CGA since Figure 2. The camera is defined by a focal length, a transformation from the origin, and bounds on the image plane rays and scene objects are both elements of the algebra. More background on CGA can be found in [13,16,22].

Camera Model and Ray Casting
A pinhole camera model is used with the geometry shown in Fig. 2. It is defined by a rotor R MV (where MV indicates model view) incorporating rotation and translation that takes the camera from the origin to its pose in space, a focal length f and two bounds x max and y max on the size of the image plane.
We take (i, j) = (0, 0) to be at the bottom left hand corner of the image. For an image of width w and height h, the world coordinates of the point P ij at the centre of pixel (i, j) are given by: We then generate the ray from the camera centre, L ij , that passes through P ij , via the expression where X 0 is the origin transformed by the model-view rotor R MV to the position of the camera.

Ray Geometries for Basic Objects
Initially we will start with some basic objects representable as blades in CGA. The ray-tracer will thus initially concentrate on rendering planes, spheres and circles/discs, an example of which is shown in Figs. 1 and 3.

Ray-Object Intersections
In order to compute intersections between blades, the meet operator (∨) is used. We will, for the proposes of this paper, always take the meet with respect to the full 5D space rather than to the join of the blades. Thus, if X is an r-grade blade, Y is an s-grade blade, and the number of basis vectors in the algebra is n, then : where Z m indicates the m-grade component of the multivector Z, and I 5 represents the 5D pseudoscalar of the algebra [22].

Planes.
A plane is a 4-blade and a ray is a 3-blade so the meet gives a 2-blade. If the meet itself is 0, the line lies in the plane. If the meet squared is 0, there is no finite intersection. Otherwise, the intersection point, X, of a line L with a plane Φ, satisfies the following: L ∨ Φ = λX ∧ n ∞ [22], where λ is a scalar. When extracting the 3D intersection point x, we need to account for the sign and magnitude of the line in our extraction, we can do this via the constant of proportionality λ: Therefore, x can be extracted from the e i e and e iē coefficients (for i ∈ {1, 2, 3}) by dividing by λ, the eē coefficient.

Spheres.
Spheres are also 4-blades and so once again, taking the meet with a ray gives a 2-blade, F . With spheres, there can be zero, one or two points of intersection corresponding to the cases where F 2 < 0, F 2 = 0 and F 2 > 0 respectively.  If F = A ∧ B (with A and B null vectors) and F 2 ≥ 0, the points can be extracted from the point pair/blade, F , by the following formula [22]: where σ a , σ b are scalar constants. If we define F = A ∧ B with F oriented in the same direction as our ray L, then A is the point closest to the origin of the ray, P 0 , as long as our sphere is 'in front of' the ray source. To ensure the alignment we can pre-normalise our sphere S via the following expression: For any given sphere, its dual can square to a positive or negative number; however, by carrying out the normalisation in Eq. (6), we ensure that all spheres, S, satisfy S * · n ∞ = −1. If the meet of a ray L and a normalised sphere S, is then formed from L ∨ S, as in Eq. (3), the resulting bivector will be ordered as A ∧ B where A is the point that the ray hits first in its orientation.
In Fig. 4, the meet of the ray (direction as shown) from a point P 0 with the smaller sphere, would result in the point pair A 2 ∧ B 2 . For the larger sphere in Fig. 4, the point pair resulting from the meet will be A 1 ∧ B 1 . We extract the points from the point pair and form the distance between these points and P 0 (via taking the inner product). We then see that for the smaller sphere the distance of the first point is less than that of the second point, whereas for the larger sphere, the distance of the first point is larger than that of the second point-which will therefore lead us to label the larger sphere as being 'behind' the point P 0 . This allows us to perform bounces only with spheres that are in front of the ray origin point P 0 . • If Y 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 Y 2 < 0, the ray does not lie in the plane of the circle and passes through the circle disc but does not intersect. • If Y 2 > 0 the line does not lie in the plane of the circle and passes outside the circle disc without intersecting. • If Y 2 = 0 (and Y = 0) the ray intersects the circumference of the circle. Figure 5 shows an example of each case along with a geometric interpretation of the form of the meet. If Y = 0 and there is an intersection, the plane containing the circle is formed by taking the wedge product between the circle and n ∞ , C ∧ n ∞ , and the intersection point is then extracted from the ray and this plane.
If Y = 0, so that the ray lies in the plane of the circle, we need to work in 2D, so that our 'meet' will result from taking the 2-part of the geometric product and dualising (with respect to the plane of the circle) to give a bivector. If the bivector has negative square there are two intersections at points A and B, so the bivector is A ∧ B. If the bivector has positive square, there is no intersection. If the bivector squares to zero there is one intersection at A, and the bivector is a ∧ n 0 , where A = F (a). In both cases, the intersection points are easily extracted.

Extracting Normals and Reflecting Rays
Extracting the normal to the surface of an object at a ray intersection point, X, and the reflection of that ray at X, are two fundamental building blocks in our ray tracer.
For a plane Φ which intersects with a ray, we can compute the reflection L of an incident ray L (we assume Φ and L have been normalised such that Φ 2 = −1, L 2 = 1) with the plane by simple sandwiching: L = ΦLΦ. The resulting line is oriented correctly, passes through the intersection point and L 2 = 1. For the case of a sphere S, we use the following formula from [22]: where X is the first point of intersection. Here, SLS is an example of an inversion, where the incoming ray/line, L, is inverted in the sphere to give a circle which passes through the two points of intersection and the origin of the sphere (note we only have a meaningful reflection if there are two points of intersection). The tangent line to this circle at the first point of intersection, X, is the reflected ray, L . Figure 6 illustrates this geometrical construction. Note that one can also form the tangent plane at X and reflect L in this plane; this is performed by the following formula: Vol. 31 (2021) Exploring Novel Surface Representations Page 7 of 33 16 Figure 5. The 5D meet of a ray, L (in red) and a circle C (black) results in a vector Y whose dual is the sphere S (grey), the properties of this sphere vary with the relative positioning of C and L. In the top case the ray passes outside of the circle, S passes orthogonally through both the circle and the ray and squares to a negative number as we would expect from a standard CGA sphere. In the middle case the ray hits the perimeter of the circle and the meet squares to zero, in this case the dual sphere is the special case of zero radius, it represents the intersection point itself.
In the bottom case the ray passes inside the circle. Here the sphere squares to a positive scalar implying it is now an imaginary sphere and in fact the circle passes through the sphere's antipodal points For a circle/disc C, we first form the plane C ∧ n ∞ = Φ in which it lies. If the ray intersects the disc (see Sect. 4.1), the reflected ray can then be found using the same formula as for the plane reflection case: L ∝ (C ∧ n ∞ ) L (C ∧ n ∞ ). 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 paper. Firstly, consider the reflected ray, L , and the incident ray, L, both normalised such that they square to 1. The normal line to the surface, N , can be simply found: The tangent line, T , (in the plane containing incident and reflected rays) can similarly be found from the sum

Ray Tracing Evolved Circles
We will now turn to an interesting class of surface that arises from the direct interpolation of CGA circles [15], examples of which are shown in Fig. 8. In order to generate such a surface, a direct interpolation is first performed between two boundary circles, C 1 and C 2 both of which are normalised such in the plane containing the incident and reflected rays can be found from L T ∝ L + L that C 2 1 = C 2 2 = 1. Our interpolation is of the form: where we take α moving between 0 and 1, which moves us from C 2 to C 1 . 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 S. This projector has only scalar and 4-vector parts and its construction is detailed in [15] and outlined in the following. Consider a quantity Σ = Σ 0 + Σ 4 . We then define the quantity 2 4 , and with this the principal square root [14] of the scalar + 4-vector, Σ, can be found as: With this square root we can then form: where S − is S with the sign of the 4-vector part reversed and 1 k = S − S. We then construct kS by reversing the sign of the 4-vector part, (kS = kS − 0 − kS − 4 ), and use this to produce the following expression for the projector S and interpolated circle C α : 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 1vector Y . If Y = 0 the ray lies in the plane of the circle and if Y = 0 and Y 2 = 0 there is one intersection. Therefore (in the case where the meet is not zero) to find the intersection point between our interpolated surface and a ray L, we need to find a value of α for which: The system must also be tested for the case of Y = 0; if an α exists such that (C α ∨ L) = 0, the ray may intersect C α once, twice or not at all. Figure 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 (C α ∨ L) 2 = 0, 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 (C α ∨ L) 2  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 C 1 and C 2 , form intermediate circles, C α , and then calculate the bounding sphere which is given by 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 [18] which is summarised as follows: 1. Ensure both spheres S 1 , S 2 are normalised according to Eq. (6) 2. Construct the line L joining the centres of both spheres L = ( Intersect the line with the first sphere to produce a point pair L ∨ S 1 = F 1 ∝ A 1 ∧ B 1 and extract A 1 using Eq. (5) 4. Intersect the line with the second sphere to produce a point pair if so S 2 encloses S 1 and so S 2 is the bounding sphere 7. If neither original sphere is enclosed by the other, the new bounding sphere is given by 1 2 (A 1 + B 2 )I 5 This iterative algorithm for the most part performs perfectly satisfactorily. When compared against specially constructed test cases for which the intersection points are known it produces negligible error. The main downside to this solution is that it is not mathematically guaranteed to give correct results especially in the case of small numbers of intermediary objects. In practice we can pre-compute large numbers of intermediary objects before rendering, allowing us to get good approximations to the function of interest. Having said that, the more intermediate objects that are created, the more computationally expensive the process is, as the root finder has to evaluate our function at each one for each ray.

'Closed Form' Solution for the Intersection of a Ray and an Evolved
Circle Surface. In this section we will use C α to be the interpolated circle; however we emphasise that the process outlined here will also hold for the intersection of rays with other evolved objects. The intersection of the ray, L, and the surface, SC α (see Eq. (12)) occurs when: Writing Σ = −C αC α , this can be rewritten as: The denominator of this expression is never infinite (other than in the uninteresting case of Σ = 0) and so does not contribute roots. Thus 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: ] is a scalar and so distributing the grade selection operators gives us: Expanding this leads to: To make progress on solving this we recall that Σ = −C αC α and that for our linear interpolation of circles we have defined C α as: As C α has terms in α of order 1 we would expect Σ to have terms of order 2.
Continuing on this train of thought one might suspect that it is possible to re-write Eq. (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 Eq. (14) we need to rearrange: We can then square both sides of the equation, eliminating the square root in the process: Expanding this out will give a polynomial in α-it turns out that this is a scalar polynomial due to the fact that Σ 4 C α has only trivector components. This polynomial can then be solved with any numerical polynomial solver such as finding the eigenvalues of the companion matrix [20]. 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 0 ≤ α ≤ 1 and some will be spurious roots corresponding to S = 0. To filter out the valid roots we simply take all roots between 0 and 1 and evaluate (L ∨ [SC α ]) 2 at these positions, selecting the roots for which |(L ∨ [SC α ]) 2 | < for some small threshold where > 0 (in our experiments = 10 −6 works satisfactorily).
An interesting point to note here is that we could extend this intersection finding method to C being higher order functions of α, so long as (L ∨ [SC α ]) 2 = 0. Generating such higher order splines through geometric primitives is described in Sect. 8.

A Comment on Rendering Speed.
For this paper our raytracer was implemented in Python with the Clifford Library [1]. 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 paper.

Analytic Form for Normals
Given the α for which the ray intersects the surface, we have both the interpolated circle, C α , and the point of intersection X. Using the result from [22], which is also used in Eq. (7), we extract a tangential line L C in the plane of the circle at X: 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 L C : some future work remains to understand how these two tangent vectors are related to the directions of principal curvature. Clearly dC(α) dα ≡Ċ α 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.Ċ · C = 0, and that the geometric product is minus itself under reversion, i.e.Ċ C = −CĊ (note that here, and in what follows, we will drop the α subscript on C). This follows from the fact that C 2 = C · C = 1 (our circles are all normalised), so that: Since C ·Ċ = −C ·Ċ and they are both scalars, this tells usĊ · C = 0. Using the fact thatC = −C, we see that CĊ = −CĊ = −(CĊ)˜. 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 [12]). Let us call this bivector, Ω C : Vol. 31 (2021) Exploring Novel Surface Representations Page 15 of 33 16 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, X, have a constant 'length' due to normalisation, so as with the circles, we can differentiate wrt α to see that X andẊ are orthogonal, i.e. X ·Ẋ = 0. If we were to define the 'velocity',Ẋ to be the inner product with the angular velocity bivector given in Eq. (19): the condition X ·Ẋ = 0 is satisfied since X · (X · B) = (X ∧ X) · B = 0. Thus, given an X 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 X according to this rule, generating a quantity we call X(α), then X(α) should lie on C α , for all α, i.e., Differentiating this (and again dropping the subscript α for clarity) and usinġ X = X · (CĊ), giveṡ We now expand this expression using standard expansion results (a · (A r ∧ B s ) = (a · A r ) ∧ B s + (−1) r A r ∧ (a · B s )): The first term on the right hand side of the first line of this expansion is zero as (CĊ) ∧ C = CĊC 5 = −C 2Ċ 5 = −Ċ 5 = 0. Since X lies on C and so X ∧ C = 0, we see that XC = CX which means that givingẊ ∧ C = −X ∧Ċ 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: The fact that lines L C and L T are perpendicular can be verified by showing that the quantity L T L C has only a bivector part (see [22] for a discussion of when intersecting lines are orthogonal-if two lines L 1 and L 2 intersect at a point, then L 1 L 2 4 = 0. In addition, if they are orthogonal, L 1 L 2 0 = 0 ). If this is the case, L T L C will reverse to minus itself.
To show this, we need to consider the reverse of (Ẋ ∧ X ∧ n ∞ )((X · C) ∧ n ∞ ). We will need the facts that that XC = CX,ẊC = −CẊ, CĊ = −ĊC, XẊ = −ẊX andC = −C. We have shown all of these identities earlier in this section. We also need an additional fact, which is thatẊ anticommutes with C. To see this we use another standard result (a ∧ (A r · B s ) = (a · A r ) · B s + (−1) r A r · (a ∧ B s )): The first term on the RHS of this equation is zero as (CĊ) · C = CĊC 1 = −C 2Ċ 1 = 0, and the second term on the RHS is also zero as X lies on C so X ∧ C = 0. ThusẊ · C = 0 andẊ therefore anticommutes with C as required.
We are now in a position to expand out (Ẋ ∧ X ∧ n ∞ )((X · C) ∧ n ∞ ): In the above we have used the facts thatẊ ∧ X =ẊX, XC = X · C and X 2 = 0. Note that the term Xn ∞ X in the final line of Eq. (24) can be written as 2(X ·n ∞ )X (from the standard reflection formula and the fact that C 2 = 0). Reversing the final line of Eq. (24) and using the commutation and anticommutation relations discussed, it is easy to show that the reverse of L T L C is indeed minus itself, implying it has only a bivector part, as required, meaning the lines are orthogonal. Note that this result relies crucially on the fact thatẊ and C anticommute, which is a good indication thatẊ lies in the right direction.
Given these two orthogonal tangent lines L C and L T , we can construct the plane tangent to the surface at X by computing the join of the two lines. Or, we can bypass the plane entirely and compute the surface normal line directly as:

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 C α is any pure-grade multivector which is a function of a scalar parameter α. Let the projection of C α onto the blade manifold be given by: where S is our blade projector. The differential of this with respect to α is given by: Thus, any closed form expression for the derivative on the manifold will first require a closed form for the derivative of the projector ∂S ∂α . Recall from Eq. (13) that we can write the projector S as a function of √ Σ, where Σ = −C αCα , and so Thus to find an expression for ∂S ∂α 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: where we are using the fact that ∂Σ ∂α g = ∂ Σ g ∂α .

∂ ∂α
Thus the derivative of √ Σ is a function only of Σ and ∂Σ ∂α which in turn can be written in terms of C α and ∂C α ∂α :

Closed Form Derivative of the Projector
With our square root derivative in place we can proceed to finding the derivative of the projector S. Recall, S 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: Now consider C α to be an interpolated circle of the form C α = C α = αC 1 + (1 − α)C 2 . The derivative of this with respect to α is a constant: This derivative is the final piece required for Eq. (26), giving us a completely closed form for ∂Cα ∂α ≡Ċ. An important point to note here is that this blade projection derivative is grade-agnostic and so can be used for objects other than just evolved circles.

Ray Tracing Evolved Point Pairs
We will return to the actual ray tracing of circles later (see Figs. 16,17), but first we turn our attention to point pairs. Due to the mathematical similarities between circles and point-pairs in CGA [13], as well as the practical desire to represent ribbon-like surfaces, we can apply similar ray-tracing methods to surfaces formed from the interpolation of point-pair bivectors representing line segments. If P 1 and P 2 are point-pairs which represent a line segment, we form a surface via: where again, S is a scalar plus 4-vector which maps the interpolated bivector onto a 2-blade. Figure 10 gives examples of such surfaces.

Closed Form Solution for the Intersection of a Ray and an Evolved Point-Pair Surface
To find the intersection point of a ray and these surfaces we again form the meet of a ray L and the form of the evolved point-pair P α . 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 S in terms of Σ where as before, Σ = −P αP α , 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 Σ 4 P α has only bivector parts): which can again be solved with a fast numerical polynomial solver. This intersection equation for linearly interpolated point pairs is of order 6, implying there are up to 6 potential hitting points. Again the same process can be used to filter the roots as was done for the roots of the circle intersection equation.

Bounding Sphere and Normal Calculation
We saw earlier that the meet will be zero if the ray hits anywhere along the carrier line of the point-pair L C = P ∧ n ∞ . Assuming the carrier line and ray do meet, the point of intersection can be extracted via the method outlined in the last section of [15]. Given that the carrier line of P α (for some α) and the ray intersect at a point X, we can then check if the intersection point is within the bounding sphere S = P 1 ∧ P 2 of the surface by ensuring: Since the endpoints of all interpolated point-pairs will lie on the surface of S (see [15]), the above condition ensures there is an intersection with the line segment and not just the carrier line of the point-pair. To find the normal to the point-pair surface we can simply use exactly the same argument, and in fact the same code, as we did before for the evolved circles case but this time extracting L C as: Figure 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, Fig. 12 shows two examples of this. As proved in [15] 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 can be found by looking for a point at which: Re-arranging gives an expression for α: Vol. 31 (2021) Exploring Novel Surface Representations Page 21 of 33 16 Figure 11. Ray tracing evolved point-pairs. Left: the scene to be rendered, in blue is a representation of the point-pair surface to be rendered, the camera frustum is shown in black, the camera axis is shown in red. Right: the resultant rendered surface of interpolated point-pairs 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 P 1 ∧ P 2 , allows the ray-tracer to detect these cases 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 L and a set of three points A, B, C 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 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 [13,22].

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 [3], which in the cubic case and with specific first order endpoint conditions are known as Hermite curves.

Linear Interpolation as a Linear Bézier Curve
The simplest form of Bézier curve is simply a linear interpolation between two vectors. If we replace the vectors with k-blades and couple with the projection to the blade manifold we have the exact same linear interpolation, although this time with α going in the other direction. Adopting a notation of C 0 as the first object and C 1 as the second: In Sects. 5.2 and 6, our analysis to extract surface normals was based on having an expression for the derivative of the pure grade multivector as a

Quadratic Bézier Curve
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 14 shows examples of orders 1, 2 and 3 Bézier interpolation through circles, along with the control objects used to generate the surface.

N th Order Bézier Curve
More generally we can say that an N th order multivector Bézier curve is of the form are known as the Bernstein polynomials. The derivative of our N th order Bézier curve is: 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 w i 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, C i , (where we use the notation C for points as we will see shortly that these can be replaced by objects) and associated tangent vectors, V i , at each end of the curve, for α ∈ [0, 1]: The derivative of the curve is: Cubic Hermite curves can be converted to cubic Bezier curves and viceversa. 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: Let us now evaluate this at the endpoint of the nth curve in a piece-wise spline where α = 1 and of curve (n + 1) where α = 0. First we note that on both curves at these points, S = 1 because the curve passes through a blade control object which requires no projection. Additionally we see that by definition of the Hermite spline at that point, the derivative in pure-grade space is shared across both curves as is the control point: where V n,n+1 and C n,n+1 are the derivative (or 'tangent') and control objects respectively of the curve that are shared between segments n and n + 1. Thus for the derivative to evaluate to the same on either side of the boundary, we only need to check that ∂S ∂α is the same either side of the boundary. Considering the equations in Sect. 6.1 we can see that ∂S ∂α is a function only of C and ∂C ∂α 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 [21]. The Kochanek-Bartels (KB) spline is an interpolating spline with three scalar design parameters t, b, c known as tension, bias and continuity respectively. For given control objects C i , C i+1 the corresponding tangents V i , V i+1 can be calculated using the control objects in the spline C i−1 and C i+2 which lie previous to, and after, the curve in the order of the spline: Setting all three scalar parameters to a value of 0 produces the commonly used Catmull-Rom spline [7]. Figure 15 shows an example of a KB spline of multivector geometric primitives.  Figure 15. A Kochanek-Bartels spline of evolved circles meshed, textured, and rendered with smooth shading

Blinn-Phong Lighting Model
To complete our discussion of the ray-tracer we need to describe its lighting model. We employ the Blinn-Phong lighting model, this simple model is well studied in the graphics literature, for more in depth information see [4]. Under this model the intensity value for each pixel for colour channel λ is given by the following expression: Note that traditionally the terms N , L i and H are standard 3D vectors. Since the language of our ray-tracer is CGA, these terms in Eq. (34) actually represent trivector lines which pass through the ray-object intersection point. As all these lines are normalised such that L 2 = 1 the inner products between them will simply give the cosine of the angles between them as in the traditional 3D vector case [12]. Other terms will be described in the following subsections.
The normal line is required across multiple terms and is crucial to the lighting model. This is given in CGA as N = 2 where L is the normalised reflected ray and L is the normalised incident ray. Note that L is necessarily calculated for tracing reflected rays. The line L i specifying the direction to the ith point light source is given by forming the wedge product between the point of intersection of the incident ray and the object, the point position of the ith light source and n ∞ and again normalising the result such L 2 i = 1. The normalised half way line H can be found by where V is the normalised line from the intersection point to the viewer.

Diffuse Term
The diffuse term is given by: Vol. 31 (2021) Exploring Novel Surface Representations Page 27 of 33 16 The line L i again represents the normalised line pointing towards specifying the the ith light source. The intensity of channel λ of the ith light source is labelled I piλ where the p subscript denotes that it is a point light source. I piλ along with the material channel reflection coefficient c λ (which effectively controls the colour of the material) and k d the diffuse property of an object material, are all scalar parameters of the model.

Specular Term
This term requires the computation of the CGA line that corresponds to the half-way vector, H, defined as the vector halfway between the vectors to the viewer, V, and the light source, so that H ∝ L i − V . The normalised half way line H can be found by where V is the normalised line from the intersection point to the viewer. k s and q are scalar parameters of the object material.

Shadow Attenuation.
In order to determine whether an object is in shadow, the ray cast from the point of intersection to the light source is used for intersection tests with objects in the scene. If this yields an intersection, the shadow attenuation constant, S i , is used, otherwise it is omitted.

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 3 shows an example of simple objects, spheres, planes and disks being rendered. Figure 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 Fig. 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 Figure 16. Left: a scene composed of only an evolved surface in blue and a camera. Right: the rendering of the scene from the camera do this we will begin by producing a set of evenly spaced points on the first object C 1 and then transform these points along a small step in α to give a second set of points. Continuing in this way we can cover the surface entirely. An appropriate transformation for this task needs to preserve the relative spacing of the points on the objects in order to produce a good quality mesh. TRS (Translation Rotation Scaling) rotors have this property and can map circles to circles, spheres to spheres and point pairs to point pairs (these quantities are sometimes known as rounds). A TRS rotor that takes one object C 1 , to another, C 2 , can be calculated with the following process: • Calculate T 1 and T 2 the translation rotors that bring C 1 and C 2 respectively to the origin • Apply T 1 and T 2 to C 1 and C 2 respectively, bringing them to the origin and producing C 1 and C 2 • Calculate the rotation rotor R 12 between the blades C 1 ∧n ∞ and C 2 ∧n ∞ (if C 1 and C 2 are spheres then we do not need a rotation rotor so set R 12 = 1) • Calculate the difference in scale between the objects by extracting their relative sizes • Use the scale to generate a dilation rotor D 12 that scales C 1 to the same size as C 2 • Compose the final TRS rotor Z 12 that takes C 1 to C 2 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 N evenly spaced points on the unit circle in the e 1 , e 2 plane by N successive rotations about the origin of a point X lying initially at X 0 = F (e 1 ) yielding X n for n ∈ 0, . . . , N i.e., for a fixed rotor R θ : where θ is chosen so that N + 1 uniformly spaced points cover the whole circle. With the TRS rotor Z 01 that maps from the unit circle at the origin to the first object C 1 it is possible to transform our points to the first object: 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 Z α to refer to the TRS rotor that maps C α to C α−δα and so can write the nth point at α as: By doing this we have effectively constructed a mapping from a coordinate system in the 2D plane of α and n to the surface manifold. This is useful as it lets us generate the mesh in the 2D plane of α and n 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 Fig. 18 shows a surface of evolved circles meshed and rendered using flat shading with ganja.js [9], Fig. 15 shows a tubular KB spline surface, meshed, textured, and shaded with a smooth vertex normal interpolation scheme.

Summary and Conclusions
In this paper 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.