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