Schneider P., Eberly D.H. Geometric Tools for Computer Graphics

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

864 Appendix A Numerical Methods

Table A.2 Comparison of operation counts for converting between representations of rotations.

Conversion A M D F C

Axis-angle to matrix 13 15 2

Matrix to axis-angle 8 7 1 2

Axis-angle to quaternion 1 5 1 2

Quaternion to axis-angle 4 2

Quaternion to matrix 12 12

Matrix to quaternion (τ>0) 6 5 1 1 1

Matrix to quaternion (τ ≤0) 6 5 1 1 3

x = λ(r

− r

), y = λ(r

− r

), and z = λ(r

− r

) require 3 additions and 3

multiplications. The total cost is 6A + 5M + 1D + 1F +1C.

If τ ≤ 0, the maximum of the diagonal entries of the rotation matrix must be

found. This requires two comparisons; call this cost 2C. For the sake of argument,

suppose that r

is the maximum. The calculation x =

√

− r

+ 1/2re-

quires 3 additions, 1 multiplication, and 1 function call. The expression λ = 1/(4x)

requires 1multiplication and 1 division. The terms w =λ(r

−r

), y =λ(r

and z = λ(r

) require 3 additions and 3 multiplications. The total cost is 6A +

5M +1D +1F + 3C.

Table A.2 is a summary of the costs of converting among the various rotation

representations.

Transformation Time

The transformation of v by a rotation matrix is the product u = Rv and requires 9

multiplications and 6 additions for a total of 15 operations.

If v = (v

, v

) and if v = v

i + v

j + v

k is the corresponding quaternion

with zero w-component, then the rotate vector u = (u

, u

) is computed as u =

i +u

j +u

k =qvq

∗

. Applying the general formula for quaternion multiplication

directly, the product p = qv requires 16 multiplications and 12 additions. The prod-

uct pq

∗

also uses the same number of operations. The total operation count is 56.

However, since v has no w-term, p only requires 12 multiplications and 8 additions—

one term is theoretically zero, so no need to compute it. We also know that u has

no w-term, so the product pq

∗

only requires 12 multiplications and 9 additions.

Using these optimizations, the total operation count is 41. Observe that conversion

from quaternion q to rotation matrix R requires 12 multiplications and 12 additions.

Transforming v by R takes 15 operations. Therefore, the process of converting to ro-

tation and multiplying uses 39 operations, 2 less than calculating qvq

∗

. Purists who

A.4 Representations of 3D Rotations 865

Table A.3

Comparison of operation counts for transforming one vector.

Representation A M Comments

Rotation matrix 6 9

Axis-angle 12 18

Quaternion 24 32 Using generic quaternion multiplies

Quaternion 17 24 Using specialized quaternion multiplies

Quaternion 18 21 Convert to matrix, then multiply

implement quaternion libraries and only use quaternions will sadly lose a lot of cycles

when transforming large sets of vertices.

The formula for transforming v using an axis-angle pair is

Rv =v +(sin θ)ˆa ×v +(1 − cos θ)ˆa ×( ˆa ×v)

As indicated earlier, sin θ and 1 − cos θ should be precomputed and stored in addi-

tion to the axis and angle, a total of 6 ﬂoats. The cross product ˆa ×ˆv uses 6 multi-

plications and 3 additions. So does ˆa × ( ˆa ×ˆv), assuming the cross product in the

parentheses was computed ﬁrst and stored in temporary memory. Multiplying the

cross products by a scalar requires 6 multiplications. Adding three vectors requires 6

additions. Therefore, we need to use 18 multiplications and 12 additions, for a total

of 30 operations.

Therefore, the rotational formulation yields the fastest transforming. The quater-

nion formulation yields the slowest transforming for a single vector. But keep in mind

that a batch transform of n vectors requires converting the quaternion to a rotation

matrix only once at a cost of 24 operations. The total operations for transforming

by quaternion are 24 + 15n. The axis-angle formulation uses 30n, so the quaternion

transformation is faster for two or more vectors. Table A.3 is a summary of the opera-

tion counts for transforming a single vector. Table A.4 is a summary of the operation

counts for transforming n vectors.

Composition

The product of two rotation matrices requires 27 multiplications and 18 additions,

for a total cost of 18A + 27M.

The product of two quaternions requires 16 multiplications and 12 additions, for

a total cost of 12A + 16M, clearly outperforming matrix multiplication. Moreover,

renormalizing a quaternion to adjust for ﬂoating-point errors is cheaper than renor-

malizing a rotation matrix using Gram-Schmidt orthonormalization.

866 Appendix A Numerical Methods

Table A.4 Comparison of operation counts for transforming n vectors.

Representation A M Comments

Rotation matrix 6n 9n

Axis-angle 12n 18n

Quaternion 24n 32n Using generic quaternion multiplies

Quaternion 17n 24n Using specialized quaternion multiplies

Quaternion 12 + 6n 12 +9n Convert to matrix, then multiply

Table A.5 Comparison of operation counts for composition.

Representation A M D F

Rotation matrix 18 27

Quaternion 12 16

Axis-angle (convert to matrix) 52 64 1 6

Axis-angle (convert to quaternion) 14 30 2 6

Composition of two axis-angle pairs is unthinkable in an application that requires

computational efﬁciency. One way to do the composition is to convert to matrices,

multiply the matrices, then extract the axis-angle pair. The two conversions from

axis-angle to matrix cost 26A + 30M + 4F , the matrix product costs 18A + 27M,

and the conversion from matrix to axis-angle costs 8A + 7M + 1D + 2F . The total

cost is 52A + 64M + 1D +6F .

Another way to do the composition of two axis-angle pairs is to convert to quater-

nions, multiply the quaternions, then extract the axis-angle pair. The two conversions

from axis-angle to quaternion cost 2A + 10M + 2D +4F , the quaternion product

costs 12A +16M, and the conversion from quaternion to axis-angle costs 4M + 2F .

The total cost is 14A + 30M + 2D +6F . Table A.5 is a summary of the operation

counts for composing two rotations.

Interpolation

This section discusses how to interpolate rotations for each of the three representa-

tion schemes.

A.4 Representations of 3D Rotations 867

Table A.6

Operation counts for quaternion interpolation.

Term AMDF

= p · q 34

θ = cos

−1

) 1

1 − t 1

(1 − t)θ 1

tθ 1

sin(θ) 1

sin((1 − t)θ) 1

sin(tθ) 1

= 1/ sin(θ) 1

= a

sin((1 − t)θ) 1

= a

sin(tθ) 1

p + a

q 48

Total 8 16 1 4

Quaternion Interpolation

Quaternions are quite amenable to interpolation. The standard operation that is used

is spherical linear interpolation, affectionately known as slerp. Given quaternions p

and q with acute angle θ between them, slerp is deﬁned as s(t; p, q) = p(p

∗

for

t ∈[0, 1]. Note that s(0; p, q) =p and s(1; p, q) =q. An equivalent deﬁnition of slerp

that is more amenable to calculation is

s(t; p, q) =

sin((1 − t)θ)p + sin(tθ )q

sin(θ)

If p and q are thought of as points on a unit circle, the formula above is a parameter-

ization of the shortest arc between them. If a particle travels on that curve according

to the parameterization, it does so with constant speed. Thus, any uniform sampling

of t in [0, 1]produces equally spaced points on the arc.

We assume that only p, q, and t are speciﬁed. Moreover, since q and −q rep-

resent the same rotation, you can replace q by −q if necessary to guarantee that

the angle between p and q treated as 4-tuples is acute. That is, p · q ≥ 0. As 4-

tuples, p and q are unit length. The dot product is therefore p · q = cos(θ). Table

A.6 shows the operation counts. Any term shown on the left that includes an already

computed term has only its additional operations counted to avoid double-counting

operations.

868 Appendix A Numerical Methods

Rotation Matrix Interpolation

The absence of a meaningful interpolation formula that directly applies to rotation

matrices is used as an argument for the superiority of quaternions over rotation

matrices. However, rotations can be interpolated directly in a way equivalent to what

slerp produces. If P and Q are rotations corresponding to quaternions p and q, the

slerp of the matrices is

S(t; P, Q) = P(P

the same formula that deﬁnes slerp for quaternions. The technical problem is to

deﬁne what is meant by R

for a rotation R and real-valued t. If the rotation has axis

ˆa and angle θ, then R

has the same rotation axis, but the angle of rotation is tθ.The

procedure for computing the slerp of the rotation matrices is

1. Compute R = P

2. Extract an axis ˆa and an angle θ from R.

3. Compute R

by converting the axis-angle pair ˆa, tθ.

4. Compute S(t; P, Q) = PR

This algorithm requires an axis-angle extraction that involves an inverse trigonomet-

ric function call and a square root operation, a couple of trigonometric evaluations

(for tθ), and a conversion back to a rotation matrix. This is quite a bit more expen-

sive than computing the slerp for quaternions, which requires three trigonometric

function calls. The quaternion interpolation is therefore more efﬁcient, but a purist

wishing to avoid quaternions in an application has, indeed, a method for interpolat-

ing rotation matrices.

Table A.7 shows the operation counts and uses the same format and rules as the

table for quaternion interpolation. Both the quaternion and rotation matrix interpo-

lation use 1 division and 4 function evaluations. However, the number of additions

and multiplications in the rotation matrix interpolation is excessive compared to that

of quaternion interpolation.

Axis-Angle Interpolation

There is no obvious and natural way to produce the same interpolation that occurs

with quaternions and rotation matrices. The only choice is to convert to one of the

other representations, interpolate in that form, then convert the interpolated result

back to axis-angle form. A very expensive proposition, just as in composition of

rotations.

A.5 Root Finding 869

Table A.7

Operation counts for rotation matrix interpolation.

Term A M D F

R = P

Q 18 27

= 0.5(Tra ce (R) − 1) 41

θ = cos

−1

) 1



d = (r

− r

, r

− r

, r

− r

) 3

= 1/|



d| 2311

ˆa =a



d 3

tθ 1

= sin(tθ ) 1

= 1 − cos(tθ ) 11

Matrix S,nocost

= I + a

S + a

18 27

Total 58 77 1 4

A.5 Root Finding

Given a continuous function



F : D ⊂R

→ R

, the problem is to ﬁnd an x (or ﬁnd

a set of points) for which



F(x) = 0.

A.5.1 Methods in One Dimension

Given a continuous function f :[a, b] → R, the ﬁrst question is whether or not

f(r)=0 for some r ∈[a, b]. I f f(a)f(b)<0, then there is at least one root. However,

there may be multiple roots. If a root r is computed, other analyses are required to

locate others. For example, if f is a polynomial and r is a root, the function can be

factored as f(t)= (t − r)

g(t),wherep ≥ 1 and g is a polynomial with degree(g) =

degree(f ) − p. The root-ﬁnding process is now continued with function g on [a, b].

If f(a)f(b)>0, there is no guarantee that f hasarooton[a, b].Forproblemsof

this type, a root-bounding preprocessing step can be used. The interval is partitioned

into t

= a + i(b − a)/n for 0 ≤ i ≤ n.Iff(t

)f (t

i+1

)<0 for some i, then that

subinterval is bisected to locate a root. A reasonable choice of n will be related to

what information the application knows about its function f .

870 Appendix A Numerical Methods

Finally, it might be necessary to ﬁnd roots of f : R → R where the domain of f

is not a bounded interval. Roots of f can be sought in the interval [−1, 1]. Any root

t of f outside this interval can be computed as t = 1/r,wherer ∈ [−1, 1]is a root of

g(r) = f(1/r).

Bisection

Bisection is the process of ﬁnding a root to a continuous function f :[a, b] → R by

bracketing a root with an interval, then successively bisecting the interval to narrow in

on the root. Suppose that initially f(a)f(b)<0. Since f is continuous, there must

bearootr ∈ (a, b). The midpoint of the interval is m = (a + b)/2. The function

value f (m) is computed and compared to the function values at the end points. If

f (a)f (m) < 0, then the subinterval (a, m) brackets a root and the bisection process

is repeated on that subinterval. If f (m)f (b) < 0, then the subinterval (m, b) brackets

a root and the bisection process is repeated instead on that subinterval. If f (m) = 0

or is zero within a speciﬁed tolerance, the process terminates. A stopping condition

might also be based on the length of the current subinterval—that is, if the length

becomes small enough, terminate the algorithm. If a root exists on [a, b], bisection is

guaranteed to ﬁnd it. However, the rate of convergence is slow.

Newton’s Method

Given a differentiable function f : R → R, an initial guess is chosen about where f

is zero, (x

, f(x

)). The tangent line to the graph at this point is used to update the

estimate to a (hopefully) better one. The tangent line is y − f(x

) = f



)(x − x

)

and intersects the x-axis at (0, x

),so−f(x

) =f



)(x

−x

). Assuming f



) =

0, solving for x

yields

= x

−

f(x

)



)

The next point in the iteration is (x

, f(x

)), and the process is repeated until a

stopping condition is met, typically one based on closeness of the function value to

zero. Unlike bisection, the iterations are not guaranteed to converge, but if there is

convergence, it is at a faster rate. Success depends a lot on the initial guess for x

Polynomial Roots

A polynomial of degree n is f(t)=



i=0

,wherea

= 0. While standard root

ﬁnders may be applied to polynomials, a better approach takes advantage of the

A.5 Root Finding 871

nature of such functions. For 2 ≤ n ≤ 4, there are closed-form equations for the

roots of the polynomial. Direct application of the formulas is possible, but numer-

ical problems tend to occur, particularly when the polynomial has a root of mul-

tiplicity larger than 1. For example, the roots of a quadratic f(t)= at

+ b

+ c

are t = (−b ±

√

− 4ac)/(2a).Ifb

− 4ac = 0, the quadratic has a double root

t =−b/(2a). However, numerical round-off errors might cause b

− 4ac =−.<0

for very small .. Another condition that leads to numerical problems is if a is nearly

zero. If so, it is possible to solve g(t) = t

f(1/t) = ct

+ bt + a = 0 and get t =

(−b ±

√

− 4ac)/(2c). But the problem still exists if c is also nearly zero. Similar

problems occur with the formulas for cubic and quartic polynomials.

An approach based on iteration schemes is to attempt to bracket the roots in a

way that each bracketing interval contains exactly one root. For each such interval,

bisection can be applied to ﬁnd the root. A hybrid scheme is also possible that mixes

bisection steps with Newton steps; the bisection step is used only when the Newton

step generates an iterate outside the current bracketing interval. The hope is that the

Newton iterates converge quickly to the root, but if they appear not to, bisection

attempts to generate better initial guesses for the Newton iteration.

Bounding Roots by Derivative Sequences

A simple approach to the bracketing problem is to partition R into intervals, the

polynomial f(t) being monotone on each interval. If it can be determined where

the derivative of the polynomial is zero, this set provides the partition. If d

and

i+1

are consecutive values for which f



) =f



i+1

) =0, then either f



(t) > 0on

, d

i+1

) or f



(t) < 0on(d

, d

i+1

). In either case, f can have at most one root on the

interval. The existence of this root is guaranteed by the condition f(d

)f (d

i+1

)<0

or f(d

) = 0orf(d

i+1

) = 0.

Solving f



(t) =0 requires the same techniques as solving f(t)=0. The difference

is that degree(f



) =degree(f ) − 1. A recursive implementation is warranted for this

problem, the base case being the constant polynomial that is either never zero or

identically zero on the real line.

If f



(t) = 0 for t ∈ (−∞, d

), it is possible that f has a root on the semi-inﬁnite

interval (−∞, d

]. Bisection does not help locate a root because the interval is un-

bounded. However, it is possible to determine the largest bounded interval that con-

tains the roots of a polynomial. The construction relies on the concepts of spectral

radius and norm of a matrix (see Horn and Johnson 1985). Given a square matrix

A, the spectral radius, denoted ρ(A), is the maximum of the absolute values of the

eigenvalues for the matrix. A matrix norm of A, denoted A, is a scalar-valued func-

tion that must satisfy the ﬁve conditions: A≥0, A=0 if and only if A = 0,

cA=|c|A for any scalar c, A + B≤A+B, and AB≤AB||.The

relationship between the spectral radius and any matrix norm is ρ(A) ≤A.Given

f(t)=



i=0

,wherea

= 1, the companion matrix is

872 Appendix A Numerical Methods

A =







−a

n−1

−a

n−2

··· −a

−a

10··· 00

01··· 00

00··· 10







The characteristic polynomial is f(t)= det(A − tI), so the roots of f are the

eigenvalues of A. The spectral norm therefore provides a bound for the roots. Since

there are lots of matrix norms to choose from, there are many possible bounds. One

such bound is Cauchy’s bound:

|t|≤max{|a

|,1+|a

|, ...,1+|a

n−1

|} = 1 + max{|a

|, ..., |a

n−1

Another bound that can be obtained is the Carmichael and Mason bound:

|t|≤







1 +

n−1



i=0

If a

= 0, then f(0) = 0, so the roots of f are bounded away from zero. It is possible

to construct lower bounds by using g(t) = [t

f(1/t)]/a

. The roots of g(t) are the

reciprocal roots of f(t). Cauchy’s bound applied to g(t), then taking reciprocals, is

|t|≥

1 + max{1, |a

|, ..., |a

n−1

The Carmichael and Mason bound is

|t|≥



1 +



n−1

i=0

These bounds are used in the recursive call to determine where f(t)is monotone.

The polynomial can be factored f(t)=t

g(t),wherep ≥0 and g is a polynomial for

which g(0) = 0. If p = 0, then f =g and f is processed for 0 <a≤|t |≤b,where

a and b are bounds computed from the previously mentioned inequalities. If p>0,

then g is processed on the intervals obtained by using the bounds from the same

inequalities.

Bounding Roots by Sturm Sequences

Consider a polynomial f(t)deﬁned on interval [a, b]. A Sturm sequence for f isaset

of polynomials f

(t),0≤i ≤m, such that degree(f

i+1

)>degree(f

) and the number

of distinct real roots for f in [a, b]isN = s(a) − s(b),wheres(a) is the number

of sign changes of f

(a), ..., f

(a) and s(b) is the number of sign changes of

A.5 Root Finding 873

Table A.8

Signs of the Sturm polynomials for t

+ 3t

− 1 at various t values.

t Sign f

(t) Sign f

(t) Sign changes

−∞ − + − + 3

−3 −+−+ 3

−2 + 0 −+ 2

−1 +−−+ 2

0 − 0 ++ 1

+1 ++++ 0

+∞ + + + + 0

Table A.9 Signs of the Sturm polynomials for (t − 1)

at various t values.

t Sign f

(t) Sign f

(t) Sign changes

−∞ − + 01

0 −+01

+∞ + + 00

(b), ..., f

(b). The total number of real-valued roots of f on R is s(−∞) −s(∞).

It is not always the case that m = degree(f ).

The classic Sturm sequence is f

(t) =f(t), f

(t) =f



(t), and f

(t) =−Remain-

der (f

i−2

i−1

) for i ≥ 2. The polynomials are generated by this method until the

remainder term is a constant. An instructive example from the article by D. G. Hook

and P. R. McAree in Glassner (1990) is f(t)= t

+ 3t

− 1. The Sturm sequence is

(t) =t

+3t

−1, f

(t) =3t

+6t, f

(t) =2t +1, and f

=9/4. Table A.8 lists the

signs of the Sturm polynomials for various t values. Letting N(a, b) denote the num-

ber of real-valued roots on the interval (a, b), the table shows that N(−∞, −3) = 0,

N(−3, −2) =1, N(−2, −1) =0, N(−1, 0) =1, N(0, 1) =1, and N(1, ∞) =0. More-

over, the number of negative real roots is N(−∞,0) =2, the number of positive real

roots is N(0, ∞) = 1, and the total number of real roots is N(−∞, ∞) = 3.

The next example shows that the number of polynomials in the Sturm sequence is

not necessarily the degree(f ) +1. The function f(t)=(t −1)

has a Sturm sequence

(t) = (t − 1)

, f

(t) = 3(t − 1)

, and f

(t) ≡ 0 since f

exactly divides f

with no

remainder. Table A.9 lists sign changes for f at various t values. The total number of

real roots is N(−∞, ∞) = 1.