Dorst L., Fontijne D., Mann S. Geometric Algebra for Computer Science. An Object Oriented Approach to Geometry

Подождите немного. Документ загружается.

506 IMPLEMENTATION ISSUES CHAPTER 18

18.2 WHO SHOULD READ WHAT

The implementation description is quite detailed and not everyone will want to read all

of it. We envision three types of audiences:

•

If you are new to geometr ic algebra, the ﬁrst two parts of the book may leave you

wondering how this is ever going to work in an actual implementation (with those

weird metrics, as a graded algebra, with so many diverse products, etc.). The fol-

lowing chapters should remove the fuzziness and magic. First, you will want to read

Chapter 20 on linear operations. This should give you the comfortable feeling that

geometric algebr a is fully consistent, since it has the structure of the linear algebra

of selected matrices. After that, you may want to read the ﬁrst part of Chapter 19 on

basis blades as it provides another viewpoint for understanding the basic products.

Then give a cursory look to the rest of the chapters, but be sure to read the bench-

marks that compare geometric algebra to traditional methods in Section 22.5.

•

If you are considering using geometric algebra in one of your programs, you might

want some background infor mation that helps you pick a particular implementa-

tion, especially the appropriate representation of the elements of geometry. In that

case we recommend skipping Chapter 19 on basis blades, but reading all of the linear

Chapter 20. Then give a cursory look at the nonlinear algorithms and efﬁciency

(Chapter 21 and 22). Finally, read the ray-tracing Chapter 23 in detail, as it will

show you actual geomet ric algebra code that might be similar to what you’re going

to write yourself, regardless of what implementation you pick.

•

If you are a hard-core coder and want to write your own implementation, or simply

want to know all the details, we suggest you read all chapters in full detail and the

proper order.

18.3 ALTERNATIVE IMPLEMENTATION APPROACHES

In the coming chapters, we will concentrate on one speciﬁc implementation approach,

where multivectors are represented as a sum of basis blades. It is what we found worked

best for low-dimensional geometric algebras useful to computer science. It is also the most

commonly used approach by others. But before we explore the detailed consequences, we

use the remainder of this chapter to mention some different approaches, just to widen

your ﬁeld of view. Some of these methods affect only the middle (linear) implementation

level, while others are so radically different that they affect every level.

18.3.1 ISOMORPHIC MATRIX ALGEBRAS

The linearity and associativ ity of geometric algebra suggests representing the whole

algebra with its geometric product by a single matrix algebra: all elements become

SECTION 18.3 ALTERNATIVE IMPLEMENTATION APPROACHES 507

represented as 2

× 2

matrices, and the geometric product is represented as the matrix

product. This is structurally pretty and quite well known among mathematicians as a

desirable representation technique, but it is computationally very expensive. It is in fact a

representation of Clifford algebra rather than geometric algebra (see Section 7.7.2), and

as such rather wasteful: in a true geometric usage of properly constructed elements such

as blades and versors, the corresponding matrices tend to be sparse, but this implemen-

tation does not make use of that property. Also, it cannot implement the contraction

and outer product as easily as the geometric product.

We will use this idea brieﬂy to implement inversion of general multivectors in Section 21.2.

18.3.2 IRREDUCIBLE MATRIX IMPLEMENTATIONS

The 2

× 2

matrices do not form the smallest matrix representation of the algebra, if we

follow common practice in mathematics and allow the ﬁeld of the linear mappings they

represent to include not only

R (the real numbers), but also C (the complex numbers)

and

H (the quaternions).

Such linear matrix representations have long been known for Clifford algebras of spaces

p,q

with arbitrary signature (p, q) (which means p + q spatial dimensions, of which p

basis vectors have a positive square, and q have a negative square; see Appendix A). Again,

each element of the algebra is represented as a matrix, and the matr ix product of two

such representatives is precisely the representation of the element that is their geometric

product. We repeat part of the table of such representations in Table 18.1 (see [49]), and

offer some structural exercises at the end of this chapter to familiarize yourself with them.

The structural advantage of such representations is that the geometric product becomes

a simple matrix multiply (although the elements of the matrix may be reals, complex

numbers, or quaternions). However, this is mostly a mathematical curiosity that will not

Table18.1: Matrix representations of Clifford algebras of signatures (p,q). Notation: R(n) are

n×n real matrices,

C(n) are n×n complex-valued matrices, and H(n) are n×n quaternion-valued

matrices. The notation

R(n) is used for ordered pairs of n × n real matrices (which you may

think of as a block-diagonal 2n× 2n matrix containing two real n × n real matrices on its diagonal

and zeros elsewhere), and similarly for the other number systems.

q =0 q =1 q =2 q =3

p =0 R(1) C(1) H(1)

H(1)

p =1

R(1)

R(2) C(2) H(2)

p =2 R(2)

R(2) R(4) C(4)

p =3 C(2) R(4)

R(4)

R(8)

p =4 H(2) C(4) R(8) C(8)

508 IMPLEMENTATION ISSUES CHAPTER 18

save processing time in practice, since approximately the same number of operations

have to be performed as in the basis-of-blades representation that is our reference

implementation. (Only if special hardware were present to handle complex number

and quaternion multiplications more efﬁciently than arbitrary multiplications, could

these methods become more efﬁcient than the multiplication of real valued matrices.)

Moreover, the matr ix representations have the disturbing disadvantage that they only

work for the geometric product. As long as they are used for the Clifford algebras for which

they were developed, this is not a problem—but it makes this representation cumbersome

for geometric algebra in general, where derived products such as the outer product and

contractions are also important. Tr ue, a similar representation of outer products can be

easily established, but one would have to switch between representations to perform one

product or the other, as both may be needed in a single application. And unfortunately,

since the contraction inner products are not associative, they are not isomorphic to matrix

algebras, so they cannot be implemented in the same framework. (Actually, this nonasso-

ciativity of the contraction will require special measures in any implementation scheme,

including ours.)

18.3.3 FACTORED REPRESENTATIONS

The mostly multiplicative structure of geometric algebra (with versors as geometric

products of vectors, and blades as outer products of vectors) has recently been explored

by the authors to produce factorized representations. This seems a viable implementation

method for high-dimensional algebras, but more research is needed.

A k-blade can be stored as a list of k vectors. The outer product of these vectors is the value

of the blade. Likewise, a versor can be stored as a list of vectors whose geometric product

is the value of the versor. The storage requirement of blades and versors becomes O(n

)

compared to O(2

) for the basis-of-blades method. Multivectors that are not blades or

versors can only be represented as a sum of multiple blades or versors, and thus require

more storage.

This multiplicative representation is radically different from the usual additive represen-

tation. As a consequence, the implementation of products and operations is also very dif-

ferent. What may be a difﬁcult problem in one implementation approach can be trivial in

the other, and v ice versa. For example, addition is simple in the basis of blades approach,

but one of the hardest problems in a factored representation. The

meet and join are

trivial with factored representation (given that you already have a linear algebra library

like LAPACK), while they lead to a rather involved algorithm in the sum-of-basis-blades

approach.

We do not discuss this approach further because it is best suited for high-dimensional

(n>10) geometric algebras, which are of less interest in this book. The details will appear

in Daniel Fontijne’s Ph.D. thesis.

SECTION 18.4 STRUCTURAL EXERCISES 509

18.4 STRUCTURAL EXERCISES

1. In a 2-D Euclidean geometric algebra on an orthonormal basis {e

, e

}, the general

element X can be written as:

X = x

+ x

According to Table 18.1, this algebra should be representable as the mat rix algebra

R(2) (i.e., we should be able to ﬁnd a 2 × 2 matrix representing the general element

so that the geometric product of elements is represented as the matrix product).

Show that the following works:

[[ x]] =



+ x

− x

+ x

− x



This is not unique; some permutations of the same principles work as well, but not

all. Why must the scalar part always be on the diagonal?

2. For a 3-D Euclidean vector space, show that you can generate a matrix algebra from

the following representation of a vector x = x e

+ y e

+ z e

[[ x]] =

⎡

⎢

⎣

⎡

⎢

⎣

z 0 x −y

0 zy x

xy−z 0

−yx 0 −z

⎤

⎥

⎦

⎤

⎥

⎦

Compute the representation of a general element of this algebra. (This represen-

tation is not in Table 18.1 because it is not the smallest matrix algebra; see the

following exercise).

3. For a 3-D Euclidean vector space, the matrix representation

C(2) from Table 18.1

may be generated by:

[[ x]] =



zx+ iy

x − iy −z



where i is the complex imaginary (so i

= −1). Verify this and compute the repre-

sentation of a general element of this algebra.

BASIS BLADES AND

OPERATIONS

All geometric algebra implementations that represent multivectors as a weighted sum of

basis blades will, at some point, have to compute products of basis blades. Among the var-

ious products, the capability to compute the geometric product (also for non-Euclidean

metrics) is clearly the minimum requirement, as the outer product and various ﬂavors of

the inner products can be derived from it using the (anti-)symmetry or grade selection

techniques of Chapter 6.

This chapter descr ibes how to implement this geometric product of basis blades through

a convenient representation of basis blades. The algorithms described are quite simple,

yet intricate enough that we need (pseudo)code to describe them precisely. Instead

of introducing our own pseudocode language, we opted to simply use Java. The code

in this chapter therefore corresponds well to our reference implementation at our web

site http://www.geometricalgebra.net, although it was sometimes slightly polished for

presentation.

Conceptually, our basis blades representation uses a list of booleans and applies Boolean

logic to the elements of that list. But in practice, the list of booleans is implemented

using bitmaps and the bitwise boolean operators. To correspond more directly to the

actual implementation, that is how we present the representation and its operators in

this chapter.

511

512 BASIS BLADES AND OPERATIONS CHAPTER 19

19.1 REPRESENTING UNIT BASIS BLADES

WITH BITMAPS

When we want to represent the geometric algebra of an n-dimensional space, we start

withasetofn independent basis vectors {e

}

i=1

. They need not be orthonormal or even

orthogonal; the geometric product of the space will essentially deﬁne their relevant geo-

metrical relationships. A basis blade is a nonzero outer product of a number of these basis

vectors (or a unit scalar). Hence, a basis blade either contains a speciﬁc basis vector or not.

Thus, each basis blade in an algebra can be represented by a list of boolean values, where

each boolean indicates the presence or absence of a basis vector. T he list of booleans is nat-

urally implemented using a bitmap, which is why we call this the bitmap representation

of basis blades.

The bits in the bitmap are assigned in the logical order: bit 0 stands for e

, bit 1 stands

for e

, bit 2 stands for e

, and so on. Table 19.1 shows how this works out. The unit scalar

does not contain any basis vectors, so it is represented by 0

(in this chapter, the postﬁx,

subscript b indicates a binary number). A blade such as e

∧e

contains e

and e

,soitis

represented by 001

+ 100

= 101

Visually, the bits are stored in reversed order relative to how we would write them as basis

blades. In a basis blade, e

is the left-most basis vector, while in bitmap representation

it is the right-most bit. While this is slightly inconvenient for reasoning about them, it

is more consistent for implementation. If an extra dimension is required, the next most

signiﬁcant bit in the bitmap is used and previous results are automatically absorbed.

Table 19.1: The bitmap representation of basis blades. Note how the representation scales

up with the dimension of the algebra.

Basis Blade Bitmap Representation

1 0

∧ e

100

∧ e

101

∧ e

110

∧ e

111

1000

... ...

SECTION 19.2 THE OUTER PRODUCT OF BASIS BLADES 513

Table 19.2: Bitwise boolean operators used in Java code examples.

Symbol Operation

& Bitwise boolean “and”

∧ Bitwise boolean “exclusive or”

>>> Unsigned bitwise “shift right”

In our reference implementation, we use 32-bit integers as our bitmaps. This allows us to

represent all basis blades of the geometric algebra of a 32-D space. This should be more

than enough for practical usage in computer science applications. Further, this basis-of-

blades method of implementing geometric algebra is only practical up to around 10-D

spaces, since the number of basis blades required to represent an arbitrary multivector

grows quickly, as 2

. The products have to process combinations of all basis blades from

both arguments, so you typically will run out of processing time before you run out of

memor y. But should you ever need more than 32 dimensions, extending the range of the

bitmaps is straightforward.

The reference implementation stores basis blades in a class named

BasisBlade. The class

has only two member variables,

bitmap (an integer) and scale (a ﬂoating point number).

scale allows the basis blade to have an arbitrary sign or scale; this is what we called the

weight in Part I. The name

ScaledBasisBlade might have been a more accurate class

description, but we found that a bit too verbose.

In Table 19.2, we show the symbols used for the bitwise boolean operations. These sym-

bols are identical to the ones Java uses. Note that in code listings, ∧ may appear as ^.

19.2 THE OUTER PRODUCT OF BASIS BLADES

We now show how to compute the products on the basis blades using the bitmap repre-

sentation. To compute the outer product of two basis blades, we ﬁrst check whether they

are dependent (i.e., whether the two blades have a common basis vector factor). If they do,

they have a bit in common, so dependence of blades is checked simply with a binary and.

If the blades are dependent, then the outcome of the outer product is 0 and the algorithm

described below does not need to be executed.

If the blades are independent, computing their outer product is done (up to a scale factor)

by taking the bitwise exclusive or of the bitmaps. For example,

∧ e

) ∧ e

= e

∧ e

514 BASIS BLADES AND OPERATIONS CHAPTER 19

is equivalent to

110

∧ 001

= 111

The hardest part in implementing the outer product is computing the correct sign for the

result. The bitmap representation assumes that the basis vectors are in canonical order

(e.g., using just the bitmap you cannot represent e

∧e

, only e

∧e

). We should employ

the

scale member variable to represent e

∧ e

as −1.0 ∗ (e

∧ e

When you compute the sign of the outer product result by hand, you count how many

basis vectors have to swap positions to get the result into the canonical order. Each swap

of position causes a ﬂip of sign. In our example, all that is required is that e

and e

swap

positions, so we get a negative sign.

To compute this sign algorithmically, we have to compute, for each 1-bit in the ﬁrst

operand, the number of less signiﬁcant 1-bits in the second operand. This is implemented

somewhat like a convolution, in that the bitmaps are slid over each other bit by bit. For

each position, the number of matching 1-bits counted. Figure 19.1 shows the implementa-

tion of this idea (the

bitCount() function counts the number of nonzero bits in a word).

The implementation of the outer product itself is in Figure 19.2; this function actually also

implements the geometric product, as described in the next section.

// Arguments ’a’ and ’b’ are both bitmaps representing

// basis blades.

double canonicalReorderingSign(int a, int b) {

// Count the number of basis vector swaps required to

// get ’a’ and ’b’ into canonical order.

a = a >>> 1;

int sum=0;

while (a != 0) {

// the function bitCount() counts the number of

// 1-bits in the argument

sum = sum + subspace.util.Bits.bitCount(a & b);

a = a >>> 1;

}

// even number of swaps — > return 1

// odd number of swaps — > return — 1

return ((sum & 1) == 0) ? 1.0 : — 1.0;

}

Figure 19.1: This function computes the sign change due to the reordering of two basis

blades into canonical order.

SECTION 19.3 THE GEOMETRIC PRODUCT OF BASIS BLADES IN AN ORTHOGONAL METRIC 515

BasisBlade gp_op(BasisBlade a, BasisBlade b, boolean outer) {

// if outer product: check for independence

if (outer && ((a.bitmap & b.bitmap) != 0))

return new BasisBlade(0.0);

// compute the bitmap:

int bitmap = a.bitmap ^ b.bitmap;

// compute the sign change due to reordering:

double sign = canonicalReorderingSign(a.bitmap, b.bitmap);

// return result:

return new BasisBlade(bitmap, sign * a.scale * b.scale);

}

Figure 19.2: This function will compute either the geometric product (in Euclidean metric

only) or the outer product, depending on the value of argument outer.

19.3 THE GEOMETRIC PRODUCT OF BASIS BLADES IN

AN ORTHOGONAL METRIC

The implementation of the geometric product is similar to that of the outer product as

long as we stay in an orthogonal metric with respect to the orthonormal basis {e

}

i=1

. Such

ametrichase

· e

= 0 for i = j, and e

· e

= m

; therefore, its metric matrix is diagonal:

· e

= m

(19.1)

where δ

is the Kronecker delta function, returning 1 when i equals j, and zero otherwise.

A particular example is the Euclidean metric, for which m

= 1 for all i, so that the metric

matrix is the identity matrix. In general, the m

can have any real value, but −1, 0, and 1

will be the most prevalent in many applications.

There are now two cases:

•

For blades consisting of different orthogonal factors, we have the usual equivalence

of the outer product and the geometric product:

∧ e

∧···∧e

= e

···e

•

However, when two factors are in common between the multiplicands, the outer

product produces a zero result, but the geometric product does not. Instead, it

516 BASIS BLADES AND OPERATIONS CHAPTER 19

annihilates the dependent-basis vectors, effectively replacing them by metric

factors. For example,

∧ e

)(e

∧ e

) = m

∧ e

In the bitmap representation, these two cases are easily merged. Both results may be com-

puted as the bitw ise exclusive or operation (i.e., 011

∧ 110

= 101

). In this sense, the

geometric product acts as a “spatial exclusive or” on basic blades.

As for the outer product, the result of the geometric product “exclusive or” should be

given the correct sign to distinguish between the outcomes e

∧ e

and e

∧ e

. We there-

fore also have to perform reordering techniques to establish the sign of the result on the

standard basis. Those sign changes are identical to those for the outer product: we count

the number of basis vector swaps required to get the result into canonical order and use

this to determine the sign.

In a Euclidean metric, all diagonal metric factors are 1, so this bit pattern and sign compu-

tation is all there is to the geometric product. The function that computes the geometric

product in this simple Euclidean case is therefore virtually the same as the function for

computing the outer product (see Figure 19.2), the only difference being the dependence

check required for the outer product. That close and convenient similarity is due to the

use of the orthonormal basis in the representation.

When the metric is not Euclidean but still diagonal, we need to incorporate the metric

coefﬁcients of the annihilated basis vectors into the scale of the resulting blade. Which

basis vectors are annihilated is determined with a bitwise and. We do not show the result-

ing code, which simply invokes the Euclidean code with this extension; you can consult

the reference implementation on our web site for the details.

19.4 THE GEOMETRIC PRODUCT OF BASIS BLADES IN

NONORTHOGONAL METRICS

In practical geometric algebra we naturally encounter nonorthonormal bases. For exam-

ple, in the conformal model we may want to represent o and ∞ directly as basis

vectors. This leads to a nondiagonal multiplication table (and therefore nondiagonal

metric matrix):

o e

∞

o 0 000−1

0 100 0

0 010 0

0 001 0

∞

−1000 0