Schneider P., Eberly D.H. Geometric Tools for Computer Graphics
Подождите немного. Документ загружается.
844 Appendix A Numerical Methods
B
´
ezout matrix for d and e, the 4 ×4 matrix M = [m
ij
] with
m
ij
=
min(4,7−i−j)
k=max(4−j ,4−i)
w
k,7−i−j−k
for 0 ≤ i ≤ 3 and 0 ≤ j ≤ 3, with w
i,j
= d
i
e
j
− d
j
e
i
for 0 ≤ i ≤ 4 and 0 ≤ j ≤ 4. In
expanded form,
M =
w
4,3
w
4,2
w
4,1
w
4,0
w
4,2
w
3,2
+ w
4,1
w
3,1
+ w
4,0
w
3,0
w
4,1
w
3,1
+ w
4,0
w
2,1
+ w
3,0
w
2,0
w
4,0
w
3,0
w
2,0
w
1,0
Thedegreeofw
i,j
is 8 −i −j . The B
´
ezout determinant det(M(z)) is a polynomial of
degree 16 in z. For each solution ¯z to det(M(z)) = 0, we need to find corresponding
values ¯x and ¯y.
Using the B
´
ezout method hides the intermediate polynomials that were con-
veniently used in the previous cases to compute the other variables. Let us find
them explicitly. The elimination process may be applied directly to d(y) = d
0
+
d
1
y +d
2
y
2
+d
3
y
3
+d
4
y
4
and e(y) =e
0
+e
1
y +e
2
y
2
+e
3
y
3
+e
4
y
4
. Define f(y)=
e
4
d(y) − d
4
e(y) = f
0
+f
1
y +f
2
y
2
+f
3
y
3
. The coefficients are f
i
=e
4
d
i
−e
i
d
4
for
all i. Define g(y) = f
3
d(y) − d
4
yf (y) = g
0
+ g
1
y + g
2
y
2
+ g
3
y
3
. The coefficients
are g
0
=f
3
d
0
, g
1
=f
3
d
1
−f
0
d
4
, g
2
=f
3
d
2
−f
1
d
4
, and g
3
=f
3
d
3
−f
2
d
4
.Nowf(y)
and g(y) are cubic polynomials. The process is repeated. Define h(y) = g
3
f(y)−
f
3
g(y) =h
0
+h
1
y +h
2
y
2
,whereh
i
=g
3
f
i
−f
3
g
i
for all i. Define m(y) =h
2
f(y)−
f
3
yh(y) =m
0
+m
1
y +m
2
y
2
,wherem
0
=h
2
f
0
, m
1
=h
2
f
1
−h
0
f
3
, and m
2
=h
2
f
2
−
h
1
f
3
.Nowh(y) and m(y) are quadratic polynomials. As we saw earlier, if the poly-
nomials have a common solution, it must be ¯y = (h
2
m
0
− h
0
m
2
)/(h
1
m
2
− h
2
m
1
).
Because the d
i
and e
i
coefficients depend on ¯z, the values h
i
and m
i
depend on ¯z.
Thus, given a value for ¯z, we compute a corresponding value ¯y.Tocompute ¯x for
a specified pair ( ¯y, ¯z), F(x, ¯y, ¯z) = a
2
x
2
+ a
1
x + a
0
= 0 and G(x, ¯y, ¯z) = b
2
x
2
+
b
1
x + b
0
= 0 are two quadratic equations in the unknown x, so a common solution
is ¯x = (a
2
b
0
− a
0
b
2
)/(a
1
b
2
− a
2
b
1
).
Example Let F(x, y, z) = (x − 1)
2
+y
2
+z
2
−4 (sphere), G(x, y, z) = x
2
+4y
2
−4z (para-
boloid), and H(x, y, z) = x
2
+ 4(y − 1)
2
+ z
2
− 4 (ellipsoid). The polynomial d(y)
with coefficients dependent on z that represents D(y, z) has coefficients d
0
=−9 +
40z −10z
2
−8z
3
−z
4
, d
1
=0, d
2
=−34 +24z +6z
2
, d
3
=0, and d
4
=−9. The poly-
nomial e(y) with coefficients dependent on z that represents E(y, z) has coefficients
e
0
=−9 −4z
2
, e
1
= 80, e
2
= 98, e
3
= 48, and e
4
=−9. The w terms for the matrix
M are
A.2 Systems of Polynomials 845
w
1,0
= 720 −3200z + 800z
2
+ 640z
3
+ 80z
4
w
2,0
=−576 +3704z − 898z
2
− 880z
3
− 122z
4
w
3,0
= 432 −1920z + 480z
2
+ 384z
3
+ 48z
4
w
4,0
= 360z − 54z
2
− 72z
3
− 9z
4
w
2,1
=−2720 +1920z + 480z
2
w
3,1
= 0
w
3,2
= 1632 −1152z − 288z
2
w
4,1
=−720
w
4,2
= 576 +216z + 54z
2
w
4,3
=−432
The determinant of M is degree 16 and has coefficients µ
i
for 0 ≤ i ≤16 given by
µ
0
= 801868087296 µ
9
= 2899639296
µ
1
=−4288520650752 µ
10
=−105691392
µ
2
= 4852953907200 µ
11
=−211071744
µ
3
=−779593973760 µ
12
= 4082400
µ
4
=−1115385790464 µ
13
= 13856832
µ
5
= 307850969088 µ
14
= 2624400
µ
6
= 109063397376 µ
15
= 209952
µ
7
=−34894540800 µ
16
= 6561
µ
8
=−3305131776
The real-valued roots to det(M)(z) = 0 are 0.258255, 1.46202, 1.60199, and 1.63258.
Observe that the coefficients of the polynomial are quite large. For numerical stability,
for each z you should evaluate the determinant using Gaussian elimination rather
than actually computing the µ
i
as a preprocess.
The cubic polynomial f(y)has coefficients that depend on z:
f
0
=−360z + 54z
2
+ 72z
3
+ 9z
4
f
1
= 720
f
2
=−576 −216z − 54z
2
f
3
= 432
846 Appendix A Numerical Methods
The cubic polynomial g(y) has coefficients that depend on z:
g
0
=−3888 + 17280z − 4320z
2
− 3456z
3
+ 432z
4
g
1
=−3240z + 486z
2
+ 648z
3
+ 81z
4
g
2
=−8208 + 10368z + 2592z
2
g
3
=−5184 −1944z − 486z
2
The quadratic polynomial h(y) has coefficients that depend on z:
h
0
= 1679616 −5598720z + 2286144z
2
+ 1189728z
3
− 26244z
4
− 52488z
5
− 4374z
6
h
1
=−3732480 −559872z
2
− 279936z
3
− 34992z
4
h
2
= 6531840 −2239488z −139968z
2
+ 209952z
3
+ 26244z
4
The quadratic polynomial m(y) has coefficients that depend on z:
m
0
=−2351462400z + 1158935040z
2
+ 399748608z
3
− 185597568z
4
− 28343520z
5
+ 15274008z
6
+ 3779136z
7
+ 236196z
8
m
1
= 3977330688 +806215680z − 1088391168z
2
− 362797056z
3
+ 30233088z
4
+ 22674816z
5
+ 1889568z
6
m
2
=−2149908480 −120932352z + 453496320z
2
− 37791360z
4
− 17006112z
5
− 1417176z
6
The ¯z-values are {0.258225, 1.46202, 1.60199, 1.63258}. The corresponding
¯y-values are computed from ¯y = (h
2
(¯z)m
0
(¯z) − h
0
(¯z)m
2
(¯z))/(h
1
(¯z)m
2
(¯z)−
h
2
(¯z)m
1
(¯z)). The ordered set of such values is {0.137429, 0.336959, 1.18963, 1.14945}.
The corresponding ¯x-values are computed from ¯x = (a
2
( ¯y, ¯z)b
0
( ¯y, ¯z) − a
0
( ¯y, ¯z)b
2
( ¯y, ¯z))/(a
1
( ¯y, ¯z)b
2
( ¯y, ¯z) − a
2
( ¯y, ¯z)b
1
( ¯y, ¯z)). The ordered set of such values is
{−0.97854, 2.32248, 0.864348, 1.11596}. As it turns out, only the ( ¯x, ¯y, ¯z) triples
(−0.97854, 0.137429, 0.258225) and (1.11596, 1.14945, 1.53258) are true intersec-
tion points. The other two are extraneous solutions (the H -values are 7.18718 and
−0.542698, respectively), an occurrence that is to be expected since the elimination
process increases the degrees of the polynomials, thereby potentially introducing
roots that are not related to the intersection problem.
A.3 Matrix Decompositions 847
A.3 Matrix Decompositions
We present a collection of various matrix factorizations that arise in many computer
graphics applications. An excellent reference for numerical methods relating to ma-
trices is Golub and Van Loan (1993b). An excellent reference for matrix analysis is
Horn and Johnson (1985).
A.3.1 Euler Angle Factorization
Rotations about the coordinate axes are easy to define and work with. Rotation about
the x-axis by angle θ is
R
x
(θ) =
10 0
0cosθ − sin θ
0 sin θ cos θ
where θ>0 indicates a counterclockwise rotation in the plane x =0. The observer
is assumed to be positioned on the side of the plane with x>0 and looking at the
origin.
Rotation about the y-axis by angle θ is
R
y
(θ) =
cos θ 0 sin θ
010
− sin θ 0cosθ
where θ>0 indicates a counterclockwise rotation in the plane y =0. The observer
is assumed to be positioned on the side of the plane with y>0 and looking at the
origin.
Rotation about the z-axis by angle θ is
R
z
(θ) =
cos θ − sin θ 0
sin θ cos θ 0
001
where θ>0 indicates a counterclockwise rotation in the plane z = 0. The observer
is assumed to be positioned on the side of the plane with z>0 and looking at the
origin.
Rotation by an angle θ about an arbitrary axis containing the origin and having
unit-length direction ˆu = (u
x
, u
y
, u
z
) is given by
R
ˆu
(θ) = I + (sin θ)S +(1 − cos θ)S
2
848 Appendix A Numerical Methods
where I is the identity matrix,
S =
0 −u
z
u
y
u
z
0 −u
x
−u
y
u
x
0
and θ>0 indicates a counterclockwise rotation in the plane ˆu · (x, y, z) = 0. The
observer is assumed to be positioned on the side of the plane to which ˆu points and
is looking at the origin.
Factoring Rotation Matrices
A common problem is to factor a rotation matrix as a product of rotations about the
coordinate axes. The form of the factorization depends on the needs of the applica-
tion and what ordering is specified. For example, you might want to factor a rotation
as R = R
x
(θ
x
)R
y
(θ
y
)R
z
(θ
z
) for some angles θ
x
, θ
y
, and θ
z
. The ordering is xyz.Five
other possibilities are xzy, yxz, yzx, zxy, and zyx. You might also envision factor-
izations such as xyx—these are not discussed here. The following discussion uses the
notation c
a
= cos(θ
a
) and s
a
= sin(θ
a
) for a = x, y, z.
Factor as R
x
R
y
R
z
Setting R = [r
ij
] for 0 ≤ i ≤ 2 and 0 ≤ j ≤ 2, formally multiplying R
x
(θ
x
)R
y
(θ
y
)
R
z
(θ
z
), and equating yields
r
00
r
01
r
02
r
10
r
11
r
12
r
20
r
21
r
22
=
c
y
c
z
−c
y
s
z
s
y
c
z
s
x
s
y
+ c
x
s
z
c
x
c
z
− s
x
s
y
s
z
−c
y
s
x
−c
x
c
z
s
y
+ s
x
s
z
c
z
s
x
+ c
x
s
y
s
z
c
x
c
y
From this we have s
y
= r
02
,soθ
y
= sin
−1
(r
02
).Ifθ
y
∈ (−π/2, π/2), then c
y
= 0 and
c
y
(s
x
, c
x
) = (−r
12
, r
22
), in which case θ
x
= atan2(−r
12
, r
22
). Similarly, c
y
(s
z
, c
z
) =
(−r
01
, r
00
), in which case θ
z
= atan2(−r
01
, r
00
).
If θ
y
= π/2, then s
y
= 1 and c
y
= 0. In this case
r
10
r
11
r
20
r
21
=
c
z
s
x
+ c
x
s
z
c
x
c
z
− s
x
s
z
−c
x
c
z
+ s
x
s
z
c
z
s
x
+ c
x
s
z
=
sin(θ
z
+ θ
x
) cos(θ
z
+ θ
x
)
− cos(θ
z
+ θ
x
) sin(θ
z
+ θ
x
)
Therefore, θ
z
+θ
x
=atan2(r
10
, r
11
). There is one degree of freedom, so the factoriza-
tion is not unique. One choice is θ
z
= 0 and θ
x
= atan2(r
10
, r
11
).
A.3 Matrix Decompositions 849
If θ
y
=−π/2, then s
y
=−1 and c
y
= 0. In this case
r
10
r
11
r
20
r
21
=
−c
z
s
x
+ c
x
s
z
c
x
c
z
+ s
x
s
z
c
x
c
z
+ s
x
s
z
c
z
s
x
− c
x
s
z
=
sin(θ
z
− θ
x
) cos(θ
z
− θ
x
)
cos(θ
z
− θ
x
) − sin(θ
z
− θ
x
)
Therefore, θ
z
−θ
x
=atan2(r
10
, r
11
). There is one degree of freedom, so the factoriza-
tion is not unique. One choice is θ
z
= 0 and θ
x
=−atan2(r
10
, r
11
).
Pseudocode for the factorization is
thetaY = asin(r02);
if (thetaY < PI/2) {
if (thetaY > -PI / 2) {
thetaX = atan2(-r12, r22);
thetaZ = atan2(-r01, r00);
} else {
// not a unique solution (thetaX - thetaZ constant)
thetaX = -atan2(r10, r11);
thetaZ = 0;
}
} else {
// not a unique solution (thetaX + thetaZ constant)
thetaX = atan2(r10, r11);
thetaZ = 0;
}
Factor as R
x
R
z
R
y
Setting R = [r
ij
] for 0 ≤ i ≤ 2 and 0 ≤ j ≤ 2, formally multiplying R
x
(θ
x
)R
z
(θ
z
)
R
y
(θ
y
), and equating yields
r
00
r
01
r
02
r
10
r
11
r
12
r
20
r
21
r
22
=
c
y
c
z
−s
z
c
z
s
y
s
x
s
y
+ c
x
c
y
s
z
c
x
c
z
−c
y
s
x
+ c
x
s
y
s
z
−c
x
s
y
+ c
y
s
x
s
z
c
z
s
x
c
x
c
y
+ s
x
s
y
s
z
Analysis similar to the xyz case leads to the pseudocode
thetaZ = asin(-r01);
if (thetaZ < PI / 2) {
if (thetaZ > -PI / 2) {
thetaX = atan2(r21, r11);
thetaY = atan2(r02, r00);
} else {
// not a unique solution (thetaX + thetaY constant)
thetaX = atan2(-r20, r22);
850 Appendix A Numerical Methods
thetaY = 0;
}
} else {
// not a unique solution (thetaX - thetaY constant)
thetaX = atan2(r20, r22);
thetaY = 0;
}
Factor as R
y
R
x
R
z
Setting R = [r
ij
] for 0 ≤ i ≤ 2 and 0 ≤ j ≤ 2, formally multiplying R
y
(θ
y
)R
x
(θ
x
)
R
z
(θ
z
), and equating yields
r
00
r
01
r
02
r
10
r
11
r
12
r
20
r
21
r
22
=
c
y
c
z
+ s
x
s
y
s
z
c
z
s
x
s
y
− c
y
s
z
c
x
s
y
c
x
s
z
c
x
c
z
−s
x
−c
z
s
y
+ c
y
s
x
s
z
c
y
c
z
s
x
+ s
y
s
z
c
x
c
y
Analysis similar to the xyz case leads to the pseudocode
thetaX = asin(-r12);
if (thetaX < PI / 2) {
if (thetaX > -PI / 2) {
thetaY = atan2(r02, r22);
thetaZ = atan2(r10, r11);
} else {
// not a unique solution (thetaY + thetaZ constant)
thetaY = atan2(-r01, r00);
thetaZ = 0;
}
} else {
// not a unique solution (thetaY - thetaZ constant)
thetaY = atan2(r01, r00);
thetaZ = 0;
}
Factor as R
y
R
z
R
x
Setting R = [r
ij
] for 0 ≤ i ≤ 2 and 0 ≤ j ≤ 2, formally multiplying R
y
(θ
y
)R
z
(θ
z
)
R
x
(θ
x
), and equating yields
r
00
r
01
r
02
r
10
r
11
r
12
r
20
r
21
r
22
=
c
y
c
z
s
x
s
y
− c
x
c
y
s
z
c
x
s
y
+ c
y
s
x
s
z
s
z
c
x
c
z
−c
z
s
x
−c
z
s
y
c
y
s
x
+ c
x
s
y
s
z
c
x
c
y
− s
x
s
y
s
z
A.3 Matrix Decompositions 851
Analysis similar to the xyz case leads to the pseudocode
thetaZ = asin(r10);
if (thetaZ < PI / 2) {
if (thetaZ > -PI / 2) {
thetaY = atan2(-r20, r00);
thetaX = atan2(-r12, r11);
} else {
// not a unique solution (thetaX - thetaY constant)
thetaY = -atan2(r21, r22);
thetaX = 0;
}
} else {
// not a unique solution (thetaX + thetaY constant)
thetaY = atan2(r21, r22);
thetaX = 0;
}
Factor as R
z
R
x
R
y
Setting R = [r
ij
] for 0 ≤ i ≤ 2 and 0 ≤ j ≤ 2, formally multiplying R
z
(θ
z
)R
x
(θ
x
)
R
y
(θ
y
), and equating yields
r
00
r
01
r
02
r
10
r
11
r
12
r
20
r
21
r
22
=
c
y
c
z
− s
x
s
y
s
z
−c
x
s
z
c
z
s
y
+ c
y
s
x
s
z
c
z
s
x
s
y
+ c
y
s
z
c
x
c
z
−c
y
c
z
s
x
+ s
y
s
z
−c
x
s
y
s
x
c
x
c
y
Analysis similar to the xyz case leads to the pseudocode
thetaX = asin(r21);
if (thetaX < PI / 2) {
if (thetaX > -PI / 2) {
thetaZ = atan2(-r01, r11);
thetaY = atan2(-r20, r22);
} else {
// not a unique solution (thetaY - thetaZ constant)
thetaZ = -atan2(r02, r00);
thetaY = 0;
}
} else {
// not a unique solution (thetaY + thetaZ constant)
thetaZ = atan2(r02, r00);
thetaY = 0;
}
852 Appendix A Numerical Methods
Factor as R
z
R
y
R
x
Setting R = [r
ij
] for 0 ≤ i ≤ 2 and 0 ≤ j ≤ 2, formally multiplying R
z
(θ
z
)R
y
(θ
y
)
R
x
(θ
x
), and equating yields
r
00
r
01
r
02
r
10
r
11
r
12
r
20
r
21
r
22
=
c
y
c
z
c
z
s
x
s
y
− c
x
s
z
c
x
c
z
s
y
+ s
x
s
z
c
y
s
z
c
x
c
z
+ s
x
s
y
s
z
−c
z
s
x
+ c
x
s
y
s
z
−s
y
c
y
s
x
c
x
c
y
Analysis similar to the xyz case leads to the pseudocode
thetaY = asin(-r20);
if (thetaY < PI / 2) {
if (thetaY > -PI / 2) {
thetaZ = atan2(r10, r00);
thetaX = atan2(r21, r22);
} else {
// not a unique solution (thetaX + thetaZ constant)
thetaZ = atan2(-r01, -r02);
thetaX = 0;
}
} else {
// not a unique solution (thetaX - thetaZ constant)
thetaZ = -atan2(r01, r02);
thetaX = 0;
}
A.3.2 QR Decomposition
Givenaninvertiblen × n matrix A, we wish to decompose it into A = QR,whereQ
is an orthogonal matrix and R is upper triangular. The factorization is just an appli-
cation of the Gram-Schmidt orthonormalization algorithm applied to the columns
of A. Let those columns be denoted as a
i
for 1 ≤ i ≤ n. The columns are linearly in-
dependent since A is assumed to be invertible. The Gram-Schmidt process constructs
from this set of vectors an orthonormal set ˆq
i
,1≤ i ≤ n. That is, each vector is unit
length, and the vectors are mutually perpendicular.
The first step is simple; just normalize ˆq
1
by
ˆq
1
=
a
1
a
1
We can project a
2
onto the orthogonal complement of ˆq
1
and represent a
2
= c
1
ˆq
1
+
p
2
,where ˆq
1
·p
2
=0. Dotting this equation with ˆq
1
leads to c
1
=ˆq
1
·a
2
. Rewritten, we
A.3 Matrix Decompositions 853
have p
2
=a
2
− ( ˆq
1
·a
2
) ˆq
1
. The vectors ˆq
1
and p
2
are perpendicular by the construc-
tion, but p
2
is not necessarily unit length. Thus, define ˆq
2
to be the normalized p
2
:
ˆq
2
=
a
2
− ( ˆq
1
·a
2
) ˆq
1
a
2
− ( ˆq
1
·a
2
) ˆq
1
A similar construction is applied to a
3
. We can project a
3
onto the orthogonal
complement of the space spanned by ˆq
1
and ˆq
2
and represent a
3
= c
1
ˆq
1
+ c
2
ˆq
2
+p
3
,
where ˆq
1
·p
3
= 0 and ˆq
2
·p
3
= 0. Dotting this equation with the ˆq
i
leads to c
i
=
ˆq
i
a
3
.Rewritten,wehave p
3
=a
3
− ( ˆq
1
·a
3
) ˆq
1
− ( ˆq
2
·a
3
) ˆq
2
. The next vector in the
orthonormal set is the normalized p
3
:
ˆq
3
=
a
3
− ( ˆq
1
·a
3
) ˆq
1
− ( ˆq
2
·a
3
) ˆq
2
a
3
− ( ˆq
1
·a
3
) ˆq
1
− ( ˆq
2
·a
3
) ˆq
2
In general for i ≥ 2,
ˆq
i
=
a
i
−
i−1
j=1
( ˆq
j
·a
i
) ˆq
j
a
i
−
i−1
j=1
( ˆq
j
·a
i
) ˆq
j
Let Q be the n ×n matrix whose columns are the vectors ˆq
i
. The upper triangular
matrix in the factorization is
R =
ˆq
1
·a
1
ˆq
1
·a
2
··· ˆq
1
·a
n
0 ˆq
2
·a
2
··· ˆq
2
·a
n
.
.
.
.
.
.
.
.
.
00··· ˆq
n
·a
n
A.3.3 Eigendecomposition
Given an n × n matrix A,aneigensystem is of the form Ax = λx or (A − λI)x =
0.
It is required that there be solutions x =
0. For this to happen, the matrix A − λI
must be noninvertible. This is the case when det(A − λI) = 0, a polynomial in λ
of degree n called the characteristic polynomial for A. For each root λ, the matrix
A − λI is computed, and the system (A − λI)x =
0 is solved for nonzero solutions.
While standard linear algebra textbooks show numerous examples for doing this
symbolically, most applications require a robust numerical method for doing so.
In particular, if n ≥ 5, there are no closed formulas for roots to polynomials, so
numerical methods must be applied. A good reference on solving eigensystems is
Press et al. (1988).
Most of the applications in graphics that require eigensystems have symmetric
matrices. The numerical methods are quite good for these since there is guaranteed