Geometric Algebra: A powerful tool for solving geometric problems in visual computing Leandro A. F. Fernandes, and Manuel M. Oliveira Instituto de Informática, UFRGS Porto Alegre, RS, Brazil Email: {laffernandes,oliveira}@inf.ufrgs.br

Abstract—Geometric problems in visual computing (computer graphics, computer vision, and image processing) are typically modeled and solved using linear algebra (LA). Thus, vectors are used to represent directions and points in space, while matrices are used to model transformations. LA, however, presents some well-known limitations for performing geometric computations. As a result, one often needs to aggregate different formalisms (e.g., quaternions and Plücker coordinates) to obtain complete solutions. Unfortunately, such extensions are not fully compatible among themselves, and one has to get used to jumping back and forth between formalisms, filling in the gaps between them. Geometric algebra (GA), on the other hand, is a mathematical framework that naturally generalizes and integrates useful formalisms such as complex numbers, quaternions and Plücker coordinates into a highlevel specification language for geometric operations. Due to its consistent structure, GA equations are often universal and generally applicable. They extend the same solution to higher dimensions and to all kinds of geometric elements, without having to handle special cases, as it happens in conventional techniques. This tutorial aims at introducing the fundamental concepts of GA as a powerful mathematical tool to describe and solve geometric problems in visual computing. Keywords-geometric algebra; geometric computing; applied mathematics; Plücker coordinates; quaternion; subspaces

I. M OTIVATION Visual computing practitioners often use linear algebra to represent geometric elements and transformations [1], [2]. For instance, vectors represent directions and, by introducing an origin to the underlying space, they can also be used to represents points. But the semantic difference between a direction vector and a point vector is not encoded in the vector type itself, and one must be aware that point vectors should be transformed differently than direction vectors [3]. In linear algebra, geometric elements of higher dimensionality (e.g., lines and planes) must be assembled from points, directions and scalar values. This way, a line in 3-D can be represented by a point and a direction vector, or by the intersection of two planes. A plane, on the other hand, is usually described by a normal vector and a distance from the origin. Such a representation of geometric elements prevents their use as computing primitives. As a result, geometric problems must be solved considering the particular type of geometric object being handled. For instance, the intersection between

two planes is computed using a different solution than the intersection between a line and a plane, or any other combination of elements. Also, each possible output must be handled on a case-by-case basis (e.g., a plane-plane intersection can result in an empty set, a straight line, or a plane). Alternatively, one can use Plücker coordinates to represent points, lines and planes as elementary types in order to operate on them [4]. However, defining a matrix for transforming geometric elements represented by Plücker coordinates is not straightforward. It is well known that, under projective geometry, matrices can be used to encode a sequence of projective transformations into a single representation. However, it is also known that matrices are not the best way to handle rotations, because they are not suitable for interpolation, and they encode the rotation axis and rotation angle in a costly way (i.e., they require several coefficients). In such a case, quaternions offer an alternative and efficient solution [5], [6]. But this extension is not automatically well-connected to other transformations in linear algebra, or to Plücker coordinates. Also, quaternions are defined only in 3-D. Geometric algebra [7]–[10] provides an alternative to conventional geometric techniques based on linear algebra. With geometric algebra one can compute directly on geometric elements (e.g., direction, points, lines, planes, circles, and spheres) as primitives. Also, it naturally generalizes and integrates useful frameworks such as Plücker coordinates, quaternions and complex numbers into a high-level specification language for geometric operations. Because of its consistent structure, geometric algebra equations are often universal and generally applicable, extending the same solution to higher dimensions and to all kinds of geometric elements. Geometric algebra, however, is still mostly unknown in the computer science community. Only in the past few years it became accessible in specialized literature [8]–[10]. This tutorial aims at providing a gentle introduction to the fundamental concepts of geometric algebra and discussing its great potential as a tool for representing and solving problems in computer graphics, computer vision, and image processing. This tutorial is organized as follows: Sections II to VI

present theoretical concepts of geometric algebra. In these sections, geometric primitives, transformations, and algebraic operations are presented in a more abstract level. For sake of simplicity, we kept the mathematical formalisms to a minimum, and illustrate the geometrical intuition with examples whenever possible. Section VII discusses how one can assume a model of geometry in order to give practical geometric interpretation to the elements discussed in previous sections. Geometric algebra implementations are enumerated in Section VIII. Section IX concludes the tutorial with some observations and directions for research opportunities involving geometric algebra. II. O RIENTED S UBSPACES AS P RIMITIVES By definition, a vector space Rn consists of a set of elements called vectors. By assuming a basis {eei }ni=1 for Rn , an arbitrary vector (i.e., a 1-D subspace) is written as a linear combination of the basis vectors. Fig. 1a illustrates a vector a = α1 e 1 + α2 e 2 + α3 e 3 in R3 . In this graphical representation, the 1-D subspace is the oriented straight line that passes through the origin and supports the arrow. The arrow’s length represents the weight of the subspace, while the arrow’s tip indicates its orientation. In geometric algebra, a 2-D subspace is spanned as the outer product of two vectors. The outer product is formally defined in Section II-B. Fig. 1b shows a 2-D subspace computed as: (1) C 2 = a ∧ b , where ∧ (the wedge symbol) denotes the outer product, and a and b are two vectors (1-D subspaces). In Fig. 1b, C 2 is the oriented supporting plane for the disk. The disk’s radius and the curved arrow are used to illustrate, respectively, the weight and orientation of the subspace. Note that the orientation of C 2 respects the order of the terms in the outer product, i.e., from a to b . Using the outer product one can span k-D oriented subspaces, for 0 ≤ k ≤ n. In geometric algebra, such subspaces are called k-blades (the terms blade and subspace are used interchangeably), and k is said to be the grade of the blade (the terms grade and dimensionality are also used interchangeably). An arbitrary blade B k presents the following properties: attitude The equivalence class α B k , for any α ∈ R. weight The value of α in B k = α J k , where J k is a reference blade with the same attitude as B k . orientation The sign of the weight relative to J k . direction The combination of attitude and orientation. A. A basis for k-D subspaces In contrast to linear algebra, geometric algebra does not represent k-D subspaces as a collection of k vectors.

(a)

(b) R3 .

Figure 1. Graphical representation of subspace in (a) A 1-D subspace a is drawn as an arrow, where the arrow’s length and the arrow’s tip represent, respectively, the weight and the orientation of a . (b) A 2-D subspace C 2 , spanned from vectors a and b , is drawn as a disk. Its weight is represented by the radius of the disk, while the orientation is the curved arrow. The basis vectors {ee i }3i=1 are not shown in this image for sake of clarity.

In fact, k-blades are “concrete” primitives by themselves, spanned from k vectors. Therefore, one needs a new basis for representing k-blades, because the basis vectors from support only 1-D subspaces. In geometric Rn , by definition,  algebra, Rn is the multivector space built  nfrom a vector space Rn . The 2n basis elements of R are defined as the k-combinations of vectors from the set {eei }ni=1 n n n (i.e., k=0 k = 2 ). A linear combination of such basis elements is called a multivector. As an example, the basis of  R3 is presented in Fig. 2 (top) and its multivector structure is shown at the bottom. By combining zero elements out of the set {eei }3i=1 , one takes the basis blade “1” for 0-D subspaces, i.e., values in R. As one would expect, the basis blades for 1-D subspaces are 2the3 three basis vectors defining R , the three basis blades are R3 . For 2-D subspaces in computed by combining two elements out of the set of three basis vectors. Finally, the outer product of all basis vectors define basis blade for 3-D subspaces. Note that blades 3 the R3 are a scaled version of the whole 3-D space. Such in blades are often called pseudoscalar, and the unit positive pseudoscalar is typically denoted by: I n = e 1 ∧ e 2 ∧ · · · ∧ e n .

(2)

The order of the vectors defining each basis blade (e.g., use e 1 ∧ e 2 or e 2 ∧ e 1 for the fifth basis blade in Fig. 2) is a matter of convention. Fortunately, converting between different conventions is easy, and relies on applying the properties of the outer product. B. The outer product Formally, the outer product is a mapping: r s r+s ∧: Rn × Rn → Rn .

(3)

Scalars  R = 0 R3

Vector Space  R3 = 1 R3

e2 ∧ e 3,



e1 ∧ e 2 ∧ e 3

Bivector Space 2 3 R

⎪ ⎪ ⎪ ⎩

e1 ∧ e 3,

⎧ ⎪ ⎪ ⎪ ⎨

e1 ∧ e 2,

e3,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

e2,

⎪ ⎪ ⎪ ⎪ ⎩

e1 ,

⎧ ⎪ ⎪ ⎪ ⎪ ⎨



1,

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩



Trivector Space 3 3 R

M = γ1 + γ2 e 1 + γ3 e 2 + γ4 e 3 + γ5 e 1 ∧ e 2 + γ6 e 1 ∧ e 3 + γ7 e 2 ∧ e 3 + γ8 e 1 ∧ e 2 ∧ e 3 

Figure 2. Basis blades for multivector space R3 (top). A linear combination of such basis elements (eight altogether) is called multivector (bottom).  3 2 3 Multivectors in R can be used to represent real values (0-D subspaces in R), vectors (1-D subspaces in R3 ), 2-blades (2-D subspaces in R ), 3 3 R ). Also, as commented in Section VI, some multivectors having mixed dimensionality encode orthogonal and the whole space (3-D subspaces in transformations, such as reflections and rotations.

It is defined from a small set of properties: antisymmetry distributivity associativity scalars commute

b ∧ a , thus c ∧ c = 0 a ∧ b = −b b + c) = a ∧ b + a ∧ c a ∧ (b b ∧ c ) = (aa ∧ b ) ∧ c a ∧ (b a ∧ (β b ) = β (aa ∧ b )

a

θ

Table I presents a step-by-step example where the properties of the outer product are used to build a 2-blade C2 2 3 in R from vectors a and b in R3 . This example is equivalent to the one depicted in Fig. 1b. It is important to note that the same set of properties are used while spanning  subspaces from two blades having arbitrary grades in Rn . III. M ETRIC AND S OME I NNER P RODUCTS The concept of subspace and the construction of blades from the outer product is independent of any metric proper ties a vector space Rn , and hence a multivector space Rn , might have. However, in order to compare the weight or the angle between two subspaces having different attitudes, one needs products which depend on the metric of Rn , extending the usual vector inner product from vector algebra to higher dimensional subspaces. The vector inner product multiplies two arbitrary vectors and returns a scalar value that characterizes their metric relation: a · b = Q (aa, b ) .

(4)

It is denoted by · (the center dot symbol) and presents the following properties: symmetry a · b = b · a b + c) = a · b + a · c distributivity a · (b scalars commute a · (β b ) = β (aa · b ) In Eq. 4, Q is a scalar-valued function defining a metric on the vector space Rn . A practical way to implement Q is

θ

b

(a)

(b)

Figure 3. (a) It is well known that, in a Euclidean vector space, the vector inner product can be used to compute the smallest angle between two a · b ) / (a a b b)). (b) In the same way, the scalar vectors: θ = cos−1 ((a product introduced in Section III-A allows one to measure the smallest angle between the having   same dimensionality:  two arbitrary  subspaces  θ = cos−1 A 2 ∗ B 2 / A 2  B 2  .

by using a metric matrix: ⎡ μ11 μ12 ⎢ μ21 μ22 ⎢ M=⎢ . .. ⎣ .. . μn1 μn2

··· ··· .. .

⎤ μ1n μ2n ⎥ ⎥ .. ⎥ , . ⎦

(5)

· · · μnn

where μij = Q(ee i , e j ), for 1 ≤ i, j ≤ n. This way, M is a symmetric positive-definite matrix encoding the inner product of pairs of basis vectors {eei }ni=1 . In Section VII we show how one can give meaningful geometric interpretations for subspaces and related operations just by assuming a metric for Rn . For instance, by letting M be an identity matrix one is assuming Euclidean metric (i.e., e i ·eej is one for i = j, and zero otherwise), and the inner product becomes the well known vector dot product from linear algebra: b cos θ. a · b = aa  b

(6)

The vector dot product is typically used to measure the smallest angle between two vectors, as depicted in Fig. 3a.

Table I S TEP - BY- STEP ALGEBRAIC MANIPULATION OF THE OUTER PRODUCT OF TWO ARBITRARY VECTORS IN R3 . Step 1.

Algebraic Manipulation C 2 = a ∧ b

2.

C 2 = (α1 e 1 + α2 e 2 + α3 e 3 ) ∧ (β1 e 1 + β2 e 2 + β3 e 3 )

3.

C 2 =

4.

C 2 =

Comments The resulting blade C 2 is computed as the outer product of vectors a = α1 e1 + α2 e2 + α3 e3 and b = β1 e1 + β2 e2 + β3 e3 . Simple substitution of a and b .

α1 e 1 ∧ (β1 e 1 + β2 e 2 + β3 e 3 ) + α2 e 2 ∧ (β1 e 1 + β2 e 2 + β3 e 3 ) + α3 e 3 ∧ (β1 e 1 + β2 e 2 + β3 e 3 )

Distributivity over addition, from right to left.

α1 β1 e 1 ∧ e 1 + α1 β2 e 1 ∧ e 2 + α1 β3 e 1 ∧ e 3 + α2 β1 e 2 ∧ e 1 + α2 β2 e 2 ∧ e 2 + α2 β3 e 2 ∧ e 3 + α3 β1 e 3 ∧ e 1 + α3 β2 e 3 ∧ e 2 + α3 β3 e 3 ∧ e 3

Distributivity over addition, from left to right, and scalar commute.

5.

C 2 =

α1 β2 e 1 ∧ e 2 + α1 β3 e 1 ∧ e 3 − α2 β1 e 1 ∧ e 2 + α2 β3 e 2 ∧ e 3 − α3 β1 e 1 ∧ e 3 − α3 β2 e 2 ∧ e 3

According to antisymmetry, e j ∧ e i = −eei ∧ e j for i = j, and ek ∧ ek = 0. By making i n (Eq. 23), and the left contraction is zero when r > s (Eq. 25). Note that one can apply Eqs. 23 and 24 onto the definition of the geometric product of vectors (Eq. 19) in order to retrieve, respectively, its outer and inner parts alone. V. D UAL R EPRESENTATION OF S UBSPACES By treating subspaces as sets, one can use a Venn diagram (Fig. 7a) to show that an arbitrary k-blade A k embedded in a n-D space I n induces an (n − k)-blade. The induced subspace is the orthogonal complement of A k with respect to the whole space. Indeed, the number of basis blades in

(a)

(b)

Figure 7. Dual representation of subspaces. (a) By using a Venn diagram, a blade A k can be seen as a set (the circle) contained in another set representing the entire n-D space (II n , the interior of the surrounding square). This way, the orthogonal complement of A k is the (n − k)-D subspace included in I n and disjoint from A k . (b) Example of dual relationship in a Euclidean vector space R3 . Here, vector d is the dual representation of the blade B 2 .

 each graded part of a multivector space Rn (Section II-A) suggests the existence of a relation between k-blades and (n − k)-blades. symmetry in the number of basis  Such  elements of k Rn and n−k Rn is related to the  equivalence  n  ). on the combinations of basis vectors (i.e., nk = n−k In geometric algebra, the dual representation of a subspace A k is its (n − k)-D orthogonal complement, with the same absolute weight of A k , and a well-defined orientation. The dual of A k is computed as: A ∗k = A k  I −1 n ,

(26)

where ∗ (superscripted asterisk) denotes the taking the dual operation, and I n is the unit pseudoscalar of the n-D space (Eq. 2). Optionally, one can replace the left contraction in Eq. 26 by a geometric product because they are equivalent in this context, since A k ⊆ I −1 n . However, the use of the left contraction provides Eq. 26 with a clear geometric interpretation: returning the portion of I −1 n that is orthogonal to A k . It is important to emphasize that A ∗k is an (n− k)-blade, as suggests the outcome of A k contracted on I −1 n in Eq. 26. Therefore, the dualization defines a mapping: k n−k Rn → Rn . (27) ∗ : Fig. 7b shows a blade B 2 and its dual representation d = B ∗2 in a Euclidean 3-D space. Table II presents the step-by-step algebraic manipulation of this dualization operation regarding an arbitrary  2-blade B 2 = γ5 e 1 ∧ e 2 + γ6 e 1 ∧ e 3 + γ7 e 2 ∧ e 3 in 2 R3 . The result is d = γ5 e 3 − γ6 e 2 + γ7 e 1 . Note that the coefficients of both direct and dual representation of B 2 are the same, up to a sign in γ6 . In this way, a vector is naturally associated with a 2-blade in 3-D space. Visual computing practitioners

Table II 2 3 S TEP - BY- STEP ALGEBRAIC MANIPULATION FOR COMPUTING THE DUAL OF AN ARBITRARY 2- BLADE IN R . Step 1.

Algebraic Manipulation d = B ∗2

2.

d = B 2  I −1 n

Simple substitution of B ∗2 according to Eq. 26.

3.

d = (γ5 e 1 ∧ e 2 + γ6 e 1 ∧ e 3 + γ7 e 2 ∧ e 3 )  I −1 n

Simple substitution of B 2 .

4.

d= + +

γ5 (ee1 ∧ e 2 )  I −1 n γ6 (ee1 ∧ e 3 )  I −1 n γ7 (ee2 ∧ e 3 )  I −1 n

Comments The resulting vector d is the dual representation 2 3 of the 2-blade B2 = γ5 e1 ∧ e2 + γ6 e1 ∧ e3 + γ7 e2 ∧ e3 in R .

Distributivity over addition, from right to left.

5.

d = −γ5 (ee1 ∧ e 2 )  (ee1 ∧ e 2 ∧ e 3 ) − γ6 (ee1 ∧ e3 )  (ee1 ∧ e2 ∧ e3 ) − γ7 (ee2 ∧ e 3 )  (ee1 ∧ e 2 ∧ e 3 )

= I n = −ee1 ∧ e 2 ∧ e 3 . In turn, the scalar Simple substitution of I −1 n values −1 were commuted to the left side of the terms.

6.

d = γ5 e 3 − γ6 e 2 + γ7 e 1

By solving the left contractions, the result is a vector having the same coefficients as B 2 , up to a sign in γ6 .

are familiarized with such observation. In vector algebra, the outcome of a cross product is precisely the normal vector of the plane induced by the input vectors. Algebraically, it means that: ∗ (28) a × b ≡ (aa ∧ b ) The equivalence of the cross product (from vector algebra) and the dual of the outer product assuming a right-handed coordinate system can be confirmed by replacing the outcome of a ∧ b (Table I) and B ∗2 (Table II) into Eq. 28: a ×b =

B ∗2

= γ5 e 3 − γ6 e 2 + γ7 e 1 = (α2 β3 − α3 β2 ) e 1 + (α3 β1 − α1 β3 ) e 2

D −∗ n−k = D n−k  I n .

(31)

This way, the dual representation of a blade can be correctly mapped back to its direct representation: −∗    A ∗k = A k  I −1 n  I n = A k I −1 n I n

(32)

= A k . I −1 n I n

= 1 ensures that the relative orientation In Eq. 32, of A k will be preserved for all n. (29)

+ (α1 β2 − α2 β1 ) e 3 . It is important to note that the use of the cross product in vector algebra is restricted to 3-D spaces because the cross product is defined only in 3-D. On the other hand, the right-hand side of Eq. 28 applies to spaces having arbitrary dimensionality. At first sight, the dual of the dual representation of a blade should result in the direct representation of the blade. A ∗k )∗ = A k does However, as demonstrated in Eq. 30, (A not hold in the general case: ∗    −1 A ∗k = A k  I −1 n  I n −1 = A k I −1 n I n

−1 n(n−1)/2 . In order to correct this issue, it is I −1 n I n = (−1) necessary to define the taking the undual operation as:

(30)

= (−1)n(n−1)/2 A k . Eq. 30 shows that the successive application of two dualization operations may change the orientation of the resulting blade according to the value of n. It is because

A. Dual relationship of products The duality of subspaces allows the definition of important dual relationships between the products of geometric algebra. For instance, the dual of the outer product can be replaced by a left contraction: ∗    A r ∧ B s = A r ∧ B s  I −1 n   −1 (33) = A r  B s  I n = A r  B∗s . Conversely, the dual of a left contraction can be written as an outer product: ∗    A r  B s = A r  B s  I −1 n   −1 (34) = A r ∧ B s  I n = A r ∧ B ∗s . The use of the dualization allows the algebraic manipulation of expressions, replacing spanning by contraction, and vice versa, whenever it is convenient. The proof of the dual relationships depicted in Eqs. 33 and 34, is presented in Dorst et al.’s book [9].

a

VI. O RTHOGONAL T RANSFORMATIONS AS V ERSORS In geometric algebra, the reflection of an arbitrary vector a with respect to an invertible vector v is obtained using a sandwiching construction involving the geometric product: v a v −1 . a = −v

(35)

A sequence of reflection operations results in an orthogonal transformation applied to a . For instance, an even number of reflections gives a rotation while an odd number of reflections represents a reflection followed by a rotation, or a rotation followed by a reflection. Fig. 8 shows how the double reflection of a in vectors p and q is equivalent to the rotation of a in the p ∧ q plane. In such a case, the rotation angle φ is twice the angle between p and q . The sense of rotation is given by the orientation of the 2-blade (i.e., from p to q ) and by the sign of φ. The associativity of the geometric product allows the grouping of successive k reflection vectors applied to a :   −1 −1 v k · · · v 2 v 1 ) a v −1 a  = (−1)k (v 1 v2 · · · vk (36) = (−1)k V a V −1 , where V is a k-versor, which is defined as the geometric product of k invertible vectors. As for the inverse of blades, the inverse of V is computed using Eq. 18. As pointed out in Section IV-A, versors may have mixed grade. A versor made up of an even number of vectors is called an even versor, and its multivector is a linear combination of basis blades having even grades. Likewise, odd versor is computed from an odd number of vectors, and it is represented by a multivector obtained as a linear combination of odd-graded basis elements. The evenness or oddness of a versor is important for defining the sign change in Eq. 36. Note that (−1)k is positive for all even versors, and negative otherwise. The versor product is distributive over addition. Also, versors preserve the structure of the outer product. Therefore, the versor product can be extended from vectors (Eq. 36) to any other element represented by a multivector:  V X V −1 for even versors  X = (37)  V −1 for odd versors V X  denotes the grade involution of the basis blades In Eq. 37, X in the multivector X . The grade involution of a k-blade is  k = (−1)k A k . It exhibits a + − + − + − · · · pattern A over the values of k. As in Eq. 36, even versors in Eq. 37 do not change the sign of the transformed element. However, for odd versors, one should be careful about the sign of the resulting element. In such a case, the grade involution guarantee the proper signs of the various grade parts of X . As orthogonal transformations, versors preserve both the symmetric inner product (Eq. 4) and the outer product. This way, the structure preservation of versors holds for

p

q (p a p -1 ) q -1

/2 q − p a p -1

Figure 8. Let a be a vector in a Euclidean space. By reflecting a in vector p and, in turn, in vector q , one gets the same result as by rotating a in the p ∧ q plane by φ radians, where φ/2 is the angle from p to q .

the geometric product, and hence to all other products in geometric algebra:     (38) V (A ◦ B ) V −1 = V A V −1 ◦ V B V −1 . In Eq. 38, the ◦ symbol represents any product of geometric algebra, and, as a consequence, any operation defined from the products (e.g., inversion, projection, and duality). Unit versors encoding rotations are called rotors. As in Fig. 8, a rotor R is build from an even number of unit vectors. In contrast to other versors, the inverse of a rotor is  R. its reverse (Eq. 10), i.e., R −1 = R An alternative (and more practical) way to define rotors is to use the exponential of 2-blades. Under Euclidean metric, the rotor R encoding a rotation of φ radians in the unit plane B 2 (the case depicted in Fig. 8) is given by:       φ φ φ R = exp − B 2 = cos − sin B 2 . (39) 2 2 2 By using the exponential form one can easily define a rotation in an arbitrary plane without being concerned about the handedness of the space. This is because the sense of rotation is related only to the given angle and to the rotation plane. Also, the exponential form allows the definition of rotors directly from the 2-blades, instead of from pairs of vectors spanning them. In Euclidean 3-D spaces, rotors are closely related to unit quaternions (see [9], [10] for details). In fact, rotors are the generalization of quaternions to n-D spaces. In Section VII-C we show that the exponential of blades can be used to define rotors encoding rotation, translation, and isotropic scaling. For this, we need to define the exponential of k-blades in an arbitrary metric space: 2 3   A k A k A k + + + ··· exp A k = 1 + 1! 2! 3! ⎧ sin α for A 2k < 0 ⎪ ⎨cos α + α A k for A 2k = 0 , = 1 + A k ⎪ ⎩ α 2 cosh α + sinh α A k for A k > 0 (40)

where α =

 A 2k ). abs(A VII. M ODELS OF G EOMETRY

So far, we have presented the most fundamental concepts of geometric algebra. In this section, we present three geometric algebra models that can be used in daily computations. Selecting a model implies in assuming a metric for Rn (see Section III) and provides a practical geometric interpretation to blades and versors. Supplementary Material B presents a quick reference for defining geometric primitives and transformations in each model, from parameters that are typically used with linear algebra. A. The Euclidean vector space model As the name suggests, in the Euclidean vector space model one assumes Euclidean metric for Rn . This way, k-blades are geometrically interpreted as k-D Euclidean subspaces, i.e., oriented flats (e.g., straight lines, planes, and their higher-dimension counterparts) that pass through the origin of the vector space. In this model, versors encode reflections and rotations. We have used the Euclidean model to illustrate all examples given so far (Figs. 1, 3–6, 7b, and 8). Euclidean subspaces in Rn are important because they represent the solution set to any homogeneous system of linear equations with n variables. For instance, consider the following system: =0 2 e1 − 3 e 2 (41) e 1 − 2 e2 + 3 e3 = 0 Each equation of the system is associated with a plane that passes through the origin of R3 . As depicted in Fig. 9, the vectors f 1 = 2 e 1 − 3 e 2 and f 2 = e 1 − 2 e 2 + 3 e 3 are the normal vectors (i.e., the dual representation) of such planes. The solution set is the intersection of the planes, which can be computed using the outer product: −∗ (ff 1 ∧ f 2 ) = 9 e 1 + 6 e 2 + e 3 .

(42)

Here, the 2-blade spanned as the outer product of f 1 and f 2 is the dual of the solution. The final solution is obtained by taking its undual. The undual operation is defined in Eq. 31. The resulting subspace will be zero when the system has no solution. Note that the technique presented in Eq. 42 can also be used to solve underdetermined systems (i.e. systems with more variables than the number of homogeneous linear equations plus one). In such a case, the result is a subspace whose dimensionality is higher than one. Visual computing practitioners often use techniques like the Cramer’s rule to solve systems of linear equations. However, the products of geometric algebra give a clear geometrical interpretation.

Figure 9. The outer product can be used to solve homogeneous systems of linear equations. Here, f 1 and f 2 are the normal vectors of the planes related to the equations of the system presented in Eq. 41. The solution set is the vector defined by the intersection of the planes. It is compute as (ff 1 ∧ f 2 )−∗ .

B. The homogeneous model The homogeneous (or projective) model in geometric algebra is similar to the use of homogeneous coordinates in linear algebra. It assumes a vector space Rd+1 with basis {ee0 , e 1 , e 2 , · · · , e d } and Euclidean metric. In this model, the d-D working space (i.e., the space where the geometric primitives reside) is embedded in Rd+1 in such a way that the extra basis vector e0 is interpreted as the origin of the working space. In Fig. 10a, the plane parallel to e1 ∧ e2 is the homogeneous representation of the 2-D working space in Fig. 10b. In the homogeneous model, vectors are geometrically interpreted as points. A proper point is a vector defining a finite location (α1 , α2 , · · · , αd ) in the working space. Such a location is given by the intersection of the 1-blade with the working space (see e 0 , a , and b in Fig. 10). Unit proper points are written in the form: p = e 0 + α1 e 1 + α2 e 2 + · · · + αd e d .

(43)

Note that the coefficient assigned to e 0 in Eq. 43 is equal to one. A general proper point γ p is a weighted version of a unit point, and it is interpreted as having the same location. When a vector is parallel to the working space (as e 1 , e 2 , and c in Fig. 10) it is called an improper point, or a point at infinity. Improper points can be seen as directions, because they are in the purely directional space Rd of the total space Rd+1 . Unlike proper points, directions have the coefficient of e 0 equal to zero: d = β1 e 1 + β 2 e 2 + · · · + β d e d .

(44)

Higher dimensional oriented flat subspaces, like straight lines and planes, are spanned as the outer product of proper and improper points. For instance, the line in Fig. 10 is defined as L 2 = a ∧ b . Optionally, one can create a line from a proper point and a direction by using the same construction. The result will be a 2-blade to be used as a computing primitive.

In the homogeneous model, 3-blades are geometrically interpreted as planes. As one would expect, they are defined in d-D working spaces (for d ≥ 3) as the outer product of: (i) three proper points; (ii) two proper points and one direction; or (iii) one proper point and two directions; as far the vectors are linearly independent. As one can see, the definition of k-flats (for 0 ≤ k < d) in the homogeneous model is straightforward. It is based on the outer product of (k + 1) vectors. Blades spanned exclusively from improper points are geometrically interpreted as k-flats at infinity. Unfortunately, versors in the homogeneous model are not as powerful as in the conformal model (Section VII-C). Here, they can be used only to encode rotations around the origin. In such a case, a rotation is defined as in Eq. 39, and B 2 must be a 2-blade spanned from direction vectors. Fortunately, these rotations can be applied to any blade by using Eq. 37. In the Supplementary Material B we present formulas to perform translations, general rotations in a plane at an arbitrary location, and rigid body motion in the geometrical primitives of the homogeneous model. C. The conformal model The conformal model outperforms the homogeneous model. Unlike the homogeneous model, in the conformal model blades can be geometrically interpreted not only as directions and flats, but also as rounds (e.g., point pairs, circles, spheres, and their higher-dimension counterparts) and tangent subspaces; while versors can combine reflections, rotations, translations, and isotropic scaling. The total vector space Rd+2 of the conformal model has basis {oo, e 1 , e 2 , · · · , e d , ∞ }, where the d-D working space is enhanced with two extra dimensions: o , a null vector interpreted as the origin point (pronounced “no”); and ∞ , a null vector interpreted as the point at infinity (pronounced “ni”). They are null vectors due to the special metric assumed in the conformal model. Eq. 45 shows the multiplication table for the vector inner product (Eq. 4) of the basis vectors: · o e1 e2 .. . ed ∞

o e1 0 0 0 1 0 0 .. .. . . 0 0 −1 0

e2 0 0 1 .. . 0 0

· · · ed ··· 0 ··· 0 ··· 0 . .. . .. ··· 1 ··· 0

∞ −1 0 0 .. .

(45)

0 0

Note that o · o = ∞ · ∞ = 0, while o · ∞ = −1. One consequence of this definition is that the inner product of two unit finite points (i.e., points at a finite distance from the origin) is given in terms of the square of the Euclidean distance between them: p ·q = −

d 1  2 (αi − βi ) , 2 i=1

(46)

e2 a e0

e1

c b LX2\

(a)

(b)

Figure 10. In (a), the plane parallel to e 1 ∧ e 2 is the homogeneous representation of the 2-D working space in (b). In the homogeneous model, the geometric interpretation of blades is given by their intersection with the working space. For instance, vectors e 0 , a , and b in (a) are interpreted as proper points in (b), while vectors e 1 , e 2 , and c in (a) are interpreted as improper points, or directions, in (b). The straight line in (b) is defined by the intersection of a 2-blade with the working space in (a). In such a case, L 2 = a ∧ b .

where (α1 , α2 , · · · , αd ) and (β1 , β2 , · · · , βd ) are the location of points p and q , respectively. This way, finite points are also null vectors (i.e., p · p = 0). Have the inner product related to the Euclidean distance of points is an interesting feature of the conformal model. That allows the definition of a coordinate-free solution, because the comparison of points is independent of their location relative to an assumed origin. In the homogeneous model, the outcome of the inner product of proper points depends on how far they are from e 0 , preventing the use of versors to encode all kinds of orthogonal transformations. Unit finite points are written in the form: p = o + α1 e 1 + α2 e 2 + · · · + αd e d +

d 1   2 α ∞ , (47) 2 i=1 i

while general finite points are weighted points (γ p ) having the same location. Fig. 11a shows that the set of all unit finite points in a 2-D working space defines a paraboloid in the ∞ -direction of the total space. In this example, the total 4-D vector space is presented as a 3-D homogeneous space, where o is treated as the homogeneous coordinate. The working space shown in Fig. 11b corresponds to the plane at the bottom of Fig. 11a. Note that the paraboloid touches the working space at its origin, and that the location of finite points (e.g., o , a , b , and c ) is given by their orthogonal projection onto the working space. From the outer product of two, three, and four finite points one builds, respectively, point pairs, circles, and spheres. So, the construction of k-spheres (for 0 ≤ k < d) is straightforward. It is achieved from the outer product of (k + 2) points. Fig. 11d shows a circle defined as C 3 = a ∧ b ∧ c . Note in Fig. 11c that such circle is a cross section of the paraboloid, orthogonally projected onto the working space. Supplementary Material B includes formulas for defining kspheres from their usual center-radius parameterization.

As commented in Section VII-B, improper points in the homogeneous model are characterized as points at infinity in such a way that each direction is one of these points. In the conformal model, however, the working space is “closed”. It means that ∞ is the unique point at infinity, with a well defined location that one can approach from any direction. So, ∞ is common to all flat subspaces, because they stretch to infinity. This way, straight lines and planes are built as the outer product of ∞ with, respectively, two and three finite points. For instance, the line passing through a to b in Fig. 11e is computed as L 3 = a ∧ b ∧ ∞ . It is important to note that all the equations for construction of flat subspaces in the conformal model are backward compatible with the ones from the homogeneous model, but including ∞ . In Fig. 11e, the cross section defined by L 3 is a parabola whose orthogonal projection in the working space is a straight line (Fig. 11f). Now, note how similar C 3 and L 3 are (Figs. 11c and 11e). Both are 3-blades which define cross sections in the paraboloid. In fact, both can be interpreted as circles, where L 3 is a circle with infinite radius. In order to be interpreted as a direction, a blade must have only directional properties and no locational aspects. The location of a blade is defined in terms of the assumed origin point o . Therefore, directions (also called free blades) are built as A k ∧ ∞ , where A k ⊂ (ee1 ∧ e 2 ∧ · · · ∧ e d ). This is the natural extension of blades, which are interpreted as directions in the homogeneous model, to the conformal model. The final type of conformal blade is the tangent subspace. As the name suggests, such primitives are tangent to something. In such a case, they encode the subspace tangent to rounds or flats at a given location. Therefore, tangent subspaces have a point-like interpretation, and also a direction information assigned to them. For a given round (or flat) X k passing thought the point p , the tangent subspace at the  k . The general equation for location of p is T k−1 = p  X  k ∞ )), where p  (A tangents subspaces is T k+1 = p ∧ (−p A k ⊂ (ee1 ∧ e 2 ∧ · · · ∧ e d ) defines the direction. In computer graphics applications, 3-D geometrical models are often represented by triangular meshes. Usually, for each mesh vertex, one stores not only its position, but also its tangent space composed of: the tangent, the binormal, and the normal vector. As pointed out by Dorst et al. [9], the use of tangent blades is an elegant alternative to represent vertices in a mesh, because they encode both the positional as the tangential information in a single primitive element. Also, tangent directions can be used to encode rays, which are usually represented by a point and a direction vector. So far, we had presented the geometrical interpretation of blades as directions, flats, rounds, or tangent subspaces. Now we show how to represent orthogonal transformations as versors (Section VI). We start with reflections. In the

e2 a c

e1

o

b

(a)

(b) e2 a c

e1

o

CX3\

(c)

b

(d) e2 a e1

o

b

(e)

LX3\

(f)

Figure 11. The conformal representation of a 2-D working space is shown on the left. In such a representation, the basis vectors e 1 , e 2 , and ∞ are seen as a homogeneous space having the basis vector o as homogeneous coordinate. The 2-D working space (on the right) corresponds to the plane at the bottom of the images on the left. Points on the paraboloid are interpreted as finite points in the working space (see o , a , b , and c in (a) and (b)). As depicted in (c) and (d), the circle defined by a , b , and c is computed as C 3 = a ∧ b ∧ c . In (e) and (f), the straight line that passes thought a and b is defined as L 3 = a ∧ b ∧ ∞ .

conformal model, reflection versors are computed as: m = M ∗d+1 ,

(48)

where M d+1 is the (d − 1)-flat or the (d − 1)-sphere that acts as a mirror. For instance, in a 3-D working space (where d = 3) one can use planes and spheres as mirrors. A translation over vector t = δ1 e 1 + δ2 e 2 + · · · + δd e d is represented by the rotor:   1 1 (49) T = exp − t ∞ = 1 − t ∞ , 2 2 where the exponential function is given by Eq. 40. Rotations around the origin are defined using Eq. 39, where B 2 is a 2-blade that resides in the Euclidean part of the conformal space. A general rigid body motion is constructed by performing a rotation around the origin, followed

by a translation. In linear algebra, the matrix encoding a motion is hard to interpolate. However, in geometric algebra one interpolates such motions in closed form by using the rotor logarithm (see [9], [10] for details). Finally, rotors encoding isotropic scaling relative to the origin are given by: γ  γ γ S = exp o ∧ ∞ = cosh + sinh o ∧ ∞ . (50) 2 2 2

models with degenerated signatures. In contrast to Gaigen 2, GATL allows the compile-time optimization of computations without the need of profiling and code re-generation. Also, GATL can be accessed using a MATLAB command line interface through a set of wrapper functions. The development of a geometric algebra library whose efficiency is compared or superior to linear algebra implementations like BLAS and LAPACK is an open problem.

Here, the scaling factor is exp(γ). So, to obtain a scaling by a factor of 5, one needs to use γ = log(5). We want to emphasize that pure versors (Eq. 39, 48, 49, and 50) can be combined into a single transformation by using the associativity of the geometric product (Section IV-B); and that versors are applied to any element by using the sandwiching construction presented in Eq. 37. Also, versors always transform blades in a consistent way (e.g., translations do not affect directions, while rotations and scalings do), because the semantic differences of geometric elements are encoded in the blade types.

IX. C ONCLUSION

VIII. G EOMETRIC A LGEBRA I MPLEMENTATIONS There are a few available implementations of geometric algebra. Ashdown [12] implemented a package for Maple. It enables the user to perform calculations in a geometric algebra of arbitrary dimensions with non-degenerated signatures. Unfortunately, it can not be used with the conformal model presented in Section VII-C, because such model has degenerated signature (note in Eq. 45 the presence of null vectors o and ∞ ). Leopardi and collaborators [13] implemented GluCat, a C++ library of template classes that model the universal Clifford algebras (or geometric algebra) over the field of real numbers, with arbitrary dimension and arbitrary non-degenerated signature. As Ashdown’s code, GluCat can not be used with the conformal model. CLUCalc [14] is a software tool developed by Perwass for 3-D visualizations and scientific calculations. It provides pleasant visualizations of geometric algebra structures and transformations. CLUCalc interprets a script language called CLUScript, which has been designed to make mathematical calculations very intuitive. Gaigen 2 [15] is a stand-alone application for generating geometric algebra libraries for a target programming language (e.g., C++, Java) from a high level specification code. Gaigen 2 was designed by Fontijne and, as pointed out in [16], efficient code is achieved, for a given application, after some profiling and code re-generation be performed. Gaigen 2 is reported in the literature as the most efficient implementations of geometric algebra. We have generated the images of this tutorial using MATLAB and our own geometric algebra library: Geometric Algebra Template Library (GATL). GATL is implemented using C++ meta-programming with template. In contrast to GluCat, GATL can handle geometric algebra

This tutorial focused on the fundamental concepts of geometric algebra. In such a mathematical framework, multivectors encode subspaces and orthogonal transformations. Also, geometric operations like spanning of subspaces and metric comparisons are all derived from the geometric product, leading to the development of elegant and generally applicable solutions for geometric problems. Geometric Algebra is a powerful tool for solving geometric problems in visual computing. For instance, Perwass [17] created a model of geometric algebra where 2-D conic sections are encoded by 5-blades. He used such model to develop an image processing algorithm for distinguishing between a number of different image structures, like: corners, crossings, y-junctions, t-junctions, lines, and line segments. Bayro-Corrochano and Flores [18] applied rotors (Section VI) on color edge detection. In computer vision, Faugeras and Mourrain [19] used the formalism of the Grassmann-Cayley algebra (the foundations of geometric algebra) to model point and line correspondences between n images. In their work, Faugeras and Mourrain have pointed out that Grassmann-Cayley algebra is the simplest way to make both geometric and algebraic statements in a very synthetic and effective way. Unfortunately, they revert equations to matrix-based notation in the final derivations. Bayro-Corrochano et al. [20] used the compact notation of geometric algebra to analyze the projective structure of n uncalibrated cameras, but without reverting equations to matrix-notation. Lasenby et al. [21] investigated the problem of estimating motion from a pair of images captured from a single camera undergoing translation and rotation. Regarding a monocular calibrated camera, Rosenhahn and Sommer [22] used the conformal model to define 2-D/3-D pose estimation of different corresponding entities, such as: 2-D point/3-D point, 2-D line/3-D point, 2-D line/3-D line, 2-D conic/3-D circle, and 2-D conic/3-D sphere. Hildenbrand et al. [23] applied the conformal model (Section VII-C) on standard tasks of computer animations and robotics, like inverse kinematics of a human-arm-like robots. Zaharia and Dorst [24] emphasize the elegance of the geometric algebra approach applied on modeling and visualization of 3-D polygonal mesh surfaces; Bayro-Corrochano and Sommer [25] discussed the use of geometric approach in object modeling and collision avoidance. Jourdan et al. [26]

performed automatic tessellation of quadric surfaces using Grassmann-Cayley algebra. Finally, Lasenby et al. [27] used geometric algebra for modeling higher dimensional fractals. We expect to see more applications of geometric algebra to visual computing in the near future. This should lead to elegant solutions that immediately generalize to subspaces of arbitrary dimensions. As such, the adoption of geometric algebra as a language to describe and solve problems should contribute to improve software development productivity and to reduce program errors. These are two important goals in software engineering.

[14] C. Perwass, “CLUCalc,” 2009. [Online]. Available: http: //www.clucalc.info/ [15] D. Fontijne, “Gaigen 2,” 2009. [Online]. Available: http: //staff.science.uva.nl/~fontijne/gaigen2.html [16] ——, “Gaigen 2: a geometric algebra implementation generator,” in Proc. of the 5th Int. Conf. Generative Progr. Component Eng. ACM Press, 2006, pp. 141–150. [17] C. B. U. Perwass, “Analysis of local image structure using intersections of conics,” Instituts für Informatik und Praktische Mathematik der Universität Kiel, Germany, Tech. Rep. Nr. 0403, July 2004.

ACKNOWLEDGMENT This work was sponsored by CNPq-Brazil, grants 142627/2007-0 and 305613/2007-3. R EFERENCES [1] J. D. Foley, A. van Dam, S. K. Feiner, and J. F. Hugues, Computer graphics principles and practice in C, 2nd ed. Addison-Wesley, 1995. [2] R. I. Hartley and A. Zisserman, Multiple view geometry in computer vision. Cambridge University Press, 2000. [3] R. N. Goldman, “Illicit expressions in vector algebra,” ACM Trans. Graph., vol. 4, no. 3, pp. 223–243, July 1985. [4] J. Stolfi, Oriented projective geometry. Professional, Inc., 1991.

Academic Press

[5] W. R. Hamilton, “On a new species of imaginary quantities connected with the theory of quaternions,” in Proc. of the Royal Irish Acad., vol. 2, 1844, pp. 424–434. [6] K. Shoemake, “Animating rotation with quaternion curves,” in Proc. of the 12th ACM SIGGRAPH, vol. 19, no. 3, July 1985, pp. 245–254. [7] D. Hestenes, New foundations for classical mechanics, 2nd ed. Reidel Publishing Company, 1999. [8] G. Sommer, Geometric computing with Clifford algebras: theoretical foundations and applications in computer vision and robotics. Springer-Verlag, 2001. [9] L. Dorst, D. Fontijine, and S. Mann, Geometric algebra for computer science: an object oriented approach to geometry. Morgan Kaufmann Publishers, 2007. [10] C. Perwass, Geometric algebra with applications in engineering. Springer Publishing Company, 2009. [11] C. Doran, J. Lasenby, L. Dorst, D. Hestenes, S. Mann, A. Naeve, and A. Rockwood, “Geometric algebra,” Course 53 at SIGGRAPH, 2001. [12] M. Ashdown, “GA package for Maple,” 2004. [Online]. Available: http://www.mrao.cam.ac.uk/~maja1/software/GA/ [13] P. C. Leopardi, “GluCat,” 2009. [Online]. Available: http://glucat.sourceforge.net/

[18] E. Bayro-Corrochano and S. Flores, Applications of geometric algebra in computer science and engineering. Birkhäuser, 2002, ch. Color edge detection using rotors, pp. 333–339. [19] O. Faugeras and B. Mourrain, “On the geometry and algebra of the point and line correspondences between N images,” in Proc. of the 5th Int. Conf. Comp. Vis., 1995, pp. 951–956. [20] E. Bayro-Corrochano, J. Lasenby, and G. Sommer, “Geometric algebra: a framework for computing point and line correspondences and projective structure using n-uncalibrated cameras,” in Proc. of the 13th Int. Conf. Pattern Recognit., 1996, pp. 334–338. [21] J. Lasenby, W. J. Fitzgerald, A. N. Lasenby, and C. J. L. Doran, “New geometric methods for computer vision: an application to structure and motion estimation,” Int. J. Comput. Vis., vol. 26, no. 3, pp. 191–213, Feb./March 1998. [22] B. Rosenhahn and G. Sommer, “Pose estimation in conformal geometric algebra part II: real-time pose estimation using extended feature concepts,” J. Math. Imaging Vis., vol. 22, no. 1, pp. 49–70, Jan. 2005. [23] D. Hildenbrand, E. Bayro-Corrochano, and J. Zamora, “Advanced geometric approach for graphics and visual guided robot object manipulation,” in Proc. of the Int. Conf. Robot. Autom., April 2005, pp. 4727–4732. [24] M. D. Zaharia and L. Dorst, “Modeling and visualization of 3D polygonal mesh surfaces using geometric algebra,” Comput. Graph., vol. 28, no. 4, pp. 519–526, Aug. 2004. [25] E. Bayro-Corrochano and G. Sommer, “Object modelling and collision avoidance using Clifford algebra,” Lect. Notes Comput. Sci., vol. 970, pp. 699–704, 1995. [26] F. Jourdan, G. Hégron, and P. Macé, “Automatic tessellation of quadric surfaces using Grassmann-Cayley algebra,” in Proc. Int. Conf. Comput. Vis. Graph., 2004, pp. 674–682. [27] J. Lasenby, R. Lasenby, A. Lasenby, and R. Wareham, “Higher dimensional fractals in geometric algebra,” Cambridge University Engineering Department, Tech. Rep. CUED/F-INFENG/TR.556, 2006.