Decomposition Models

Lecture 19 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 classification 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 classified as follows: 1. Exhaustive enumeration 2. Space subdivision 3. Cell decomposition 19.2 Exhaustive enumeration 19.2.1 Definition 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 fixed 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 different resolution, and cube location and orientation. 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, filtering 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, filtering 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 (first 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). Figure 19.1: Exhaustive enumeration. Figure 19.2: Various connectivities of cells in exhaustive enumeration models 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 fixed space (size, location and orientation of primary cube) and resolution. There do not exist different 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 possible. However, for high resolution the algorithm slows down. 1 Closure means that an operation such as Boolean results in an object that can be represented by the same type of data structure. 1 2 4 1 3 2 3 4 full empty partially full Figure 19.3: Quadtree representation. 19.3 Space subdivision 19.3.1 Motivation and definitions 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 quadrants) that can be represented as an 8-ary tree (or 4-ary tree) for which efficient algorithms 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 fields such as medical applications and sonar imaging. 19.3.2 Construction of octrees To create an octree, we apply a classification 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 classification 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 sufficient 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-infinite 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, classification 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, difference): A low resolution tree could be an effective 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]: Figure 19.4: An octree model. Figure 19.5: Classification of an quadtree node with respect to a concave object 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 representation. • Validity: if no special connectivity is required, all octrees are valid. • Unambiguity and uniqueness: for a fixed resolution there is only one compacted 2 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. 2 Algorithms 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. 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äntylä [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]. 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, finite element analysis methods, scientific visualization of scalar and vector fields. • 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 different cell types, eg. box-like, tetrahedra-like, etc. Figure 19.7: A cell decomposition (finite 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. Figure 19.8: Cell data structure for a 2D model Volume Face Edge Vertex Figure 19.9: Cell data structure for a 3D model 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 approximations; • Validity: they require an intersection test for verification; • Unambiguity: they provide an unambiguous representation; • Nonuniqueness: Similarly to the Constructive Solid Geometry method, the same object can be represented at different resolutions or with different 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 Representation; • Applicability: finite element meshing, multimaterial non-homogeneous objects, visualization of fields, etc. 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 affects drag, volume affects 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 evaluation of single, double and triple integrals of the form φcurve = Z curve f (P)dL, φsurf ace = Z surf ace 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 f dV = XZ i (19.3) f dV celli 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 Z A∪B A−B f dV f dV = = Z ZA A f dV + f dV − Z ZB f dV − A∩B f dV. Z f dV A∩B (19.4) (19.5) Consequently we need to compute integrals over primitives (which can be evaluated analytiR cally) and integrals involving intersections of primitives A∩B f dV which can be approximated using a ray casting, ray classification 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 (P dx + Qdy) = ZZ ∂Q ∂P − ∂x ∂y R ! dA (19.6) 2. Divergence Theorem (also called Gauss’ Theorem) The flux of vector field F flowing 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 ∇FdV (19.7) where n is the outward unit normal vector and ∇·F = ∂F ∂F ∂F ·i+ ·j+ · k. ∂x ∂y ∂z (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 defined by r = (x(t), y(t)) , to ≤ t ≤ t1 (19.9) • Length L = = Z s1 s0 Z t1 t0 ds = q Z t1 t0 q ṙ(t) · ṙ(t)dt ẋ2 (t) + ẏ 2 (t)dt (19.10) • “Centroid” R s1 rds s0 ds rc = (xc , yc ) = Rs0s1 1 L = • “Moments of inertia” Ixx = Iyy = Ixy = 19.6.2 s1 Z s0 Z s1 s0 Z s1 s0 Z t1 q r(t) ẋ2 (t) + ẏ 2 (t)dt t0 Z y 2 ds = t1 t0 Z t1 x2 ds = xyds = (19.11) t0 Z t1 t0 q y 2 (t) ẋ2 (t) + ẏ 2 (t)dt q x2 (t) ẋ2 (t) + ẏ 2 (t)dt q x(t)y(t) ẋ2 (t) + ẏ 2 (t)dt (19.12) (19.13) (19.14) 3D curves Let a 3D curve be defined by r = (x(t), y(t), z(t)) , to ≤ t ≤ t1 • Length L = = s1 Z s0 Z t1 t0 ds = q Z t1 t0 q (19.15) ṙ(t) · ṙ(t)dt ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt • “Centroid” R s1 rds s0 ds rc = (xc , yc , zc ) = Rs0s1 1 L = Z t1 t0 • “Moments of inertia” Ixx = Iyy = Izz = Ixy = Iyz = Ixz = Z s1 s0 Z s1 s0 Z s1 s0 Z s1 s0 Z s1 s0 Z s1 s0 2 2 (y + z )ds = (x2 + z 2 )ds = 2 2 (x + y )ds = xyds = yzds = xzds = Z t1 t0 Z t1 t0 Z t1 t0 Z (19.16) (19.17) q r(t) ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt t1 t0 Z t1 t0 Z t1 t0  y 2 (t) + z 2 (t)  q x2 (t) + z 2 (t)  2 2 q x (t) + y (t) q ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt (19.19) q ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt (19.20) x(t)y(t) ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt q y(t)z(t) ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt q (19.18) x(t)z(t) ẋ2 (t) + ẏ 2 (t) + ż 2 (t)dt (19.21) (19.22) (19.23) 19.7 Integral properties of surface patches 19.7.1 Planar regions Let us consider a planar region as in Figure 19.10 C:r(t)











































































































































 A





























































































































 Figure 19.10: Planar region A • Area A= ZZ dA A (19.24) Using Green’s theorem with Q(x, y) = x and P (x, y) = −y, then 1 1 (Qx − Py ) = (1 + 1) = 1 2 2 (19.25) where subscripts x, y denote partial derivatives. We can rewrite equation (19.24) using Green’s theorem as A = ZZ A 1 (Qx − Py )dxdy 2 1 = (P dx + Qdy) 2 C I 1 (−ydx + xdy) = 2 C I 1 = (xẏ − y ẋ)dt 2 C I (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, numerical integration methods may be used more easily (see Appendix). • Centroid RR rdA , A dA rc = (xc , yc ) = RRA where A is the shaded area, and ZZ Let Q(x, y) = x2 , P (x, y) 2 A ZZ rdA = xdA, A ZZ  (19.27) ydA . A = 0, then Qx − P y = x − 0 = x (19.28) Therefore ZZ A xdA = I P dx + Qdy = C I x2 dy C (19.29) where C is the complete boundary of A. Thus, RR xdA A dA RRA xc = 1 = xdA A A I 1 = x2 ẏdt A C ZZ (19.30) 2 Let Q(x, y) = 0, P (x, y) = − y2 , then Qx − P y = 0 + y = y (19.31) Similarly, yc = RR ydA A dA = 1 A ZZ RRA = − 1 A ydA A I y 2 ẋdt (19.32) C • Moments of inertia 1. Ixx = A y 2 dA 3 Let Q(x, y) = 0, P (x, y) = − y3 , then RR Qx − P y = y 2 (19.33) Using Green’s Theorem, Ixx = ZZ = − A I (Qx − Py )dxdy = C 3 y ẋdt I C −y 3 dx (19.34) 2. Iyy = A x2 dA 3 Let Q(x, y) = x3 , P (x, y) = 0, then RR Qx − P y = x2 (19.35) Using Green’s Theorem, Iyy = = ZZ I A (Qx − Py )dxdy = 1 4 I C 3 x3 dy 3 (19.36) x ẏdt C 3. Ixy = A xydA 2 Let Q(x, y) = x2y , P (x, y) = 0, then RR Qx − Py = xy + 0 = xy (19.37) Using Green’s Theorem, Ixy = = ZZ 1 2 A I (Qx − Py )dxdy x2 y ẏdt C (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) ∈ A, where A is a given parametric domain, as in Figure 19.11. r(u,v) v u Figure 19.11: Curved surface patch • Area A = = ZZ dA ZZA |ru × rv |dudv = A ZZ √ A EG − F 2 dudv (19.39) where E, F and G are the first fundamental form coefficients E = ru ·ru , F = ru ·rv , G = rv · rv . (see Chapter 2). • “Centroid” RR rdA A dA rc = (xc , yc , zc ) = RRA 1 = A • “Moments of inertia” Ixx = Iyy = Izz = Ixy = Ixz = Iyz = ZZ A √ [x(u, v), y(u, v), z(u, v)] EG − F 2 dudv √ [y 2 (u, v) + z 2 (u, v)] EG − F 2 dudv ZZA √ [x2 (u, v) + z 2 (u, v)] EG − F 2 dudv ZZA √ [x2 (u, v) + y 2 (u, v)] EG − F 2 dudv ZZA √ [x(u, v)y(u, v)] EG − F 2 dudv ZZA √ [x(u, v)z(u, v)] EG − F 2 dudv ZZA √ [y(u, v)z(u, v)] EG − F 2 dudv ZZ A (19.40) (19.41) (19.42) (19.43) (19.44) (19.45) (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 dV V (19.47) Choose r = xi + yj + zk (19.48) then ∇·r= ∂r ∂r ∂r ·i+ ·j+ ·k=3 ∂x ∂y ∂z (19.49) Using the Divergence (or Gauss’) Theorem, V = ZZZ V 1 dV = 3 ZZZ V ∇ · rdV 1 ZZ 1 ZZ r · ndA = r · n|ru × rv |dudv 3 A 3 A ZZ 1 = r · (ru × rv )dudv 3 A = given that n= ru × r v |ru × rv | (19.50) (19.51) • Centroid RRR rdV V dV rc = (xc , yc , zc ) = RRRV Choose 1 r = x2 i 2 (19.52) (19.53) then ∂r ∂r ∂r ·i+ ·j+ ·k=x ∂x ∂y ∂z ZZZ ZZZ 1 1 = xdV = ∇ · rdV V V V V ZZ 1 2 1 x (i · n)dA = V A2 1 ZZ 1 2 x (i · (ru × rv ))dudv = V A2 ∇·r = xc (19.54) (19.55) Similarly, expressions are obtained for yc , zc : yc zc 1 1 2 = y (j · (ru × rv ))dudv V A2 ZZ 1 2 1 z (k · (ru × rv ))dudv = V A2 ZZ (19.56) (19.57) • Moments of inertia Ixx = ZZZ V (y 2 + z 2 )dV (19.58) Choose r = (y 2 + z 2 )xi (19.59) then ∇·r = ∂r ∂r ∂r ·i+ ·j+ · k = y2 + z2 ∂x ∂y ∂z (19.60) Thus Ixx = = = ZZZ 2 V ZZ ZZA A 2 (y + z )dV = ZZZ V ∇ · rdV (y 2 + z 2 )x(i · n)dA (y 2 + z 2 )x(i · (ru × rv ))dudv (19.61) (x2 + z 2 )y(j · (ru × rv ))dudv (19.62) (x2 + y 2 )z(k · (ru × rv ))dudv (19.63) Similarly, Iyy = Izz = Ixy = Ixz = Iyz = 19.9 ZZ ZZA ZZA ZZA xyz(k · (ru × rv ))dudv (19.64) xzy(j · (ru × rv ))dudv (19.65) A yzx(i · (ru × rv ))dudv (19.66) ZZA Example: solid of revolution r(t)=x(t)i+y(t)j y y end1 end2 θ z x(t1) x(t2) x Figure 19.12: Solid of revolution Let a solid of revolution be defined by r(t, θ) = x(t)i + y(t) sin θj + y(t) cos θk, t1 ≤ t ≤ t2 , 0 ≤ θ ≤ 2π. (19.67) and we assume that the two end caps are closed off with planar disks (see Figure 19.12). • Surface area of surface of revolution (with end caps) rt = (ẋ, ẏ sin θ, ẏ cos θ) rθ = (0, y cos θ, −y sin θ) (19.68) (19.69) i j k ẋ ẏ sin θ ẏ cos θ 0 y cos θ −y sin θ rt × r θ = = (−y ẏ sin2 θ − y ẏ cos2 θ)i + ẋy sin θj + ẋy cos θk = −y ẏi + ẋy sin θj + ẋy cos θk |rt × rθ | = q A = ZZ q y 2 ẏ 2 + ẋ2 y 2 = y ẋ2 + ẏ 2 dA = A = 2π Z Z t2 t1 t2 t1 Z 2π 0 q (19.70) (19.71) q y ẋ2 + ẏ 2 dθdt y ẋ2 + ẏ 2 dt (19.72) • Volume rt × rθ = (−y ẏ, ẋy sin θ, ẋy cos θ) (19.73) r · (rt × rθ ) = (x, y sin θ, y cos θ) · (−y ẏ, ẋy sin θ, ẋy cos θ) = −xy ẏ + ẋy 2 sin2 θ + ẋy 2 cos2 θ = −xy ẏ + ẋy 2 (19.74) 1 1 t2 2π 1 r · (rt × rθ )dtdθ − = xdxdy + xdxdy 3 t1 0 3 Aend1 3 Aend2 Z Z π π 1 t2 2π (−xy ẏ + ẋy 2 )dθdt − x(t1 )y 2 (t1 ) + x(t2 )y 2 (t2 ) = 3 t1 0 3 3 Z t2 π 2π π (19.75) = (−xy ẏ + ẋy 2 )dt − x(t1 )y 2 (t1 ) + x(t2 )y 2 (t2 ) 3 t1 3 3 Z Z V ZZ ZZ Using integration by parts, − Z t2 t1 xy ẏdt = −[xy 2 ]tt21 + 2 Z t2 t1 (xy)0 ydt 2 = −x(t2 )y (t2 ) + x(t1 )y (t1 ) + Z t2 t1 (ẋy 2 + xy ẏ)dt (19.76) Thus −2 Z t2 t1 xy ẏdt = −x(t2 )y 2 (t2 ) + x(t1 )y 2 (t1 ) + Z t2 t1 ẋy 2 dt (19.77) The volume, therefore, is Z Z 2π t2 2 π 2π t2 π V = ẋy dt − (−xy ẏ)dt − x(t1 )y 2 (t1 ) + x(t2 )y 2 (t2 ) 3 t1 3 t1 3 3 = π Z t2 ẋy 2 dt (19.78) t1 (corroborating the obvious formula from elementary calculus) • Centroid 1. i · (rt × rθ ) = −y ẏ 1 Z t2 Z 2π 1 2 π 2 π 2 xc = x (−y ẏ)dtdθ − x (t1 )y 2 (t1 ) + x (t2 )y 2 (t2 ) V t1 0 2 2V 2V Z π 2 π t2 (x (t2 )y 2 (t2 ) − x2 (t1 )y 2 (t1 )) (19.79) (−x2 y ẏ)dt + = V t1 2V Integrate by parts − Z t2 t1 2 x y ẏdt = −[x2 y 2 ]tt21 2 + Z 2 t2 t1 (x2 y)0 ydt 2 2 = −x (t2 )y (t2 ) + x (t1 )y (t1 ) + Thus −2 and Z t2 t1 Z t2 t1 x2 y ẏdt = −x2 (t2 )y 2 (t2 ) + x2 (t1 )y 2 (t1 ) + 2 (2xẋy + x2 ẏ)ydt (19.80) Z t2 t1 π Z t2 xẋy 2 dt xc = V t1 (corroborating the obvious formula from elementary calculus) 2. j · (rt × rθ ) = ẋy sin θ ZZ 1 2 1 y ẋy sin θdtdθ yc = V A2 Z t2 Z 2π 1 = y 3 ẋ sin θdθdt 2V t1 0 Z 1 t2 3 y ẋ[− cos θ]2π = 0 dt 2V t1 = 0 3. k · (rt × rθ ) = ẋy cos θ 1 1 2 z ẋy cos θdtdθ V A2 1 Z t2 Z 2π = ẋyz 2 cos θdθdt 2V t1 0 Z 1 t2 ẋyz 2 [sin θ]2π = 0 dt 2V t1 = 0 zc = xẋy 2 dt (19.81) (19.82) (19.83) ZZ (19.84) 19.10 Appendix: Review of numerical integration methods 19.10.1 Trapezoidal rule of integration Here we compute ab 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. R • Area in each subinterval [xi , xi+1 ] Z xi+1 f (x)dx ≈ (f (xi ) + f (xi+1 )) xi h ∆x = (fi + fi+1 ) 2 2 (19.85) where ∆x = xi+1 − xi . • [a, b] subdivided into n-subintervals Z b a f (x)dx = n X h h (fi + fi+1 ) = (f1 + 2f2 + 2f3 + · · · + 2fn + fn+1 ) 2 i=1 2 f(x) x1=a x2 x3 x4 ... xn xn+1=b Figure 19.13: Trapezoidal rule of integration • Error in Trapezoidal rule [3] Local error Global error h3 00 f (ξ), xi < ξ < xi+1 12 h2 − (b − a)f 00 (ξ), a < ξ < b 12 − (19.86) 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 f2 = c f3 = ah2 + bh + c (19.88) (19.89) (19.90) If we solve for a, b, and c, we obtain y f2 f3 f1 −h 0 h x Figure 19.14: Simpson’s rule f1 − 2f2 + f3 2h2 f3 − f 1 b = 2h c = f2 a = (19.91) (19.92) (19.93) Thus, Z h −h a 3h b [x ]−h + [x2 ]h−h + c[x]h−h 3 2 2a 3 h + 2ch = 3 h = [f1 + 4f2 + f3 ] 3 (ax2 + bx + c)dx = (19.94) • Area in [a, b] subdivided into n-subintervals Z b a f (x)dx = h (f1 + 4f2 + 2f3 + 4f4 + 2f5 + · · · + 2fn−1 + 4fn + fn+1 ) 3 (19.95) • Error in Simpson’s rule [3] Local error Global error 19.10.3 h5 (4) f (ξ), xi < ξ < xi+1 90 b − a 4 (4) − h f (ξ), a < ξ < b 180 − 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 ) = C1 h2 , where C1 is constant. Then, True value = T0,1 + C1 h2 (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 C 1 . Subtracting (19.97) from four times (19.98) yields 1 True value ≈ T0,2 ≡ T1,1 + (T1,1 − T0,1 ) 3 Similarly, we can obtain for ∆x = (19.99) h 4 C1 2 h 16 (19.100) 1 T1,2 = T2,1 + (T2,1 − T1,1 ) 3 (19.101) True value ≈ T2,1 + From (19.98) and (19.100), we obtain We can make a further improvement by using T0,2 and T1,2 and setting up the relations True value = T0,2 + C2 h4 True value ≈ T1,2 + C2 (19.102) h 2 !4 (19.103) and hence T0,3 = T1,2 + 1 (T1,2 − T0,2 ) 15 (19.104) T0,1 T1,1 T2,1 T0,2 T1,2 T0,3 We can arrange procedures in a matrix form as We can also summarize the rule as Improved value = more accurate +  1 (more accurate − less accurate) 2n − 1  where n is the exponent on h in the error term O(hn ). R1 Examples: 0 cos(x)dx = sin(1) = 0.841470985 T0,1 = T1,1 = T2,1 = T0,2 = T1,2 = T0,3 = 1 (cos(0) + cos(1)) = 0.770151153 2 1 1 T0,1 + cos(0.5) = 0.823866857 2 2 1 1 T1,1 + (cos(0.25) + cos(0.75)) = 0.837083751 2 4 1 T1,1 + (T1,1 − T0,1 ) = 0.841651966 3 1 T2,1 + (T2,1 − T1,1 ) = 0.841489382 3 1 T1,2 + (T1,2 − T0,2 ) = 0.841479221 16 which is accurate to 5 significant digits. 19.10.4 Double integrals ZZ A f (x, y)dA (19.105) 1. Over a rectangular domain y d c b a I= ZZ A f (x, y)dA = Z b a Z d c ! f (x, y)dy dx = x Z d c Z b a ! f (x, y)dx dy (19.106) We will use the trapezoid rule in x and y directions. Let b−a n d−c = n hx = (19.107) hy (19.108) we start with y = c, y = c : I1 = = Z b a f (x, c)dx hx (f1 (c) + 2f2 (c) + · · · + 2fn (c) + fn+1 (c)) 2 (19.109) where fi (c) = f (xi , c). Similarly, y = c + h y : I2 = = Z b f (x, c + hy )dy a hx (f1 (c + hy ) + 2f2 (c + hy ) + · · · + 2fn (c + hy ) + fn+1 (c + hy )) 2 (19.110) ··· y = d : In+1 = = Z b a f (x, d)dy hx (f1 (d) + 2f2 (d) + · · · + 2fn (d) + fn+1 (d)) 2 (19.111) We now sum I1 , I2 , · · · , In+1 in y direction in terms of the trapezoidal rule I= hy (I1 + 2I2 + · · · + 2In + In+1 ) 2 (19.112) 2. Over a curved boundary domain, see Figure 19.15 First we need to find 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 bi − a i n d−c = n hix = (19.113) hy (19.114) y an+1 bn+1             d c a2 a1 b1 b2 C(t)=(x(t),y(t)) x Figure 19.15: A curved boundary domain We start with y = c: y = c : I1 = = y = c + h y : I2 = = Z b1 a1 h1x 2 Z f (x, c)dx (f1 (c) + 2f2 (c) + · · · + 2fn (c) + fn+1 (c)) (19.115) b2 a2 h2x 2 f (x, c + hy )dx (f1 (c + hy ) + 2f2 (c + hy ) + · · · + 2fn (c + hy ) + fn+1 (c + hy )) (19.116) ··· y = d : In+1 = = Z bn+1 an+1 hn+1 x 2 f (x, d)dx (f1 (d) + 2f2 (d) + · · · + 2fn (d) + fn+1 (d)) (19.117) These formulae can be extended to curved domain boundaries bounding multiply connected domains. Bibliography [1] L. Bardis and N. M. Patrikalakis. Topological structures for generalized boundary representations. 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 Å. Björck. Numerical Methods. Prentice-Hall, Inc., Englewood Cliffs, 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. Stolfi. 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 Cliffs, New Jersey, 1976. [8] C. M. Hoffmann. 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äntylä. 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. [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, Scientific Visualization of Physical Phenomena, pp. 211–232. Tokyo: Springer-Verlag, 1991.