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

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

252 THE VECTOR SPACE MODEL: THE ALGEBRA OF DIRECTIONS CHAPTER 10

bivectors is automatically measured in a plane perpendicular to the common

line of the two planes, and this plane and the angle are found as the bivector

angle log(V/U).

•

Vector and Bivector (2-Blade). When we have a unit bivector U and a unit vec-

tor v, we previously deﬁned the cosine of their angle through the contraction, as

in Figure 3.2(b). With the geometric product, we can proceed slightly differently,

deﬁning the full bivector angle. We try to determine a unit vector w in the plane U

and perpendicular to v, such that it can rotate v into the plane over an angle α as

, after which that rotated version of v and w together span U. The sketch of

Figure 10.2 shows that this mouthful is the algebraic demand

U = v e

This fully determines both w and α as aspects of the geometric product vU. It is sim-

plest to make a rotor out of the element vUby undualization to show its bivector:

vUI

= e

(π/2+α)

Taking the logarithm retrieves this bivector, which gives all parameters in one oper-

ation (though it is a pity not to get the actual bivector angle α wI

Angular relationships between three directions can also be deﬁned. They are much

more involved, because there are various standard ways of characterizing the parame-

ters of the spherical triangle, as depicted in Figure 10.3. These are its vertices (represented

Figure 10.2: The angle between a vector and a bivector (see text).

SECTION 10.2 ANGULAR RELATIONSHIPS 253

Figure 10.3: A spherical triangle and its characterizing parameters.

by three vectors), its sides (the angles between the vectors), its angles (the angles between

the bivectors containing its sides), and its altitudes (angle between vector and plane bivec-

tor). The various combinations of these quantities provide the laws of spherical geometry.

Geometric algebra again permits compact and computationally complete speciﬁcation of

the relationships.

The bivectors containing the vertices can be deﬁned through

b ≡ e

c ≡ e

a ≡ e

(10.10)

where

c are the unit vectors pointing at the vertices of the spherical triangle. The

weights A, B, C of the bivectors A, B, C are the edge lengths, measured as angles on the unit

sphere. The bivectors can be used to deﬁne the angles of the spherical triangle, through

A ≡−e

B ≡−e

C ≡−e

(10.11)

with

C now denoting the unit bivectors and I

the unit pseudoscalar of R

3,0

.The

weights of the vectors a, b, c thus deﬁned are the internal angles a, b, c of the spherical

triangle in Figure 10.3.

Multiplication of the unit vector expressions leads to

= 1,

254 THE VECTOR SPACE MODEL: THE ALGEBRA OF DIRECTIONS CHAPTER 10

and multiplication of the unit bivector expressions gives

= −1.

We have met the former in the guise of the multiplication of rotors in Section 7.3.4.

That should give you the conﬁdence that these expressions can indeed lead to the cor-

rect expressions for sine and cosine formulas, when split in their different grades. To

be speciﬁc, you can work out that the scalar part of the unit vector expression gives the

identity

cos C = cos A cos B + sin A sin B cos c

known as the cosine law for sides. The minus sign in the bivector expression leads to subtle

differences between the laws for the sides and for the angles. The scalar part is now

cos c = −cos a cos b + sin a sin b cos C.

which is known as the cosine law for angles.

For the full details on these and other laws of spherical trigonometry, we refer to Appendix

A of [29], highly recommended reading for anybody using spherical trigonometry. When

we get more versed in geometric algebra, we should probably not reduce the geometry

of these spherical triangles to the set of classical scalar relationships, learning instead to

compute directly with the rotors, vectors, and bivectors involved.

10.2.3 ROTATION GROUPS AND CRYSTALLOGRAPHY

A famous old puzzle in the mathematics of physics has been the classiﬁcation of crys-

tallographic groups: how many ways are there to make crystals in space? This is clearly

a geometric problem, and one would hope that geometric algebra could help out. The

vector space model, which provides the algebra of directions, can be used to solve a sub-

problem: how can one select the local symmetry of reﬂections and rotations to combine

well, in the sense that multiple reﬂections and rotations generate only a ﬁnite number of

elements? Such a set of operators is called a point group.

We use a unit vector a (with a

= 1) to denote a local plane of reﬂection. When used as

the normal vector of this plane, it generates the reﬂection of a vector x as −axa. This

transformation is insensitive to the sign of a. However, we can in principle distinguish the

planes a and −a by their orientation, and for accurate classiﬁcation of the point groups

we should use such oriented planes.

In a 3-D crystal, there are several symmetry planes for reﬂection. In the algebra, they

are each denoted by their vector versors. As the reﬂections combine as operators, new

symmetries appear as the products of these versors. To form a crystal, the combined set

of operators should form a ﬁnite set, so that a single point (atom) transformed by all

SECTION 10.2 ANGULAR RELATIONSHIPS 255

operators of the crystal is only equivalent to a ﬁnite number of other atoms around the

same point.

An even number of reﬂections generates a rotation, which is represented by the geometric

product of the corresponding normal vectors to give an operator like R ≡ ab. A six-fold

symmetry at the local point implies that applying the rotation operator x → R x



R six

times gives back the original x. Two possible conditions would satisfy this demand: R

= 1

and R

= −1. The latter is more speciﬁc, and uses the possibility of geometric algebra to

represent oriented rotations (see also Section 7.2.3). Since we need to reﬂect other oper-

ators (not just points), we use the −1. That gives the most accurate classiﬁcation of the

symmetries of the crystal. In general, R

= 1 for a p-fold symmetry.

To use a particular example in 2-D, suppose we have two generating vectors a and b

representing reﬂecting planes. These induce symmetries in the plane. To be a four-fold

symmetry, we have the three demands:

= 1, b

= 1, (ab)

= 1.

Clearly the generating vectors of this point group should be unit vectors with a relative

angle of π/4 to satisfy these demands. The rotation and reﬂection operators that are pos-

sible by the combinations of these vectors are listed in Table 10.1. There are 16 symmetry

operators, which together form the group 2H

of [30]. Each of these elements is a local

symmetry V

of the crystal. If you start with a single location vector x and generate all

elements x



x V

−1

, you ﬁnd 8 equivalent locations for atoms in the symmetry group

of this crystal. This is demonstrated in the programming example in Section 10.7.2.

The deﬁning relations on two vectors is sufﬁcient to determine all points groups in 2-D.

In 3-D, three generating vectors are required. If we have the demands in the form

(ab)

= (bc)

= (ca)

= 1,

Table10.1: The operators of the point group 2H

, with deﬁning relations a

= b

= (ab)

= 1.

The terms positive and negative are relative (but consistent), and based on the ordering of a

and b. These generate the symmetries of the crystal.

Positive Negative Positive Negative

Rotations Rotations Reﬂections Reﬂections

1 −1 a −a

ab −ab aba −aba

(ab)

(ba)

−bab bab

−ba ba −b b

256 THE VECTOR SPACE MODEL: THE ALGEBRA OF DIRECTIONS CHAPTER 10

then the rotors generated by this can be written as

(ab) = e



π/p

, (bc) = e



π/q

, (ca) = e



π/r

in which I is the pseudoscalar of the 3-D space, and a



and so on are the poles (unit axes)

of the rotations. Careful analysis performed by Coxeter (in a different, non-GA represen-

tation) reveals that a necessary condition for ﬁnitely representable solutions is:

> 1.

The integer solutions to this equation then generate the point groups in 3-D space. The

complete classiﬁcation may be found in [30].

Geometric algebra thus classiﬁes the groups by sets of vectors in relative positions; these

not only represent the reﬂection planes (which is how more classical solutions also char-

acterize the symmetry), but they can immediately be used to generate the operators that

perform all the operations in the group simply by the geometric product.

The extension to crystallographic groups imposes the additional condition that the point

groups should remain closed under translation. We cannot treat those with the same sim-

plicity within the vector space model. Our reference [30] develops the crystallographic

groups further, using the conformal model described in Chapter 13 to include transla-

tions with the same ease.

10.3 COMPUTING WITH 3-D ROTORS

The Euclidean space R

3,0

is a structural model of orientations in 3-D. In this model, rota-

tions are represented by rotors. Those are structure-preserving and easily interpolated.

We have already treated many of their structural properties when we introduced them in

the general setting of n-dimensional algebras. Even though they did not then necessarily

represent 3-D rotations, those provided a good and immediate illustration. We revisit and

extend that material to provide the practical tools for 3-D rotors: how to ﬁnd them from

frames, how to determine their bivectors, and how to interpolate them.

10.3.1 DETERMINING A ROTOR FROM A ROTATION PLANE

AND AN ANGLE

The most natural way to ﬁnd a rotor in any geometrical problem is to know the rota-

tion plane represented by the 2-blade I and the desired rotation angle ␾ (measured to be

positive in the orientation of I). Their product gives the geometric quantity I␾, which

does not depend on the choice of the orientation of the plane. We called this the bivector

angle in Section 7.2.2 (that is actually a misnomer in more than three dimensions, where

SECTION 10.3 COMPUTING WITH 3-D ROTORS 257

it should then more properly be called a 2-blade angle, but that is rather a mouthful). The

rotor in terms of the bivector angle is

R = e

−I␾/2

We have seen how this is essentially a quaternion in Section 7.3.5, and in the programming

exercise of Sections 7.10.3 and 7.10.4 showed how to convert between rotation matrices,

quaternions, and rotors. In Section 7.3.4, we gave a pleasant visualization of the

composition of rotations as the addition of half-angle arcs on a unit sphere.

If you have been able to determine the cosine c of the rotation angle in some way

(for instance, through a scalar product), you can use this fairly directly to determine the

rotor as

R =



(1 + c)/2 −



(1 − c)/2 I. (10.12)

The square roots are expressions of the cosine and sine of the half angle in terms of the

cosine. When you take the standard positive values of the square roots, the orientation of

I denotes the direction of rotation.

10.3.2 DETERMINING A ROTOR FROM A FRAME

ROTATION IN 3-D

A rotor can be found when you know how enough vectors rotate. In any number of dimen-

sions, there is a cheap way to make the unique rotor turning unit vector a to unit vector b

in their common plane. This of course involves the half angle between them, which may

seem expensive to compute. However, a unit vector in that half-angle direction is easily

found by adding the vectors a and b and normalizing. The rotor is then the geometric

ratio of this vector with a. That gives

R =

1 + ba



2(1+ a · b)

(10.13)

You should realize that this is one of many rotations that can turn a into b. It is the simplest

in the sense that it only involves their common plane. The formula is unstable when a and

b are close to opposite; when they are truly opposite, there is no unique rotation plane to

turn one into the other.

In general, you would need to know the image of several vectors before you can determine

the exact relative rotor. In 3-D Euclidean space, there is a compact formula to retrieve the

rotor from the three vectors of a frame {e

} (which does not even have to be orthogonal)

and their images {f

}.Itis

R ∼ 1 +



i=1

which needs to be properly scaled to become a rotor (for which of course R



R = 1).

Note that it uses the reciprocal frame of the vectors f

(see Section 3.8). It does not give

258 THE VECTOR SPACE MODEL: THE ALGEBRA OF DIRECTIONS CHAPTER 10

the correct result for rotations over π (returning zero instead), and is unstable near

rotations close to π. These consequences were already explored in the programming

exercise of Section 7.10.3.

We give a rather advanced derivation of this formula, invoking vector derivatives. The

only reason for doing this is to show the general n-dimensional pattern, should you ever

need it, before we home in on the 3-D-only case.



R e



R e

= R ∂

(



R a)

= n − 2R



0



R

+ 2



R

+ 4



R

+ ···



= 3 − 4 R



R

= 3 − 4 R (



R −R

)

= 4 R

R − 1

and the result follows, since this shows that R is proportional to 1 +



in 3-D.

10.3.3 THE LOGARITHM OF A 3-D ROTOR

As we saw in Section 7.4.3, general rotors are the exponentials of bivectors. In the 3-D

vector space model of

3,0

, it is fairly straightforward to retrieve that bivector from the

rotor. In such a 3-D space, the rotors are exponentials of 2-blades, since all bivectors are

2-blades in 3-D.

An explicit principal logarithm (as explained in Section 7.4.4) is easy to give for a rotor

R = exp(−I␾/2). It can be determined by writing out the expression for the rotor in its

grades and reassembling those using standard trigonometric functions on the grade parts:

−I␾/2 = log(R)

= log



exp(−I␾/2)



= log



cos(␾/2) − I sin(␾/2)



R

R



atan



R



R





. (10.14)

There are some special cases that should be borne in mind. Obviously, when the scalar

R

equals zero, the division in the atan is ill-deﬁned, even though there is still a well-

deﬁned rotation plane, and therefore a well-deﬁned logarithm. In an implementation,

you should use the

atan2 function to provide numerical stability and additionally make

(10.14) valid for the second and third quadrant. In two other cases, namely when the

rotor equals 1 or −1, the rotation plane is ambiguous. For the identity rotation, the small

angle makes the behavior around R = 1 still numerically stable, and the logarithm is

virtually zero. But the ambiguity at R = −1 cannot be resolved without making arbitrary

choices. These cases are recognizable in the pseudocode of Figure 13.5 (where the rotation

logarithm is presented as part of the logarithm of the general rigid body motion).

SECTION 10.3 COMPUTING WITH 3-D ROTORS 259

This formula (10.14) can be extended from a rotor to a nonunit versor; there is then an

additional term of log(R).

10.3.4 ROTATION INTERPOLATION

When you want to interpolate between two known orientations, you can do this by divid-

ing the rotor between the extreme poses in equal amounts. This requires being able to take

the n th root of a rotor through using its logarithm. For the Euclidean rotors in the 3-D

vector space model of

3,0

, this can be done explicitly.

Let us suppose we know the initial and ﬁnal rotor in the interpolation sequence as R

and

, respectively. Then the total rotation that needs to be performed is characterized by the

rotor R ≡ R

. To perform this total rotation in n steps, one needs to apply the rotor

r ≡ R

1/n

, n times. In the geometric algebra formulation, this is simply

R = R

= e

−I␾/2

⇒ r = e

−I␾/(2n)

= e

log(R)/n

This formula requires the bivector corresponding to the rotor, which we derived in (10.14)

as its principal logarithm. The rotor that needs to be applied to the original element after

k applications is r

, and of course we should have r

= R

. That gives the simple

program of programming exercise 10.7.1. See also Figure 10.4.

An alternative is to give the resulting rotor not in multiplicative form, but to use the

trigonometric functions to give a closed expression of the interpolation going on between

[X]

Figure 10.4: The interpolation of rotations illustrated on a bivector X. The poses R

[X] and

[X] are interpolated by performing the rotor R

in eight equal steps.

260 THE VECTOR SPACE MODEL: THE ALGEBRA OF DIRECTIONS CHAPTER 10

and R

. This is the way it is found in the quaternion literature, where this is known as

slerp interpolation (for spherical-linear interpolation). To derive the formula, note that the

requirement R

= cos(␾/2) − I sin( ␾/2) gives IR



− R

cos(␾/2)



/ sin(␾/2).

Then the linear interpolation is achieved by rotation over a fraction λ␾/2 of the angle,

from R

towards R

. This is the rotor e

−λI␾/2

, which may be expressed as



cos(λ␾/2) − I sin(λ␾/2)



sin(␾/2) cos (λ␾/2) − R

cos(␾/2) sin (λ␾/2) + R

sin(λ␾/2)

sin(␾/2)

sin((1 − λ)␾/2)

sin(␾/2)

sin(λ␾/2)

sin(␾/2)

(10.15)

This is the linear interpolation formula for rotations (see [56]). It is valid in n-dimensional

space.

In our software, we prefer not to use this explicit form, but instead the structurally simpler

formulation in terms of the incremental rotor r. In the bivector formulation, one can

easily design more sophisticated interpolation algorithms that interpolate between several

rotors, such as bivector splines.

10.4 APPLICATION: ESTIMATION IN THE VECTOR

SPACE MODEL

For applications in synthetically generated computer graphics, the generative techniques

for rotations sufﬁce, but in the analytical ﬁelds of computer vision and robotics you need

to determine rotations based on noisy data.

10.4.1 NOISY ROTOR ESTIMATION

If your data is noisy, you should of course not establish the rotation based on three frame

vectors as in Section 10.3.2, but instead use a rotor estimation technique. You could try to

use the frame estimation on triplets of rotors and average their bivectors, but it is better

to do a proper estimation minimizing a well-deﬁned cost function.

We encountered such a technique in Section 8.7.2, when we used the rotation estimation

problem to give an example of multivector differentiation. Its outcome is a useful result

that you can easily implement even without understanding the details of its derivation.

10.4.2 EXTERNAL CAMERA CALIBRATION

When you have a set of cameras of unknown relative positions and attitudes observ-

ing one scene, you can only integrate their views if you have calibrated the setup (i.e.,

estimated the parameters of their relative geometry). You can collect the data for such

SECTION 10.4 APPLICATION: ESTIMATION IN THE VECTOR SPACE MODEL 261

a calibration by moving a spherical marker around in the scene and synchronizing the

cameras to record their various views of it at the corresponding times, as illustrated

in Figure 10.7. The observed data is of course inherently noisy, and you need to do

some processing to determine the “best” estimate for their relative poses. We describe

an algorithm for this, taken from [38], culminating in the programming exercise of

Section 10.7.3. Our source assumes that the cameras have been calibrated internally.

This involves a determination of the parameters of their optics and internal geome-

try, so that we can interpret a pixel in an image as coming from a well-determined

spatial direction relative to its optical axis. Geometrically, it turns the camera into a

measurement instrument for spatial directions.

LetusconsiderM + 1 cameras. We arbitrarily take one of them as our reference camera

number 0, and characterize the location of the optical center of camera j relative to the

center of camera 0 by the translation vector t

and its orientation by the rotor R

.We

mostly follow the notation of [38] for easy reference, which uses bold capitals for vectors

in the world, and lowercase for vectors within the cameras. We will simplify the situation

by assuming that the marker is visible in all cameras at all times (our reference deals with

occlusions; this is not hard but leads to extra administration).

The marker is shown N times at different locations X

in the real world. Relative to camera

j, it is seen at the location X

given implicitly by

= t

+ R



However, all that camera j can do is register that it sees the center of the marker in its

image, and (using the internal calibration) know that it should be somewhere along the

rayindirectionx

from its optical center. The scaling factor along this ray is σ

;ifwe

would know its value, the camera would be a 3-D sensor that would measure σ

= X

and then R

and t

could be used to compute the true location of the measured points.

But the only data we have are the x

for all the cameras. All other parameters must be

estimated. This is the external calibration problem.

The estimation of all parameters is done well when the reconstructed 3-D points are not

too different from their actual locations. When we measure this deviation as the sum of

the squared differences, it implies that we want to minimize the scalar quantity

Γ=



j=1



i=1



− t

− R





Now partial differentiation with respect to the various parameters can be used to derive

partial solutions, assuming other quantities are known. This employs the geometric dif-

ferentiation techniques from Chapter 8. The results are simple to interpret geometrically

and are the estimators you would probably have proposed even without being able to

show that they are the optimal solutions.