Schneider P., Eberly D.H. Geometric Tools for Computer Graphics
Подождите немного. Документ загружается.
824 Chapter 13 Computational Geometry Topics
13.12.3 Volume of a Polyhedron
The discussion of this section is the direct extension of the ideas for area of a polygon
in 2D. Consider a tetrahedron V
0
, V
1
, V
2
, V
3
. Setting V
i
= (x
i
, y
i
, z
i
), methods of
basic algebra and trigonometry can be used to show that the signed volume of the
tetrahedron is
Volume(V
0
, V
1
, V
2
, V
3
) =
1
6
det
1111
x
0
x
1
x
2
x
3
y
0
y
1
y
2
y
3
z
0
z
1
z
2
z
3
(13.10)
Assuming the points are not all coplanar, half of the permutations produce a positive
value and half produce the negative of that value.
Let V =(x, y, z) be an arbitrary point in space. The following algebraic identity
is true:
Volume(V
0
, V
1
, V
2
, V
3
) = Volume(V , V
0
, V
1
, V
2
) + Volume(V , V
1
, V
2
, V
3
)
+ Volume(V , V
2
, V
3
, V
0
) + Volume(V , V
3
, V
0
, V
1
)
(13.11)
The identity can be verified by expanding the determinants on the right-hand side of
the equation and performing algebraic operations to show that the result is the same
as the determinant on the left-hand side. The geometric motivation is the same as in
2D. Some of the signed volume terms are positive; some are negative. Whenever two
tetrahedra that share V overlap and have opposite sign volumes, there is cancellation.
The volume of a simple polyhedron P whose faces are all triangles can be com-
puted by using the extension of Equation 13.3 to three dimensions. The faces are
assumed to be counterclockwise oriented as viewed from the outside of the polyhe-
dron. The point V is arbitrary; in practice it is chosen to be the zero vector
0.
Volume(P ) =
face F
Volume(
0, F .V
0
, F .V
1
, F .V
2
) (13.12)
The analytic construction using methods of calculus is presented below. Consider
a polyhedron, denoted S, that encloses a simply connected region R in space. Let
the n faces of the polyhedron be named S
i
for 0 ≤ i<n. Let each face S
i
have unit-
length outer normal ˆn
i
and vertices P
i,j
for 0 ≤ j < m(i), counterclockwise ordered
when viewed from outside, where it is emphasized that the total number of vertices,
m(i), depends on the face i. A formula for the volume enclosed by the polyhedron
is obtained by an application of the Divergence Theorem, the direct generalization
of Green’s Theorem from 2D to 3D, to a simply connected region R with boundary
surface S
R
∇·Fdxdydz=
S
F ·ˆndσ
13.12 Area and Volume Measurements 825
where F(x, y, z) =(F
1
(x, y, z), F
2
(x, y, z), F
3
(x, y, z)) is a differentiable vector field,
∇·F = ∂F
1
/∂x + ∂F
2
/∂y + ∂F
3
/∂z is the divergence of the vector field, ˆn is an
outward-pointing normal vector to the surface S, and dσ represents an infinitesimal
element of the surface. This formula can be found in many standard texts on calculus
(for example, Finney and Thomas 1996).
The volume formulation of the Divergence Theorem arises when we choose F =
(x, y, z)/3. In this case,
∇·F ≡1 and
Volume(R) =
R
dx dy dz =
1
3
S
ˆn · (x, y, z) dσ
Since the polyhedron is a disjoint union S =
n−1
i=0
S
i
, the integral on the right
becomes a sum of integrals, each integral related to a single polygonal face of the
polyhedron
S
ˆn · (x, y, z) dσ =
n−1
i=0
S
i
ˆn
i
· (x, y, z) dσ
where ˆn
i
is the unit-length outer normal to S
i
. The plane of S
i
is ˆn
i
·(x, y, z) =c
i
for
some constant c
i
. Any point on the polygon determines the constant, in particular
c =ˆn
i
· P
0,i
. The integral further reduces to
n−1
i=0
S
i
ˆn
i
· (x, y, z) dσ =
n−1
i=0
S
i
c
i
dσ =
n−1
i=0
c
i
Area(S
i
)
Substituting the formula from Equation 13.9, the volume formula is
Volume(R) =
1
6
n−1
i=0
( ˆn
i
· P
0,i
) ˆn
i
·
m(i)−1
j=0
(P
i,j
× P
i,j +1
)
(13.13)
Appendix
ANumerical
Methods
A.1 Solving Linear Systems
The general form of an m × n system of linear equations is
a
1,1
x
1
+ a
1,2
x
2
+···+a
1,n
x
n
= c
1
a
2,1
x
1
+ a
2,2
x
2
+···+a
2,n
x
n
= c
2
.
.
.
a
m,1
x
1
+ a
m,2
x
2
+···+a
m,n
x
n
= c
n
Frequently, linear systems are written in matrix form:
Ax =c
a
1,1
a
1,2
··· a
1,n
a
2,1
a
2,2
··· a
2,n
.
.
.
a
m,1
a
m,2
··· a
m,n
x
1
x
2
.
.
.
x
n
=
c
1
c
2
.
.
.
c
n
827
828 Appendix A Numerical Methods
The matrix A isreferredtoasthecoefficient matrix, and the matrix
a
1,1
a
1,2
··· a
1,n
c
1
a
2,1
a
2,2
··· a
2,n
c
2
.
.
.
a
m,1
a
m,2
··· a
m,n
c
n
is known as the augmented matrix.
A.1.1 Special Case: Solving a Triangular System
In order to introduce a general method for solving square linear systems, let’s look at
a special case: that of a triangular system. Such a system is one in which a
i,j
= 0 for
i>j; it is called triangular because of its characteristic shape:
a
1,1
x
1
+ a
1,2
x
2
+····································+a
1,n
x
n
= c
1
a
2,2
x
2
+····································+a
2,n
x
n
= c
2
.
.
.
a
n−2,n−2
x
n−2
+ a
n−2,n−1
x
n−1
+ a
n−2,n
x
n
= c
n−2
a
n−1,n−1
x
n−1
+ a
n−1,n
x
n
= c
n−1
a
n,n
x
n
= c
n
We can solve this by a technique called back substitution. If we look at the last of
the equations, we see that it has just one unknown (x
n
), which we can solve trivially:
x
n
=
c
n
a
n,n
Now, since we have a solution for x
n
, we can substitute it into the second-to-last
equation and again trivially solve for x
n−1
:
x
n−1
=
c
n−1
− a
n−1,n
c
n
a
n,n
a
n−1,n−1
and so on backward, until we have reached the first equation and solve for the first
unknown, x
1
.Ingeneral,x
k
is computed by substituting the previously computed
values of x
n
, x
n−1
, ···, x
k+1
in the kth equation:
x
k
=
c
k
−
n
m=k+1
a
k,m
x
m
a
k,k
A.1 Solving Linear Systems 829
A.1.2 Gaussian Elimination
Gaussian elimination is perhaps the most widely used method for solving general
linear systems. It is based on the back-substitution technique we just covered: first,
the system is converted step by step into triangular form, and then the previously
described back-substitution technique is applied. In fact, we’ve already seen a (trivial)
example of this algorithm in Section 2.4.2.
Gauss observed that a system of linear equations can be subjected to certain
transformations and still have the same set of solutions. These modifications are
1. exchanging two rows
2. multiplying one row by a nonzero factor
3. replacing two rows with their sum
We saw transformations 2 and 3 in the earlier example, but you may be wondering
what the point of exchanging rows might be. In the world of infinite-precision math-
ematics, it is irrelevant, but in the world of computer implementation, it is quite a
beneficial transformation.
Consider a simple system of two equations in two unknowns, taken from Johnson
and Riess (1982):
(1) 0.0001x
1
+ x
2
= 1
(2) x
1
+ x
2
= 2
If we want to eliminate the x
1
from (2), we’d need to multiply (1) by
−10000 ×(1): −x
1
− 10000x
2
=−10000
1 × (2): x
1
+ x
2
= 2
Sum: −9999x
2
=−9998
yielding us a new system:
(1) 0.0001x
1
+ x
2
= 1
(2) −9999x
2
=−9998
If we complete the solution by computing
x
2
=
9998
9999
≈ 0.99989999
830 Appendix A Numerical Methods
and substitute this back into (1):
0.0001x
1
+ x
2
= 1
0.0001x
1
+
9998
9999
= 1
0.0001x
1
= 1 −
9998
9999
0.0001x
1
=
1
9999
x
1
=
10000
9999
x
1
≈ 1.00010001
One “problem” with this is that the above calculation is assuming infinite preci-
sion. Instead, consider what happens if we calculated with, say, six digits and stored
only three. Then, if we again multiply (1) by −10,000, we would get
−10000 ×(1): x
1
− 10000x
2
=−10000
1 × (2): x
1
+ x
2
= 2
Sum: −10000x
2
=−10000
because of round-off, and x
2
= 1. If we substitute this back into (1), we’d get
0.0001x
1
+ x
2
= 1
0.0001x
1
+ 1 = 1
0.0001x
1
= 0
x
1
= 0
which is quite wrong. Of course, this example was contrived to show a worst-case
scenario, but truncation and round-off errors will tend to compound themselves in
each successive step.
Another precision-related problem can occur in the division step of the back
substitution. Recall that the basic back-substitution step is
x
k
=
c
k
−
n
m=k+1
a
k,m
x
m
a
k,k
That is, the coefficient appears in the denominator. As division in a computer is most
accurate when the divisor has as large an absolute value as possible, it would be best
A.1 Solving Linear Systems 831
if the x
k
with the largest coefficient (among all x
k,n
) were the one subjected to the
division.
Here’s an example: Suppose we have the following system in three unknowns:
(1) x
1
− 3x
2
+ 2x
3
= 6
(2) 2x
1
− 4x
2
+ 2x
3
= 18
(3) −3x
1
+ 8x
2
+ 9x
3
=−9
Thinking now of how this might be implemented in a computer program, we’d want
to “rearrange” the system so that the largest coefficient of x
1
comes first
(1) −3x
1
+ 8x
2
+ 9x
3
=−9
(2) 2x
1
− 4x
2
+ 2x
3
= 18
(3) x
1
− 3x
2
− 2x
3
= 6
and then multiply (1) by 2 and (2) by 3 and add, yielding a new (2), and multiply (1)
by 1 and (3) by 3 and add, yielding a new (3)
2 ×(1): −6x
1
+ 16x
2
+ 18x
3
=−18
3 × (2): 6x
1
− 12x
2
+ 6x
3
= 54
Sum: 4x
2
+ 24x
3
= 36
1 × (1): −3x
1
+ 8x
2
+ 9x
3
=−9
3 × (3): 3x
1
− 9x
2
− 6x
3
= 18
Sum: −1x
2
+ 3x
3
= 9
which gives us a new system:
(1) x
1
− 3x
2
+ 2x
3
= 6
(2) 4x
2
+ 24x
3
= 36
(3) −1x
2
+ 3x
3
= 9
We can then multiply (3) by 4 and add it to (2), yielding a new (3)
1 × (2) 4x
2
+ 24x
3
= 36
4 ×(3) −4x
2
+ 12x
3
= 36
Sum: 36x
3
= 72
which we can then solve by back substitution for (1, −3, −2).
832 Appendix A Numerical Methods
The pseudocode of the basic algorithm is
// j indexes columns (pivots)
for (
j = 1 to n − 1)do
// Pivot step
find
l such that a
l,j
has the largest value among (a
j,j
, a
j+1,j
, ···, a
n,j
)
exchange rows l and j
// Elimination step
//
k indexes rows
for (
k = 1 to n)do
// Form multiplier
m =−
a
k,j
a
j,j
// Multiply row j by m and add to row k
for (i = j + 1 to n)do
a
k,i
= a
k,i
+ ma
k,j
// Multiply and add constant for row j
c
k
= c
k
+ mc
j
// Back substitute
x
n
=
c
n
a
n,h
for (k = n − 1 downto 1)do
x
k
=
c
k
−
n
m=k+1
a
k,m
x
m
a
k,k
A.2 Systems of Polynomials
The last section showed how to solve systems of linear equations. Given n equations in
m unknowns,
m
j=0
a
ij
x
j
=b
i
for 0 ≤i<n, let the system be represented in matrix
form by A x =
b,whereA = [a
ij
]isn × m, x = [x
j
]ism × 1, and
b = [b
i
]isn × 1.
The n × (m + 1) augmented matrix [A|
b] is constructed and row-reduced to [E|c], a
matrix that has the following properties:
Thefirstnonzeroentryineachrowis1.
If the first nonzero entry in row r is in column c, then all other entries in column
c are 0.
All zero rows occur last in the matrix.
If the first nonzero entries in rows 1 through r occur in columns c
1
through c
r
,
then c
1
<...<c
r
.
If there is a row whose first m entries are zero, but the last entry is not zero, then
the system of equations has no solution. If there is no such row, let ρ = rank([E|c])
A.2 Systems of Polynomials 833
denote the number of nonzero rows of the augmented matrix. If ρ = m, the system
has exactly one solution. In this case E = I
m
, the m × m identity matrix, and the
solution is x =c.Ifρ<m, the system has infinitely many solutions, the solution set
having dimension m −ρ. In this case, the zero rows can be omitted to obtain the ρ ×
(m +1) matrix [I
ρ
|F|c
+
], w he re I
ρ
is the ρ ×ρ identity matrix, F is ρ ×(m −ρ), and
c
+
consists of the first ρ entries of c.Letx be partitioned into its first ρ components
x
+
and its remaining m − ρ components x
−
. The general solution to the system is
x
+
=c
+
− Fx
−
, where the x
−
are the free parameters in the system.
Generic numerical linear system solvers for square systems (n = m) use row re-
duction methods so that (1) the order of time for the algorithm is small, in this case
O(n
3
), and (2) the calculations are robust in the presence of a floating-point num-
ber system. It is possible to solve a linear system using cofactor expansions, but the
order of time for the algorithm is O(n!), which makes this an expensive method for
large n.However,n = 3 for many computer graphics applications. The overhead for
a generic row reduction solver normally uses more cycles than a simple cofactor ex-
pansion, and the matrix of coefficients for the application is usually not singular (or
nearly singular) so that robustness is not an issue, so for this size system the cofactor
expansion is a better choice.
Systems of polynomial equations also arise regularly in computer graphics ap-
plications. For example, determining the intersection points of two circles in 2D is
equivalent to solving two quadratic equations in two unknowns. Determining if two
ellipsoids in 3D intersect is equivalent to showing that a system of three quadratic
equations in three unknowns does not have any real-valued solutions. Computing
the intersection points between a line and a polynomial patch involves setting up and
solving systems of polynomial equations. A method for solving such systems involves
eliminating variables in much the same way that you do for linear systems. How-
ever, the formal calculations have a flavor of cofactor expansions rather than row
reductions.
A.2.1 Linear Equations in One Formal Variable
To motivate the general idea, consider a single equation a
0
+ a
1
x = 0 in the variable
x.Ifa
1
= 0, there is a unique solution x =−a
0
/a
1
.Ifa
1
= 0 and a
0
= 0, there are no
solutions. If a
0
= a
1
= 0, any x is a solution.
Now consider two equations in the same variable, a
0
+a
1
x =0 and b
0
+b
1
x =0,
where a
1
=0 and b
1
=0. The first equation is multiplied by b
1
, the second equation is
multiplied by a
1
, and the two equations are subtracted to obtain a
0
b
1
−a
1
b
0
=0. This
is a necessary condition that a value x be a solution to both equations. If the condition
is satisfied, then solving the first equation yields x =−a
0
/a
1
. In terms of the row
reduction method for linear systems discussed in the last section, n = 2, m = 1, and
the augmented matrix with its reduction steps is
a
1
−a
0
b
1
−b
0
∼
a
1
b
1
−a
0
b
1
a
1
b
1
−a
1
b
0
∼
a
1
b
1
−a
0
b
1
0 a
0
b
1
− a
1
b
0
∼
1 −a
0
/a
1
0 a
0
b
1
− a
1
b
0