13.472J/1.128J/2.158J/16.940J
COMPUTATIONAL GEOMETRY
Lecture 19
Prof. N. M. Patrikalakis
Copyrightc 2003MassachusettsInstituteofTechnology
Contents
Decomposition models 2
19.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
19.2 Exhaustive enumeration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
19.2.1 De nition and construction methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
19.2.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
19.2.3 Properties of exhaustive enumeration methods . . . . . . . . . . . . . . . . . . . . . . . 5
19.3 Space subdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
19.3.1 Motivation and de nitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
19.3.2 Construction of octrees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
19.3.3 Algorithms for octrees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
19.3.4 Properties of octrees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
19.3.5 Binary space subdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
19.4 Cell decompositions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
19.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
19.4.2 Cell tuple data structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
19.4.3 Properties of cell decompositions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Integral properties of geometric models 14
19.5 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
19.6 Integral properties of curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
19.6.1 Planar curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
19.6.2 3D curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
19.7 Integral properties of surface patches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
19.7.1 Planar regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
19.7.2 Curved surface patch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
19.8 Solids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
19.9 Example: solid of revolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
19.10Appendix: Review of numerical integration methods . . . . . . . . . . . . . . . . . . . . . 25
19.10.1 Trapezoidal rule of integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
19.10.2 Simpson’s rule of integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
19.10.3 Romberg integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
19.10.4 Double integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
Bibliography 31
1
Decomposition models
19.1 Introduction
Decomposition models are representations of solids via combinations (unions) of basic special
building blocks glued together. Alternatively, decomposition models may be considered to
represent solids in terms of a subdivision of space (see also Lecture 1 for more details on the
classi cation of these models). Various types of decomposition models are created by:
various building blocks
various combination methods used to create the model.
In order of increasing complexity, decomposition models are classi ed as follows:
1. Exhaustive enumeration
2. Space subdivision
3. Cell decomposition
2
19.2 Exhaustive enumeration
19.2.1 De nition and construction methods
Exhaustive enumeration is a representation by means of nonoverlapping cubes of uniform size
and orientation, see Figure 19.1. An object is represented by a three dimensional Boolean array.
Each cell represents a cubic volume of space. If a cell intersects with the region of interest it
has a true value. Otherwise, the value is false. This can be pictured as a box divided into 3D
cubical pixels, with 0 assigned if empty and 1 assigned if full. This representation involves:
A regular subdivision of 3D space within a cube of given size which is partitioned and
oriented in a certain way within a global coordinate system. The subdivision is made up
of sub-cubes (3D pixels) of given size. Reference and access to each sub-cube is made by
three integer indices i; j; k.
For xed space of interest we need a 3-D array, Cijk of binary data:
Cijk =
( 1 if the sub-cube i; j; k intersects the solid
0 if the sub-cube i; j; k is empty (19.1)
Construction of exhaustive enumeration models requires an alternate representation or
measurements (eg. digital tomograghy, medical scanning, sonar data, acoustic tomography
data, etc). Usually the primary data type for such construction is a B-Rep or a CSG model
or another exhaustive enumeration model at di erent resolution, and cube location and orien-
tation.
Operations on exhaustive enumeration models are easy. Boolean operations for example
(especially for models within the same cube at the same resolution) are direct. Similarly
visualization and integral computations are very easy. However, for higher quality rendering,
ltering methods to estimate accurate surface normals may be involved [16].
The binary matrix (19.1) typically represents a valid solid. However disconnected cells or
cells with low degree of connectivity as in Figure 19.2 are undesirable. For the results of Boolean
operations, ltering may be needed to maintain connectivity of cells. Strict connectivity occurs
when each full cell has at least one full neighbor across a face.
19.2.2 Applications
Applications of exhaustive enumeration methods include:
Underwater environment representation.
Finite element meshing ( rst step in an algorithm to build such a mesh).
Medical 3D data representation.
Preprocessing representation for speeding up operations on other representations (eg.
approximating integral properties such as volume, center of gravity, moments of inertia).
3
Figure 19.1: Exhaustive enumeration.
Figure 19.2: Various connectivities of cells in exhaustive enumeration models
4
19.2.3 Properties of exhaustive enumeration methods
Properties of exhaustive enumeration methods include:
Expressive power: these methods are an approximation scheme and do not form a primary
representation, especially within CAD/CAM applications.
Unambiguity and uniqueness for xed space (size, location and orientation of primary
cube) and resolution. There do not exist di erent representations for the same object
under these conditions (which is not true of many other representation methods such as
B-rep or CSG described before).
Memory intensive: eg. for a linear resolution of 256, 2563 integer elements for Cijk in
equation 19.1 leading to 16M bits and this is a bare minimum.
Closure 1 of operations (eg. Booleans).
Computational ease for algorithms: VLSI implementation for volume rendering is possi-
ble. However, for high resolution the algorithm slows down.
1Closure means that an operation such as Boolean results in an object that can be represented by the same
type of data structure.
5
1 2
3
4
empty
1
2 3
4
full
partially full
Figure 19.3: Quadtree representation.
19.3 Space subdivision
19.3.1 Motivation and de nitions
Some of the motivations behind space subdivision methods include:
Exhaustive enumeration is memory intensive and typically has low accuracy.
Smaller memory requirements are possible, if adaptive subdivision is used;
Octree/quadtree representations lead to a recursive subdivision into 8 octants (or 4 quad-
rants) that can be represented as an 8-ary tree (or 4-ary tree) for which e cient algo-
rithms are also known.
In an octree representation a solid region is represented by hierarchically decomposing a usually
cubic volume of space into successively smaller cubes (8 of them). Hierarchical division and
cube orientation usually follows the spatial coordinate system. An example of quadtree, the
two dimensional analogue, is shown Figure 19.3.
This is a trivial example. The method can continue to many more levels for a much more
complex model. Some tolerance for the minimum size block is required. In addition, this
very concise representation would become very large if the coordinate system was changed; for
example, rotated 45 degrees.
This method leads to a quick way to compute the area and other integral properties of
a region. It is often used in data analysis in elds such as medical applications and sonar
imaging.
6
19.3.2 Construction of octrees
To create an octree, we apply a classi cation procedure to a given solid (represented using the
Boundary Representation, Constructive Solid Geometry, or Exhaustive Enumeration methods,
etc.) and decide if a given node of the octree is:
Exterior to solid (white);
Interior to solid (black);
Partially interior to solid (grey).
The classi cation procedure is used recursively. It is based on Boolean solid operations,
especially intersection. Figure 19.4 provides an simple example of octree representation.
In general, the decision if a given node of the octree is white, black, or grey is not an
easy task. For the simple case of a convex solid object, it is su cient to classify the eight
vertices of the given node of the octree (which is a cube) with respect to the solid. This
can be accomplished by for example casting a half-in nite ray from the point intersecting
the solid’s surfaces in a number of (multiplicity one) intersection points. If the number of
such intersection points is even/odd, the point is outside/in (or on the surface of the solid).
However, for a concave solid object, classi cation of the six faces of the cube with respect to
the solid is necessary, see Figure 19.5 for an illustration in the 2-D case. This requires surface
intersections with a planar patch. The memory and processing computation required for a 3-D
object is on the order of the surface area of the object [12] [9]. Depending on the object and
the resolution, this can still represent a large storage requirement.
19.3.3 Algorithms for octrees
Various algorithms for octrees are developed in Meagher [12] and are summarized here:
1. Tree generation or conversion from other representation methods were discussed above
in Section 19.3.2.
2. Set operators (union, intersection, di erence): A low resolution tree could be an e ective
preprocessor for a B-rep model in processes like interference checking.
3. Geometric transformations (translation, rotation, scaling).
4. Analysis procedures (integral, volume properties, connected components).
5. Rendering [16].
As an example we consider set or Boolean operations. Set operations lead to simple tree
traversal
Intersection(Tree A, Tree B) = Tree C
Trees are traversed in synchronous fashion and a case analysis for the types of nodes is
performed. We use the terms \black"= in-solid,\white" = out-of-solid. At each level of
subdivision there are three cases [11]:
7
Figure 19.4: An octree model.
Figure 19.5: Classi cation of an quadtree node with respect to a concave object
8
1. If nodes n1 and n2 are both leaves, then the resulting node n3 is black if n1 and n2 are
both black; otherwise n3 is white.
2. Either n1 or n2 is a leaf. If the leaf node is black, then the subtree of the non{leaf node
is copied to the resultant octree; otherwise the resulting node is white.
3. If nodes n1 and n2 are non-leaves, then the algoritm considers recursively their children
as above.
The complexity of such an intersection algorithm is proportional to the size of the smaller
tree.
Not all algorithms are in the form of a simple tree traversal. Some algorithms may require at
worst, a traversal up to the root and back down to a neighbor. Examples of such algorithms are
surface rendering (shading), transparency rendering, and extraction of connected components.
19.3.4 Properties of octrees
Some of the properties of octrees include:
Expressive power: they are an approximation scheme and do not form a primary repre-
sentation.
Validity: if no special connectivity is required, all octrees are valid.
Unambiguity and uniqueness: for a xed resolution there is only one compacted2 octree;
Memory requirements: not as large as exhaustive enumeration models, yet typically much
larger than Boundary Representation and Constructive Solid Geometry models;
Closure of operations: for example for Boolean operations and transformations;
Computational ease: octrees are somewhat more complex than exhaustive enumeration.
2Algorithms such as set operations can create octrees with unnecessary nodes (eg. an internal nodes whose
children are all black). Such nodes can be removed with a relatively simple tree traversal algorithm.
9
19.3.5 Binary space subdivision
Figure 19.6: Binary subdivision tree.
Beyond octrees, an alternative type of tree representation involves dividing nodes into 2
rather than 8 components, see M antyl a [11] and Figure 19.6. Subdivisions are performed in the
X, Y, and Z coordinate directions sequentially. Binary trees are typically somewhat smaller
than octrees and they can be converted to linear arrays containing special symbols [11].
10
19.4 Cell decompositions
19.4.1 Motivation
The motivation for cell decomposition methods is:
Use of elements other than cubes, see Figure 19.7 for an example.
Applications such as design of inhomogeneous (eg. composites) and functionally graded
materials, nite element analysis methods, scienti c visualization of scalar and vector
elds.
Cells are parametrized instances of a generic cell type, eg. a cell bounded by quadratic
curves and surfaces.
Cells are homeomorphic to spheres.
Cells meet at a vertex, edge, face otherwise the representation is invalid.
Cells are disjoint and non-overlapping.
Cells may belong to di erent cell types, eg. box-like, tetrahedra-like, etc.
Figure 19.7: A cell decomposition ( nite element mesh).
19.4.2 Cell tuple data structure
A cell decomposition can be represented using the cell-tuple data structure [2] which applies for
n-D models, see also [1] for a review and summary of other related data structures such as the
Quad-edge structure [6] for 2D models and the Facet-edge pair structure [4, 5] for 3-D models.
Figures 19.8 and 19.9 present 2D and 3D examples. This data structure can be mapped into
a relational database or a graph structure.
11
Figure 19.8: Cell data structure for a 2D model
Volume
Face
Edge
Vertex
Figure 19.9: Cell data structure for a 3D model
12
19.4.3 Properties of cell decompositions
The properties of cell decomposition methods are:
Expressive power: they are very general and accurate, not necessarily requiring approxi-
mations;
Validity: they require an intersection test for veri cation;
Unambiguity: they provide an unambiguous representation;
Nonuniqueness: Similarly to the Constructive Solid Geometry method, the same object
can be represented at di erent resolutions or with di erent types of mesh (eg. hexahedral,
tetrahedral, etc.);
Generation: It typically is done by conversion from other representations with or without
geometric approximation;
Conciseness: memory utilization is less than octrees, yet more than Boundary Represen-
tation;
Applicability: nite element meshing, multimaterial non-homogeneous objects, visualiza-
tion of elds, etc.
13
Integral properties of geometric models
19.5 Introduction
One of the important advantages of using a CAD model for representing and designing an
object is that we can easily compute the integral properties of such models such as edge curves,
faces and volumes. Integral properties include length, area, centroid, moment of inertia, and
volume. These are very useful in preliminary design. For example, surface area a ects drag,
volume a ects the carrying capacity of a vehicle, centroids are useful in hydrostatic balance,
moments of inertia are used in dynamics and in hydrostatic stability calculation (for ships).
Computation of the integral properties of curves, surface patches and solids involves eval-
uation of single, double and triple integrals of the form
curve =
Z
curve
f(P)dL; surface =
Z
surface
f(P)dS; solid =
Z
solid
f(P)dV (19.2)
where is the required property, P is a point and f is a real-valued function, which depends on
the type of property required. We have studied three classes of solid representation methods in
the previous chapters, namely Decomposition methods, Constructive Solid Geometry (CSG)
methods, and Boundary Representation (B-rep) methods.
For decomposition methods, the integral over the solid reduces to a sum of integrals
Z
solid
fdV = X
i
Z
celli
fdV (19.3)
where celli is the i th cell which is either full or partially full. For the case of exhaustive
spatial enumeration the cells are constant-sized cubes and for the octree decompositions they
are variable-sized cubes [10, 13]. For high resolution models it is enough to consider all the
celli to be full entirely and for those cases the resulting integrals are elementary and can be
computed using simple analytic forms.
As we have studied in the previous chapters, CSG is a tree whose nodes represent the
Boolean operators and the leaves are the primitive solids. Therefore the computation of integral
properties of CSG solids consists of applying the following formula recursively [10]:
Z
A[B
fdV =
Z
A
fdV +
Z
B
fdV
Z
A\B
fdV (19.4)
Z
A B
fdV =
Z
A
fdV
Z
A\B
fdV: (19.5)
14
Consequently we need to compute integrals over primitives (which can be evaluated analyti-
cally) and integrals involving intersections of primitives RA\B fdV which can be approximated
using a ray casting, ray classi cation and integral approximation method.
Boundary representation, which is the most generally used representation today, represents
the object in terms of their boundary elements (e.g. vertices, edges, faces). For evaluating the
integral properties for B-rep solids, the following theorems from vector calculus are useful [7]:
1. Green’s Theorem
If C is a piecewise smooth, simple closed curve that bounds a region R, and if P(x; y)
and Q(x; y) are continuous functions which have continuous partial derivatives along C
and throughout R, then
I
C
(Pdx + Qdy) =
ZZ
R
@Q
@x
@P
@y
!
dA (19.6)
2. Divergence Theorem (also called Gauss’ Theorem)
The ux of vector eld F owing outward through a closed surface S equals the integral
of the divergence of F over the region R bounded by S;
ZZ
S
F ndA =
ZZZ
R
rFdV (19.7)
where n is the outward unit normal vector and
r F = @F@x i + @F@y j + @F@z k: (19.8)
where i, j, k are the unit coordinate vectors.
In the sequel we will apply these theorems to compute the integral properties of geometric
models represented by the B-rep method (in 1-3 dimensions).
19.6 Integral properties of curves
19.6.1 Planar curves
Let a planar curve be de ned by
r = (x(t); y(t)); to t t1 (19.9)
Length
L =
Z s1
s0
ds =
Z t1
t0
q
_r(t) _r(t)dt (19.10)
=
Z t1
t0
q
_x2(t) + _y2(t)dt
15
\Centroid"
rc = (xc; yc) =
Rs1
s0 rdsR
s1
s0 ds
(19.11)
= 1L
Z t1
t0
r(t)
q
_x2(t) + _y2(t)dt
\Moments of inertia"
Ixx =
Z s1
s0
y2ds =
Z t1
t0
y2(t)
q
_x2(t) + _y2(t)dt (19.12)
Iyy =
Z s1
s0
x2ds =
Z t1
t0
x2(t)
q
_x2(t) + _y2(t)dt (19.13)
Ixy =
Z s1
s0
xyds =
Z t1
t0
x(t)y(t)
q
_x2(t) + _y2(t)dt (19.14)
19.6.2 3D curves
Let a 3D curve be de ned by
r = (x(t); y(t); z(t)) ; to t t1 (19.15)
Length
L =
Z s1
s0
ds =
Z t1
t0
q
_r(t) _r(t)dt (19.16)
=
Z t1
t0
q
_x2(t) + _y2(t) + _z2(t)dt
\Centroid"
rc = (xc; yc; zc) =
Rs1
s0 rdsR
s1
s0 ds
(19.17)
= 1L
Z t1
t0
r(t)
q
_x2(t) + _y2(t) + _z2(t)dt
\Moments of inertia"
Ixx =
Z s1
s0
(y2 + z2)ds =
Z t1
t0
y2(t) + z2(t)
q
_x2(t) + _y2(t) + _z2(t)dt (19.18)
Iyy =
Z s1
s0
(x2 + z2)ds =
Z t1
t0
x2(t) + z2(t)
q
_x2(t) + _y2(t) + _z2(t)dt (19.19)
Izz =
Z s1
s0
(x2 + y2)ds =
Z t1
t0
x2(t) + y2(t)
q
_x2(t) + _y2(t) + _z2(t)dt (19.20)
Ixy =
Z s1
s0
xyds =
Z t1
t0
x(t)y(t)
q
_x2(t) + _y2(t) + _z2(t)dt (19.21)
Iyz =
Z s1
s0
yzds =
Z t1
t0
y(t)z(t)
q
_x2(t) + _y2(t) + _z2(t)dt (19.22)
Ixz =
Z s1
s0
xzds =
Z t1
t0
x(t)z(t)
q
_x2(t) + _y2(t) + _z2(t)dt (19.23)
16
19.7 Integral properties of surface patches
19.7.1 Planar regions
Let us consider a planar region as in Figure 19.10
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a1a0a2a0a1a0a1a0a2a0a1a0a1a0a2a0
C:r(t)
A
Figure 19.10: Planar region A
Area
A =
ZZ
A
dA (19.24)
Using Green’s theorem with Q(x; y) = x and P(x; y) = y, then
1
2(Qx Py) =
1
2(1 + 1) = 1 (19.25)
where subscripts x, y denote partial derivatives. We can rewrite equation (19.24) using
Green’s theorem as
A =
ZZ
A
1
2(Qx Py)dxdy
= 12
I
C
(Pdx + Qdy)
= 12
I
C
( ydx + xdy)
= 12
I
C
(x _y y _x)dt (19.26)
If x(t); y(t) are piecewise polynomial functions, the above integral can be evaluated from
a symbolic/analytic integration formula but this is typically tedious. By contrast, nu-
merical integration methods may be used more easily (see Appendix).
Centroid
rc = (xc; yc) =
RR
A rdARR
A dA
;
17
where A is the shaded area, and
ZZ
A
rdA =
ZZ
A
xdA;
ZZ
A
ydA
: (19.27)
Let Q(x; y) = x22 ; P(x; y) = 0, then
Qx Py = x 0 = x (19.28)
Therefore
ZZ
A
xdA =
I
C
Pdx + Qdy =
I
C
x2dy (19.29)
where C is the complete boundary of A. Thus,
xc =
RR
A xdARR
A dA
= 1A
ZZ
A
xdA
= 1A
I
C
x2 _ydt (19.30)
Let Q(x; y) = 0; P(x; y) = y22 , then
Qx Py = 0 + y = y (19.31)
Similarly,
yc =
RR
A ydARR
A dA
= 1A
ZZ
A
ydA
= 1A
I
C
y2 _xdt (19.32)
Moments of inertia
1. Ixx = RRA y2dA
Let Q(x; y) = 0; P(x; y) = y33 , then
Qx Py = y2 (19.33)
Using Green’s Theorem,
Ixx =
ZZ
A
(Qx Py)dxdy =
I
C
y3dx
=
I
C
y3 _xdt (19.34)
18
2. Iyy = RRA x2dA
Let Q(x; y) = x33 ; P(x; y) = 0, then
Qx Py = x2 (19.35)
Using Green’s Theorem,
Iyy =
ZZ
A
(Qx Py)dxdy = 14
I
C
x3
3 dy
=
I
C
x3 _ydt (19.36)
3. Ixy = RRA xydA
Let Q(x; y) = x2y2 ; P(x; y) = 0, then
Qx Py = xy + 0 = xy (19.37)
Using Green’s Theorem,
Ixy =
ZZ
A
(Qx Py)dxdy
= 12
I
C
x2y _ydt (19.38)
If x(t); y(t) are piecewise polynomial functions, the above integrals can be evaluated
from a symbolic/analytic integration formula but this is typically tedious. By contrast,
numerical integration methods may be used more easily (see Appendix).
19.7.2 Curved surface patch
Let us consider a curved surface patch r = r(u; v), with (u; v) 2 A, where A is a given
parametric domain, as in Figure 19.11.
r(u,v)
v u
Figure 19.11: Curved surface patch
19
Area
A =
ZZ
A
dA
=
ZZ
A
jru rvjdudv =
ZZ
A
pEG F 2dudv (19.39)
where E; F and G are the rst fundamental form coe cients E = ru ru; F = ru rv; G =
rv rv. (see Chapter 2).
\Centroid"
rc = (xc; yc; zc) =
RR
A rdARR
A dA
= 1A
ZZ
A
[x(u; v); y(u; v); z(u; v)]pEG F 2dudv (19.40)
\Moments of inertia"
Ixx =
ZZ
A
[y2(u; v) + z2(u; v)]pEG F 2dudv (19.41)
Iyy =
ZZ
A
[x2(u; v) + z2(u; v)]pEG F 2dudv (19.42)
Izz =
ZZ
A
[x2(u; v) + y2(u; v)]pEG F 2dudv (19.43)
Ixy =
ZZ
A
[x(u; v)y(u; v)]pEG F 2dudv (19.44)
Ixz =
ZZ
A
[x(u; v)z(u; v)]pEG F 2dudv (19.45)
Iyz =
ZZ
A
[y(u; v)z(u; v)]pEG F 2dudv (19.46)
Integrals 19.39-19.46 may be evaluated numerically as in the Appendix.
19.8 Solids
For solids described by the B-rep method it is convenient to transform volume integrals into
surface integrals by means of the divergence theorem.
Volume
V =
ZZZ
V
dV (19.47)
Choose
r = xi + yj + zk (19.48)
then
r r = @r@x i + @r@y j + @r@z k = 3 (19.49)
20
Using the Divergence (or Gauss’) Theorem,
V =
ZZZ
V
dV = 13
ZZZ
V
r rdV
= 13
ZZ
A
r ndA = 13
ZZ
A
r njru rvjdudv
= 13
ZZ
A
r (ru rv)dudv (19.50)
given that
n = ru rvjr
u rvj
(19.51)
Centroid
rc = (xc; yc; zc) =
RRR
V rdVRRR
V dV
(19.52)
Choose
r = 12x2i (19.53)
then
r r = @r@x i + @r@y j + @r@z k = x (19.54)
xc = 1V
ZZZ
V
xdV = 1V
ZZZ
V
r rdV
= 1V
ZZ
A
1
2x
2(i n)dA
= 1V
ZZ
A
1
2x
2(i (ru rv))dudv (19.55)
Similarly, expressions are obtained for yc; zc:
yc = 1V
ZZ
A
1
2y
2(j (ru rv))dudv (19.56)
zc = 1V
ZZ
A
1
2z
2(k (ru rv))dudv (19.57)
Moments of inertia
Ixx =
ZZZ
V
(y2 + z2)dV (19.58)
Choose
r = (y2 + z2)xi (19.59)
21
then
r r = @r@x i + @r@y j + @r@z k = y2 + z2 (19.60)
Thus
Ixx =
ZZZ
V
(y2 + z2)dV =
ZZZ
V
r rdV
=
ZZ
A
(y2 + z2)x(i n)dA
=
ZZ
A
(y2 + z2)x(i (ru rv))dudv (19.61)
Similarly,
Iyy =
ZZ
A
(x2 + z2)y(j (ru rv))dudv (19.62)
Izz =
ZZ
A
(x2 + y2)z(k (ru rv))dudv (19.63)
Ixy =
ZZ
A
xyz(k (ru rv))dudv (19.64)
Ixz =
ZZ
A
xzy(j (ru rv))dudv (19.65)
Iyz =
ZZ
A
yzx(i (ru rv))dudv (19.66)
19.9 Example: solid of revolution
r(t)=x(t)i+y(t)j
x
yy
z
θ
end1 end2
x(t1) x(t2)
Figure 19.12: Solid of revolution
Let a solid of revolution be de ned by
r(t; ) = x(t)i + y(t) sin j + y(t) cos k; (19.67)
t1 t t2; 0 2 :
and we assume that the two end caps are closed o with planar disks (see Figure 19.12).
22
Surface area of surface of revolution (with end caps)
rt = (_x; _y sin ; _y cos ) (19.68)
r = (0; y cos ; y sin ) (19.69)
rt r =
i j k
_x _y sin _y cos
0 y cos y sin
= ( y _y sin2 y _y cos2 )i + _xy sin j + _xy cos k
= y _yi + _xy sin j + _xy cos k (19.70)
jrt r j =
q
y2 _y2 + _x2y2 = y
q
_x2 + _y2 (19.71)
A =
ZZ
A
dA =
Z t2
t1
Z 2
0
y
q
_x2 + _y2d dt
= 2
Z t2
t1
y
q
_x2 + _y2dt (19.72)
Volume
rt r = ( y _y; _xy sin ; _xy cos ) (19.73)
r (rt r ) = (x; y sin ; y cos ) ( y _y; _xy sin ; _xy cos )
= xy _y + _xy2 sin2 + _xy2 cos2
= xy _y + _xy2 (19.74)
V = 13
Z t2
t1
Z 2
0
r (rt r )dtd 13
ZZ
Aend1
xdxdy + 13
ZZ
Aend2
xdxdy
= 13
Z t2
t1
Z 2
0
( xy _y + _xy2)d dt 3x(t1)y2(t1) + 3x(t2)y2(t2)
= 2 3
Z t2
t1
( xy _y + _xy2)dt 3x(t1)y2(t1) + 3x(t2)y2(t2) (19.75)
Using integration by parts,
Z t2
t1
xy _ydt = [xy2]t2t1 +
Z t2
t1
(xy)0ydt
= x(t2)y2(t2) + x(t1)y2(t1) +
Z t2
t1
( _xy2 + xy _y)dt (19.76)
Thus
2
Z t2
t1
xy _ydt = x(t2)y2(t2) + x(t1)y2(t1) +
Z t2
t1
_xy2dt (19.77)
23
The volume, therefore, is
V = 2 3
Z t2
t1
_xy2dt 2 3
Z t2
t1
( xy _y)dt 3x(t1)y2(t1) + 3x(t2)y2(t2)
=
Z t2
t1
_xy2dt (19.78)
(corroborating the obvious formula from elementary calculus)
Centroid
1. i (rt r ) = y _y
xc = 1V
Z t2
t1
Z 2
0
1
2x
2( y _y)dtd
2V x
2(t1)y2(t1) +
2V x
2(t2)y2(t2)
= V
Z t2
t1
( x2y _y)dt + 2V (x2(t2)y2(t2) x2(t1)y2(t1)) (19.79)
Integrate by parts
Z t2
t1
x2y _ydt = [x2y2]t2t1 +
Z t2
t1
(x2y)0ydt
= x2(t2)y2(t2) + x2(t1)y2(t1) +
Z t2
t1
(2x _xy + x2 _y)ydt (19.80)
Thus
2
Z t2
t1
x2y _ydt = x2(t2)y2(t2) + x2(t1)y2(t1) + 2
Z t2
t1
x _xy2dt (19.81)
and
xc = V
Z t2
t1
x _xy2dt (19.82)
(corroborating the obvious formula from elementary calculus)
2. j (rt r ) = _xy sin
yc = 1V
ZZ
A
1
2y
2 _xy sin dtd
= 12V
Z t2
t1
Z 2
0
y3 _x sin d dt
= 12V
Z t2
t1
y3 _x[ cos ]2 0 dt
= 0 (19.83)
3. k (rt r ) = _xy cos
zc = 1V
ZZ
A
1
2z
2 _xy cos dtd
= 12V
Z t2
t1
Z 2
0
_xyz2 cos d dt
= 12V
Z t2
t1
_xyz2[sin ]2 0 dt
= 0 (19.84)
24
19.10 Appendix: Review of numerical integration meth-
ods
19.10.1 Trapezoidal rule of integration
Here we compute Rba f(x)dx as in Figure 19.13 by decomposing the integral into n subintervals
[3]. The area over each subinterval is approximated by the area of the trapezoid. The integral
is then obtained by the sum of areas over all the subintervals.
Area in each subinterval [xi; xi+1]
Z xi+1
xi
f(x)dx (f(xi) + f(xi+1)) x2 = h2 (fi + fi+1) (19.85)
where x = xi+1 xi.
[a; b] subdivided into n-subintervals
Z b
a
f(x)dx =
nX
i=1
h
2(fi + fi+1) =
h
2(f1 + 2f2 + 2f3 + + 2fn + fn+1) (19.86)
x1=a x2 x3 x4 ... xn xn+1=b
f(x)
Figure 19.13: Trapezoidal rule of integration
Error in Trapezoidal rule [3]
Local error h
3
12f
00( ); xi < < xi+1
Global error h
2
12(b a)f
00( ); a < < b
25
19.10.2 Simpson’s rule of integration
For the trapezoidal rule, we approximate the curve by n straight line segments, while in
Simpson’s rule of integration we approximate it with a parabola in a piecewise manner [3].
y = ax2 + bx + c: (19.87)
Therefore, this rule requires the use of three base points.
Area in two subintervals [ h; h]
Let the three points be x = h; 0; h, and the corresponding y values be f1; f2, and f3,
see Figure 19.14. Then,
f1 = ah2 bh + c (19.88)
f2 = c (19.89)
f3 = ah2 + bh + c (19.90)
If we solve for a; b; and c, we obtain
?h 0 x
y
f2
f1 f3
h
Figure 19.14: Simpson’s rule
a = f1 2f2 + f32h2 (19.91)
b = f3 f12h (19.92)
c = f2 (19.93)
Thus,
Z h
h
(ax2 + bx + c)dx = a3[x3]h h + b2[x2]h h + c[x]h h
= 2a3 h3 + 2ch
= h3[f1 + 4f2 + f3] (19.94)
Area in [a; b] subdivided into n-subintervals
Z b
a
f(x)dx = h3(f1 + 4f2 + 2f3 + 4f4 + 2f5 + + 2fn 1 + 4fn + fn+1) (19.95)
26
Error in Simpson’s rule [3]
Local error h
5
90f
(4)( ); xi < < xi+1
Global error b a180 h4f(4)( ); a < < b
19.10.3 Romberg integration
Let us compute the integration of f(x) using the trapezoidal rule over the interval [a; b] with
x = h. If we denote the output of the trapezoidal rule as T0;1, then
True value = T0;1 + O(h2) (19.96)
Now let us assume that O(h2) = C1h2, where C1 is constant. Then,
True value = T0;1 + C1h2 (19.97)
If we double the number of subintervals such that x = h2, then,
True value T1;1 + C1
h
2
!2
(19.98)
There are two unknowns in equation (19.97) and (19.98), True value and constant C1. Sub-
tracting (19.97) from four times (19.98) yields
True value T0;2 T1;1 + 13(T1;1 T0;1) (19.99)
Similarly, we can obtain for x = h4
True value T2;1 + C116h2 (19.100)
From (19.98) and (19.100), we obtain
T1;2 = T2;1 + 13(T2;1 T1;1) (19.101)
We can make a further improvement by using T0;2 and T1;2 and setting up the relations
True value = T0;2 + C2h4 (19.102)
True value T1;2 + C2
h
2
!4
(19.103)
and hence
T0;3 = T1;2 + 115(T1;2 T0;2) (19.104)
27
T0;1 T0;2 T0;3
T1;1 T1;2
T2;1
We can arrange procedures in a matrix form as
We can also summarize the rule as
Improved value = more accurate +
1
2n 1
(more accurate less accurate)
where n is the exponent on h in the error term O(hn).
Examples: R10 cos(x)dx = sin(1) = 0:841470985
T0;1 = 12(cos(0) + cos(1)) = 0:770151153
T1;1 = 12T0;1 + 12 cos(0:5) = 0:823866857
T2;1 = 12T1;1 + 14(cos(0:25) + cos(0:75)) = 0:837083751
T0;2 = T1;1 + 13(T1;1 T0;1) = 0:841651966
T1;2 = T2;1 + 13(T2;1 T1;1) = 0:841489382
T0;3 = T1;2 + 116(T1;2 T0;2) = 0:841479221
which is accurate to 5 signi cant digits.
19.10.4 Double integrals
ZZ
A
f(x; y)dA (19.105)
1. Over a rectangular domain
x
y
a b
c
d
I =
ZZ
A
f(x; y)dA =
Z b
a
Z d
c
f(x; y)dy
!
dx =
Z d
c
Z b
a
f(x; y)dx
!
dy (19.106)
28
We will use the trapezoid rule in x and y directions. Let
hx = b an (19.107)
hy = d cn (19.108)
we start with y = c,
y = c : I1 =
Z b
a
f(x; c)dx
= hx2 (f1(c) + 2f2(c) + + 2fn(c) + fn+1(c)) (19.109)
where fi(c) = f(xi; c). Similarly,
y = c + hy : I2 =
Z b
a
f(x; c + hy)dy
= hx2 (f1(c + hy) + 2f2(c + hy) + + 2fn(c + hy) + fn+1(c + hy))
(19.110)
y = d : In+1 =
Z b
a
f(x; d)dy
= hx2 (f1(d) + 2f2(d) + + 2fn(d) + fn+1(d)) (19.111)
We now sum I1; I2; ; In+1 in y direction in terms of the trapezoidal rule
I = hy2 (I1 + 2I2 + + 2In + In+1) (19.112)
2. Over a curved boundary domain, see Figure 19.15
First we need to nd ai; bi. We equate y(t) with c+(i 1)hy, where c = min y(t). Then
we solve for t, leading to two (or more) intersections. Plugging the resulting t into x(t)
yields ai and bi.
Similar to the rectangular domain case, we use the trapezoidal rule in both x and y
directions.
Let
hix = bi ain (19.113)
hy = d cn (19.114)
29
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a5a3a4a3a4a3a4a3a4a3a5a3
x
y
a1a2
an+1
b1 b2
bn+1
C(t)=(x(t),y(t))
d
c
Figure 19.15: A curved boundary domain
We start with y = c:
y = c : I1 =
Z b1
a1
f(x; c)dx
= h
1
x
2 (f1(c) + 2f2(c) + + 2fn(c) + fn+1(c)) (19.115)
y = c + hy : I2 =
Z b2
a2
f(x; c + hy)dx
= h
2
x
2 (f1(c + hy) + 2f2(c + hy) + + 2fn(c + hy) + fn+1(c + hy))
(19.116)
y = d : In+1 =
Z bn+1
an+1
f(x; d)dx
= h
n+1
x
2 (f1(d) + 2f2(d) + + 2fn(d) + fn+1(d)) (19.117)
These formulae can be extended to curved domain boundaries bounding multiply con-
nected domains.
30
Bibliography
[1] L. Bardis and N. M. Patrikalakis. Topological structures for generalized boundary repre-
sentations. MITSG 94-22, MIT Sea Grant College Program, Cambridge, Massachusetts,
September 1994.
[2] E. Brisson. Representing geometric structures in d dimensions: Topology and order.
Discrete and Computational Geometry, 9:387{426, 1993.
[3] G. Dahlquist and A. Bj orck. Numerical Methods. Prentice-Hall, Inc., Englewood Cli s,
NJ, 1974.
[4] D. P. Dobkin and M. J. Laszlo. Primitives for the manipulation of three-dimensional
subdivisions. In Proceedings of the Third ACM Symposium on Computational Geometry,
pp. 86{99, Waterloo, Canada, June 1987.
[5] D. P. Dobkin and M. J. Laszlo. Primitives for the manipulation of three-dimensional
subdivisions. Algorithmica, 4:3{32, 1989.
[6] L. Guibas and J. Stol . Primitives for the manipulation of general subdivisions and the
computation of Voronoi diagrams. ACM Transactions on Graphics, 4(2):74{123, April
1985.
[7] F. B. Hildebrand. Advanced Calculus for Applications. Prentice-Hall, Inc., Englewood
Cli s, New Jersey, 1976.
[8] C. M. Ho mann. Geometric and Solid Modeling: An Introduction. Morgan Kaufmann
Publishers, Inc., San Mateo, California, 1989.
[9] G. M. Hunter and K. Steiglitz. Operations on images using quad trees. IEEE Transactions
on Pattern Analysis and Machine Intelligence, 1(2):145{153, 1979.
[10] Y. T. Lee and A. A. G. Requicha. Algorithms for computing the volume and other integral
properties of solid objects, I: Known methods and open issues. Communications of the
ACM, 25(9):635{641, September 1982.
[11] M. M antyl a. An Introduction to Solid Modeling. Computer Science Press, Rockville,
Maryland, 1988.
[12] D. Meagher. Geometric modeling using octtree encoding. Computer Graphics and Image
Processing, 19:129{147, June 1982.
31
[13] M. E. Mortenson. Geometric Modeling. John Wiley and Sons, New York, 1985.
[14] A. A. G. Requicha. Representations of solid objects - theory, methods and systems. ACM
Computing Surveys, 12(4):437{464, December 1980.
[15] V. Shapiro. Solid modeling. In G. Farin et al., editor, Handbook of Computer Aided
Geometric Design, Chapter 20, pp. 473{518. Elsevier, Amsterdam, 2002.
[16] R. Yagel, D. Cohen, and A. Kaufman. Context sensitive normal estimation for volume
imaging. In N. M. Patrikalakis, editor, Scienti c Visualization of Physical Phenomena,
pp. 211{232. Tokyo: Springer-Verlag, 1991.
32