c Copyright by Holun Cheng, 2002 ALGORITHMS FOR SMOOTH AND DEFORMABLE SURFACES IN 3D BY HOLUN CHENG B.S., Chinese University of Hong Kong, 1994 M.Phil., Hong Kong University of Science and Technology, 1996 THESIS Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computer Science in the Graduate College of the University of Illinois at Urbana-Champaign, 2002 Urbana, Illinois Abstract In this thesis, we study the skin surface as a new paradigm for the deformable surfaces. The skin surface handles deformation and topology changes robustly, supported by the un- derlying structure of Delaunay triangulations and alpha shapes. The surface serves as a deformable manifold in various disciplines, such as computer graphics, molecular modeling, and mechanical engineering. We develop an algorithm and software for the construction and visualization of the skin surface in 3D in various ways, namely, a parametric representation, static and dynamic triangulations. The triangulation algorithm is guaranteed to terminate with a high quality triangle mesh. In our investigation, geometric properties of the skin serve as the foundation of our proofs and insights for the algorithms. The proofs can be extended to the meshing of other low degree surfaces, such as NURBS. The surfaces created by the software bring stability in finite element methods and visualization of molecular structures to scientists. iii To Ngan-chung Koon and Kam-chung Cheng. iv Acknowledgments I absolutely cannot do this thesis alone. There are many people I would like to thank heartily. First, Prof. Herbert Edelsbrunner gave me all these ideas, guidance, supports. However, I would like to thank especially his patient on me. Without his help, this thesis will not exist. Also, I would like to thanks my defense committee,namely, Herbert, Prof. Jeff Erickson, Prof. Robert Haber, Prof. Michael Garland, Prof. John Hart, and Prof. Mehdi T. Harandi. I deeply understand forming the committee and agree with a time for the defense is very difficult for them. Moreover, staffs in computer science, namely Barb Cicone and Erna Amerman keep helping me, especially when I am not in Illinois. My parents, Kam-chung Cheng and Ngai-chung Koon, together with my sister, Elaine Ho-ling Cheng gave me incredible support through out my life. My parents are the more wonderful parents in the world and I deeply thank them for every bit of me they gave me. Moreover, my brothers and sisters in Christ help me a lot emotionally, listen to my complains,and , pray for me, etc. Wai-yip Chang, Man-yee Lo and Wing-yee Cheung called me many times from Hong Kong regularly (I hope I did not make them bankrupt). In North Coralina, I would like to thank specially Diana and William Sun, Annisa and Ben Lee, Ruth and Brian Howard, Connie Fan, and Sannie Chong and her mother. They are like my relatives here even though I just knew them for less than two years. Their (daily!) prayer supports lift me up from time to time. Also, Jasmine Ku, Cynthia and Billy Saelim care me a lot and always be my 'complain-listeners'. Lynda and Dennis Hui, Sophia and Season Wong, Stella and Dacon Chow also gave me a lot of support. Kathy Lee delivered my machine from Hong Kong, by person! Once, Cissy Ho in Charlotte even drove two hours v to deliver a Nanking beef noodle for me. Frank Fong is my great helper... in computer games. They are all incredibly kind brothers and sisters in Christ. My friend Ka-kit Wong solves my catering problem recently. Yusu Wang Pang, and Andrew Ban are my friends here in Duke Computer Science Department who went to coffee with me from time to time, not to mention that Dr. Xiaobai Sun's great help in the last minute. Also, I would like to thank Damrong Guoy, King-shan Lui for helping me in Champaign when I was not there. Especially thanks Rev. Mike McQueen for his prayer in Illinois. There are too many people to thank. At last, I thank my God in the highest who made and saved me. May all the glory belongs to Him. This research was supported by NSF under grant CCR-0096229 and the Army under grant DAAG 55-98-1-0177. Portion of this thesis have been previously published [21, 22, 23, 24, 25]. vi Table of Contents Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Deformable Surface Modeling: A Brief Review . . . . . . . . . . . . . . . . . 1 1.2 A New Approach: The Skin . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Chapter 2 Geometric Fundamentals . . . . . . . . . . . . . . . . . . . . . . . 9 2.1 Weighted Points and Unions of Balls . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Simplicial Complexes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Weighted Voronoi Diagrams and Delaunay Triangulations . . . . . . . . . . . 11 2.4 Alpha Shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5 Orthospheres and the Size of a Simplex . . . . . . . . . . . . . . . . . . . . . 15 Chapter 3 Skin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1 Sphere Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1.1 Sphere Arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.1.2 Convex and Affine Hulls of Spheres . . . . . . . . . . . . . . . . . . . 20 3.1.3 Orthogonality and Coaxal Systems . . . . . . . . . . . . . . . . . . . 21 3.2 Skin Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 Skin Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.1 Complementarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.2 Envelopes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.3 Mixed Cells and Centers . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.4 Skin Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 Complementarity and Sandwiching Spheres . . . . . . . . . . . . . . . . . . . 28 3.4.1 Complementarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4.2 Sandwiching Spheres . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Chapter 4 Skin Surfaces in B´ezier Form . . . . . . . . . . . . . . . . . . . . 31 4.1 B´ezier Curves and Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2 Sphere Patches in B´ezier Form . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2.1 Stereographic Projection . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2.2 Trimming Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.3 Hyperboloids of Two Sheets . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.1 Parametrization by Rational Projection . . . . . . . . . . . . . . . . 37 4.3.2 Trimming Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 vii 4.4 Hyperboloids of One Sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.4.1 Stack of Two-dimensional Projections . . . . . . . . . . . . . . . . . . 41 4.4.2 Trimming Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Chapter 5 Skin Triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1 Restricted Delaunay Triangulation . . . . . . . . . . . . . . . . . . . . . . . 48 5.2 Curvature and Normals of Skin Surface . . . . . . . . . . . . . . . . . . . . . 49 5.2.1 Maximum Curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2.2 Curvature Continuity . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.2.3 Curvature Variation . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.2.4 Normal Variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.3 Triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.3.1 Distance Claims . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.3.2 Normal Lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.4 Closed Ball Property . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.4.1 Voronoi Edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.4.2 Voronoi Polygon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.4.3 Voronoi Polyhedron . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Chapter 6 Static Skin Triangulation Algorithm . . . . . . . . . . . . . . . . 64 6.1 Curvature Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.2 Generating a Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 6.3 Surface Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.4 Refinement Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.4.1 Edge Flipping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.4.2 Edge Contraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.4.3 Point Insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.5 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.5.1 Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.5.2 Refinement Termination . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.5.3 Close Dual Vertices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Chapter 7 Dynamic Skin Algorithm . . . . . . . . . . . . . . . . . . . . . . . 78 7.1 Maintaining Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 7.2 Vertex Movement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 7.2.1 Parametrization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.2.2 Crossing Mixed Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 7.3 Event Scheduling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 7.4 Topology Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.1 Growth Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.2 Hot Spots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.4.3 Time for Change . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 7.4.4 Special Cap Construction . . . . . . . . . . . . . . . . . . . . . . . . 89 viii 7.4.5 Special Waist Construction . . . . . . . . . . . . . . . . . . . . . . . 90 7.4.6 Special Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 7.4.7 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.5 Scheduling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 7.5.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 7.5.2 Early Warnings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.5.3 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 7.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 Chapter 8 Measurements of Skin Curves . . . . . . . . . . . . . . . . . . . . 105 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 8.2 Union of Disks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 8.2.1 Area and Perimeter of Two Disks . . . . . . . . . . . . . . . . . . . . 108 8.2.2 Area and Perimeter of Union of Disks . . . . . . . . . . . . . . . . . . 109 8.2.3 Notations and Decomposition of Derivatives Computation . . . . . . 110 8.2.4 Slope-preserving Motion . . . . . . . . . . . . . . . . . . . . . . . . . 111 8.2.5 Distance-preserving Motion . . . . . . . . . . . . . . . . . . . . . . . 112 8.2.6 Assembly of Equations . . . . . . . . . . . . . . . . . . . . . . . . . . 114 8.3 2D Skin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 8.3.1 Partition of Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 8.3.2 Centers of Triangles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 8.3.3 Area and Perimeter of 2D Skin . . . . . . . . . . . . . . . . . . . . . 121 8.3.4 Derivatives of Area and Length . . . . . . . . . . . . . . . . . . . . . 127 8.3.5 Simplex Size Changes . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 8.3.6 Boundary Motions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 8.3.7 Derivatives from Boundary Motions of Hyperbola . . . . . . . . . . . 131 8.3.8 Assembly of Equations of Area Derivative of Skin . . . . . . . . . . . 133 8.3.9 Theorems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 8.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 8.5 Table of Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 Chapter 9 Results and Future Applications . . . . . . . . . . . . . . . . . . 142 9.1 Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 9.1.1 Static Skin and Refinement . . . . . . . . . . . . . . . . . . . . . . . 142 9.1.2 Dynamic Triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . 144 9.2 Future Applications of Skin Surface . . . . . . . . . . . . . . . . . . . . . . 147 9.2.1 Molecular Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 9.2.2 Computer Animation . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 Chapter 10 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 Vita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 ix List of Tables 7.1 The four types of topological changes . . . . . . . . . . . . . . . . . . . . . . 86 7.2 Conditions (I)-(V). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 9.1 Performance of static skin algorithm. . . . . . . . . . . . . . . . . . . . . . . 143 9.2 The relationship between Q and the minimum angle . . . . . . . . . . . . . . 144 9.3 Time comparison for static and dynamic skin algorithms. . . . . . . . . . . . 147 x List of Figures 1.1 Modeling molecules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Models of deformable surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Growing a skin surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1 Union of disks and its growth . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Four kinds of simplices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Union of spheres and its Voronoi diagram . . . . . . . . . . . . . . . . . . . . 12 2.4 The Delaunay triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5 Alpha complex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.6 Orthogonal circles in R2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.7 Example of two orthospheres in B : . . . . . . . . . . . . . . . . . . . . . . . 17 3.1 Convex hull of 2 spheres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Coaxal system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 A simple skin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4 Shrinking operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.5 Envelope of a shrunken coaxal system. . . . . . . . . . . . . . . . . . . . . . 25 3.6 Voronoi mixed cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.7 The decomposition of skin surface . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1 B´ezier curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2 The Sphere patch in B´ezier form . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.3 From Poincar´e model to Klein disk . . . . . . . . . . . . . . . . . . . . . . . 37 4.4 Klein disk to hyperboloid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.5 Rotation of a line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.1 Restricted Delaunay triangulation of a surface . . . . . . . . . . . . . . . . . 48 5.2 Restricted Delaunay triangulation of two circles in R2 . . . . . . . . . . . . 49 5.3 Sandwiching spheres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.4 Normal variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.5 Long distance claim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.6 Triangle normal lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.7 Voronoi polygon lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 6.1 Non-flippable edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6.2 Edge flipping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 xi 6.3 Edge contraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.4 Circumcenter lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 7.1 Integral lines of the circle and hyperbola . . . . . . . . . . . . . . . . . . . . 81 7.2 Topology change . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 7.3 Special flipping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 7.4 Top view of a waist . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 7.5 Voronoi polyhedron in special sampling . . . . . . . . . . . . . . . . . . . . . 97 7.6 Chart of edge contraction scheduling . . . . . . . . . . . . . . . . . . . . . . 100 7.7 Chart of point insertion scheduling . . . . . . . . . . . . . . . . . . . . . . . 101 7.8 Graph of the dependence between time and norm . . . . . . . . . . . . . . . 102 8.1 Two grids of a beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 8.2 Distance from a center to the edge . . . . . . . . . . . . . . . . . . . . . . . 107 8.3 Perimeter and area of a clipped disk . . . . . . . . . . . . . . . . . . . . . . 108 8.4 Rotational motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 8.5 Distance-preserving motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 8.6 Multi-disk area derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 8.7 Skin decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 8.8 Motion of boundaries in the 1-mixed cell. . . . . . . . . . . . . . . . . . . . . 119 8.9 Boundaries relationship of the hyperbola. . . . . . . . . . . . . . . . . . . . . 119 8.10 Computing the center of a triangle . . . . . . . . . . . . . . . . . . . . . . . 120 8.11 The first quadrant of a hyperbola. . . . . . . . . . . . . . . . . . . . . . . . . 123 8.12 Four different possible configurations of the first quadrant. . . . . . . . . . . 123 8.13 The length and area of a hyperbola . . . . . . . . . . . . . . . . . . . . . . . 124 8.14 The eight cases of a disk intersecting a triangle . . . . . . . . . . . . . . . . . 126 8.15 General movement of vertices in a triangle . . . . . . . . . . . . . . . . . . . 129 8.16 Naming of points in the triangle mixed cell . . . . . . . . . . . . . . . . . . . 135 9.1 The skin software. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 9.2 The backbone of a molecule. . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 9.3 Two skin surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 9.4 cdk2 molecule. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 9.5 Two views of Gramicidin A. . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 9.6 Two different values of C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 9.7 Transitions of the dynamic skin triangulation. . . . . . . . . . . . . . . . . . 149 9.8 Four skin surfaces with different weighted surface area and volume. . . . . . 150 9.9 Two views of Gramicidin A with its tunnel. . . . . . . . . . . . . . . . . . . 150 9.10 A human head . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 9.11 A sequence of morphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 xii Chapter 1 Introduction Many geometric structures in scientific and engineering computations can be modeled as steady state representations. However, many application areas, including molecular bi- ology, engineering, and computer graphics, require computations on dynamic structures [8, 50, 60, 67, 68, 70, 77, 83, 85, 88, 94]. The coordination of spatial movement is essential, and changes of topology need to be handled robustly. Practitioners use various geometric models to approach these problems. One of the modeling techniques of those structures uses deformable surfaces that represent boundaries of structures, which move according to time or, modification by optimization algorithms or users. We differentiate the study into three levels of concepts. We refer the top level as to the deformable phenomena which are going to be modeled. The middle level is the higher level model, or shape representations of the phenomena. For each representation, it has a lower level representation or implementation. There are many ways to represent a deformable surface in the middle level, and this thesis investigates one specific representation, which is referred as the 'skin'. First, we give some reviews of previous models. 1.1 Deformable Surface Modeling: A Brief Review Deformable surfaces, have many kinds of shape representations [12]. Some current repre- sentations of deformable surfaces are meshes, parametric surfaces [93], and implicit surfaces 1 [17, 92]. Because of the diversity of needs, various application areas use different representa- tions accordingly. This section reviews some major representations used in several disciplines and the challenges involved. Biomolecular Sciences. Biology and chemistry research study inter-molecular behaviors by investigating properties about surfaces of molecules [15, 28, 29]. For example, in order to analyze the ligand docking problem, researchers study the surface electric potentials and geometric properties to interpret how molecules bind to each other. These analyses are facil- itated by representing the surfaces with the van der Waals, molecular, and solvent accessible surface models of molecules [66]. Modeling each atom with a ball of van der Waals radius, the van der Waals surface is the boundary of the union of balls. The molecular and solvent accessible surface models use a probe to define the surfaces. The probe is a sphere with an adjustable radius that rolls over the van der Waals surface. The molecular surface [78] is the boundary of the volume which the probe cannot penetrate. The solvent accessible model [67] is the trace of the probe center when the probe is rolling on the van der Waals surface. probe movement Van der Waals Model Solvent Accessible Model Molecular Model Figure 1.1: Modeling molecules The disadvantage of these models is the lack of smoothness of the surface. In the van der 2 Waals and solvent accessible models, no tangent continuity is maintained at the intersections of the balls. Even when the molecular surface is mostly smooth, it can still have sharp spurs between atoms that are close to each other. Since sharp edges exist at intersections among spheres, triangle meshes of the surface with a good quality cannot be guaranteed [64, 65]. This leads to ill-conditioned numerical computations, for example, in finite element applications. Mechanical and Materials Engineering. Physical simulations with deformable surfaces arise in many engineering research situations [8, 49, 81, 82, 84, 86], for example, the modeling of the boundary between solid and liquid portions of metal during solidification [77]. Such simulations are facilitated by representing the surface with a triangulation or mesh. This process is called meshing. The simulation involves numerical computations which require the triangles in the mesh to be well-shaped. This usually means the triangles need to avoid having large and small angles. In order to achieve this requirement, refinement algorithms are applied to maintain the mesh quality. However, it is a challenge to construct a correct algorithm that terminates, gives desired triangle qualities, and preserves topology. Some of the existing refinement algorithms fail for poorly conditioned surfaces [72]. One of the difficulties is to prove that the algorithm terminates. Heuristics are added in some algorithms in order to stop the infinite loops, but the results can be unsatisfying, with poor triangle quality. Another challenge is preserving the surface topology during refinement. Undesirable holes or extraneous components may be created in the mesh due to the inability of the algorithm to cover all cases. In addition to meshes, engineers adopt the raster model to represent shapes in designing structures, which is referred to as the structural optimization [84]. Elements of grids or voxels are used to represent structures in R2 and R3 respectively. The grid contains values ranging from 0 (empty) to 1 (solid). During the optimization process, finite element application compute the fraction of area, lengthen of the perimeter, and their derivatives within a circular 3 window. With this local information the boundary of the shape is modified. Because of this locality, this method is more efficient than shapes represented by level sets which are zero sets of a function. The disadvantage of the representation is that it has no exact boundary. In situations involving the exact location of the boundary, such as pressure load computation, this problem can be acute if precise interpretation of the boundary geometry is an issue. Other than the raster model, the notion of using level sets has been proposed before. One difficulty of using the model is preserving the possibility of introducing solid material to avoid regions during the optimization process. Santosa investigated level sets in inverse problems involving obstacles [82]. The "bubble method" [49] and other method based on topological derivatives [86] also work with direct surface representations. However, in these methods, a topology change requires a discrete change in the surface representation, and does not result from a variation of continuous design parameters. Computer Graphics. Deforming objects occur often in the field of animation [50, 60, 70, 85], and parametric surfaces are popular in modeling these objects [93]. Animation artists model creatures with surfaces, such as Non-Uniform Rational B-Splines (NURBS) [76]. These parametric surfaces are widely supported by most graphics hardware, as well as software libraries [95]. However, modeling is complicated with NURBS if the animation involves topology changes, for example, poking a hole on an object surface. One solution to reduce the complication is trimming curves. More complicated patches, such as a rectangular surface with a hole, can be represented with these trimming curves. This reduces the number of patches used. Disadvantages of clipping the trimming curves in the parametric domain include expensive computations in solving high degree equations and numerical instability. In addition, controlling the deformation of some simple patches of NURBS is a complicated process [61, 89], and researchers are looking for more intuitive manipulation techniques to enhance user interaction. [52, 58, 69, 74]. Other than parametric surfaces, another representation of shapes is the implicit surface 4 [17]. Implicit surfaces, or level sets, can be defined by the zero-set of a function [12], e.g. metaballs modeled by Gaussian basis functions [16]. Turk and O'Brien showed the use of implicit surfaces for deformation by interpolating different shapes [92]. During morphing, some of the intermediate shapes are not tangent continuous. This usually happens at the time of topology changes when the curvatures and surface normals become undefined. Therefore, it is imperative to know the exact time and location for topology changes in order to extract the portion of the surface with undefined curvatures. Hart investigated such critical moments using Morse theory [54, 55]. For general implicit surfaces, it is still a challenge to acquire such information. Moreover, the manipulation of the implicits is often neither interactive nor intuitive [91, 92] and parametrizing, as well as triangulating implicit surfaces for display and computational purposes, remain challenges [12]. Summary. Deformable surfaces modeled by different shape representations are useful in many fields, such as molecular modeling, engineering and computer graphics. Although some representations have been introduced, each of them has some deficiencies. The sur- faces of molecules lack smoothness. The quality of meshes always depends on the surface properties, such as curvatures. However, not all refinement algorithms give high-quality and faithful meshes. The raster model cannot give a exact boundary. Parametric surfaces favor triangulating and parametrizing, but not morphing and ray tracing. Implicit surfaces get the opposite pros and cons. By far, it seems that the search of a perfect surface is still in progress. 1.2 A New Approach: The Skin This thesis develops a new paradigm, namely, the Skin. The skin surface, defined by Edels- brunner [37, 40], is an implicit surface that favors morphing and ray tracing, and at the same time, is capable of being parametrized and triangulated. It has the following useful 5 properties to tackle the difficulties mentioned in the previous section. Surface Smoothness. The skin surface is a C1-continuous two dimensional manifold. Although non-smooth points are created during the time of topology changes, such moments are still predictable. Special treatment to deal with these degenerate portions is covered. Capability of Deformation. The skin can deform freely with smooth transitions even when topology changes occur [25]. Other surface deformations are often constrained by the topology of the initial and final shapes; for example, two different shapes can morph into each other only if they have the same genus. Decomposable Into Quadratic Patches. The skin is decomposed into quadratic patches and each patch is clipped within a mixed-cell. This quality favors ray-tracing by the partition of space and the low degree of the surface. Mesh Quality. The skin surface can be meshed with good quality triangles while growing. Our algorithm guarantees a minimum angle of 21 and termination [24]. Although the skin is not a perfect surface, it supplements previous studies of shape representations. See Figure 1.2. 1.3 Outline The main goal of this thesis is to investigate and implement the skin surface. During the construction of the software, many difficulties were encountered, and solving them led to the discovery of more properties. This thesis describes three versions of skin software. The first version constructs a NURBS representation of the skin surface to give an interactive visualization [23]. This version relies on the hardware triangulation of the parametric patches. The second one constructs a 6 Deformable Surfaces represented by Unions of Polygonal NURBS Skin Spheres modeled (or implemented) by NURBS Static Mesh Dynamic Mesh Figure 1.2: Models of deformable surfaces triangular mesh using a static approach that generates a mesh and then applies a refinement algorithm on the surface. This version serves as a testing platform of the refinement process for the next version. The last version constructs a triangular mesh following a dynamic approach which constructs a skin mesh through growing the spheres, as shown in Figure 1.3. This version maintains the mesh quality as well as the correctness of the topology during growing. In this thesis, we first present geometric fundamentals such as Delaunay triangulations and alpha shapes in Chapter 2. The definition of the analytical skin surface with sphere algebra is introduced in Chapter 3. The first implementation is presented with parametric surfaces in Chapter 4. The second method with piecewise-linear approximation of the skin is presented afterwards with the proof of its topological correctness in Chapters 5 and 6. The focus of this work is the final implementation, namely, the dynamic skin algorithm in Chapter 7. Chapter 8 discusses the measurements of two-dimensional skins, namely, the volume, area, and their derivatives. These measurements are applied in the topology optimization of structural designs. The description of the skin software is in Chapter 9 and 7 Figure 1.3: Growing a skin surface the conclusion is in Chapter 10. 8 Chapter 2 Geometric Fundamentals In this chapter, geometric background is developed before defining the skin surface. The concepts of simplicial complex [51, 75], Voronoi diagrams [7, 46], Delaunay triangulation and alpha shape [44] are introduced [36]. Generally these definitions apply to any dimensions but are illustrated in R2, and R3 in this work. General position is assumed and the specific assumption is stated before use. 2.1 Weighted Points and Unions of Balls Let bi = (zi, wi) denote a weighted point in Rd × R, where zi = (zi,1, zi,2, zi,3, ..., zi,d) Rd is called the center and wi R is called the weight. Define B as a set of weighted points in Rd: B = {bi | i = 1..n}. We can view B as a collection of balls in Rd, with zi as the centers and ri = wi as the radii. (The notations wi and ri will be used interchangeably, exploiting the relationship wi = r2i.) Points without weights are often referred to as unweighted points. Unweighted points are a special case of weighted points with all wi = 0. We define the growth of these balls by increasing their radii, with a parameter , where 2 R. Growing means increasing the radius of each ball as ri( ) = r2i + 2, with 2 9 continuously increasing in R. So the new set of growing balls is defined as B( ) = {bi( ) | bi( ) = (zi, ri( )2), bi B}. Figure 2.1: Union of disks and its growth The union of balls B( ) = {x | x B( )} is the underlying space of balls in B( ). Figure 2.1 shows an example in R2 and how it grows. If r2i + 2 < 0, the radius ri( ) is imaginary and bi( ) is empty. bi( ) will not contribute to B( ) until 2 r2i. 2.2 Simplicial Complexes Define the affine hull and convex hull of a point set S Rd as aff(S) and conv(S) respectively: aff(S) = {x Rd | x = izi, zi S, i = 1}, conv(S) = {x Rd | x = izi, zi S, i = 1, i 0}. A set T Rd is affinely independent if every point x T is not in aff(T - {x}). A necessary condition for T to be affinely independent is that card(T ) d + 1. A simplex = conv(T ) is 10 the convex hull of an affinely independent point set T . is called a k-simplex if its dimension is k = dim( ) = card(T ) - 1. In Figure 2.2, we show four kinds of simplices in R3 as an example. The four kinds of simplices, from left to right, are a 0-simplex or vertex, an 1- simplex or edge, a 2-simplex or triangle, and a 3-simplex or tetrahedron. For any subset S T , = conv(S) is called a face of and is a coface of . To be consistent, we define the empty set to be the (-1)-simplex and a face of every simplex. A finite set K of simplices is a simplicial complex if: i. if K and is a face of then K, and ii. if , K then = is a face of both and . Figure 2.2: Four kinds of simplices. 2.3 Weighted Voronoi Diagrams and Delaunay Triangulations Recall that B is a finite set of weighted points. Define the weighted squared distance of a point x Rd from bi B as: i(x) = x - zi 2 - wi. This is also called the power of x from bi. Let x, y be the dot product of two vectors x and y. For any two weighted points bi, bj B, any point x at equal weighted distance from bi 11 and bj satisfies the following equations: i(x) = j(x) x - zi 2 - r2i = x - zj 2 - r2j x 2 - 2 x, zi + zi 2 - r2i = x 2 - 2 x, zj + zj 2 - r2j 2 x, zi - zj - zi 2 + zj 2 + r2i - r2j = 0. The last equation shows that the set is a hyperplane, and the set of points x with i(x) j(x) is a linear halfspace of Rd. We define the Voronoi cell i of bi B as the intersection of these halfspaces: i = {x Rd | i(x) j(x), bj B}. Figure 2.3: Union of spheres and its Voronoi diagram Each Voronoi cell i is a convex polyhedron in Rd formed by the intersection of finitely many halfspaces. The set { i | bi B} partitions Rd. Also note that there are situations where i = . The corresponding weighted point bi is called redundant. 12 The Voronoi cell of a set X B is defined as X = i. bi X The geometric interpretation of X is that it is the set of points in Rd for which the weighted square distance is simultaneously minimized by the weighted points bi X. The cell X can be and usually is empty. The Voronoi complex of B is the set of all non-empty Voronoi cells VB = { X | X = , X B}. The Voronoi complexes are called the Voronoi diagrams of weighted points, and also called power diagrams. Figure 2.3 shows an example of the decomposition in R2. For the following discussion, the first assumption for general position is there are no d+1 points in B such that they are all orthogonal to a common sphere. Two spheres bi and bj are orthogonal to each other zi - zj 2 - r2i - r2j = 0. Define X as the convex hull of the centers of the spheres in X: X = conv({zi | bi X}). The Delaunay triangulation of B is DB = { X | X VB}, which is a simplicial complex in Rd if the points are in general position; see Figure 2.4. Before continuing, we describe the scenario when the radii of the spheres in B grow with . If 2 is increased, the Delaunay triangulation does not change. Recall that the hyperplane of common distance from two weighted points bi( ) and bj( ) is 2 x, zi - zj - zi 2 + zj 2 + ri( )2 - rj( )2 = 0. Note that ri( )2 - rj( 2) = r2i + 2 - r2j - 2 = r2i - r2j. The parameter thus cancels in the equation, which implies that the hyperplane remains fixed when changes. Hence, 13 Figure 2.4: The Delaunay triangulation the Voronoi diagram and Delaunay triangulation are the same for all values of , i.e. DB = DB( ), for all 2 R. 2.4 Alpha Shape The alpha complex of B( ) is K( ) = { X DB | B( ) X = }. K( ) is a subcomplex of DB and the underlying space | K( ) | is called the alpha shape. B( ) grows continuously as 2 increases from - to (figure 2.5). For every simplex X DB, there is a corresponding value such that the union of spheres `just touches' X. We call this value of the birth time X of X: B( X) X = , and B( ) X = for all 2 < 2X A filter is a sequence of simplices 1, 2, ..., m such that Ki = { j | 1 j i} is a simplicial complex for all 0 i m. The sequence of prefixes Ki is called a filtration. With 14 Figure 2.5: Alpha complex: When the spheres increase their radii, the alpha complex grows in discrete steps. the birth time, we define the filter of the alpha complex as follows. For any i = X and j = Y , i < j if: 1. 2X < 2Y , or 2. 2X = 2Y and dim( i) < dim( j). Let m be the number of simplices in DB, The alpha filtration is = K0 K1 ... Km = DB, where Ki - Ki-1 = { i}, for all 1 i m. 2.5 Orthospheres and the Size of a Simplex Define i,j = zi - zj 2-r2i-r2j. Two spheres bi and bj are orthogonal to each other, denoted as bi bj, if i,j = 0. If both r2i and r2j are positive, the two orthogonal spheres intersect in a circle and form right angles at the intersections. For each simplex X DB, one or more spheres are orthogonal to all elements in X. For example, if X = {bi}, there is an infinite 15 Figure 2.6: Orthogonal circles in R2: The dashed circles are orthogonal to the shaded circle. number of orthogonal spheres as shown in figure 2.6. To define the set of orthospheres B of B, pick the orthosphere with the smallest squared radius for each X. The center of the smallest orthosphere of X is z X = aff( X) aff( X), and the radius is r X = z X - zi 2 - r2i, for all bi X. Formally, B = {b 2 X = (z X , r X ) | X DB}. Figure 2.7 shows two orthospheres of two simplices. 16 Figure 2.7: Example of two orthospheres in B : One orthosphere of an edge with card(X) = 2 (dashed circle) and one for a triangle with card(X) = 3 (the dotted circle) Size of a Simplex. As the set B can be grown to B( ), define the set B ( ) as the growing of B according to : B ( ) = {b X( ) = (z X, w X( )) | X DB}, where w 2 2 X ( ) = r X ( ) = z X - zi 2 - r2i - 2 = r X - 2. Define the size of a simplex 2 2 X at as sX ( ) = -r X ( ) = 2 - r X . Note that just sX = sX(0) will be used instead of sX( ) if is not involved in the discussion. 17 Chapter 3 Skin The skin surface, defined by Edelsbrunner [37, 40], is introduced in this chapter. This is an implicit surface [17] with similar flexibility of modeling as metaballs [16]. The skin surface has a better control on the surface and is decomposable into quadratic algebraic surfaces while metaballs are defined by exponential functions. First the sphere algebra is defined with arithmetic operations including addition and multiplication of spheres. These operations facilitate the definition of convex and affine hulls of spheres. Then, we define the skin. It is hard to illustrate the skin surface without mentioning the complementary part of it. For this reason, the set of orthogonal spheres and complementarity is discussed. The goal is to give a clear picture of how skins in R2 and R3 look, although the skin can be generalized to any dimension. 3.1 Sphere Algebra The arithmetic operations on spheres are addition, scalar multiplication, and shrinking. With addition and scalar multiplication, the convex and affine hulls of a set of spheres can be defined. Then, with the shrinking operation, the skin is defined. 18 3.1.1 Sphere Arithmetic Given two spheres, bi and bj, and a real number c, addition , scalar multiplication, and power of spheres are defined as follows: (zi, wi) + (zj, wj) = (zi + zj, wi + wj + 2 zi, zj ) c · (zi, wi) = (c · zi, c · (wi - (1 - c) zi 2)) (zi, wi)c = (zi, c · wi). The first two equations are the standard operations on vectors in Rd+1 under the paraboloid lifting map (zi, wi) (zi, zi 2 + wi). Summation of Spheres Lemma. If b n n n 0 = i=1 bi, then z0 = i=1 zi and w0 = i=1 wi + 2 i =l zi, zl . Proof. The case holds for n = 1. For n = 2, it holds also directly from the addition operation defined previously. Assume it is true for n = k, and n n n i=1 bi = i=1 zi, i=1 wi + 2 k+1 i =l zi, zl . For n = k + 1, the sphere b0 = i=1 bi is equal to the sum of first n spheres plus b n+1 k+1. By the addition operation, the center z0 is i=1 zi. and the weight of b0 is k k w0 = wk+1 + wi + 2 zi, zl + 2 zk+1, zi i=1 i =l i=1 k+1 k+1 = wi + 2 zi, zl . i=1 i=1,l =i 19 3.1.2 Convex and Affine Hulls of Spheres With the above operations, define the convex and affine hulls of a set of spheres, X = {b1, b2, ..., bn} as follows: aff(X) = { ibi | i = 1} bi X i conv(X) = { ibi | i = 1, i 0}. bi X i A sphere bi in aff(X) or conv(X) is called an affine combination or convex combination of X respectively. For any sphere bj aff(X) (or in conv(X)), its center and weight are zj = izi i wj = ( iwi + ( 2i - i) zi 2) + 2 izi, lzl i i =l 2 = iwi + izi - i zi 2. i i i Zero-set Lemma. For any set of spheres X = {b n 1, b2, ..., bn} and any i with i i = 1, the sphere b0 = i ibi aff(X) is the zero-set of i i(x). Proof. Recall that for a weighted point bi = (zi, wi), the squared weighted distance function is i(x) = x - zi 2 - wi. The affine combination of the functions is i( x - zi 2 - wi) = i( x 2 - 2 x, zi + zi 2 - wi) i i = x 2 - 2 x, izi + i( zi 2 - wi) i i 2 2 = x - izi - izi + i( zi 2 - wi) i i i The zero set of i i i(x) is exactly the sphere ibi. 20 Figure 3.1: Convex hull of 2 spheres: Two circles (left) and their convex hull (right). An example in R2 is shown in Figure 3.1. We simplify the illustration by using just two circles, bi and bj. If bi intersects bj, the point(s) of intersection satisfies both i(x) = 0 and j(x) = 0, which implies that t1 i(x) + t2 j(x) = 0 for any t1 + t2 = 1. This means all spheres in aff({bi, bj}) pass through the intersection points. 3.1.3 Orthogonality and Coaxal Systems For any two spheres, bi and bj, recall that bi bj if i,j = 0. If a sphere b0 is orthogonal to all bi in X, then it is orthogonal to every sphere bk aff(X). To show this, we check 0,k. 0,k = z0 - zk 2 - w0 - wk 2 2 = z 0 - izi - w0 - iwi - izi + i zi 2 i i i i 2 2 = z 0 2 + izi - 2 z0, izi - w0 - iwi - izi + i zi 2 i i i i i = i( z0 2 + zi 2 - 2 z0, zi - w0 - wi) i = i 0,i i = 0 21 0,k = 0 implies that b0 is orthogonal to every bk in aff(X). Moreover, if two sets of spheres X1 and X2 are such that every bi X1 is orthogonal to every bj X2, then every bl aff(X1) is orthogonal to every bm aff(X2). Such a configuration is referred to as a coaxal system. Figure 3.2 shows an example in R2 with card(X1) = card(X2) = 2. Figure 3.2: Coaxal system: Every solid circle is orthogonal to every dashed circle. 3.2 Skin Definition 1 Define the square root operation on a sphere b 2 i as bi = bi = (zi, wi/2) and B = { bi | bi B} for a set of spheres. The operation shrinks every bi in B by dividing its radius by 2. We define the body of B as body(B) = conv(B), where the operation of union of spheres is defined as B = {x}. b B x b 22 Define (P ) is the boundary of a set of point P . The skin surface as the envelope of the skin body, which is: skin(B) = ( conv(B)). Figure 3.3 gives an example of a simple case with card(B) = 2. Figure 3.3: A simple skin: The convex hull of two spheres (left) and the shrunken conv(B) (right). 3.3 Skin Decomposition The skin surface is the boundary of an infinite set of spheres. In this section, we give a finite description of the skin's shape. First the complementary property of the shrunken affine hull of spheres is studied with less than or equal to four spheres. Then the mixed cell concept is introduced. Finally, the shape of the skin surface can be described as a whole. 23 3.3.1 Complementarity Given two orthogonal balls bi and bj, the distance between their centers is zi - zj 2 = r2i + r2j r2i + r2j - ( 2ri/2 - 2rj/2)2 (r = i + rj)2 . 2 This implies the sum of the shrunken radii is at most the distance between their centers. Furthermore, ri = rj implies the distance between centers equals the sum of shrunken radii. This means that after shrinking, the two orthogonal spheres will touch each other only if they are of the same radius; otherwise, they are disjoint. See Figure 3.4. Figure 3.4: Shrinking operation: After shrinking two orthogonal spheres, they will touch each other only if they have the same radius If X1 and X2 are orthogonal sets of spheres, then aff(X1) and aff(X2) have the same envelopes. See Figure 3.5. 3.3.2 Envelopes We first study skins of up to four weighted points in R3. For card(X) = 1, 2 the shape of the envelope of conv(X) can be illustrated directly. The rest of the cases need the help of the complementary part. Note that X is a subset of B with card(B) > 3. 24 Figure 3.5: Envelope of a shrunken coaxal system. card(X) = 1, 2. The case of 1 sphere is simple and the case of 2 spheres is discussed here. If X just contains 1 sphere, conv({bi}) = bi and the shrunken envelope is a sphere. The envelope of aff(B) with card(B) = 2 is a hyperbola. We simplify the illustration by putting B = {bi, bj} with zi = (ai, 0) and zj = (aj, 0) in R2. So, with t1 + t2 = 1, the affine hull is aff(B) = {((t1ai + (1 - t1)aj, 0), (t1ai +(1-t1)aj)2 -t1a2i -t2a2j +t1wi +t2wj) | t1 R}. To find the envelope of aff({bi, bj}), we put it in one higher dimension along t1. Shrinking 25 each of the radius by 1/ 2 yields f(t1, x1, x2) = (x1 - (t1ai + (1 - t1)aj))2 + x22 - ((t1ai + (1 - t1)aj)2 - t1a2i - t2a2j + t1wi + t2wj)/2 = (x1 - t1(ai - aj) - aj))2 + x22 - ((t1(ai - aj) + aj)2 - t1a2i - t2a2j + t1wi + t2wj)/2 f = -2(x t 1 - t1(ai - aj) - aj))(ai - aj) 1 -(2(t1(ai - aj) + aj)(ai - aj) - a2i + a2j + wi - wj)/2 where f(t1, x1, x2) = 0 is a circle at t1. The envelope of aff(B) is the projection of the silhouette of f-1(0) as viewed along the t1 direction. It is the set of points for f = 0. We t1 simplify the computation by putting a constant c = (-a2i + a2j + wi - wj)/2(ai - aj) and d = (t1a2i + (1 - t1)a2j - t1wi - (1 - t1)wj)/2: f 0 = t1 = -2(x1 - t1(ai - aj) - aj)) - (t1(ai - aj) + aj) - c t1 = (2x1 + c - aj)/(ai - aj) f(t1, x1, x2) = f((2x1 + c - aj)/(ai - aj), x1, x2) = (x1 - 2x1 - c)2 + x22 - (2x1 + c)2/2 + d = -x21 + x22 + c2/2 + d. The terms of x1 is at most degree 1 in the term d because the terms of t1 in d and the terms of x1 in t1 are both degree 1. Expanding the terms, the degree of the x21 is -1 and x22 is 1. This shows the envelope, f(t1, x1, x2) = 0, is a hyperbola. For higher dimensions, more terms of x23, x24, ..., x2d are added but their coefficients will be still equal to 1 in f(t1, x1, x2, ..., xd). Hence, the envelope of two spheres is a hyperboloid with x1 as the rotational axis in any 26 dimension. card(X) = 3, 4. When card(X) = 4, X is a tetrahedron. There will be a sphere, bi, which is orthogonal to all four spheres in X. Recall that the envelopes of the affine hulls of two sets of spheres are the same. So the envelope of 4 spheres is the shrunken sphere bi. When card(X) = 3, X is a triangle. There exists a smallest orthosphere b X which is orthogonal to all spheres in X. Also, there should be at least a tetrahedron Y which is a coface of X in the Delaunay triangulation. b Y will also be orthogonal to all spheres in X. So, skin(X) = skin({b Y , b X}) and which is a hyperboloid. 3.3.3 Mixed Cells and Centers For each Delaunay simplex X in DB, define the mixed cell as the Minkowski sum of X with its dual Voronoi cell X, scaled by 1/2: µX = ( X + X)/2. In R3, the dimension of µX is always 3. All the mixed cells partition R3. Figure 3.6 shows examples of the four different types of mixed cells corresponding to different cardinalities of X. Furthermore, define the center of each mixed cell as z X. Figure 3.6: Voronoi mixed cells: From left to right: a typical Voronoi polyhedron, a Voronoi polygon times a Delaunay edge, a Voronoi edge times a Delaunay triangle, a Delaunay tetrahedron. 27 3.3.4 Skin Surface We can now describe the overall picture of the skin surface. The mixed cells partition R3. For each mixed cell µX, the surface of skin(X) is clipped within µX. Each patch is either a clipped sphere or hyperboloid in R3. All these patches are composed to form the skin surface. See Figure 3.7. 3.4 Complementarity and Sandwiching Spheres In closing, this section states a rather special and important property of skin surfaces : complementarity. 3.4.1 Complementarity Recall that for every B there is an orthogonal sphere set B = {b X = (z X, (r X)2) | X DB}. The mixed cells decomposed the skin into finitely many patches. Within a cell µX, only spheres in conv(X) defined by the Delaunay simplex X and the Voronoi cell X are relevant. In other words, within µX the skin of B looks the same as the envelope of aff(X), and symmetrically, the skin of B looks the same as the envelope of aff({b Y | Y X}). This is true for every mixed cell and together they cover the entire R3. Thus, let Z be the centers of B, the two skins are complementary to each other within conv(Z), i.e. (body(B) body(B )) conv(Z) = conv(Z) (body(B) body(B )) conv(Z) = skin(B) conv(Z) = skin(B ) conv(Z). 28 3.4.2 Sandwiching Spheres Within each mixed cell, there are two sets of spheres of conv(B) and conv(B ). There is a unique pair of spheres from each of the sets which touches the skin surface due to the complementarity in section 3.3.1. Sandwiching Property. For every point x on the skin of B, there are two unique spheres Sx conv(B) and Tx conv(B ) that pass through x. These two spheres Sx and Tx are externally tangent, and have equal radius. The skin surface stays outside both spheres, and is thus tangent to them at x. Sx and Tx are referred to as the sandwiching spheres at x because they squeeze the surface flat in a neighborhood of x. The two spheres are equally large from the complementarity in Section 3.3.1, which is important for further proofs on curvature properties. 29 (a) card(X) = 1 (b) card(X) = 2 (c) card(X) = 3 (d) card(X) = 4 (e) The Skin Surface Figure 3.7: The decomposition of the skin surface into four types of patches and the skin surface. 30 Chapter 4 Skin Surfaces in B´ezier Form The skin surface can be decomposed into patches of spheres and hyperboloids. For display purposes, all patches can be represented in bi-quadratic B´ezier form. This overcomes one of the difficulties of implicit surfaces, namely, parametrizing and triangulating the implicits [12]. The motivation for this representation is that B´ezier surfaces can be easily translated into NURBS patches [79, 93], which are supported by graphical systems like OpenGL [95]. Moreover, it is not necessary to fit the surface with general high degree B-Spline patches, as in [35, 71] since we know the patches only consist of spherical and hyperbolic patches. In this chapter, the B´ezier representation of spheres, one sheet hyperboloids, and two sheet hyperboloids are derived. The trimming curves that clip the patches are also computed within their mixed cells. First, definitions and basic properties of B´ezier curves and surfaces are introduced. Then the B´ezier representation of spheres and hyperboloids, as well as the trimming curves are derived. Without loss of generality, we can treat each of the spheres as if their centers are the origin. The common tip of the asymptotic double cones of any hyperboloid is translated to the origin. The rotational axis is also aligned with the x3-axis in R3. 31 4.1 B´ezier Curves and Surfaces Rational B´ezier curves and surfaces are defined by control points and weights. Control points in R3 govern the general shape of the curve or surface. Also any B´ezier curve or surface is always inside the convex hull of its control points. The weights control the contribution of each control point to the curve (or surface). The curve will be closer to a control point if its weight is large, as shown in Figure 4.1. (1) (2) (1) (2) (1) (1) (1) (1) Figure 4.1: B´ezier curve: An example of a B´ezier curve with 4 control points (left). The curve is attracted to the control points when the weights increase (right). The B´ezier basis of order n , also called Bernstein basis, is the set of polynomials n Bn,i(t) = ti(1 - t)n-i i with t [0, 1] and i = 0, 1, ..., n. These form a basis for the vector space of degree-n polynomials. Rational B´ezier Curve. A degree-n rational B´ezier curve in R3 is defined in terms of n + 1 control points, ci R3, and weights, wi R, with i = 0..n. The curve is denoted Q = {Q(t) | t [0, 1]}. Q(t) is defined as follows: n Q(t) = i=0 ciwiBn,i(t) n . i=0 wiBn,i(t) Rational B´ezier Surface. A degree n × n rational B´ezier surface is defined by a matrix of control points, ci,j, and weights, wi,j, with 0 i, j n. The surface is denoted by 32 Q = {Q(u, v) | u, v [0, 1]}. Q(u, v) is defined as follows: n n Q(u, v) = i=0 j=0 ci,jwi,jBn,i(u)Bm,j(v) n n . i=0 j=0 wi,jBn,i(u)Bn,j(v) Power Basis. There are other representation of curve. One of them is known as the power basis. Any curve Q(t) can be represented in the form: n piti, i=0 for some control points pi. For example, the parabola y = x2 in R2 can be represented as Q(t) = 2i=0 piti with p0 = (0, 0), p1 = (1, 0) and p2 = (0, 1). Given a curve in power basis as Q(t) = 3i=0 piti, we can compute its B´ezier coefficients with a simple matrix transformation [76]: c0 1 0 0 0 p0 c 1 1 0 0 p 1 = 3 1 c 1 2 1 0 p 2 3 3 2 c3 1 1 1 1 p3 . 4.2 Sphere Patches in B´ezier Form Spheres and hyperboloids are degree-2 surfaces and can be represented as bi-quadratic B´ezier patches. In this section, the B´ezier form of the sphere is first derived by using stereographic projection. The projection helps in parametrizing the sphere surface. Then the clipping of B´ezier patch by planes through the trimming curve is derived. A unit sphere is used with equation x21 + x22 + x23 = 1. 33 4.2.1 Stereographic Projection First, the line through the north pole (0, 0, 1) of the unit sphere and a point (u, v, -1) on the plane x3 = -1 is defined by the equation x1 0 u x = t 0 + (1 - t) v . 2 x3 1 0 For each pair u, v, the line intersects the sphere at a unique point other than the north pole. The point of intersection is: 4u x1 = , u2 + v2 + 4 4v x2 = , u2 + v2 + 4 u2 + v2 - 4 x3 = . u2 + v2 + 4 This stereographic projection gives us a parameterization of the sphere in u and v. Consider the portion of the sphere with u, v [-2, 2]. This part of the surface is a sphere patch bounded by four arcs , as shown in Figure 4.2. This patch covers the lower half of the sphere. The other half of the sphere can be completed by reflecting the patch through the x1x2-plane. Normalize the parameters u, v and translate the domain from u, v [-2, 2] to u , v [0, 1]. By putting u = u+2 and v = v+2, the new parametrization of the point (x 4 4 1, x2, x3) on the sphere is: 4u - 2 x1 = , 4u 2 - 4u + 4v 2 - 4v + 3 4v - 2 x2 = , 4u 2 - 4u + 4v 2 - 4v + 3 4u 2 - 4u + 4v 2 - 4v + 1 x3 = . 4u 2 - 4u + 4v 2 - 4v + 3 34 Figure 4.2: The sphere patch in B´ezier form By equating this with the B´ezier form, (x1, x2, x3) = Q(u , v ), we solve for the nine control points and weights as follows: (-2, 2, 1, 1), (0, 2, -1, 1), (2, 2, 1, 1), 3 3 3 3 3 3 3 (-2, 0, -1, 1), (0, 0, 3, -1), (2, 0, -1, 1), 3 3 3 (-2, -2, 1, 1), (0, -2, -1, 1), (2, -2, 1, 1). 3 3 3 3 3 3 3 Here, the fourth coordinates are the weights. 4.2.2 Trimming Curves Each B´ezier patch is clipped within a mixed cell, which is the intersection of finitely many linear half-spaces. We write each half-space in the form x1a1 + x2a2 + x3a3 d, 35 with a1, a2, a3, d R. Substituting Q(u, v) into the above equation, and assuming a3 > d, we have 4ua1 + 4va2 + (u2 + v2 - 4)a3 d(u2 + v2 + 4) 4ua1 + 4va2 + (u2 + v2)(a3 - d) 4(d + a3) 2a 2 2a 2 a2 u + 1 + u + 2 4 1 + a22 + a23 - d2 . a3 - d a3 - d (a3 - d)2 Normalizing to [0, 1], the new inequality in u and v is: a 2 a 2 a2 v + 2 - a3 + d + u + 1 - a3 + d 1 + a22 + a23 - d2 . 2(a3 - d) 2(a3 - d) 4(a3 - d)2 This shows that the trimming curve in the (u , v )-domain is a circle. If a3 < d, the previous inequality changes sign. If a3 = d, ua1 + va2 2a3 which boundary is the following line in the (u , v ) domain: 2u a1 + 2v a2 = a1 + a2 + a3. 4.3 Hyperboloids of Two Sheets Similar to the sphere case, we first describe the parametrization of a hyperboloid in the (u , v )-domain. Then we compute the control points and weights of the B´ezier surface. Finally, we give the trimming curve formulas. Hyperboloids, unlike spheres, are not bounded. To avoid difficulties with unbounded patches, two halfspaces x3 h1 and x3 -h2 are used to clip the hyperboloid within the mixed cell. 36 4.3.1 Parametrization by Rational Projection The hyperboloid with equation -(x21 + x22) + x23 = 1, is rotationally symmetric about the x3-axis and its closest approach to the origin is 1. We develop a rational polynomial parametrization by projecting one model of the hyperbolic plane to another, namely, from the Poincar´e model to the Klein model of hyperbolic space, and then to the hyperboloid model [19]. Figure 4.3: From Poincar´e model (left) to Klein model (right) Poincar´e to Klein Disk. The Poincar´e model of the hyperbolic plane is the unit disk with geodesics that are circular arcs meeting the bounding circle in two right angles. The Klein model is also the unit disk but its geodesics are straight line segments that start and end on the bounding circle, see Figure 4.3. The two models agree at infinity, which is represented by the boundary circle. The map H : B2 B2 that changes the Poincar´e model into the Klein model is 2(u, v) (a, b) = . u2 + v2 + 1 37 This transformation maps a point with norm d = u2 + v2 to one with Euclidean norm u2 + v2 = 2d . The image of the disk is itself, H(B2) = B2. d2+1 x3 x1 1 Figure 4.4: Klein disk to hyperboloid From Klein Disk to Hyperboloid. The Klein model relates to the hyperboloid model of the hyperbolic plane via projection from a point. The hyperboloid model consists of the southern sheet of the canonical hyperboloid. The Klein model is the unit disk in the plane x3 = -1, and the projection is from the origin, see Figure 4.4. The map G : B2 R3 from the Klein model to the hyperboloid model is defined by (x1, x2, x3) = G(a, b) where: a x1 = 1 - a2 - b2 b x2 = 1 - a2 - b2 1 x3 = - 1 - a2 - b2 Rational Parametric Forms. The hyperboloid surface is the composition of these two maps: 2u x1 = , u2 + v2 - 1 2v x2 = , u2 + v2 - 1 u2 + v2 + 1 x3 = . u2 + v2 - 1 38 The domain of u, v is [- , ]. If u2 + v2 - 1 < 0, Q(u, v) maps to the lower sheet of the hyperboloid and its complementary region can be obtained by mirroring the lower sheet. Height Bound. There are two kinds of halfspaces in the mixed cell that clip the hy- perboloid. The clipping planes of the first type are parallel to the x1x2-plane and cut the hyperboloid horizontally. Two of them serve as upper and lower bounds: x3 h1 and x3 -h2. First pick the lower sheet of the hyperboloid, which satisfies u2 + v2 - 1 < 0. Substituting the value of x3 into the lower bounding plane gives u2 + v2 + 1 -h u2 + v2 - 1 2 h (u2 + v2) - 2 - 1 0. h2 + 1 which is a circle in the (u, v) domain. Let l = h2-1 , and now u, v [-l, l]. The next step h2+1 is to normalize (u, v) to (u , v ) [0, 1]2 as follows: u + l u = , 2l v + l v = . 2l The parametric representation is now: 2(2lu - l) x1 = , 4l2u 2 - 4l2u + 4l2v 2 - 4l2v + 2l2 - 1 2(2lv - l) x2 = , 4l2u 2 - 4l2u + 4l2v 2 - 4l2v + 2l2 - 1 4l2u 2 - 4l2u + 4l2v 2 - 4l2v + 2l2 + 1 x3 = . 4l2u 2 - 4l2u + 4l2v 2 - 4l2v + 2l2 - 1 39 By equating the these formulas with Q(u, v), we compute the nine control points: 2l 2l 2l2 + 1 c0,0 = - , - , , 2l2 - 1 2l2 - 1 2l2 - 1 c0,1 = (2l, 0, -1), 2l 2l 2l2 + 1 c0,2 = - , , , 2l2 - 1 2l2 - 1 2l2 - 1 c1,0 = (0, 2l, -1) 2l2 - 1 c1,1 = 0, 0, , 2l2 + 1 c1,2 = (0, -2l, -1), 2l 2l 2l2 + 1 c2,0 = , - , , 2l2 - 1 2l2 - 1 2l2 - 1 c2,1 = (-2l, 0, -1), 2l 2l 2l2 + 1 c2,2 = , , . 2l2 - 1 2l2 - 1 2l2 - 1 The weights are w1,0 = w0,1 = w2,1 = w1,2 = -1, w0,0 = w2,0 = w0,2 = w2,2 = 2l2 - 1, and w1,1 = -2l2 - 1. 4.3.2 Trimming Curves Additional trimming curves are formed by intersecting the halfspaces of the form x1a1 + x2a2 d. Keeping the condition u2 + v2 - 1 < 0 gives d(u2 + v2 - 1) - 2a1u - 2a2v 0 a a a2 (u - 1 )2 + (v - 2 )2 1 + 1 + a22 . d d d2 40 Once again, algebraically it shows these trimming curves are circles in both (u, v) and (u , v ) domains. 4.4 Hyperboloids of One Sheet First align the axis of rotation of the hyperboloid with x3-axis. Then clip the hyperboloid in half with the halfspace x1 0. The reason is the parametrization here would extended to infinity if the the entire hyperboloid is used. 4.4.1 Stack of Two-dimensional Projections 3 x 1 x 2 x Figure 4.5: The line (1, u, u) For a one-sheet hyperboloid of the form x21 + x22 - x23 = 1, 41 there is a line l(u) on the hyperboloid such that l(u) = (1, u, u). As shown in Figure 4.5, the points on this line are rotated to complete the hyperboloid by multiplying the rotational matrix v2-1 - cos - sin 0 - 2v 0 v2+1 v2+1 R(v) = sin - cos 0 = 2v v2-1 0 , v2+1 v2+1 0 0 1 0 0 1 where v = tan . So the hyperboloid is the set of points {l(u)T R(v) | u, v R}. The 2 parametrization of each point is: v2 - 2uv - 1 x1 = , v2 + 1 2v + u(v2 - 1) x2 = , v2 + 1 x3 = u. However, there is a missing line of the hyperboloid when v = ± . By cutting the hyperboloid at x1 = 0, the portion of the surface that surrounds the missing line is eliminated. The hyperboloid is completed by reflecting the clipped portion with x1 = 0. The portion of the hyperboloid clipped by x1 0 has the following implicit equation in the (u, v)-domain: v2 - 2uv - 1 0. Before this curve in B´ezier form is represented, the limit of the curve must be known. Again bound the hyperboloid with two horizontal planes at the top and the bottom. Let the plane be x3 = h and x3 = -h and consider the portion of the hyperboloid between them. This 42 implies the following bounds for the parameter used to implicitly represent the curve on the hyperboloid: |u| h. To simplify the notation, consider the portion of the curve v2 -2uv -1 = 0 where u is within [-h, h]. The corresponding bounds for v are: |v| h + h2 + 1. Let bv = h + h2 + 1, so that -bv v bv. Map (u, v) [-h, h] × [-bv, bv] to (u , v ) [0, 1] × [0, 1] by putting u + h u = , 2h v + b v = v . 2bv The new parametrization of the surface is b2 x v(2v - 1)2 - 2hbv(2u - 1)(2v - 1) - 1 1 = , b2v(2v - 1)2 + 1 2b x v(2v - 1) + 2h(2u - 1)(b2v(2v - 1)2 - 1) 2 = , b2v(2v - 1)2 + 1 x3 = h(2u - 1). 43 Again, equating with the B´ezier form, the control points are -b2 hb2 c v + 2bvh + 1 v - h + 2bv 0,0 = , - , -h , b2v + 1 b2v + 1 b2 -h(b2 c v + 1 v + 1) 0,1 = , , -h , b2v - 1 b2v - 1 2b hb2 c vh + b2v - 1 v - h - 2bv 0,2 = , , -h , b2v + 1 b2v + 1 b2 -2b c v - 1 v 1,0 = , , 0 , b2v + 1 b2v + 1 b2 c v + 1 1,1 = , 0, 0 , b2v - 1 b2 2b c v - 1 v 1,2 = , , 0 , b2v + 1 b2v + 1 2b hb2 c vh + b2v - 1 v - h - 2bv 2,0 = , , h , b2v + 1 b2v + 1 b2 h(b2 c v + 1 v + 1) 2,1 = , , h , b2v - 1 b2v - 1 b2 hb2 c v - 2bvh - 1 v - h + 2bv 2,2 = , , h , b2v + 1 b2v + 1 and their weights are w2,1 = w0,1 = w1,1 = -b2v + 1, w0,2 = w1,2 = w0,0 = w2,2 = w2,0 = w1,0 = b2v + 1. 4.4.2 Trimming Curves We first discuss the trimming curves formed by the halfspace x1 0, then by general halfspaces. Recall that the surface is bounded by x1 0 to eliminate the missing line. The bounding curve in the (u, v) domain is v2 - 2uv - 1 = 0, which is a hyperbola. The curve is a degree- 2 B´ezier curve with three control points. By substituting the B´ezier curve Q(t) into the 44 implicit form, we compute the control points and weights of this curve in the uv-domain: c0 = (h + h2 + 1, -h, 1), 1 c1 = ( , 0, h2 + 1), h2 + 1 c2 = (-h + h2 + 1, h, 1), . The (u , v )-domain, they are c0 = (1, 0, 1), h h2 + 1 + h2 + 2 1 c1 = ( , , h2 + 1), 2h h2 + 1 + 2h2 + 2 2 h2 + 1 c2 = ( , 1, 1). h + h2 + 1 General clipping halfspaces are specified by inequalities of the form a1x1 + a2x2 d. Sub- stituting Q(u, v) into the inequality yields a1(v2 - 2uv - 1) + a2(2v + u(v2 - 1)) - d(v2 + 1) 0. The bounding curve is obtained by changing the inequality to equality. Solving for u, we obtain the equation d(v2 + 1) - 2a u = 2v - a1(v2 - 1) . a2(v2 - 1) - 2va1 Setting v = t, the rational representations of u and v in t, with u = f(t) and v = t = g(t), are r(t) r(t) f(t) = d(t2 + 1) - 2a2t - a1(t2 - 1), g(t) = t(a2(t2 - 1) - 2ta1), r(t) = a2(t2 - 1) - 2ta1. 45 These equations are in the power basis. They can be transformed into the B´ezier basis as described at the beginning of this chapter. 4.5 Conclusion We compute the B´ezier representation of the patches in the skin surfaces. The implemen- tation is on Silicon Graphics workstation with the OpenGL library. We did not further improve this implementation mainly because of the hardware. First, the implementation is not efficient because the patches are retriangulated for every frame. Secondly, there are gaps between patches because of the non-conformal triangulations of neighboring patches along the edges they share. This motivates us to search for a better approximation of the skin surface. 46 Chapter 5 Skin Triangulation We give the description of the skin surface triangulation in this chapter. The general idea is that if enough points are sampled on the skin surface according to the curvature, the surface can be approximated by a restricted Delaunay triangulation, which has the same topology [90]. This is referred to as the surface reconstruction problem. Several works have been done on reconstructing surfaces. Hoppe gives an algorithm for reconstructing surfaces represented by the zero set of a signed distance function [57]. Curless and Levoy also give an algorithm to reconstruct the laser-range data with estimations of error and surface normals [30]. The alpha shape software presented by Edelsbrunner reconstructs surfaces from uniform sampling [39, 44]. However, these algorithms do not analyze on the relationship between the original surface and the reconstructed triangulation. Amenta, Bern, and Dey construct algorithms based on the restricted Delaunay triangulation with a dense sampling of points [6, 4, 5]. They proved that the triangulation is homeomorphic to the surface if the points on the surface are sampled proportionally to the local feature size, which is the distance to the medial axis. The proof of the homeomorphism of the skin triangulation is similar to the previous works of Amenta, Bern and Dey. We use the local length scale instead of the local feature size. First, we introduce the restricted Delaunay triangulation in this chapter. Then, we discuss some properties of the curvatures and normals on the skin surface. Finally, we give 47 a proof of homeomorphism of the triangulation. 5.1 Restricted Delaunay Triangulation Given a surface F in R3, consider a finite set of points, V F . For each a V , the corresponding Voronoi cell a intersects F . F is partitioned into regions similar to the Voronoi diagram in R2. Finally construct a restricted Delaunay triangulation on the surface as if in R2, as shown in Figure 5.1. Formally, the restricted Delaunay triangulation is the Figure 5.1: Restricted Delaunay triangulation of a surface: Starting with a surface (left), we sample a set of points on the surface and intersect the surface with the Voronoi cells of the points (middle). Finally we construct the triangulation according to the Voronoi configuration. nerve of the collection of restricted polygons a F : DV = conv(U) | U V, F a = . a U If V is not properly sampled, the final triangulation may not have the same topology with F . Figure 5.2 shows an example of under-sampling, resulting in an incorrect triangulation. Over-sampling will improve the triangulation, but decrease the performance. In chapter 6, we describe the condition for sampling as well as the details of points placement. In this chapter, we describe the triangulations of skins, F = skin(B), of weighted point sets. 48 Figure 5.2: Restricted Delaunay triangulation of two circles in R2: When a Voronoi edge intersects the surface, there is an edge. A well sampled restricted Delaunay triangulation keeps the right topology (left). Wrong topology will result if the sampling is not done correctly(right). 5.2 Curvature and Normals of Skin Surface This section proves that the maximum curvature is continuous and satisfies a Lipschitz condition. Use this to control local density in the triangulation. This section also proves a one-sided Lipschitz condition for the normal direction. 5.2.1 Maximum Curvature Given a surface F , a point x on F and tangent vector tx, the normal curvature of F is that of a geodesic passing through x in the direction tx. The maximum curvature is the function : F R that maps x F to the maximum normal curvature at x. Define the slope of the asymptopic double cone to be a1/a2 for a hyperboloid in the form of a1x21 - a2(x22 + x23 + ... + x2n) = R2. For a hyperboloid of revolution with the slope of 1, the minimum curvature is measured within planes containing the symmetry axis (along meridians), and the maximum curvature is measured in planes normal to the axis (along latitudes). Explicit expressions for are easy to compute. If the sphere and hyperboloids are in the form of: x21 + x22 + x23 = R2, and x21 + x22 - x23 = R2, 49 the maximum curvatures are, respectively, = 1/R, and = 1/ R2 + 2x23 = (x21 + x22 + x23)-1. This implies that the maximum curvature at x is one over the distance of x from the origin. Iso-curvature Lemma. Every point x R3 belongs to exactly one hyperboloid of the form x21 + x22 - x23 = R2, and the maximum curvature of that hyperboloid at x is = 1/ x . 5.2.2 Curvature Continuity To prove that varies continuously over the skin surface, consider two infinite families of spheres that define the skin as their common envelope. For a finite set of spheres B, let S = conv(B) and T = conv(B ). The skin of B is F = ( S) = ( T ). The family of S defines F from inside, and T defines the surface from outside. For any point x F , there are unique spheres Sx S and Tx T that pass through x. We make essential use of the Sandwich Property stating that Sx and Tx have the same size in section 3.4.2. It is convenient to define (x) = 1/ (x) as the radius of maximal curvature, and for reasons that will become clear later we refer to (x) : F R as the local length scale. Curvature Sandwich Lemma. For every point x F , the local length scale (x) is the common radius of Sx and Tx. Proof. If x belongs to a sphere patch, then that patch either lies on Sx or on Tx and (x) is obviously the radius. Now suppose x belongs to a hyperboloid patch. The hyperboloid is obtained by revolving a hyperbola around one of its two symmetry axes. As indicated in Figure 5.3, the hyperbola is the common envelope of two families of circles, one centered along each of the two symmetry axes. By the Sandwich Property, Sx and Tx have equal 50 radii. Because x is halfway between the centers of Sx and Tx, that radius is equal to the distance of x from the origin. By the Iso-curvature Lemma, this distance is x = (x). a x Sx b Tx Figure 5.3: The circles Sx and Tx sandwich the hyperbola: Depending on whether we revolve the hyperbola around the vertical or the horizontal axis, we get a one-sheeted or a two-sheeted hyperboloid. The sandwiching spheres, and their common radius, vary continuously with the point x F . As stated in the Sandwich Property, each point on the surface is sandwiched by a unique pair of sphere with the same size. Since this is true also for the points at the boundary of the mixed-cell, the Curvature Sandwich Lemma thus implies that the maximum curvature varies continuously over the skin surface (except at centers, where it blows up). In fact, at every point x the local length scale (x) equals the distance from x to the center z X of the mixed cell µX that contains x. 5.2.3 Curvature Variation We now strengthen the result that (x) is continuous by showing that it varies rather slowly. In fact, we can extend its reciprocal (x) to a function defined on all of R3 and show that (x) has Lipschitz constant 1. As we have seen, within any mixed cell µX, is simply the distance to the center z = z X. By the definition of the mixed complex, this is a continuous 51 function in R3. Within µX, the triangle inequality gives the Lipschitz bound, | (x) - (y)| = | x - z - y - z | x - y . By applying this to the pieces of the line segment from x to y contained in different mixed cells, we obtained the following result. Curvature Variation Lemma. For all points x, y in space, | (x) - (y)| x - y . Note that the extension of to a function R3 R describes the length scale of all surfaces in the family defined by the growth model of deformation. 5.2.4 Normal Variation The tangent or C1-continuity of the skin surface follows from the Sandwich Property. We strengthen this result by proving a one-sided Lipschitz condition for the normal vectors. Specifically, we prove an upper bound that relates the angle between two normal vectors at points x, y to the Euclidean distance between them and to their length scales. The outward unit normal vector at x F is denoted nx, and the angle between two normals is nxny = arccos nx, ny . In proving the upper bound, we consider again the one-parameter family of skin surfaces generated by increasing squared radii with time. For points x = (x1, x2, x3) on a sphere in standard form the unit normals are nx = ±x/ x , and for points x on a hyperboloid in standard form they are nx = ±(x1, x2, -x3)/ x . In both cases, the normals are the same along a line passing through the origin, and they vary their directions as the point is rotated about the origin. The formulas imply that the normals of points x and y in two adjacent mixed cells are the same if x and y are mirror images of each other across the separating plane. This property is illustrated in Figure 5.4. Normal Variation Lemma. Let x and y be points on F with distance x - y < (x). The angle nxny between the surface normals at x and y is at most arcsin( x - y / (x)). 52 y' y p1 x Figure 5.4: Normal variation: The normal of a point y can be found by reflecting by the intersection point pi. Proof. Consider first the case where x and y belong the same mixed cell, and translate the coordinates so that the center is at the origin. Given x and the distance x - y between the two points, the angle subtended at the origin is a maximum if x 2 = x - y 2 + y 2. Hence: x - y nxny arcsin , (x) as claimed. Consider second case where x and y lie in different mixed cells. The directed line segment from x to y passes through i 1 planes h1, h2, ..., hi separating adjacent mixed cells. Let pj = xy hj be the intersection points ordered from x to y. Construct a polygonal path that starts at x whose length is x - y . It is obtained from xy by reflecting the portion after pj across the plane hj, for j = i, i - 1, ..., 1 in this order. The endpoint y of the path is contained inside the sphere with radius x - y around x, which implies that the angle between x and y subtended at the origin is < arcsin( x - y / (x)). Since ny is normal to the sphere or hyperboloid defined for the mixed cell of x that passes through y , is also the angle between nx and ny. The claim follows. The proof of the Normal Variation Lemma does not require that x and y belong to the same skin surface. The claimed inequality holds more generally far any points x, y R3 with 53 normals defined by the one-parameter family of skin surfaces. Suppose the distance between x and y is x - y < (x). Then the Normal Variation Lemma implies nxny < arcsin , which is the form used most often in this work. 5.3 Triangulation A finite set V F is an -sampling if for every point x F there is a vertex a V whose distance from x is x - a < (x). The goal of this section is to prove that the restricted Delaunay triangulation defined by an -sampling is homeomorphic to the skin surface, provided the following condition holds: 0 0 (I) where 0 = 0.279... is a root of 2 2 f( ) = 2 cos arcsin + arcsin - . 1 - 1 - Note that f( ) is defined for -1 1, and that it is non-negative for 0 3 0. The goal of this section is to prove that, for < 0.279..., the restricted Delaunay triangu- lation, DV is a triangulation of F = skin(B). Following the standard topology terminology, this means the underlying space of DV is homeomorphic to F . As shown by Edelsbrunner and Shah [47], it suffices to prove that every non-empty common intersection of restricted Voronoi polygons is a closed topological ball of the appropriate dimension, namely 3 minus the number of polygons. If this is the case, DV has the closed ball property. Formulate this property in terms of the (unrestricted) Voronoi polyhedra. Here general 54 position is assumed the intersection of k + 1 = 2, 3, 4 Voronoi polyhedra is a polygon, edge and vertex respectively. Depending on the case, the intersection with the skin surface is: * case k = 0: a closed disk, * case k = 1: empty or a closed interval, * case k = 2: empty or a single point, * case k = 3: empty. The case k = 0 corresponds to a single Voronoi polyhedron, which has non-empty intersec- tion with F because its generating point lies on F . We establish four technical lemmas in preparation of proving that DV has the closed ball property. 5.3.1 Distance Claims If two surface points lie in the same Voronoi polyhedron, then they cannot be far from each other, and if they lie on a line that is almost normal to the surface, then they cannot be close to each other. We quantify the first claim under the assumption that V is an -sampling. Short Distance Claim. If points x and y on F belong to a common Voronoi polyhedron defined by a vertex in an -sampling V F , then x - y < (2 /(1 - )) (x). Proof. Let a be the generating point of the common Voronoi polyhedron. By the - sampling assumption, x - a < (a) and y - a < (y). Using triangle inequality, x - y < ( (x) + (y)). The Curvature Variation Lemma now implies (x) (x) - x - y > (1 - ) (y) - (x), 55 and hence (1 + ) (x) > (1 - ) (y). The distance between x and y is therefore: 1 + x - y < 1 + (x) 1 - 2 = (x), 1 - as claimed. The bound on the distance is better if one of the points generates the Voronoi polyhedron. Assuming x = a, we have a - y < (y) (a) + a - y , which implies a - y < (a). 1 - This version of the Short Distance Claim is needed in the proof of the Voronoi Polyhedron Lemma below. Next quantify the second claim, which is independent of V . x y Skin surface (x) Figure 5.5: Long distance claim: y is out of the sphere touching x with radius (x). The length of xy must be larger than 2 (x) cos . Long Distance Claim. Suppose a line meets F in two points x and y and forms an angle smaller than with the surface normal at x. Then x - y > 2 (x) cos . 56 Proof. By the Sandwich Property, there are two spheres of radius (x) that both pass through x and locally sandwich the surface, see Figure 5.5. The line meets the two spheres at x and at points at distance larger than 2 cos (x) on both side. The skin surface contains no points inside either sandwiching sphere, which implies the claimed lower bound for x - y . Playing the Short and Long Distance Claims against each other, we reach contradictions proving various claims. 5.3.2 Normal Lemmas If the vertices of a short edge or a triangle with small circumcircle lie on the skin surface, then the edge or triangle lies almost flat. We quantify both claims. For an edge ab let tab = (b-a)/ b - a be the unit vector of edge ab. The first result is an immediate corollary to the Long Distance Claim: Edge Normal Lemma. The angle between and edge ab and the surface normal at its vertex a is tabna > /2 - arcsin( a - b /2 (a)). A common use of the Edge Normal Lemma is when ab belongs to the restricted Delaunay triangulation of an -sampling. Then (a) > (1- ) (x), where x is a point in the intersection of the dual Voronoi polygon with the skin surface. Hence a - b < 2 (x) < 2 (a)/(1- ). The angle between ab and the surface normal at a is then: tabna > - arcsin . 2 1 - Next consider the triangle normal lemma. Assume the angle at a inside the triangle abc is no smaller than the angles at b and c. Let Rabc be the radius of the circumcircle and let nabc be the outward unit normal vector of abc. 57 Triangle Normal Lemma. If a is a vertex of the triangle abc with greatest angle, then the angle between the normal of abc and the surface normal at a is nabcna < arcsin(2Rabc/ (a)). z ( ) a A b a c Figure 5.6: Triangle normal lemma: The dashed sandwiching spheres meet the solid sphere around a in two parallel dotted circles. Vertices b and c are placed to maximize the angle between the triangle normal and the surface normal at a. Proof. Consider the two spheres of radius (a) that locally sandwich the surface at a, as shown in Figure 5.6. The face angle at a is at least /3 and the length of the edges ab and ac is at most 2Rabc each. To compute a bound on the angle between na and nabc, assume a - c a - b and consider the sphere with radius a - b around a. It intersects the sandwiching spheres in two parallel circles. Let 2X be the distance between these two circles and not that X/ a - b = a - b /2 (a) by dropping perpendicular from z to the midpoint of ab and using similar triangles. Hence 2X = a - b 2/ (a). Since the angle at a is greater than or equal to the ones at b and c, bc is the longest edge of abc. The angle between the edge bc and the planes of the intersection circles is therefore less than a - b 2 2R arcsin arcsin abc . b - c (a) (a) This is an upper bound for the angle between the two normal vectors at a. 58 Suppose that abc belongs to the restricted Delaunay triangulation of an -sampling. Then (a) > (1 - ) (s), where x is a point of the intersection between the dual Voronoi edge and the skin surface. Hence Rabc < (x) < ( /(1 - )) (a). The angle between the two normals at a is then 2 nabcna < arcsin 1 - 5.4 Closed Ball Property It is now ready to prove the closed ball property for the restricted Delaunay triangulation, assuming V is an -sampling of F where < 0.269.... Assume general position and consider the three cases in turn: first Voronoi edges, then Voronoi polygons, and finally Voronoi polyhedra. 5.4.1 Voronoi Edge Voronoi Edge Lemma. A Voronoi edge of V intersects the skin surface in at most one point. Proof. Assume there is a Voronoi edge that intersect F in at least two points, x and y. Let abc be the dual triangle in the restricted Delaunay triangulation. The Triangle Normal Lemma gives an upper bound for the angle between the normal of abc and the surface normal at a. The Normal Variation Lemma gives an upper bound for the angle between the surface at a and x. Together they imply an upper bound for the angle xi between nabc and nx: nabcna + nanx 2 arcsin + arcsin . 1 - 59 The angle is also the angle between the Voronoi edge and nx. The Long Distance Claim implies x - y > 2 (x) cos , which because < 0 contradicts the upper bound x - y < (2 /(1 - )) (x) implies by the Short Distance Claim. 5.4.2 Voronoi Polygon Voronoi Polygon Lemma. The intersection of a Voronoi polygon of V with the skin sur- face is either empty or a closed topological interval. L L L x x x y y y Figure 5.7: Voronoi polygon lemma: A Voronoi polygon intersecting the skin in a circle to the left and two intervals to the right. Proof. Assume there is a Voronoi polygon whose intersection with the skin surface contains a topological circle or two topological intervals, as shown in Figure 5.7. Let ab be the dual edge in the restricted Delaunay triangulation, and let x be an arbitrary point of the intersection. If x lies on a circle, then let L be the line in the plane of the polygon that intersects the circle in a right angle at x. Lnx L nx for any line L in the same plane and passing through x. Choose L to minimize the angle with na. The Edge Normal Lemma implies an upper bound for the angle between L and the surface normal at a. The Normal Variation Lemma implies an upper bound on the angle between the surface normals 60 at a and x. Together these inequalities imply Lnx arcsin + arcsin . 1 - This angle is less than the upper bound for xi in the proof of the Voronoi Edge Lemma, which implies a contradiction between the two distance claims. In the case of two intervals let L be a line connecting x to the closest point y on the other interval. If y lies in the interior, then L intersects the interval in a right angle at y. In this case we get a contradiction with the same argument as above only with x and y interchanged. Otherwise, y is an endpoint of the interval and lies on a Voronoi edge. The angle between L and ny is less than that between the Voronoi edge and ny. Thus a contradiction occurs with the same argument as used in the proof of the Voronoi Edge Lemma. 5.4.3 Voronoi Polyhedron Voronoi Polyhedron Lemma. The intersection of a Voronoi polyhedra of V with the skin surface is a closed topological disk. Proof. Assume there is a Voronoi polyhedron whose intersection with the skin surface con- tains a closed 2-manifold (without boundary), a 2-manifold with boundary other than a circle, or two circles. In the first case let L be a line that intersects the 2-manifold in two points, x and y, and forms a right angle at x. A contradiction occurs between the two distance claims as before. For the rest of the proof, let a be the generating vertex of the Voronoi polyhedron and assume the intersection between this polyhedron and the skin surface is a 2-manifold with boundary, F . This 2-manifold with boundary can be different from a disk either because it is non-orientable, it contains a handle, or it has at least two boundary circles. The non-orientability of F contradicts the orientability of F . If F has a handle but only one 61 boundary circle, then homology theory gives us a pair of simple closed curves in F that intersect each other transversely exactly once. Along either one of these curves, there is a point such that the line normal to F that passes through that point meets the other curve, and hence F again. This gives a contradiction to the two distance claims. A more elaborate argument is needed for the case where there are two or more boundary circles. Then either F is connected, and in the simplest case is an annulus, or it is disconnected, and in the simplest case consist of two disks. By the remark after Short Distance Claim, the distance between a and a point y F is a - y < ( /(1 - )) (a). Let L be the normal line at a and note that it contains the line segment of length 2 (a) that connects the centers of the two spheres sandwiching the surface at a. This line segment is contained in the Voronoi polyhedron, which implies that the polyhedron is fairly tall and slim. Consider a plane that contains L and intersects at least two boundary circles of F . Such a plane exists for else a plane can be found through L that intersects no boundary circle at all. However, then L meets F in at least two points, and again a contradiction occurs to the two distance lemmas. The plane that meets two boundary circles intersects the Voronoi polyhedron in a convex polygon and F in at least two connected curves. One of the curves contains a. The second curve may be assumed to lie on one side of L. Let L be the line passing through its two endpoints, which both lie on the boundary of the convex polygon. The line L intersects the sphere with radius (a) /(1- ) around a and it does not intersect the line segment connecting the centers of the two sandwiching spheres. The angle between L and the surface normal at a is therefore naL arcsin( /(1- )). By the intermediate value theorem there is a point y on the second curve whose curve normal n y is also normal to L . Hence nan y /2 - arcsin( /(1 - )). Since the surface normal at y is also normal to the tangent line parallel to L , its angle with na is at least this large. From the Normal Variation Lemma, yield nany < arcsin( /(1 - )). Putting both inequalities together, /2 < 2 arcsin( /(1 - )). This is equivalent to > 2 - 1 = 0.414... and contradicts Condition (I). This completes the proof of the Voronoi Polyhedron Lemma 62 for the final case where F has at least two boundary circles. 5.5 Summary The three Voronoi lemmas establish that for any < 0.297..., the restricted Voronoi diagram of an -sampling V has the closed ball property. The result of [47] implies that the underlying space of the restricted Delaunay triangulation is homeomorphic to the skin surface. General Homeomorphism Theorem. The restricted Delaunay triangulation of an - sampling triangulates the skin surface, for < 0.279.... For the purpose of changing the topology of the skin surface we rely on point distribu- tions that locally violate the -sampling condition. We give a separate proof of the closed ball property in Section 7 and thus obtain a Special Homeomorphism Theorem for such distributions. 63 Chapter 6 Static Skin Triangulation Algorithm We concluded that the restricted Delaunay triangulation with -sampling of a skin surface is homeomorphic to the analytical skin in the previous chapter. Here, we discuss the implemen- tation of skin mesh construction, including analysis of the algorithm. For skin construction, we discuss construction of the restricted Delaunay triangulation and the refinement. Refine- ment improves the mesh quality by reducing the minimum angle of the triangles on the sur- face. These algorithms have been proved to give guaranteed-quality mesh [8, 13, 20, 27, 32]. Most of the works prove that their algorithms give a good quality mesh and terminate, if the surface is smooth or some bounds exist on the variation of the surface normals. For the skin case, we know the surface normal and curvature properties of the surface, and we prove our algorithm is correct and terminates based on these properties. First, we provide a description of the curvature adaptation of the algorithm. Then we discuss sampling a subset of points, V F = skin(B), according to curvature. After constructing the restricted Delaunay triangulation, a refinement algorithm is applied to ensure the triangle quality. At the end of this chapter, we analyze the algorithm to ensure correctness. 64 6.1 Curvature Adaptation The goal of our algorithm is to triangulate with edges and triangles of size roughly propor- tional to the local length scale, (x) = 1/ (x). Define the size of an edge ab as half of its length, (Rab = a - b /2) and the size of a triangle abc as the radius Rabc of the circumcircle. The length scale at the vertices of an edge or triangle are roughly but usually not exactly the same. We worry about the edges getting too short, so we compare their sizes with the maximum length scale at any of their vertices. For triangles, we worry about them growing too large, so we compare their sizes with the minimum length scale of their vertices. We define ab = max{ (a), (b)}, and abc = min{ (a), (b), (c)}. Our algorithm is formulated using two positive constants, C and Q. C controls how closely the triangulation approximates the skin surface, and Q controls the quality of the triangles. The following two inequalities are maintained as invariants, which we refer to as the Lower Size Bound and the Upper Size Bound, respectively: C Rab > · Q ab for every edge ab D. [L] Rabc < CQ · abc for every triangle abc D. [U] It is not necessary to check for long edges and small triangles explicitly. This is because an edge of size Rab CQ ab belongs to two triangles that both violate [U]. Symmetrically, a triangle of size Rabc (C/Q) abc has three edges that violate [L]. Appropriate values of C, Q will be determined in the analysis of the algorithm but we already anticipate C = 0.08, Q = 1.65 as a feasible assignment in section 6.5.1. 65 Minimum Angle Lemma. A triangle that satisfies Conditions [L] and [U] has minimum angle larger than arcsin(1/Q2). Proof. Let abc be a triangle, with bc its shortest edge and R its circumradius. abc/ bc 1 by definition of length scale. Using Conditions [L] and [U], we have R CQ Q2 abc . b - c 2(C/Q) bc 2 The minimum angle is bac, and b - c = 2R sin bac. Hence bac = arcsin( b - c /2R) arcsin(1/Q2), as claimed. The Minimum Angle Lemma suggests choosing Q as small as possible, contingent upon satisfying all constraints needed to prove the algorithm correct. For Q = 1.65 from section 6.5.1, the minimum angle is larger than 21.54... , and the maximum angle is smaller than 180 - 2 · 21.54 = 136.90... . 6.2 Generating a Sample Within each mixed cell, there is either a sphere or hyperboloid patch. We sample points on the patch within each mixed cell for the restricted Delaunay triangulation. In order to have a distribution according to the curvature, we first estimate the number of sample points needed in each patch. Then we generate random points on the surface with a probability distribution proportional to the curvature. We estimate the number of points on a patch by multiplying the weighted area by a constant c, which is the number of points per unit-weighted area. We define the weighted area of a patch to be the integral of the curvature (x) over the area. So if A is the patch, the weighted area is WA = (x)dA. For example, the weighted area of a unit sphere is simply its surface area 4 because (x) = 1 for all x on the sphere. Now we discuss how to compute weighted area on a hyperboloid. Parametrize a hyperbola 66 by the Klein model K(p) = (k1, k2) = ( p , 1 ). A hyperboloid is formed by rotating 1-p2 1-p2 points with height x3 = k1 and radius x21 + x22 = k2 around the x3 axis, or swapping k1 and k2 depending on whether it is a 1-sheeted or 2-sheeted hyperboloid respectively. Also, since the hyperboloid is infinite, a bound pmax is needed in the parametric domain. Differentiating K yields K(p) 1 p = , . p (1 - p2) 1 - p2 (1 - p2) 1 - p2 The length scale of a point on the hyperboloid is the distance to the center, so 1 + p2 (K(p))2 = K(p) 2 = . 1 - p2 Define the derivative of length scale to be K(p) 2 ds2 = p dp 1 + p2 = dp2. (1 - p2)3 For 2-sheeted hyperboloid, we rotate a point with radius k1 around x3-axis, which gives the weighted area as: pmax 2 k W 1 A(pmax) = ds 0 2(p) pmax p = 2 dp 0 (1 - p2) 1 + p2 2 1 + p -1 + p pmax = arctanh - arctanh . 2 2(1 + p2) 2(1 + p2) 0 For 1-sheeted hyperboloid, we have the radius of revolution of each point as k2 which gives 67 the weighted area: pmax 1 WA(pmax) = 2 dp 0 (1 - p2) 1 + p2 2 1 + p -1 + p pmax = arctanh + arctanh 2 2(1 + p2) 2(1 + p2) 0 Next we determine the constant c which is number of point per unit weighted distance. We compute c using only the sphere since the resulting constant will be consistent for both spheres and hyperboloids. The minimum length of an edge of the triangle on the unit sphere is 2C/Q. Each point occupies a disk with approximate radius C/Q and thus approximate area C2/Q2. Hence, we choose c = Q2 4C2 . Another way of creating the initial triangulation is meshing independently within each mixed cell, which does not require three-dimensional Delaunay triangulation. Since the patches can be parametrized, we can create a mesh for each patches and then connect the neighboring ones. However, the main difficulty is the configuration of patches, or in another word, the connectivity of clipped patches. Clipping spherical or hyperbolic patches within their mixed cells may result in very complicated or degenerated patches. Previous work has been done in a similar way without the knowledge of connectivity of the patches and suffers the problem of lacking of robustness [22]. The connectivity is computed by clipping the triangulation with the mixed cells numerically, i.e., if a certain face of the mixed cell cuts the mesh, the patch is connected with another patch in the other mixed cell and their meshes are going to be merged. The main fault is that the clipping planes of the mixed cells may miss the triangulation but actually they intersect the real surface. The combinatorial information of the connectivity of clipped patches is essential to ensure the integrity of the composed surface and left as an open problem. 68 6.3 Surface Construction With a sampled point set, we construct the restricted Delaunay triangulation. First we compute the Delaunay triangulation of the point set. Then, we keep those triangles which are either: * On the convex hull and the attached tetrahedrons have their centers in the skin body, or * Not on the convex hull but one of the attached tetrahedrons has its center in the skin body and the other has its center outside the body. We call it the in-body test which picks the triangles that belongs to the restricted Delaunay triangulation since their Voronoi edges intersecting the skin surface. The next step is to further refine the triangulation in order to satisfy Conditions [L] and [U]. 6.4 Refinement Operations Random points sampled on the skin surface usually do not result in a high-quality trian- gulation. In order to satisfy conditions [L] and [U], three refinement operations are used to improve the quality. These operations are edge flipping, edge contraction, and point in- sertion. The edge contraction operation ensures condition [L] and point insertion ensures condition [U]. Then edge flipping maintains the restricted Delaunay triangulation. First, the three operations are described. After each operation, the local neighborhoods are modified and may lead to other operations. The algorithm pushes `suspicious' edges and triangles onto three stacks: - Flip Stack (FS) for edges - Contraction Stack (CS) for edges - Insertion Stack (IS) for triangles 69 The elements in the stacks are popped, checked, and refined if necessary. The algorithm will terminate when all the stacks are cleared. The order of popping does not affect the correctness of the algorithm since there is no infinite loop (proved in section 6.5.2). For efficiency, we pop the flip stack first, then the contraction stack, and finally the insertion stack. The reason for this order is that the point insertion do not cause edge contraction, so we do not need to check the contraction stack again. The flip stack is popped first because the flipping operation improve the condition of the triangles, namely, reducing triangle sizes and increasing edge lengths. After each operation, new `suspects' are created and pushed onto the corresponding stacks. 6.4.1 Edge Flipping b b d d flip(ab) a a c c Figure 6.1: Non-flippable edge: One more triangle abd (duplicated) is created after flipping the edge ab Flip Check. Before flipping an edge, we check to see that the edge is `flippable' as well as belonging to the restricted Delaunay triangulation. An edge is not flippable if the flip could create a second copy of an already existing triangle; see Figure 6.1. This problem will occur if either a or b is degree 3. After this combinatorial test, we check whether this edge belongs to the restricted Delaunay triangulation as described below. If not, we flip the edge to restore the local Delaunay-hood. To test whether an edge ab belongs to the restricted Delaunay triangulation with neigh- boring triangle abc and abd, compute the orthocenter of the tetrahedron abcd. If the center 70 c c b b a a d d Figure 6.2: Edge flipping lies within the skin body, the edge is locally Delaunay because the Voronoi edges of the neighboring triangles intersect the surface. However, the test is expensive. First, it involves computing the orthocenter of the 4 vertices. Secondly, we check whether this center is within the skin body. This involves solving matrices and algebraic equations. We improve the prac- tical efficiency of this test by approximation. For the two configurations of triangles (before flipping and after flipping), we choose the one with the larger minimum angle among the triangles. We flip the edge if the minimum angle of the two triangles is increased after flip- ping. Since locally the triangulation is very close to a plane, this test is an approximation of the flip check. We do not have a prove for the equivalence of this approximation with the in-body check in Section 6.3, Because of the neighborhood of an edge of the mesh is close to a plane due to -sampling, the second method can approximate the in-body test and experimental results shows that the speed is improved. Post-operation. After flipping the edge, new triangles are formed. Local Delaunay-hood may not be preserved. Let triangle abc and dba be two adjacent triangles and edge ab is flipped. After the operation, we push the following: * ac, cb, ad and bd to FS * adc and dbc to IS * cd to CS 71 6.4.2 Edge Contraction An edge contraction occurs if an edge ab is too short (violating condition [L]). However, an edge ab is not contractable if it has the link conflict [34]. If no conflict exist, contract the edge by edge contraction through point deletion, since contracting an edge by snapping one vertex to another may cause folded triangles. Link Conflict Check. To preserve the topology, the link conflict check is essential [34]. An edge ab can be contracted only if link(a) link(b) = link(ab). e b b a a a d Figure 6.3: Edge contraction: (left) before contraction ab is too short and b is going to be deleted ; (middle) flip with inD() test until b has degree 3;(right) b is deleted . Edge Contraction by Point Deletion. Contract the edge ab by deleting the point b together with all edges and triangles in the star. The process is similar to deleting a point in a Delaunay triangulation in 2D [33]. The star of b is a polygon. For each vertex p on this polygon let p- and p+ be its predecessor and successor in an ordering around the polygon. The algorithm converts the star of b by creating one triangle p-pp+ at a time by flipping. For every edge bp, check the points p, p- and p+ with the following inD() test. The edge bp will be flipped if inD() returns true. Function inD(). Any three points p, p- and p+ define a circumcircle with center c. The function inD() return true if: 72 1 The projection of b on the plane of the circle is inside the circle, and 2 All other vertices on the polygon, except p, p-, and p+, are all outside of the circle. Condition 2 maintains the Delaunay property of the triangulation after deletion of the point b. Condition 1 checks if the angle of a corner of the polygon exceeds , for example, the edges da and ae in figure 6.3. We say they are concave edges. Flipping the edge pb can cause a folded triangle even if condition 1 is satisfied. Post-operation. In the flipping stage, edges and triangles are pushed to the stacks as in the normal flipping operations. After the degree of b reaches 3, delete the three adjacent triangles and fill the hole with a new triangle. Only this new final triangle is pushed onto IS. 6.4.3 Point Insertion If a triangle abc is too big, it violates condition [U]. An additional point is added to destroy this triangle. However, before adding a point, we check to see if the edges of abc can be flipped. If none of the edges can be flipped, we try to add a point. We perform these operations in this order because adding a new point may not destroy abc if the triangle has a large angle and its circumcenter is far away. However, in this case, the triangle cannot be part of the restricted Delaunay triangulation, so flipping an edge removes it. To add a new point, let L be the line of points at equal distance from a, b and c. Starting from triangle abc, walk on the triangulation until the triangle that intersects L is found. The walk is a direct path to the line since for each step we head for the triangle closer to the line L. After locating the triangle, intersect with the skin patch by solving the root of a quadratic equation and insert the intersection point x. After inserting the point, the local region of triangles may not keep the restricted De- launay property anymore. We repair this by flipping all the edges in the link of x until all 73 edges in link(x) satisfy the flip check. Post-operation. The triangles in the star of x are pushed onto IS. This is sufficient because no short edges are created by point insertion and the flipping establishes the restricted Delaunay property of the new triangulation. 6.5 Analysis We now analyze the algorithm and prove its correctness. First we prove that the algorithm terminates. Next we prove that the result is the restricted Delaunay triangulation of an -sampling, which implies that the triangulation is homeomorphic to the skin. 6.5.1 Conditions We first prove that point insertions do not generate edges that violate the Lower Size Bound [L]. That proof requires that Q is not too large. We then prove that the restricted Voronoi vertex dual to a triangle can be found near the circumcenter of that triangle. That proof requires that the vertices of the triangulation form an -sampling. Finally, we prove that the vertices indeed form an -sampling, with satisfying Condition (I), i.e. < 0 = 0.279.... The closed ball property established in Section 5 then implies that the triangulation produced is homeomorphic to the skin surface. That proof relies on the quality of the approximation, which is guaranteed by the algorithm provided CQ is not too large. For ease of reference we collect the conditions before deriving them. Q2 - 4CQ - 2 > 0, (II) 2/(1 + )2 - 4/4 > C2Q2, (III) 74 where = - 2C( + 1)/(Q + 2C). We get (II) and (III) as sufficient conditions for the proofs of the No-Short-Edge Lemma and the Sampling Lemma below. Condition (II) is equivalent to Q > 2C + 4C2 + 2. Assuming = 0 = 0.279..., we can satisfy Conditions (II) and (III) by setting C = 0.08 and Q = 1.65. In this case = 0.166.... Small improvements are possible. 6.5.2 Refinement Termination An edge contraction may cause other edge contractions, but this cannot go on forever because we eventually violate the Upper Size Bound. Similarly, a vertex insertion may cause other vertex insertions, but this cannot go on forever because we eventually violate the Lower Size Bound. It is possible that an edge contraction causes vertex insertions, but a vertex insertion cannot cause edge contractions. This is because a vertex insertion cannot create edges of size below the allowed threshold. This prevents infinite loops in spite of the algorithm's partially conflicting efforts to avoid short edges and large triangles simultaneously. Let abc be the triangle that causes the addition of the dual restricted Voronoi vertex x F . No-Short-Edge Lemma. Every edge xy created during the addition of x has size larger than (C/Q) xy. Proof. We have Rabc CQ abc. The sphere with center x that passes through a, b, c has radius X Rabc and it contains no vertex other than x inside. Every new edge xy has therefore length x - y X CQ abc. Assume without loss of generality that abc = (a). We use the Curvature Variation Lemma to derive upper bounds for the length scales at x and y: 1 (x) (a) + X + 1 x - y , CQ 1 (y) (x) + x - y + 2 x - y . CQ 75 Hence x - y max{ (x), (y)} Rxy = . 2 4 + 2/CQ Condition (II) implies C/Q < CQ/(4CQ + 2), and therefore Rxy > (C/Q) xy, as claimed. 6.5.3 Close Dual Vertices Consider the point addition triggered by the triangle abc violating the Upper Size Bound. As before, we denote the line of points at equal distance from a, b, c by L, the circumcenter of abc by z, and the point of L F closest to z by x. We prove an upper bound on the distance between x and z assuming an -sampling of F . Circumcenter Lemma. The distance between x and z is x - z < ( 2/2) abc. x a z Figure 6.4: Circumcenter lemma: Dashed sphere of radius (x) passing through a, b, c, x and bold circle with center z passing through a, b, c. Proof. Assume abc = (a) (b), (c). We have x - a (x) by assumption of - sampling and therefore abc/(1 + ) (x) by the Curvature Variation Lemma. We get an upper bound on the distance between x and z by assuming (x) is as small as possible and a,b,c lie on the sandwiching sphere with radius (x) = abc/(1 + ) passing through x. This configuration is sketched in Figure 6.4. Note that x - z / x - a /2 (x) by equality of angles formed by orthogonal sides. Therefore, x - a 2 2 x - z = (x), 2 (x) 2 76 as claimed. The relevance of the Circumcenter Lemma to the curvature adaptation algorithm should be obvious. When the triangle abc violates Condition (III), we need first to find its dual vertex in the restricted Voronoi diagram and then add this vertex to V . This vertex is the point x, and the Circumcenter Lemma gives a bound on how far from z we have to search before we are guaranteed to find x. As shown in the proof of the Voronoi Edge Lemma, each additional point y L F is too far from z to possibly belong to the Voronoi edge dual to abc. 77 Chapter 7 Dynamic Skin Algorithm The algorithm for constructing the skin triangulation dynamically is described in this chap- ter. Assuming we have a skin triangulation D0 at time t0, we grow the skin gradually to another triangulation D1 at t0. We call D0 and D1 the states of the mesh. During the growth, triangles are distorted and violate the Conditions [L] and [U]. These are the invalid triangles and refinements are applied. Once the triangle is refined and valid, we compute a time step t when the next refinement is needed. The triangle gets a certificate and is guaranteed to be valid before the period of t is passed. We call t the time step. We call this framework the Kinetic Data Structures (KDS) [2, 3, 9, 56, 59]. The Ph.D thesis of Basch is a good resource of KDS [9], and he states that KDS are simply made of two components: a proof of correctness of the attributes (Conditions [L] and [U] in this case), and a priority queue called the event queue. The events in the event loop are topology changes and refinement operations. The event loop starts with a correct triangulation (which may be empty). The algorithm processes the events according to their time. For each event of topology changes, we know when it happens by the alpha value of the corresponding simplex. At that instant, the curvature will rise to infinity, so our previous triangulation contains infinitely small triangles because of Condition [L]. Special treatment is needed in the region that is close to the position of topological changes. In these special regions, local portions of the triangulation are not maintained anymore within short time intervals as described in the previous chapter. The Upper Size 78 Bound and Lower Size Bound are maintained at the times when the special sampling begins and finishes. This special treatment is described in Section 7.4. For refinement operations, the oracle of t computation involves solving equations of high degrees. The difficulties come from the non-linear paths of the points with the changes of its local length scale. Alternatively, we estimate a time to check an edge or triangle is valid, which is described in Section 7.5. In this chapter, we first prove that the density is maintained during growing. Secondly we describe how to update the vertices on the skin triangulation followed by the overall scheduling algorithm. The refinement operations are described in the previous chapter, and we give the details of special sampling of the topology event in this chapter. At the end of this chapter, we give a method to estimate t for the checking of both edges and triangles. We prove the edge estimation is correct but not for the triangles. An alternative experimental strategy is used to attack this problem. 7.1 Maintaining Density We now show that the algorithm for curvature adaptation maintains the -sampling property of the vertex set. Recall that this means that for every point x F there is a vertex a V whose distance from x is a - x < (x). The constant is to be chosen less than 0 = 0.297.... It is interesting to see that the two Size Bounds by themselves are too weak to imply -sampling. We can put four points near each other on a sphere in such a way that all four triangles and six edges satisfy [L] and [U]. Nevertheless, the boundary of the tetrahedron is miserably inadequate approximation of the sphere surface. We argue that the algorithm cannot get to this problematic state, because of the way it would temporarily have to violate the two Size Bounds. In other words, we use continuity in time to prove the claim on sampling. In stating the result, we assume the skin surface deforms continuously with 79 time. The details of the motions are discussed in the next chapter. For now we disallow metamorphoses, which simply means the topology changes while the skin grows. Let t0 < t1 be two points in time so the topological type is constant within [t0, t1]. We write F (t) for the skin surface at time t and V0, V1 for the vertex sets at time t0 and t1 respectively. Sampling Lemma. If V0 is an -sampling of F (t0), then V1 is an -sampling of F (t1). Proof. Assume the opposite and let t [t0, t1] be the first moment in time when the skin surface is not -sampled. Then there is a point x F (t) such that no vertex lies inside the sphere with center x and radius X = (x). By minimality of t, the sphere increase the sphere while keeping its center on the surface. Vertex a remains on the sphere and we permit no vertices inside the sphere. Let y F (t) be the center when we reach the other two vertices, b and c. The radius of the new sphere is Y X because the radius can only increase from x to y. Using the Curvature Variation Lemma, we get abc (x)(1 + ) and therefore ( /(1 + ) (x) Y . Assume without loss of generality that abc = (a), and let z be the circumcenter of abc. The Upper Size Bound implies z - a = Rabc < CQ (a). Using the Circumcenter Lemma, we get an upper bound on the square distance between y and a, Y 2 = y - z 2 + z - a 2 4 < 2 4 abc + C2Q2 2abc. This implies 2 4 < + C2Q2, (1 + )2 4 which contradicts Condition (III). For example for C = 0.08, Q = 1.65 the Sampling Lemma holds for all in an interval with endpoints 0.15... and 0.98.... 80 7.2 Vertex Movement The intuition for moving vertices is taken from Morse theory, which considers structures that arise in sweeping out a smooth manifold [73]. The skin surface is a cross section at each moment during the sweep, and the manifold is the stack of cross sections in the time direction. In other words, the manifold is the graph of M : R3 R that maps a point x to the time t at which x belongs to the surface F (t). Hence F (t) = M-1(t). A metamorphosis of F corresponds to a critical point of M. Any two cross sections F (t0) and F (t1), where the interval [t0, t1] is free of critical points, are diffeomorphic, and an explicit diffeomorphism : F (t0) F (t1) can be constructed from the vector field grad(M(x)). Figure 7.1: Dotted integral lines of a solid growing circle and a solid growing hyperbola The step from time t0 to time t1 thus amounts to moving each vertex a Di along its integral curve to a = (a) Dj. In the growth model the integral lines are pieces of straight lines in spherical patches and hyperbolas in hyperbolic patches, as illustrated in Figure 7.1. To see this, note that (2x1, -2x3) are the normal vectors of the family of hyperbolas x21 - x23 = ±R2, and that (2x3, 2x1) are the normal vectors of the family 2x1x3 = ±R2 obtained by rotating the first family through an angle of /4. The three-dimensional picture is obtained by revolving the hyperbolas in Figure 7.1 about the x3-axis. The first family of hyperbolas turns into the one-parameter family of hyperboloids. The second family turns into a two-parameter family of hyperbolas each orthogonal to each of these hyperboloids. 81 7.2.1 Parametrization It is convenient to parametrize the integral curves by time so that points can be moved by evaluation. Each integral curve is decomposed by the mixed complex into segments of lines and hyperbolas. We first consider the case of a line segment inside a mixed cell constructed from a Delaunay vertex and its dual Voronoi polyhedron. After translating the center of the cell to the origin, the mixed cell is swept out by a sphere in standard form x21 + x22 + x23 = s, for s 0. We thus get integral rays that start at the origin and go to infinity, and we clip each such ray to the mixed cell. If the origin lies inside the mixed cell, then it is the source of an entire sphere of integral lines. We parametrize that sphere by longitude and latitude, [0, 2 ] and [- , ]. For each pair of angles we have a half-line , : R R3 defined by cos cos s , (s) = sin cos s . sin s The case of a mixed cell constructed from a Delaunay tetrahedron and its dual Voronoi vertex is symmetric, with integral lines ending rather than starting at the origin. If the origin lies inside the mixed cell, then it is the sink of an entire sphere of integral lines. We next consider the case of a mixed cell constructed from a Delaunay edge and its dual Voronoi polygon. We assume the hyperboloid sweeping out the mixed cell is in standard form x21 + x22 - x23 = s, for s R. The integral lines are hyperbolas with the x3-axis as one asymptote and a line through the origin in the x1x2-plane as the other. We parameterize the family with the longitudinal angle, [0, 2 ], and the minimum distance to the origin, 82 R 0. For each pair of parameters we get a hyperbola ,R : R R3 defined by ± cos u ± sin u ,R(s) = , R2 ± 2 u with u = (s + s2 + r2)/2. To check the correctness of the parametrization note that the points ,R(s) satisfy the equation of the hyperboloid and the equations of the orthogonal hyperbola. The case of a mixed cell constructed from Delaunay triangle and its dual Voronoi edge is symmetric, with integral curves oriented in the opposite direction. In either case we obtain a parametrization in time by setting s = R2 + t. Note that in all four cases of integral lines, the speed of the parametrization depends only on the distance to the center of the mixed cell: / s = 1/(2 (s) ). This is consistent with the length of the gradient of M(x) = ±x21 ± x22 ± x23 being independent of the choice of signs: grad M(x) = 2 x . 7.2.2 Crossing Mixed Cells We described the movement of vertices within specific mixed cells. The vertices within the mixed cells will hit the boundaries and move across to other mixed cells. For vertices on the sphere patches, it is trivial to intersect the vertex trajectories, which are straight lines, with the boundary planes of the mixed cells. Whenever a vertex hit a boundary plane, the vertex move to the other mixed cell which is divided by this plane. For the case of hyperboloid, we can find the intersection of the trajectory and the bound- ary plane in the same way. However, this will involve solving the equation of the trajectory. We avoid this by first finding the plane which the trajectory lies on. Then we find the intersection of the trajectory plane with the boundary planes of the mixed cells which are parallel to the rotational axis. The vertex will hit the boundary plane which has its inter- 83 section closest to the rotational axis. Currently, we use brute force to search the boundary plane which is hit by a moving point. If | B |= n, the worst computational time for searching can be O(n) for each movement. This can be improved by point-location algorithm within a mixed cell to improve the running time. 7.3 Event Scheduling Vertices move continuously along their integral lines, but updating their positions continu- ously is computationally infeasible. The common escape from this dilemma is the time-slicing method, which takes discrete time steps and advances all vertices from time t0 to time t1 without intermediate stop. There are drawbacks to time-slicing related to the difficulty of choosing the right step size. We follow an alternative approach and take different time steps at different locations. This is done by prioritizing the four types of operations that occur at discrete moments in time, which are edge flips, edge contractions, vertex insertions, and metamorphoses. Edge flipping, edge contractions and vertex insertions are discussed in pre- vious chapters and the time for checking of these operations in Section 7.5. We will introduce the detail of topology changes in the metamorphoses section. In the KDS paradigm, we should store the trajectory of the points. However, we store the coordinates of each point in this work to represent the trajectory. Coordinate updates are done lazily, moving a vertex when and only when it is used by one of the other four operations. This result in a time-warped surface with different pieces reflecting the state at different times. To bring the entire surface to the present time, we simply update all the vertex coordinates, and by assumed correctness of the prioritization this requires no other changes in the triangulation. We use a priority queue to schedule the operations in time. We denote this queue by q. Starting with nothing in space, we first throw all the metamorphoses operations into q 84 and we extract the operation according to time in q. Then for any triangle in the current triangulation, we estimate a time when a checks for refinement is required. We also push this operation onto q. So the algorithm of constructing the dynamic skin is the processing of q until it is empty or the desired time t is reached. 7.4 Topology Adaptation The way the skin surface is connected can change during deformation. This section studies when, where, and how these changes happen in the growth model. It also describes how we locally modify the general sampling strategy to avoid the computational impossibility of sampling infinitely many points accumulating at the location of infinite curvature. 7.4.1 Growth Model We recall that the growth model of deformation is defined by changing the square radius of each sphere (zi, ri) from wi = r2i at time 0 to wi + t at time t R. Computationally, this is the simplest kind of deformation of the skins because it keeps the mixed complex invariant. Each mixed cell contains a (possibly empty) sphere or hyperboloid patch of the skin surface. After normalization, the equation of each sphere or hyperboloid at time t is t x21 + x22 ± x23 = ±R2 + . 2 A metamorphosis happens when the right-hand side vanishes at time t = 2R2, and it happens at the center, but only if the center lies in the interior of its mixed cell. If the center lies outside, the portion of the sphere or hyperboloid that passes through the center is not part of and thus does not affect the skin surface. Using local considerations, we can reduce the list of metamorphoses to the four given in table 7.1. Cases k = 0, 3 correspond to an appearing or disappearing sphere. Cases k = 1, 2 85 k type of change and inverse 1 creating/annihilating a component 2 adding/removing a handle 3 closing/opening a tunnel 4 filling/starting a void Table 7.1: The four types of generic topological changes that happened during growth/shrinking. correspond to switching a hyperboloid from two sheets to one, or vice versa. In each case we can interpret the center as a critical point of the map M : R3 R whose level sets M-1(t) are the skin surfaces at time t. Cases k = 0, 3 correspond to minima and maxima, and Cases k = 1, 2 to two types of saddle points. The gradient of M vanishes at all these points and also at centers that lie on the boundary of their mixed cells. The latter centers correspond to degenerate critical points in the sense that an arbitrarily small perturbation of B suffices to turn them into regular points [45]. 7.4.2 Hot Spots Common to every metamorphosis is a local drop in length scale, which reaches zero at the moment and point of the metamorphosis. We analyze the situation in some detail. Let H be a positive real number. The hot portion of the skin surface F is the set of points with length scale H or smaller, FH = {x F | (x) H}. By the Iso-curvature Lemma, we have (x) H only if x is sufficiently close to the center of a sphere or hyperboloid. Let zX be such a center. We call the ball X = {y R3 | y - zX H} 86 the hot ball of X. We pick H to be at most half of the shortest distance between all the centers to their nearest boundaries of the mixed cells and denote the balls as R3H to be the hot portion of space. Hot Spot Lemma. FH = F R3H The Hot Spot Lemma follows directly from the Iso-curvature Lemma and does not need a separate proof. In words, a point x F belongs to the hot portion of the skin surface if and only if it belongs to the hot portion of space. In the growth model the hot portion of space is constant, while the hot portion of the skin changes as the surface moves through that portion of space. The radius of the hot ball in any mixed cell is half of the distance from the center to the nearest mixed cell boundary. This guarantees pairwise disjointness. We refer to mixed cells with distance zero as degenerate cases. Our algorithm does not handle these cases because we are not sure if there is topology changes. For example, four coplanar points lies on a circle have two triangles with centered on one of the diagonals of the convex hull. Topology occurs at the centers with no different from the case of a non-degenerate triangle. Possible solution may be clustering simplices in DB into complexes, for example, merge the two triangles into a quadrilateral with one center by more complicated data structure other than simplicial complex. However in reality, there are cases which are very close to degenerate cases, e.g. a regular pentagon. Because of numerical imperfection, a pentagon cannot be represented perfectly by floating points. This causes mixed cell with very small distance between center and its cell boundary. In cases like this, our algorithm works but with poor performance. The reason is because there are huge numbers of operations near the centers. Improving the computational time for these cases is also an open problem. 87 7.4.3 Time for Change The hot portion is more difficult to triangulate than the rest of the skin surface. One reason is metamorphosis; another is accumulation of vertices in a small region. The sphere case is relatively harmless, because the area decreases at the same rate as the density requirement increases. Indeed, a constant number of vertices suffices to shrink a sphere to an arbitrarily small size. The case of a hyperboloid that approaches its limiting double-cone is more problematic, because the number of vertices near the center grows beyond any bound. To circumvent the computational impossibility of sampling infinitely many points, we change the sampling strategy inside the hot balls. We give up on -sampling to get a sparse sampling, but we preserve the closed ball property. The triangulation algorithm remains oblivious to the changed sampling density and maintains a restricted Delaunay triangulation. 2R 2R 0 1 2R H 1 2R 0 hH hot spheres Figure 7.2: Topology change: Head on view of start, middle, end configuration generated by special sampling strategy taking a 2-sheeted to a 1-sheeted hyperboloid. Consider a two-sheeted hyperboloid and translate time such that the metamorphosis happens at time t = 0. The hyperboloid enters its hot ball at time -2H2, turns into a double- cone at time 0, and leaves the hot ball as a one-sheeted hyperboloid at time 2H2. The special sampling strategy that allows us to go through this motion depends on a parameter 0 < h < 1. Special sampling begins at time t0 = -2H2h2 when the two-sheeted hyperboloid enters the ball of radius Hh, and it ends at time t1 = 2H2h2 when the one-sheeted hyperboloid leaves that ball, as shown in Figure 7.2. At time t0, the hyperboloid intersects the boundary 88 of the hot ball in two hot circles. The shape adaptation algorithm moves these circles along their integral lines, which implies that they grow from radius R0 = H (1 - h2)/2 at time t0 to radius R1 = H (1 + h2)/2 at time t1. Simultaneously, the distance between the two circles decreases from 2R1 to 2R0. We define the hot sphere to pass through the two hot circles. At time t0 and t1, it is the boundary of the hot ball which radius is H(1 - h4)1 R 4 0R1 = . 2 In the open time interval between t0 and t1, it is cocentric and smaller than that boundary. General sampling applies outside the hot sphere and special sampling applies on and inside that sphere. 7.4.4 Special Cap Construction At time t0, we kick off special sampling by creating the double-cap as the start configuration of the metamorphosis representing the intruding portion of the two-sheeted hyperboloid, as shown in the left drawing of Figure 7.2. Consider one sheet of the hyperboloid and let a be its intersection point with the symmetry axis. Let b0, b1, ..., bl-1 be the vertices of a regular l-gon along the hot circle in this sheet. We mirror these points across the symmetry plane of the hyperboloid and get points a , b 0, b 1, ... on the other sheet. First we add the tip of the pyramid a and the base points bi to the triangulation. During the point insertion, flipping is applied also as in the normal insertion operations. The only difference is that these insertions may produce short edges. The second step is to delete any point within the hot sphere using edge contraction. Finally, we check if more edge contraction is needed. After the edge contractions, every point connected to the tip is either the base of the pyramid or a point outside the hot sphere. We flip any edge with one endpoint outside the hot sphere. The cap is formed at this point. However, because of the final flipping, there 89 may be triangles near the base of the cap that are too large. One more final operation is to refine the triangles involved. This may cause special flipping, which will be discussed in a Section 7.4.6. Section 7.4.7 will derive sufficient conditions for h and l that guarantee the above algo- rithm successfully constructs the double-cap as the start configuration of the metamorphosis. By this we mean that (1) a,a are the only vertices inside and bi, b i are the only vertices on the hot sphere, (2) the link of a in D is the regular l-gon of vertices bi, and symmetrically the link of a is the l-gon of vertices b i. Assuming C = 0.08 and Q = 1.65, we will see that h = 0.98 and l = 5 are feasible values for the two constants. For ease of reference we say the vertices and edges in the links of a, a are hot and the vertices, edges, and triangles in the stars of a, and a are very hot. 7.4.5 Special Waist Construction The end configuration of the metamorphosis is identical to the start configuration of the inverse metamorphosis. As shown in the right drawing of Figure 7.2, it consists of two rings of triangles forming a spool representing the intruding portion of the one-sheeted hyperboloid. Let u0, u1, ..., um-1 be the vertices of a regular m-gon along the waist where the hyperboloid intersects its symmetry plane. Similarly, let w0, w1, ..., wm-1 be the vertices of another regular m-gon along one of the two hot circles, rotated by /m relative to the m-gon along the waist. Finally, let w 0, w 1, ..., w m-1 be the vertices of the mirror m-gon on the other hot circle. We start by adding each ui to the triangulation, and then each wi and w i. Then we delete all other vertices within the hot sphere by edge contraction. After the contractions we will have ui connected to wi, w i, another uj, or a point out of the hot ball. We flip any edge leaving the hot ball. we flip this edge. This will again create some oversized triangles. Finally, we refine these triangles, which involves special flipping. 90 Section 7.4.7 derives sufficient conditions for h and m that guarantee the above algorithm successfully constructs the spool as the start configuration of the inverse metamorphosis. Assuming C = 0.08 and Q = 1.65, we will see that h = 0.98 and m = 40 are feasible values for the two constants. For ease of reference we again say that the vertices and edges along the two hot circles are hot and that the vertices, edges, and triangles between the two hot circles are very hot. In the forward direction we switch from the double-cap to the spool at time 0, and in the backward direction we do it the other way round. The 1 to 2 transformation is easier because we just need to meld the m-gon of the waist into a single vertex and then split that vertex into two. In the forward direction we first meld a and a into a single vertex and then expand that vertex into a regular polygon (interleaving angularly between the bi). The expansion creates two new rings of triangles between the new polygon and the polygons representing the two hot circles. This is done following the angular order of the vertices around the symmetry axis. 7.4.6 Special Sampling The main difference between special and general sampling is that the former gives up on the Lower Size Bound [L] for hot edges and on the Upper Size Bound [U] for very hot triangles. The length of hot edges is bounded from above because the Upper Size Bound applies to the incident triangles outside the hot sphere. A more detailed analysis of edge and triangle sizes including a proof of the closed ball property (in spite of special sampling) is given in Section 7.4.7. The goal of special sampling is to maintain the double-cap and the spool during the first and the second halves of the time interval. It acts primarily by modifying general sampling for points on and inside the hot sphere. As a general rule, an edge is contracted by removing an endpoint that is not hot. Cases where both endpoints are hot occur only at the end of the metamorphosis (or its inverse) and will be discussed separately. There are two ways in which 91 general sampling can intrude into the hot sphere: by adding a point inside that sphere and by flipping a hot edge. In both cases we prevent the intrusion by bisecting the endangered hot edge bc, as illustrated in Figure 7.3. Specifically, we add the midpoint q of the shorter hot circle arc that connects b with c. The addition of q may create edges that violate the Lower Size Bound. Of these we contract the ones that are not hot, always making sure we remove the endpoint that is not hot. As discussed above, we choose H small enough so that hot spheres cannot get too close to each other and every non-hot edge has at least one non-hot vertex. Infinite loops cannot occur because each iteration leaves an addition hot vertex behind. The hot circle gets denser and intrusions into the hot sphere get progressively more difficult. d b q c p a Figure 7.3: Special flipping: The hot edge bc is bisected either because the dual restricted Voronoi vertex of bcd lies inside the hot sphere or edge flipping attempts to change bc to da. Special sampling maintains the special configuration, but it does not guarantee the two Size Bounds. They must therefore be enforced by clean-up operation algorithmically at the end of the metamorphosis. The clean-up operation is correct if we maintain the closed ball property, which is initially guaranteed by special sampling. While maintaining that property might be difficult in general, we can use the insights gained from the proofs of the two Persistence Lemmas in Section 7.4.7 and add points only on the two hot circles. 92 7.4.7 Analysis This section analyzes the point configurations generated by fipling. Recall that Hh < H is the length scale threshold that triggers the start and end of special sampling. In the forward direction we start with a two-sheeted hyperboloid that enters the ball with radius Hh around the center, and we end with a one-sheeted hyperboloid that exits the same ball. In the backwards direction the events are the same in reverse order. Sizes at Transition. Refer to the double-cap shown in Figure 7.2. The l+1 points on one sheet form a regular l-sided cap. The l vertices of the base lie on the hot circle with radius R0 = H (1 - h2)/2, which lies in a plane at distance R1 = H (1 + h2)/2 from the center. Note that R20 + R21 = H2. Define b = bi and c = bi+1, with indices modulo l. Independent of the index i, the lengths of the edges of abc are ab = 2Rab = 2R1(R1 - Hh), bc = 2Rbc = 2R0 sin . l Any isosceles triangle with sides of length E and height L has circumradius E2/2L. The height of abc is La,bc = 4R2ab - R2ab. The circumradius is therefore 4R2ab/2La,bc, which is R R 1(R1 - Hh) abc = . 2R1(R1 - Hh) - R20 sin2( /l) Next refer to the spool shown in Figure 7.2. The 3m points form three parallel regular m-gons. The distance between two contiguous planes is R0, and the circumradii of the three m-gons are R1, Hh, R1. Define u = ui, v = ui+1, w = wi, x = wi+1, with indices modulo m. 93 Independent of the index i, the lengths of the edges uv and wx are uv = 2Ruv = 2Hh sin , m wx = 2Rwx = 2R1 sin . m To compute Rvw, Ruvw, and Rvwx, we consider the projection of the middle and outer m-gons onto a plane parallel to the two m-gons, as shown in Figure 7.4. The distance between the projections of w and (u + v)/2 is R1 - Hh cos( /m), and that between the projections of v and (w + x)/2 is R1 cos( /m) - Hh. We get the heights Lw,uv and Lv,wx of the two triangles by taking the distances to three dimensions, which means squaring, adding R20, and taking square roots. The length of an edge connecting the middle m-gon with one of the two outer m-gons is the root of R2uv + L2w,uv, which is 2Rvw = 2R1(R1 - Hh cos( /m)). We compute the circumradii of the two isosceles triangles again from their edges and heights. x P v w h u Figure 7.4: Top view of a waist: Portion of the spool in Figure 7.2 projected onto a plane parallel to the 2k-gons. 94 In particular, the circumradius of uvw is 4R2/2Lw,uv, and that of vwx is 4R2vw/2Lv,wx. Hence, R R 1(R1 - Hh cos( /m)) uvw = , H2 - 2R1Hh cos( /m) + H2h2 cos2( /m) R R 1(R1 - Hh cos( /m)) uvw = . R21 - 2R1Hh cos( /m) + R21h2 cos2( /m) Smooth Transition. We derive necessary and sufficient conditions for h, l, and m that guarantee a smooth transition from the general to the special sampling strategy. By this we mean that the configurations at the start of a metamorphosis is an -sampling and satisfies both Size Bounds. At the end of the metamorphosis, the Size Bounds are enforced by eliminating offending edges and triangles through edge contraction and vertex insertion. The result is a triangulation whose vertex set is an -sampling of the surface; see also the remark immediately following the proof of the Sampling Lemma. The length scale at vertices a,u,v is Hh, and that at b,c,w, and x is H. The Lower and Upper Size Bounds are therefore equivalent to Rab, Rbc, Ruv/h, Rwx, Rvw > (C/Q)H and Rabc, Ruvw, Rvwx < CQHh respectively. The inequalities for Rvw, Ruv, and Ruvw are redundant because Rab < Rvw, Rwx < Ruv/h, and Ruvw < Rvwx for all h < 1. In addition to requiring that the triangles abc, uvw, and vwx satisfy the Upper Size Bound, it is convenient to require that their radii are less than the locally allowed minimum edge length. This extra requirement implies that after adding points on and inside the hot sphere, all old points inside or on the hot sphere are too close to at least one new point and thus are deleted. It follows that all remaining old vertices lie outside the hot sphere. We thus have the following two conditions: Rab, Rbc, Rwx > HC/Q, (IV) Rabc, Rvwx < min{Q, 2/Q}CHh. (V) Conditions (I)-(V), summarized in Table 7.2, are satisfied for = 0.279, C = 0.08, 95 Q = 1.65, h = 0.98, l = 5 and m = 40. We summarize the results assuming this assignment of constants. (I) 0 0 = 0.297... (II) Q2 - 4CQ - 2 > 0 (III) 2/(1 + )2 - 4/4 > C2Q2 (IV) Rab, Rbc, Rwx > HC/Q (V) Rabc, Rvwx < min{Q, 2/Q}CHh Table 7.2: Conditions (I)-(V). Transition Lemma. The triangulation at the start of a metamorphosis satisfies the two Size Bounds and its vertex set is an -sampling of the skin portion within a mixed cell. As mentioned earlier, the same does not automatically hold for the end configurations of metamorphoses, but it can be enforced algorithmically by refinement. The purpose of bounding the size of triangles in Condition (V) by (2/Q)Ch is to guarantee that the algo- rithm given in Section 7.4 constructs the special configurations without having to search for remaining old vertices inside the hot sphere. To prove this algorithm correct, we also need to show that these configurations are part of the restricted Delaunay triangulation, which follows from the Persistence Lemmas proved below. Persistence of Triangulation. We show that the special configuration exist as subcom- plexes of the restricted Delaunay triangulation during the entire time interval of a meta- morphosis. Consider the simplices in the Delaunay complex spanned by hot vertices and their dual Voronoi polyhedra. Figure 7.5 sketches both for the double-cap before and the spool after the transition. During the first half of the time interval, the hot vertices span two pyramids, each the reflection of the other across the symmetry plane of the hyperboloid. The points are in degenerate position, which implies that the Delaunay complex 1 of the hot points contains polyhedra that are more complicated than tetrahedral. Specifically, 1 consists of the two-pyramids joined by an edge connecting their apices and a ring of prisms 96 Figure 7.5: Solid Delaunay complex and dotted Voronoi polyhedra of hot ball configuration before and after the double cone. and quadrangles around that edge. As usual, D denote the restricted Delaunay triangulation of the entire vertex set V . Persistence Lemma A. At any time in [-2H2h2, 0], the intersection of 1 and D consists of two rings of triangles forming the double-cap and of their edges and vertices. Proof. We first show that the edges, polygons, and polyhedra in 1 that do not belong to the double-cap also do not belong to D. The edges connecting the two caps have dual Voronoi polygons which lie in the symmetry plane separating the two sheets and therefore cannot intersect the hyperboloid. To see that they do not intersect any other part of F , we consider the sandwiching spheres defined for points of F inside the hot sphere. The Voronoi polygons are contained in the union of balls bounded by these spheres, else they would imply an empty sphere that intersects the hyperboloid in a patch outside the hot sphere that is large enough to contradict the -sampling property. Detailed computations of a lower bound 97 for the size of such an implied patch are omitted. Since D is a complex, it also does not contain the Delaunay polygons and polyhedra incident to the excluded edges. The base polygons of the two pyramids in 1 have their dual Voronoi edges on the symmetry line of the hyperboloid. For the same reason as above, these edges are contained in the union of balls bounded by the sandwiching spheres of points of F inside the hot sphere. We second show that the triangles abc of the double-cap belong to D. At time t0 = -2H2h2 this is true because these triangles have circumspheres that are small enough that every point of F inside these spheres would belong to edges that violate the Lower Size Bound. At any time t0 < t < 0 this is true because any violation is prevented by the algorithm before it occurs. During the second half of the time interval, the hot vertices from three convex polygons in three parallel planes. The middle polygon is a regular m-gon in the symmetry plane of the hyperboloid, and the other two are reflections of each other across that plane and are inscribed in the two hot circles. The Delaunay complex 2 of the hot points is again degenerate, consisting of the above mentioned three polygons, which form the top and bottom facets of two drum-like polyhedra. The two drums are surrounded by a ring of four-sided pyramids alternating with tetrahedra. Persistence Lemma B. At any time in (0, 2H2h2), the intersection of 2 and D consists of two rings of triangles forming the spool and their edges and vertices. Proof. The edges, polygons, and polyhedra in 2 that do not belong to the spool have their dual Voronoi polygons, edges and vertices either in the symmetry plane or the symme- try axis of the hyperboloid. For the reason mentioned in the proof of the Persistence Lemma A, these polygons, edges, and vertices are contained in the union of balls bounded by sand- wiching spheres of points of F inside the hot sphere. The corresponding edges, polygons, and polyhedra of 2 thus do not belong to D. 98 The reminder of the proof establishes that the triangle uvw and vwx belong to D. Im- mediately after time t = 0 this is true because the triangles in the double-cap belonged to D immediately before that time t = 0. At any time 0 < t < t1 this is true because any violation is prevented by the algorithm before it occurs. The two Persistence Lemmas also hold for the reverse metamorphosis, which changes a one-sheeted into a two-sheeted hyperboloid. To see this, run time backwards and exchange the arguments that establish that the two special configurations are subcomplexes of D when they are first created. These arguments are contained in the respective last paragraphs of the two proofs. Summary. The two Persistence Lemmas establish that the closed ball property of the restricted Voronoi diagram is maintained even inside the hot spheres that guide the algorithm through the various metamorphoses. Special Homeomorphism Theorem. The restricted Delaunay triangulation of the points chosen by special sampling triangulates the skin surface inside each hot sphere. Together with the General Homeomorphism Theorem this implies that we have a trian- gulation of the skin surface at all times. 7.5 Scheduling When the skin is growing, the triangles can get distorted. From time to time, checking on each triangle is necessary. We have already discussed how to check and refine the mesh when it is needed. The remaining problem is how to determine when to check. In this section, we will introduce relaxed scheduling as a method to avoid the computationally expensive root computation. We complete a prove for edge contraction scheduling, but the checking time in point insertion scheduling is estimated by experiments. 99 7.5.1 Method The key technical insight about the dynamic skin triangulation algorithm is that we can find constants C, h, l, and m and interval [Q0, Q1] such that Conditions (I) through (V) are satisfied for all Q [Q0, Q1]. Instead of fixing Q and contracting an edge uv when its size- scale ratio reaches C/Q, we suggest to contract uv at some time while its size-scale ratio is in the interval ( C , C ]. After the ratio enters this interval at the upper end it may leave Q1 Q0 again at the upper end or the edge may get contracted, but the ratio is not allowed to leave the interval at the lower end, since that would violate the Lower Size Bound for Q = Q1. The two possible scenarios are illustrated in Figure 7.6. ratio C Q0 contraction C Q1 time Figure 7.6: The curve shows the development of the size-scale ratio of an edge in time. The edge is contracted before the ratio reaches the lower bound. Vertex insertions are treated symmetrically, that is, a triangle is removed by adding its dual point on F at some moment when the size-scale ratio at the triangle is in the interval [CQ0, CQ1). The ratio can enter and leave the interval at the lower end but not the upper end. The two possible scenarios are illustrated in Figure 7.7. We call ( C , C ) and [CQ Q 0, CQ1) the lower and upper size buffers. The quality of the 1 Q0 triangulation is guaranteed because all edges satisfy the Size Bounds for Q = Q1. The correctness of the algorithm is guaranteed because edge contractions and vertex insertions 100 ratio CQ1 insertion CQ0 time Figure 7.7: The curve shows the development of the size-scale ratio of a triangle in time. The triangle is removed before its ratio can rise above the upper bound. are executed only if the Size Bounds for Q = Q0 are violated. 7.5.2 Early Warnings We say an edge matures at the time its ratio reaches C , and similarly a triangle matures Q1 at time its ratio reaches CQ1. The relaxed scheduling method depends on an early warning algorithm that reports an edge or triangle before it matures. The algorithm may err and produce false alarms, but it must not miss a maturation. False alarms cost time but do not cause any harm. Suppose an edge uv is reported at time t0. It is contracted if its ratio lies inside the lower size buffer, and it is not contracted if it lies outside that buffer, which can only be because Ruv > C . Note that in the second case, the ratio has to cross the entire uv Q0 lower size buffer before it can reach the lower end, C . The early warning algorithm is based Q1 on this observation and issues the next warning before the ratio had enough time to cross that buffer. Triangles are treated similarly. 101 7.5.3 Analysis In this section we derive lower bounds on the amount of time it takes for an edge to change its size by more than some threshold value. From these upper bounds we will derive upper bounds on the time it takes an edge to pass through the entire size buffer. norm time Figure 7.8: Graph of the dependence between time and norm. At any moment, the speed is minus the slope of the tangent line. Traveling Point. We recall that the speed of a point u F is 1 , assuming the patch 2 x that contains it is in standard form. The distance traveled by u in a small time interval is maximized if we move u straight towards the origin, which for example happens if u lies on a shrinking sphere. The amount of time it takes u to reach the origin is u 2. If we start the motion at time - u 2 we reach the origin at time 0 and get a dependence between time and distance from the origin as -t = u 2, whose graph is shown in Figure 7.8. The distance from the origin is therefore f(t) = u = -t and the speed is -df (t) = dt 1 = 1 , as required. To move from point u to point x = (1 - )u thus takes time 2 -t 2 u u 2 - x 2 = (2 - 2) u 2. By the Curvature Variation Lemma, the difference in length scale between two points u and x is at most their distance: (u) - (x) u - x . If that distance is u = (u) then 102 (x) (1 + ) (u). In other words, by traveling a distance (u) from u we can grow the length scale by at most a factor of 1 + . Shrinking Edge. Suppose uv is an edge that satisfies the Lower Size Bound for Q = Q0 at time t0, that is, Ruv > C Q uv. The edge cannot shrink faster than if both endpoints lie on 0 the same shrinking sphere. After traveling a distance (a) = (b) to points x and y, the fraction of the length is x - y R = xy 1 - . (7.1) u - v Ruv We get a lower bound on the size of xy in terms of its scale from the Lower Size Bound for uv and inequality 7.1: Rxy (1 - )Ruv C > (1 - ) Q uv 0 C = (1 - ) max{ (u), (v)} Q0 1 - C max{ (x), (y)} 1 + Q0 1 - C = 1 + Q xy. 0 If we set 1- = Q0 we get R 1+ Q xy > C xy, or in words, the edge xy is guaranteed to satisfy 1 Q1 the Lower Size Bound for Q = Q1. In other words, if we let the points travel a distance times the length scale, with = Q1-Q0 , then the edge still satisfies the relaxed version of the Q1+Q0 Lower Size Bound. The corresponding time interval in which we can be sure that the edge uv does not become unacceptably short is Q Q t = 1 - Q0 2 - 1 - Q0 · 2 Q uv. 1 + Q0 Q1 + Q0 103 For example, for Q0 = 1.65 and Q1 = 2.3, we get = 0.165 and t = 0.302 · 2uv. Point Insertion. The case of triangles is more complicated and no theoretical result for t is proved yet. The quality of a triangle is determined by triangle center. Solving the trajectory of the center is relatively more complicated then the edge contraction predication. In the other hand, we can estimate t for point insertion experimentally. Using different values of t, we observe that for a triangle abc, t = 0.15 · 2abc works for the above early warning theory. This means, if we use this value, we expect that no triangle will have size larger than CQ abc at any time. 7.6 Conclusion We develop the dynamic skin triangulation in this chapter. The algorithm maintains a quality mesh during growth of the skin surface. We would like to generalize our algorithm to more general models of deformations, where the skin can deform not only by the growth model but weighted points in B can both move and change their radius arbitrarily. We may not be able to determine the integral lines explicitly as in Section 7.2, but moving vertices along integral lines is convenient but not necessary for the algorithm fortunately. An approximation of that movement may suffice. For small time steps, the triangulation changes only a small amount and may be maintained with the methods described in this chapter. This is still an open problem and it will be interesting to see the deformation of the mixed cells and the patches in them. The solution may involve techniques related to the alpha shape morphing [25, 26, 40] and kinetic geometric data structures such as the works of Basch [9]. 104 Chapter 8 Measurements of Skin Curves 8.1 Introduction In the field of mechanical engineering, the shape and topology of a structure determine its strength and cost. The goal of structural optimization is to distribute material to minimize a cost function under one or more constraints [10]. During the optimization, the boundary of a structure is modified and new components or holes are also introduced. The problem with topology changes is called generalized or variable-topology shape optimization [81]. We assist the optimization problem by providing deformable shapes with exact and smooth boundaries. In this chapter, we will concentrate on two-dimensional designs. A structure is represented by grids in present practice [53, 81], as illustrated in Figure 8.1. Each grid element is assigned an area fraction-like design parameter ranging from 0, signifying a void element, to 1, signifying a solid element. Finite element analysis is performed on these grids to optimize a structure in structural optimization, for example, computing a lightest and strongest floor support inside an airplane. The finite element analysis modifies the shape locally according to the measurements of the portion within a small disk. As shown in Figure 8.1, a structure (upper) is optimized into a better one with a shorter perimeter (lower) by finite element analysis [53]. We will not go into the details of analysis, and we refer the reader to the works of Bendsøe, Haber, and Jog [10, 53]. Rather, we study the deformable shape representation that give the geometric and differential information for the analysis. 105 Figure 8.1: Cross sections of two beam designs represented by grids [53]. The grid representation does not work well when the exact location of the boundary is an issue, such as in the case of optimizing a structure under a load or pressures because the boundary of the structure is too fuzzy. A deformable shape with exact boundary is needed to solve this kind of problems. In order to provide a shape representation with an exact boundary and deformation, we propose representing the structure by using union of disks and bodies bounded by skin curves [26]. General level sets have difficulties such as introducing topology changes into the shapes. Variation in topology usually requires additional representation of shapes [53]. The skin and union of disks, on the other hand, provide the following information for the algorithm of finite element analysis and enable the process to deform and change topology without complications. Given a window, which is a disk in R2, with center x, we compute: * The measurements, namely, area, perimeter, and their derivative of the shapes within the window, and 106 * The nearest point p on the boundary of the shapes to x, the distance xp , and the tangent line passing through p. In this chapter, we will focus on how to compute the area, the perimeter, as well as their derivatives of the entire union of disks or skin curve in R2. Note that these qualities are also useful in molecular computations [78, 87]. The computation of these qualities clipped by a window and the nearest point p are not done yet in this work. We assume the centers of the skin and union of disks move independently. We denote the velocity of the center zi by z i and the rate of changes of the radius ri by r i = (2ri)-1dr2i/dt = (r2i) /2ri. The finite element analysis gives B = {b i = (z i, r i)} for the current configuration of B, and we can compute the area A and the perimeter P of the shapes of B and the skin of B as well as their derivatives, A and P as a feedback to the analysis. The variable used in our analysis are listed at the end of this chapter. 8.2 Union of Disks zi Hij zj Figure 8.2: Distance from a center to a Voronoi edge. To calculate the area, perimeter, and their derivatives of a union of disks, we will break down the problem locally. Each disk bi B is clipped within its Voronoi cell. We can locally 107 compute the measures within each cell, the intersection of finitely many half-planes. We first study the behavior of the case of two disks, then we generalize to arbitrary unions of disks. y x Figure 8.3: Perimeter and area of a clipped disk : For a disk in standard form, x2 + y2 = 2, we integrate its length and area from - to . We denote the distance from zi to the bisector between bi and bj as Hij, as in Figure 8.2. 8.2.1 Area and Perimeter of Two Disks We first give area and perimeter formulas for a clipped disk: {(x, y) R2 | x2 + y2 2, and y }. The perimeter P and area A are P = 2 ( - ), (8.1) A = 2( - ) + , (8.2) where = arccos and = 2 - 2 as in Figure 8.3. Given two disks bi = (zi, r2i) and bj = (zj, r2j), the center cij of the smallest disk orthogonal 108 to both is cij = zi + (1 - )zj, where = 1 - r2i-r2j and l 2 2l2 ij = zi - zj [42]. The distance between the disk center zi and ij the bisector between bi and bj is Hij = zi - cij = (1 - ) zi - zj = (1 - )lij l r2 = ij + i - r2j . (8.3) 2 2lij We put = ri and = Hij in Equations (8.1) and (8.2) and get the perimeter Pij and area Aij of the union of disks as follows: H P ij ij = 2ri - arccos , ri H A ij ij = r2i - arccos + H r2 r ij i - H2ij. i 8.2.2 Area and Perimeter of Union of Disks Before we give the equations of the measurements, we first define some geometric entities. Define bji = {x R2 | j(x) i(x) 0}, which is the portion of bi on bj's side of the bisector. Define also Aji = area(bji) and Ajki = area(bji bki). Similarly, let Pji and Pjk i be the perimeters of the arcs of bji and bji bki respectively. The perimeter and area of the union of 109 disks is PB = 2 ri - P ji - P jk i , i j j,k AB = r2i - Aji - Ajk i . i j j,k These inclusion-exclusion formulas were proved by Edelsbrunner [38]. 8.2.3 Notations and Decomposition of Derivatives Computation Given B = {(zi, r2i)} and B = {(z i, (r2i) )}, our goal is to compute the derivatives of perime- ter and area due to the disks' movements. The change of each radius affects the size of a disk inside its Voronoi cell as well as its edge. Note that each edge of the cell moves in the direction parallel to the corresponding Delaunay edge when only the radius changes. The other cause of change is the disk center movement z i, which can be decomposed into two components, one parallel to the Delaunay edge zizj and one perpendicular. The parallel component is z (z ¯z i j - zi) ij = z . z i, i lij The normal component is then z i - ¯z ij. The parallel component ¯z ij causes a slope-preserving motion of the bisector of bi and bj. The other component causes a distance-preserving motion that rotates the edge of the cell around the other center zj. For derivative computation of the slope-preserving motion, we decompose the calculation into several components, in two classes: 1. Voronoi cell boundary motion: Hij/ zi, Hij/ zj, Hij/ rj and Hij/ ri. 2. Area and perimeter changes due to boundary motion and radius changes: Aij/ Hij, Pij/ Hij, Aij/ ri and Pij/ ri. 110 For slope-preserving motion, we can compute Hij/ lijl ij instead of Hij/ ziz i+ Hij/ zjz j, where l ij = (- z i, zj - zi - z j, zi - zj )/lij = z i - z j, zi - zj /lij. In the following two subsections, we will discuss slope and distance preserving motions. 8.2.4 Slope-preserving Motion Motion of Boundaries. We compute the derivative of Hij directly by differentiating it with respect to ri and lij using Equation (8.3): H H H H ij ij ij ij = l r r l ij + i + j, (8.4) ij ri rj where Hij 1 r2 = - i - r2j , lij 2 2l2ij Hij r = i , and ri lij Hij r = - j . rj lij Hence, we have 1 r2 r r H i - r2j ir i jr j ij = - l - 2 2l2 ij + ij lij lij Derivatives of Disks. We first give the generic derivatives of area and perimeter of a disk clipped by a half-plane as in Figure 8.3. The perimeter and area are expressed in Equations (8.1) and (8.2), and = 2 - 2. The derivatives of the area and the perimeter with 111 respect to and are P = 2 - - , (8.5) P 2 = , (8.6) A = L, (8.7) A = 2 . (8.8) We did not consider as an independent variable because = arccos . Then by putting = ri, = Hij, = ij = arccos Hij/ri and = dij, we get Pij H = 2 - ij , (8.9) r ij - i dij Pij 2r = i , (8.10) Hij dij Aij = P r ij, (8.11) i Aij = 2d H ij, (8.12) ij where dij = r2i - H2ij. 8.2.5 Distance-preserving Motion The components ¯z ij and ¯z ji change the cell wall distances Hij and Hji in the slope-preserving motion. The other components cause rotational motions and change only the perimeter. Suppose we keep bj fixed and move bi only. Let t1 = ¯z ij and t2 be the perpendicular component of z ij in the counter-clockwise direction about zj. Define t2 = z i - ¯z ij if zi moves in the counter-clockwise direction and t2 = - z i - ¯z ij if zi moves in the clockwise direction, as in Figure 8.4. Denote the two endpoints of the arc of bj as x+ij and x-ij, so that x+ij is the endpoint on the right side of the Delaunay edge in the direction from zi to zj and 112 x+ij r z' j i t2 zj z t i 1 l ij x-ij Figure 8.4: Rotational motion: The rotational component of the movement of a vertex contributes to changes in the arc length. x-ij is the other endpoint on the left. Note that x+ij = x-ji and x-ij = x+ji. The perimeter of bj decreases at x+ij and increases at x-ij during counter-clockwise motion. In a two-disk configuration the effect can be ignored; the perimeter changes at the two endpoints cancel each other as in Figure 8.5. However, this motion will affect the length Pji if one of the two endpoints is covered by a third disk. The rate of perimeter change at each point is t2rj/lij. Define +ij to be 1 if x+ij exists and 0 otherwise. With the same notation for -ij, the rate of change of the perimeter due to this rotational motion is (- +ij + -ij)t2rj/lij. decreasing Pj t2 zj zi increasing Pj Figure 8.5: Distance-preserving motion does not contribute to length changes in a two-disk configuration. On the other hand, consider the upwards motion caused by t2 at zi as a relative motion of going down at zj. This affects the perimeter of bi and the rate of change is t2ri/lij. This 113 effect causes a change of ( +ij - -ij)t2ri/lij in the perimeter of bi. The overall rate of change in the perimeter in this rotational motion due to z i is, thus z Sign(z i - ¯ zij (ri - rj) i × (zj - zi)) · ( + ij - - ij ) , lij where z i × (zj - zi) is the cross product of vectors z i and zj - zi in R3, and we consider the `out-of-paper' direction to be positive and the `into-paper' direction to be negative. 8.2.6 Assembly of Equations Put the equation together now. For perimeter Pij, the derivative P ij of a clipped disk is: P P P ij ij ij = r H r i + ij [applying Equations (8.9) and (8.10)] i Hij H 2r = 2 - ij i ij - r H d i + ij [applying Equations (8.3) and (8.4)] ij dij 1 l r2 2r 1 r2 r r = 2 - ij i - r2j j i - r2j i j ij - + r - l r r . d i + ij + i - j ij 2 2lij dij 2 2l2ij lij lij The sum of derivatives of two neighboring clipped disks is: 1 (ri - rj)2 2( - ij) + - l d l ij ij ij P 1 (ri - rj)2 ij + P ji = r . i r j l ij 2( - - l ij) + d l ij ij ij r i + rj (ri - rj)2 1 - . dij l2ij By letting 1 (r i - rj)2 ij = l , d ij - ij lij 114 we simplify the sum of derivatives to 2( - ij) - ij P 2( - ij) - ij ij + P ji = r . i r j l ij r i + rj · l ij ij Perimeter Derivative of Union of Disks Theorem. The perimeter derivative of the union of disks, B, is l + P ij ij + - ij B = 2 ri i + ij(ri + rj) - 1 l 2 b ij i B i =j z +Sign(z i - ¯ zij (ri - rj) i × (zj - zi)) · ( + ij - - ij ) , lij with j P ji - j,k P jk i i = 1 - , 2 ri z ¯z i ij = z i, (zj - zi)/lij , and z i 1 (r i - rj)2 ij = l , d ij - ij lij where i is the fraction of perimeter inside the Voronoi cell of bi. The area derivative, A ij, of one clipped disk in a two disks configuration is A A A ij ij ij = H r H ij + i [applying Equations (8.12) and (8.11)] ij ri = 2dijH ij + Pijr i. 115 Since H ij + H ji = l ij and dij = dji, the sum of derivatives of A ij and A ji is: Pij A ij + A ji = r i r j l ij P ji 2dij zj 2dij ij zi 2ri i Figure 8.6: For a multi-disk configuration, the area derivatives are 2 rir i i and 2dijl ij ij due to the changes of ri and lij respectively. In a two disk configuration, the area derivative of a clipped disk is the perimeter Pij times the radius change r i. For a multi-disk configuration, we reuse i [0, 1] defined in the perimeter derivative theorem; area derivative due to radius change is 2 rir i i. Secondly, the area derivative due to change of lij is equal to the chord length between bi and bj, which is 2dij in the two-disk case. For a multi-disk configuration, define ij [0, 1] to be the fraction of chord length within the Voronoi cell. The rate of change of area due to change of lij is then l ij times the length of the chord between bi and bj, which is 2dij ij. See Figure 8.6. Area Derivative of Union of Disks Theorem. The area derivative of the union of disk, 116 B, is A B = 2 rir i i + 2dijl ij ij. bi B i =j The formula actually states that if there are decreases in the radii or distances between disk centers, the area decreases also. This related to the proof of Kneser-Poulsen conjecture with continuous movement of disk centers [14]. 8.3 2D Skin In this subsection we will derive the formulas of area, perimeter, as well as their derivatives, for a two-dimensional skin body. We will first introduce how to partition the space in order to break down the problem. Then we will study the center of a triangle. With these information, we will give the area and perimeter of a planar skin body. Then we will derive the derivatives of area and perimeter. After discussing the size and boundary changes, we will derive the formulas for the derivatives. For convenience, we define µi = µ{bi}, µij = µ{bi,bj}, and µijk = µ{bi,bj,bk} as mixed cells mentioned in Chapter 3. Also, we define ci, cij, and cijk to be the centers of the simplices zi, zizj, and zizjzk respectively. 8.3.1 Partition of Space To compute the perimeter and area of skin(B), as well as their derivatives, we partition the plane into mixed cells. Each mixed cell contains either a circle or a hyperbola; see Figure 8.7. Each mixed cell is the intersection of a collection of half-planes. When the disks move, the mixed cells change shape also because the motion of the half-planes as in Figure 8.8. All the half-planes are bounded by a line with some distances from the centers. These distances 117 zi h h ij ij ijk w wijk hji 0-mixed cell 1-mixed cell 2-mixed cell zj Figure 8.7: The skin is decomposed into circular and hyperbolic arcs. can be negative and will be defined as in Figure 8.7. Specifically, hij is negative if vertex bi is hidden by edge ij in the Delaunay triangulation. Similarly, wijk is negative if triangle ijk hides edge ij. 0-mixed Cell. Any Delaunay vertex bi determines 0-mixed cell µi in the partition. The cell is formed by shrinking the Voronoi polygon of bi by a factor of 2. For every Delaunay edge zizj, a corresponding half-plane clips bi. We define the distance from zi to the boundary of this half-plane to be hij. Note that hij = Hij/2, where Hij is defined in Section 8.2. 118 % % % ' ' ' ' ' ' & & & ( ( ( ( ( % % % ' ' ' ' ' ' & & & ( ( ( ( ( % % % ' ' ' ' ' ' & & & ( ( ( ( ( % % % ' ' ' ' ' ' & & & ( ( ( ( ( % % % & & & # # # # # # $ $ $ $ $ ! ! ! ! ! ! ! " " " " " " " # # # # # # $ $ $ $ $ ! ! ! ! ! ! ! " " " " " " " # # # # # # $ $ $ $ $ ! ! ! ! ! ! ! " " " " " " " Figure 8.8: Motion of boundaries in the 1-mixed cell. zj zj ij c ij c hji c h ijk ji wijk zi zi hij ij h zk Figure 8.9: Boundaries relationship of the hyperbola. 1-mixed Cell. A Delaunay edge zizj forms a rectangular cell µij with two kinds of bound- aries. If we align the Delaunay edge with the y-axis, the horizontal boundaries separate µij from µi and µj with a distance of hij and hji respectively. The vertical boundaries separate µij from µjim and µijk if triangle zjzizm and zizjzk exist in the Delaunay triangulation. We define the distance between cij to these two boundaries as wjim and wijk respectively. We assume wij = if there is no triangle attached. 119 2-mixed Cell. The mixed cell of a Delaunay triangle zizjzk is formed by shrinking the Delaunay triangle by a factor of one half. A circular void is clipped inside with its center at cijk. The distances from cijk to the boundaries of the neighboring 1-mixed cells, µij, µjk and µki, are wijk, wjki and wkij respectively. Formulas of the wijk will be given in Section 8.3.2. 8.3.2 Centers of Triangles y zk (p , p ) 1 2 c =(c ,c ) ijk 1 2 zi 2w z x ijk j (a ,0) 1 (b ,0) 1 Figure 8.10: Computing the center of a triangle: Aligning the edge zizj at x-axis and the center of the edge at origin. cijk = (c1, c2) = (0, 2wijk) will be the center of the triangle. We define wijk to be half of the perpendicular distance of the orthocenter cijk of triangle zizjzk to the edge zizj as in Figure 8.9. To compute wijk, we align the edge zizj to the x-axis such that the center cij of that edge is at the origin. Hence, cijk lies on the y-axis, and the only missing value is the y-coordinate of cijk. The configuration is shown in Figure 8.10. We let cijk = (c1, c2), zi = (a1, 0), zj = (b1, 0) and zk = (p1, p2) as in Figure 8.10. By Section V.18.2 of [42], putting a0 = 1(r2 (r2 (r2 2 i - a21), b0 = 12 j - b21), and p0 = 12 k - p21 - p22), we 120 obtain the followings: -a 0 0 1 a1 -a0 1 -b b 0 0 1 1 -b0 1 -p 0 p2 1 p1 -p0 1 c1 = and c2 = . a 1 0 1 a1 0 1 b b 1 0 1 1 0 1 p 1 p2 1 p1 p2 1 Simplifying these equations, we obtain: b c 0 - a0 1 = = 0 a a 0 = b0 if a1 = b1, 1 - b1 a c 1(-b0 + p0) + a0(b1 - p1) - b1p0 + p1b0 2 = p2(b1 - a1) p -a = - 0 + 1b0 + a0b1 p2 p2(b1 - a1) p -a = - 0 + 1a0 + a0b1 since a p 1 = b1 2 p2(b1 - a1) a = 0 - p0 . p2 Hence, c a w 2 0 - p0 ijk = = . (8.13) 2 2p2 8.3.3 Area and Perimeter of 2D Skin We now discuss the measurements of different kinds of mixed cells. 0-mixed Cell. A disk bi with radius ri/ 2 is clipped inside the 0-mixed cell. For every edge z j j izj we define bi = {x bi | x - zi, zj - zi > cij - zi, zj - zi /2}. bi is the portion of b j i on the µij side of the mixed cell wall that divides µi and µij. Define Pi 121 and P jk j j k j jk i to be the perimeters of bi and bi bi . Similarly, define Ai and Ai to be the areas. The perimeter and area of bi clipped in µi is P j jk i = 2 ri - Pi + Pi , j j,k r2 A i j jk i = - A + A . 2 i i j j,k 1-mixed Cell. Before we compute the perimeter and area of the hyperbola pieces in the 1-mixed cell, we first introduce the computation of the length of a general curve in an implicit form and the elliptic integral. Any curve can be represented in an implicit form x = f(y). The length of such curve whose y value within the range of [y1, y2] is equal to: y 2 dx 2 1 + dy. y dy 1 The elliptic integral of the second kind, E(u, v), is defined as: u 1 - v2t2 E(u, v) = dt. 0 1 - t2 The elliptic integral is not solvable analytically [18], but fast numerical solvers are available. To measure the hyperbola, we align the edge zizj to the y-axis and the center of the hyperbola to the origin, such that zi is above zj. Let the equation of the aligned hyperbola be x2 - y2 + R2ij = 0. The goal is to compute the measurements within the rectangle of [-wijk , w ] × [-h 1 ijk2 ji, hij]. We assume the origin is within the rectangle. The hyperbola can be divided into four quadrants by the x and y-axis. We first derive the formulas for the upper right quadrant, see Figure 8.11. Define ij = h2ij -R2ij to be the x value of the hyperbola endpoint when y = hij. We have four different configurations of the upper right quadrant with respect to different values of Rij, wijk and hij, see Figure 8.12. 122 ijk w hij Figure 8.11: The first quadrant of a hyperbola. y wijk y ijk w y wijk y wijk hij h hij h ij ij x x x x 2ij R >0, < w R >0, 2 ijk w R <0, 2 ijk R <0, 2 ijk ijk ij < ijk w ij > ijk wijk ij > ijk Figure 8.12: Four different possible configurations of the first quadrant. We will derive the formulas for the first case of the upper right quadrant. We express x in terms of y and differentiate x with respect to y : x = y2 - R2ij, dx y y = = . dy y2 - R2 x ij The length of a hyperbola piece from the tip, (0, Rij), to the point touching the mixed 123 y y h h ij ij R R ij ij x x w w ij ij Figure 8.13: The length and area of a hyperbola: For a hyperbola in standard form, x2 - y2 + R2ij = 0, we integrate its length and area from y = Rij to y = hij. cell boundary, ( ij, hij), is: hij dx P = 1 + ( )2dy R dy ij hij 2y2 - R2 y = ij dy [substituting t = ] R y2 - R2 R ij ij ij hij/Rij 1 - 2t2 = Rij dt 1 1 - t2 h ij /Rij 1 - 2t2 1 1 - 2t2 = Rij dt - dt 0 1 - t2 0 1 - t2 = Rij(E(hij/Rij, 2) - E(1, 2)). Here, E is the elliptic integral of the second kind. Since 1-2t2 1-t2 is imaginary when t [ 0.5, 1], E(1, 2) has an imaginary part. However, it is canceled out if hij Rij, which is the normal case. Hence, we have P = Real(Rij(E(hij/Rij, 2) - E(1, 2))). (8.14) 124 The area, A, from the tip Rij to height hij is: hij A = y2 - R2ijdy Rij 1 hij = y y2 - R2 y2 - R2 2 ij - R2ij ln(y + ij) Rij 1 = h h2 h2 2 ij ij - R2ij - R2ij ln(hij + ij - R2ij) + R2ij ln Rij 1 h = h ij + ij . (8.15) 2 ij ij - R2ij ln Rij We skip the derivation of other 3 cases. They can be computed by addition and subtraction of the measurements. The results of the four cases for the perimeter in the first quadrant is given here as: f1(Rij, min(hij, w2 P ijk + R2ij)) if R2ij > 0 ijk = (8.16) f1( -1Rij, min( ijk, wijk)) if R2ij < 0 and area: f 2(Rij, hij) if R2ij > 0, ijk < wijk f w2 2(Rij, ijk + R2ij) Aijk = , (8.17) +(h w2 ij - ijk + R2ij)wijk if R2ij > 0, ijk > wijk -f -1R 2( ij, min(wijk, ijk)) +hij min(wijk, ijk) if R2ij < 0. where: f1( , ) = Real( (E( / , 2) - E(1, 2)), and, (8.18) 1 + 2 - 2 f2( , ) = 2 - 2 - 2 ln . (8.19) 2 125 We name the perimeter in a quadrant, Pijk, to be the one on the zi and zk sides in µij. For example, if we align the Delaunay edge zizj to the y-axis with zi above zj, such that zk is on the right and zm is on the left, the perimeters of four quadrants in a clockwise order will be named as Pijk, Pjik, Pjim and Pijm starting from the upper right one. The same index rule is applied on the area Aijk. For the four pieces of hyperbola inside µij, we sum up the four quadrants for perimeter and area with sign correction: P = Sign(hij) × Sign(wijk) × Pijk + Sign(hji) × Sign(wjik) × Pjik + Sign(hji) × Sign(wjim) × Pjim + Sign(hij) × Sign(wijm) × Pijm, A = Sign(hij) × Sign(wijk) × Aijk + Sign(hji) × Sign(wjik) × Ajik + Sign(hji) × Sign(wjim) × Ajim + Sign(hij) × Sign(wijm) × Aijm, Figure 8.14: The eight cases of a disk intersecting a triangle. 2-mixed Cell. For a triangle, the area of the skin body inside is equal to the area of the triangle subtracting the area of the clipped disk in the triangle. All possible combinations of intersection are shown in Figure 8.14. We denote the disk bijk to be the disk in the triangle with radius Rijk > 0 and center cijk. The distance from cijk to ci is 4w2ijk + 4h2ij. 2R2ijk is the squared radius of a disk with 126 center cijk that is orthogonal to bi. With this relationship, we have: 2R2ijk + r2i = 4w2ijk + 4h2ij 4w2 R ijk + 4h2ij - r2i ijk = . (8.20) 2 We define b ij ij ijk = {x bijk | x - cijk, cij - cijk > cij - cijk 2/2}. Define Pijk and P k ij ki jk ij k ijk to be the perimeter of bijk and bijk bijk . Moreover, Aijk and Aijk as the areas. The perimeter and area of the skin body in a 2-mixed cell is: P = 2 R ij k ijk - Pijk + Pijk , i,j k A = Area(µ ij k ijk) - R2ijk + Aijk - Aijk . i,j k 8.3.4 Derivatives of Area and Length Similar to the case of union of disks, we decompose the computation of the derivatives into parts. The skin model is relatively more complicated. The approach to derive the formulas is decomposed into 3 parts: 1. Changes of R : R and R . ri zi 2. Boundary motion: The motion, (z i, r i), leads to shape change of the mixed cell, hij , hij , wijk and wijk . zi ri zi ri 3. Change of area and length due to boundary motion and R : L , A , hij hij L , A . L and A . wijk wijk R R 8.3.5 Simplex Size Changes There are three types of simplices in 2D Delaunay triangulation which form three types of mixed cells. 127 0-mixed Cell. The radius of the disk, Ri, within the mixed cell is ri/ 2 as stated in Section 8.3.3. Hence, R i = r i/ 2. 1-mixed Cell. For the 1-mixed cell, we recall that sX is the size of the simplex X with X = {bi, bj}. In order to fit in the previous area computation of hyperbola, we define 2R2ij = -s2X. The reason is first Rij is equal to 1/ 2 of sX in magnitude. Also, when s2X is positive, we want R2ij to be negative to meet the computation of measurements in Section 8.3.3. The disk with center cij , and radius -s2X, is orthogonal to bi. The distance between the two centers, cij and ci = zi, is 2hij. The orthogonality implies -s2X + r2i = (2hij)2 and hence: 2R2ij = 4h2ij - r2i, 4h2 R ij - r2i ij = . (8.21) 2 The partial derivatives of Rij with respect to hij and ri: Rij 2 = 2h h ij ij 4h2ij - r2i 2h = ij , Rij Rij r 2 = - i hij 2 4h2ij - r2i r = - i . 2Rij The derivative of Rij is: R R R ij ij ij = h r h ij + i ij ri 1 r = 2h i r . (8.22) R ijh ij - i ij 2 128 2-mixed Cell. From Equation (8.20), the partial derivatives of Rijk are: Rijk 2w = ijk , wijk Rijk Rijk 2h = ij , and hij Rijk Rijk r = - i . ri 2Rijk Hence, the derivative of Rijk is: 1 r R ir i ijk = 2w . (8.23) R ijkw ijk - 2hijh ij + ijk 2 8.3.6 Boundary Motions We can express the boundary distances of mixed cells in terms of hij and wijk. The derivative of hij is equal to (Hij) /2 and H ij is computed in Equation (8.4), so: 1 1 r2 r r h i - r2j i j ij = - l r r . (8.24) 2 2 2l2 ij + i - j ij lij lij For the derivative of wijk, we follow Equation (8.13). Recall that wijk = a0-p0 in Section 2p2 p' -((p -a )a' +(b -p )b' )/l 2 1 1 2 1 1 2 ij y 2 p' y z :k (p , p ) z : (p , p ) 1 2 k 1 2 p'1 p'1 a'2 b'2 c=(c ,c ) c=(c ,c ) 1 2 2w x 1 2 a' ijk 2w x 1 a' ijk 1 (a ,0) b' b' i z : 1 z : (a ,0) 1 1 z : (b ,0) j i z : (b ,0) 1 1 j 1 Figure 8.15: General movement of vertices in a triangle (left), and rearrangement of the movement to fix an edge on x-axis (right). 129 8.3.2 with p0 = (r2p - p21 - p22)/2 and a0 = (r2i - a21)/2. The partial derivatives are: wijk p = 1 , p1 2p2 wijk 1 a = - 0 - p0 , p2 2 2p22 wijk r = - p , rp 2p2 wijk r = i , ri 2p2 wijk a = - 1 . a1 2p2 However, we assumed that the two vertices zi and zj are on the x-axis as in Figure 8.10. In general cases, the vertices will move vertically to leave the x-axis also. Let a 2 and b 2 be their vertical velocity as shown in Figure 8.15 (left). To take this movement into account, we can still fix the edge on the x-axis but we subtract the motion a 2 and b 2 from p 2 as in Figure 8.15 (right). The final vertical motion of zk will be p 2 -((p1 -a1)a 2 +(b1 -p1)b 2)/lij. Hence, the final w ijk is : p 1 a (p w 1 0 - p0 1 - a1)a 2 + (b1 - p1)b 2 ijk = p - p 2p 1 + 2 - 2 2 2p22 lij r a r b r - p r 1 a i r 1 b j r 2p p - 1 + i - 1 + j. (8.25) 2 2p2 2p2 2p2 2p2 130 8.3.7 Derivatives from Boundary Motions of Hyperbola Before we derive the derivatives with respect to hij, wijk and ri, we first compute the deriva- tives of f1 and f2 defined in Equations (8.18) and (8.19). First, f 1 is: f f f 1 1 1( , ) = + 1 2 - 2 2 2 - 2 2 = f + 1( , ) - 2 - 2 2 - 2 f 2 2 - 2 = 1( , ) + - . (8.26) 2 - 2 Note that for any point, ( 2 - 2, ), on the hyperbola x2 - y2 + = 0, 2 2 - 2 is the distance between the upper endpoint and the origin. Secondly, the derivative of f2 is: f f f 2 2 2( , ) = + = ln + 2 - 2 . (8.27) + 2 - 2 Differentiate the measurements in the first quadrant in Equations (8.16) and (8.17), we get: f w2 1(Rij, hij) if R2ij > 0, hij < ijk + R2ij f w2 w2 P 1(Rij, ijk + R2ij) if R2ij > 0, hij > ijk + R2ij ijk = , f -1R 1( ij, ijk) if R2ij < 0, wijk > ijk f 1( -1Rij, wijk)) if R2ij < 0, wijk < ijk 131 with ijk = h2ij - R2ij. We derive each of the formulas into: f 2h2 h f 1(Rij, hij) ij - R2ij ijR ij 1(Rij, hij) = R h , R ij + ij - ij h2ij - R2ij Rij f 1(Rij, w2ijk + R2ij) 2w2 w f ijk + R2ij ijkR ij 1(Rij, w2ijk + R2ij) = R w , R ij + ijk - ij w2ijk + R2ij Rij f -1R 2w2 w f 1( ij, wijk) ijk + R2ij ijkR ij 1( -1Rij, wijk) = R w , R ij + ijk - ij w2ijk + R2ij Rij f -1R 2h2 h f 1( ij, ijk) ij - R2ij ijR ij 1( -1Rij, ijk) = R h . R ij + ij - ij h2ij - R2ij Rij For area derivatives, we have: f 2(Rij, hij) if R2ij > 0, ijk < wijk f w2 w2 A 2(Rij, ijk + R2ij) + ((hij - ijk + R2ij)wijk) if R2ij > 0, ijk > wijk ijk = , -f -1R 2( ij, ijk) + (hij ijk) if R2ij < 0, ijk < wijk -f 2( -1Rij, wijk) + (hijwijk) if R2ij < 0, ijk > wijk with the first case equals to: R f ij 2(Rij, hij) = Rij ln R ij + h2ij - R2ijh ij. hij + h2ij - R2ij 132 The second case is equal to: f 2(Rij, w2ijk + R2ij) + ((hij - w2ijk + R2ij)wijk) R = R ij ij ln R ij + wijk( w2ijk + R2ij) wijk + w2ijk + R2ij +(hij - w2ijk + R2ij)w ijk + (hij - w2ijk + R2ij) wijk R = R ij ij ln R ij + (hij - w2ijk + R2ij)w ijk + wijkh ij. wijk + w2ijk + R2ij The third case is: -f 2( -1Rij, ijk) + (hij ijk) -1R = R ij ij ln R ij - hij ijk + hij ijk + ijkh ij ijk + hij -1R = R ij ij ln R ij + ijkh ij. ijk + hij Finally, the fourth case is: -f 2( -1Rij, wijk) + (hijwijk) -1R = R ij ij ln R ij - w2ijk + R2ijw ijk + (wijkhij) wijk + w2ijk + R2ij -1R = R ij ij ln R ij + (hij - w2ijk + R2ij)w ijk + wijkh ij. wijk - w2ijk + R2ij 8.3.8 Assembly of Equations of Area Derivative of Skin Recall the definitions of P j jk i and Pi in Section 8.3.3. We define i as the fraction of perimeter in µi as: j jk i = 1 - ( Pi - Pi )/( 2 ri). j j,k 133 Also we define ij to be the fraction of chord length of bi on µi µij. With the direction from zi to zj, we define +ij to be 1 when the endpoint of the arc of bi exists on the right, and 0 if it does not. Similarly, let -ij represents the other for the left point. The distance- preserving components are Sign(z i × zizj) · z i - ¯z ij at zi and Sign(z j × zjzi) · z j - ¯z ji at zj in a counter-clockwise fashion. The rate of change in perimeter due to this effect is: 2 -(Sign(z i × zizj) · z i - ¯z ij + Sign(z j × zjzi) · z j - ¯z ji ) · ( +ij - -ij)lij · . ri For the distance-preserving motion, we put = ri/ 2 and = hij into Equations (8.5) and (8.6) and get: 2h h 2r P = 2 - arccos ij - ij r i h r i + ij i r2i/2 - h2ij r2i/2 - h2ij 2h 2 = 2 - arccos ij + (r r ih ij - hijr i. i r2i/2 - h2ij The perimeter derivative of the clipped disk, bi in µi is: + 2(r P ij + - ij ih ij - hijr i) i = ir i + 2 j r2i/2 - h2ij 2 -(Sign(z i × zizj) · z i - ¯z ij + Sign(z j × zjzi) · z j - ¯z ji ) · ( +ij - -ij)lij · . ri For area, we have derivative of bi in µi as: A i = 2 ri ir i + 2 r2i/2 - h2ij h ij j 134 1-mixed Cell. The perimeter and area derivatives of hyperbola in µij is: P = Sign(hij) × Sign(wijk) × P ijk + Sign(hji) × Sign(wjik) × P jik + Sign(hji) × Sign(wjim) × P jim + Sign(hij) × Sign(wijm) × P ijm, A = Sign(hij) × Sign(wijk) × A ijk + Sign(hji) × Sign(wjik) × A jik + Sign(hji) × Sign(wjim) × A jim + Sign(hij) × Sign(wijm) × A ijm, 2-mixed Cell. Let bijk to be the disk in µijk with center cijk and radius Rijk. We define ijk as: P ij - P k i,j ijk k ijk ijk = 1 - , 2 Rijk which is the fraction of perimeter of bijk in µijk. We also define xijk and xjik to be the zk cijk zi ijk w xijk xjik zj Figure 8.16: The naming of points in µijk. points of intersection of bijk with aff(zizj), such that xijk is closer to zi and as xjik closer to zj. Define ijk to be 1 if xijk exist inside µijk and 0 otherwise, and same to jik, as in Figure 8.16. If we consider only the slope-preserving motions of µijk, each edge zizj contributes a rate 135 of change in perimeter of bijk as follows: w 2(R P = 2 - arccos ijk R ijkw ijk - wijkR ijk) . R ijk + ijk R2ijk - w2ijk according to Equations (8.5) thru (8.6). Hence, the derivative of perimeter in µijk is: 2(R P = ijkw ijk - wijkR ijk) ijk + jik ijkR ijk + · + 2 ij R2ijk - w2ijk R ijk(Sign( ¯ z ji, zjzi ) ¯z ji + h ji - Sign( ¯z ij, zizj ) ¯z ij - h ij)( ijk - jik) . wijk Let ijk = length(xijkxjik zizj)/length(xijkxjik). The area derivative is: A = area(µijk) - ijkRijkR ijk - 2 R2ijk - w2ijk ijkw ijk, i,j where area(µijk)' is (lijw ijk + ljkw jki + lkiw kij)/4. the boundary change 8.3.9 Theorems Skin Perimeter Theorem. The length of the skin curve of B is the sum of perimeter of the three types of mixed cells. The perimeter in all 0-mixed cell is: 2 r j jk i - Pi + Pi . i j j,k The perimeter in all 1-mixed cell is: Sign(hij)×Sign(wijk)×Pijk +Sign(hji)×Sign(wjik)×Pjik + ij Sign(hji) × Sign(wjim) × Pjim + Sign(hij) × Sign(wijm) × Pijm , 136 which f1(Rij, min(hij, w2 P ijk + R2ij)) if R2ij > 0 ijk = f1( -1Rij, min( ijk, wijk)) if R2ij < 0 where ijk = h2ij - R2ij and f1( , ) = Real( (E( / , 2) - E(1, 2)). The perimeter in all 2-mixed cell is: 2 R ij k ijk - Pijk + Pijk . ijk i,j k Skin Area Theorem. The area of the skin body of B is the sum of area of the three types of mixed cells. The area in all 0-mixed cell is: r2 i - A j + A jk . 2 i i i j j,k The area in all 1-mixed cell is: Sign(hij)×Sign(wijk)×Aijk +Sign(hji)×Sign(wjik)×Ajik + ij Sign(hji) × Sign(wjim) × Ajim + Sign(hij) × Sign(wijm) × Aijm , with f 2(Rij, hij) if R2ij > 0, ijk < wijk Aijk = , f w2 w2 2(Rij, ijk + R2ij) + (hij - ijk + R2ij)wijk if R2ij > 0, ijk > wijk -f2( -1Rij, min(wijk, ijk)) + hij min(wijk, ijk) if R2ij < 0; 137 where: 1 + 2 - 2 f2( , ) = 2 - 2 - 2 ln . 2 The area in all 2-mixed cell is: Area(µ ij k ijk) - R2ijk + Aijk - Aijk . ijk i,j k Skin Perimeter Derivative Theorem. The perimeter derivative of the skin body of B is the sum of perimeter derivative of the three types of mixed cells. The perimeter derivative in all 0-mixed cell is: + 2(r ij + - ij ih ij - hijr i) ir i + 2 i j r2i/2 - h2ij 2 -(Sign(z i × zizj) · z i - ¯z ij + Sign(z j × zjzi) · z j - ¯z ji ) · ( +ij - -ij)lij · . ri The perimeter derivative in all 1-mixed cell is: Sign(hij)×Sign(wijk)×P ijk +Sign(hji)×Sign(wjik)×P jik + ij Sign(hji) × Sign(wjim) × P jim + Sign(hij) × Sign(wijm) × P ijm . with f 2h2 1(Rij,hij) ij -R2ij R h if R R ij + h2 ij - hijR ij R ij > 0, ijk < 0 ij ij -R2ij ij f1(Rij, w2ijk+R2ij)R 2w2ijk+R2ij w if R R ij + w2 ijk - wijkR ij R ij > 0, ijk > 0 P ij ijk+R2ij ij ijk = 2w2 f1( -1Rij,wijk) ijk+R2ij R w if R ij + ijk - wijkR ij ij < 0, wijk < ijk Rij w2 Rij ijk+R2ij f1( -1Rij, ijk)R 2h2ij-R2ij h if R R ij < 0, wijk > ijk ij ij + h2 ij - hijR ij ij -R2ij Rij 138 with ijk = hij - w2ijk + R2ij. The perimeter derivative in all 2-mixed cell is: 2(R ijkw ijk - wijkR ijk) ijk + jik ijkR ijk + · + 2 ijk ij R2ijk - w2ijk R ijk(Sign( ¯ z ji, zjzi ) ¯z ji + h ji - Sign( ¯z ij, zizj ) ¯z ij - h ij)( ijk - jik) . wijk Skin Area Derivative Theorem. The area derivative of the skin body of B is the sum of area derivative of the three types of mixed cells. The area derivative in all 0-mixed cell is: 2 ri ir i + 2 r2i/2 - h2ij h ij . i j The area derivative in all 1-mixed cell is: Sign(hij)×Sign(wijk)×A ijk +Sign(hji)×Sign(wjik)×A jik + ij Sign(hji) × Sign(wjim) × A jim + Sign(hij) × Sign(wijm) × A ijm , which R ln Rij R h2 ij ij + ij - R2ijh ij if R2ij > 0, ijk < wijk hij+ h2 ij -R2ij Rij ln Rij R ij + ijkw ijk + wijkh ij if R2ij > 0, ijk > wijk A wijk+ w2ijk+R2ij ijk = . -1Rij R ln R ij ij + ijkh ij if R2ij < 0, ijk < wijk ijk+hij R -1Rij ij ln R w ij + ijkw ijk + wijkh ij if R2ij < 0, ijk > wijk ijk- w2ijk+R2ij 139 with ijk = hij - w2ijk + R2ij. The area derivative of all 2-mixed cell is: l ijw ijk + ljkw jki + lkiw kij - 2 R2 . 4 ijkRijkR ijk - ijk - w2ijk ijkw ijk ijk i,j 8.4 Conclusion In this chapter, we investigate the measurements of the union of disks and skin in order to provide deformable shapes with discrete boundaries in the structural optimization process. We compute the formulas of the area, the perimeter, and their derivatives of these shapes. It would also be interesting to analyze the second derivatives, which could be useful in accelerating the global design cycle of topology optimization. The optimization process needs knowledge of these measurements within a window, which is a disk in the two dimensional case. Although we did not complete this window computa- tion, the calculation of a portion of the shapes remains the same with the formulas given in this chapter, if the entire mixed cell of that portion is contained in the window. The new consideration is the mixed cells that intersect the boundary of the window, and this involves clipping of the shapes. We also did not solve the second part of the problem, namely, given a point x, find the closest point p on the shapes and compute the their distances and tangent lines on p. This can be done by first locating the Voronoi cell which x lies by point location algorithms [31] and then search p by traversing the neighboring Voronoi cells or mixed cells. Another direction of this research is the computation of measurements in the three di- mensional case, namely, the surface area, volume, and their derivatives of the unions of ball and skin body in R3. The study in two dimensions gives a blueprint of the computations in three dimensions but we expect that the cases in R3 are more complicated. For spherical patches, the inclusion-exclusion formula is still applicable for the measurements, but not for hyperbolic patches in the skin case. The clipping of a hyperboloid with a cylinder and the computation of its measurements may be difficult, not to mention the window computation. 140 8.5 Table of Variables Variable names Description The radius of a general disk The distance between a clipping half-plane to the center of a disk Half the chord length, which is 2 - 2 arccos( / ) ri The radius of a disk bi zi The center of a disk bi lij The length of zizj, also equals to Hij + Hji l ij, z i, r i The rate of change of lij,zi and ri Hij lij/2 + (r2i - r2j)/2lij dij Half the chord length of clipped disk bi, which is r2i - H2ij i The fraction of exposed circle ij Indicator of corners in the disk union ij The fraction of chord of bi µij c Centers of simplices ¯z ij The vector component of z i in the zizj direction hij Hij/2 wijk cij - cijk /2, negative if ijk hides ij E Elliptic integral of second kind Ri The radius of bi = ri 2 Rij x2 - y2 + R2ij = 0 is the hyperbola for the 1-mixed cell Rijk The radius of disk in µijk +ijk, -ijk The indicators of corners in skin curve 141 Chapter 9 Results and Future Applications One of the major works in this thesis is the construction of the skin software, which is implemented on the PC platform, with C++ as the language and OpenGL as the graphics library [21]. In this chapter, we first show some of the software compatibilities, namely, the static and dynamic skin triangulations. Then we demonstrate some skin surfaces and describe potential applications. 9.1 Software We present two versions of the skin software, namely the static and dynamic versions, and some comparisons between the two implementations. 9.1.1 Static Skin and Refinement Figure 9.2 shows the backbone (or the stick model) of a molecule from the Protein Data Bank [1]. Disregarding the hydrogen atoms, we generate the skin surface of the molecule with van der Waals radii. Figure 9.3 (left) shows the mesh of the restricted Delaunay triangulation of a random point set, which is generated according to a probability density proportional to the local curvature of the surface. We can see the bad triangles with large and small angles. After refinement is applied, the triangle qualities are improved as in Figure 9.3 (right). We proved that the minimum angle is larger than 21 , but experimental results show the minimum 142 Figure 9.1: The skin software. angle can go up to 23 , and sometimes even 26 . The overall time for computing the surface is 1 minute and 50 seconds on a Pentium III 400MHz machine and the program is complied by Microsoft Visual C++ under release setting. The minimum angle of the triangles on the surface is 23.91 . More performance of the static skin computation is given in Table 9.1. Although the mesh quality is good, the computing time is long for molecules with more than 2000 atoms. model name computing no. of no. of vertices no. of time defining spheres on surface triangles Simple (Figure 9.7) 00:00:10 4 2376 4748 Simple protein (Figure 9.2) 00:01:44 34 13623 27242 Gramicidin A (Figure 9.5) 02:05:12 319 167832 335740 cdk2 molecule (Figure 9.4) 13:03:34 2370 591179 1182006 Table 9.1: Performance of static skin algorithm. We proved that the refinement will terminate with parameters C = 0.08 and Q = 1.65. Recall that C controls the density of triangles, and Q controls the quality. Since the minimum 143 Figure 9.2: The backbone of a molecule. angle is arcsin 1Q2, we can see that the quality of triangles improves as the parameter Q decreases. However, the minimum angle cannot be guaranteed by the theoretical proof when Q continues to decrease. Table 9.2 shows some experimental results that occur beyond the theoretical minimum angle of arcsin (1/Q2) in which the minimum angle cannot be guaranteed anymore. C 0.08 0.08 0.08 Q 1.65 1.5 1.45 min. angle 26.37 27.36 28.14 arcsin (1/Q2) 21.55 26.39 28.40 Table 9.2: The relationship between Q and the minimum angle. For the density of triangles, Figure 9.6 shows the effect of changing C. As C decreases, the density of triangles increases and the mesh approximates the analytical surface better. The trade-off is that more triangles are used, so the surface consumes more memory and degrades the performance. 9.1.2 Dynamic Triangulation Figure 9.7 shows an example of the dynamic triangulation. The skin in the demonstration is defined by four spheres. The algorithm first creates four spheres as in Figure 9.7 (a). When 144 Figure 9.3: Two skin surfaces of the above molecule. two of the spheres get closer, the tips of the two hyperboloid enter the hot sphere and start the special sampling. Figure 9.7(b) shows the special caps created inside the hot sphere with m = 5 vertices as the bases. The white dots indicate the points involved in the special sampling. When the two tips touch each other, the topology changes and a double cone is formed as in Figure 9.7(c). Then a waist is formed as in Figure 9.7(d) and the points leave the hot sphere. Hence, the special sampling ends for this topology change. Figure 9.7(e) shows another kind of change where a tunnel is filled. A special waist is formed when the tunnel is small enough to touch the hot sphere as in Figure 9.7(f). It will then be filled later. Figure 9.7(g) shows such a filled tunnel and that the inside and outside surfaces are separated from each other. A void is created inside the body. After this void vanishes, the skin is completed as in Figure 9.7(h). The time for computing the skin surface in Figure 9.3 (right) by dynamic skin algorithm is about 30 minutes. It is slower than the static skin algorithm because more time is spent on scheduling in the priority queue. On the other hand, the time for the static skin algorithm is proportional to the weighted surface area and to the weighted volume in the dynamic case. Here, the weighted surface area and volume are the integral of the curvature over 145 Figure 9.4: cdk2 molecule. Front view Side view Figure 9.5: Two views of Gramicidin A. the surface and volume respectively. However, for the weighted volume, the curvature is undefined at the positions of topology change. More specifically, the dynamic algorithm takes time proportional to the weighted volume outside the hot spheres. For objects with higher weighted volume to surface area ratio, the dynamic skin has better performance. In Figure 9.8, the four skin surfaces have different ratios and their computing time is shown in Table 9.3. 146 Figure 9.6: Two different values of parameter, (left) C = 0.08 and (right) C = 0.04. Static Skin Dynamic Skin 0.55 00:30s 01:13 0.52 00:39s 00:50 0.51 00:45s 00:42 0.5001 01:43 00:31 Table 9.3: Time comparison for static and dynamic skin algorithms. 9.2 Future Applications of Skin Surface We now illustrate some more skin surfaces we created with our software and describe two closely related applications: molecular modeling and computer graphics. 9.2.1 Molecular Modeling Chemists and biologists study how compounds react [66]. After a scientist knows which compounds can react with each other, the binding and separation of the molecules are investigated [63]. This can be facilitated by modeling and visualizing of the molecules [11]. Researchers also study how these molecules deform, for example, the protein folding problem [48]. Geometric information of molecules need a deformable model to approximate the actual behavior. The actual behavior of the molecules can be model through the geometric 147 information of its deformable model. The skin surface has the intuition of formation by spheres as well as the deformable property. By the complementarity property of the skin surface we can compute the void, pocket or tunnel of a molecule [41]; see Figure 9.9. This computation is useful in the ligand docking problem [43] [62]. 9.2.2 Computer Animation Figure 9.10 shows how skin surfaces can model some geometric objects represented by a set of spheres [6]. Objects in a scene do not need every detail when they are being viewed from far away. By choosing different values we can have different levels of detail. This property may help represent objects in multiple resolutions [80]. Moreover, the skin is a deformable surface that can handle topology changes. Figure 9.11 shows a sequence of morphing between two shapes. Each instance of the skin surface is generated independently. In the future, we can investigate how to deform one instance to another directly. This property can help in computer animation problems as well as the study of molecular motion. It is possible to approximate the electron distributions with the skin surface while a molecule is deforming during a chemical reaction to calculate electric potentials on the surface. 148 (a) (b) (c) (d) (e) (f) (g) (h) Figure 9.7: Transitions of the dynamic skin triangulation. 149 =0.55 =0.52 =0.51 =0.5001 Figure 9.8: Four skin surfaces with different weighted surface area and volume. Top view Front view Figure 9.9: Two views of Gramicidin A with its tunnel. 150 Figure 9.10: A human head modeled by skin surface with two different values. 151 Figure 9.11: A sequence of morphing from a question mark to a human figure. 152 Chapter 10 Conclusion In this thesis, we study the skin surface as a new paradigm for the deformable surfaces. The skin surface is an implicit surface that can deform, and at the same time, process the ability to be decomposed, parametrized and triangulated. The deformation ability gives the skin freedom to model deformable structures such as molecules. We also construct two pieces of skin software that give a NURBS representation and (static and dynamic) piecewise-linear approximation of the skin surface. This work also contributes to the following areas: 1. Computational geometry 2. Implicit/parametric surfaces 3. Meshing 4. Structural optimization The skin surface captures the advantages of both implicit and explicit (parametric) sur- faces. The implicits, or level-sets, are computationally expensive in surface meshing. Also, computing a parametrization on the surface, expecially on implicits, cannot be effective without topology information. Parametric surfaces often model objects with a large number of high degree patches in order to maintain curvature and surface normal continuity. In addition, the degree and number of patches are often the limiting factors of most algorithms, such as ray-tracing. The skin overcomes these disadvantages with the integration of low- degree patches in space partitions, parametrization, and meshing algorithms as introduced in this thesis. 153 The skin has a high degree of freedom to model objects and it is computationally efficient with partitions of space. However, it is a smooth surface with inability to model objects with sharp edges or corners. This problem is not serious in R2 but becomes acute in R3. This is a critical drawback to modeling general objects and trade-off for smoothness. Further investigation is expected in several aspects. First, we investigate the preservations of the mesh integrity during growing, but not for general movements. General movements are much harder than the growth model because deformation of the mixed complex is involved and there will be more types of topology changes other than creations of simplices. Secondly, all meshing and parametrizing algorithms are derived on the skin surface with shrinking an infinite collection of spheres by a factor of one half. The maximum curvature continuity will be destroyed by factors other than 1/2 and the proof of homeomorphism needs to be fixed by using a local feature size other than the maximum curvature. New triangulation algorithms need to be investigated. Third, the skin can be applied to more fields. We model some molecules with the skin, but numerical computations are yet to be carried out. The measurements of skin are useful information in the topology optimization. We investigate the two-dimensional cases. Three dimensional skin surfaces and bodies are more complicated but definitely computable. The volume, surface area and their derivatives of the skin body are useful in three dimensional structural designs. Other than scientific and engineering disciplines, it would be interesting to see the skin surfaces in fine arts such as sculpting. 154 References [1] The Protein Data Bank. URL: http://www.rcsb.org/pdb/ . [2] Agarwal, P., Basch, J., Guibas, L., Zhang, L., and Hershberger, J. De- formable Free Space Tiling for Kinetic Collision Detection. New Directions in Alg. and Computational Robotics (2001), 83­96. [3] Agarwal, P., Eppstein, D., Guibas, L., and Henzinger, M. Parametric and Kinetic Minimum Spanning Tree. Proc. 39th Symp. Foundations of Comp. Sc., IEEE (1998), 596­605. [4] Amenta, N., and Bern, M. Surface Reconstruction by Voronoi Filtering. Discrete and Comp. Geometry 22 (1999), 481­504. [5] Amenta, N., Bern, M., and Kamvysselis, M. A New Voronoi-based Surface Reconstruction Algorithm. Proc. SIGGRAPH 1998 (1998), 415­421. [6] Amenta, N., Choi, S., Dey, T., and Leekha, N. A Simple Algorithm for Home- omorphic Surface Reconstruction. Proc. 16th Ann. Sympos. Comput. Geom. (2000), 213­222. [7] Aurenhammer, F. Voronoi Diagrams - a Study of a Fundamental Geometric Data Structure. ACM Comput. Surveys 23 (1991), 345­405. [8] Baker, T. J. Mesh Movement and Metamorphosis. Proc. 10th international Meshing Roundtable (2001), 387­396. 155 [9] Basch, J. Kinetic Data Structures. Ph.D. Thesis, Computer Science Dept., Stanford U., URL: http://graphics.stanford.edu/ jbasch/publications/phdthesis/ . [10] Bendsøe, M. Optimization of Structure Topology, Shape and Material. Springer- Verlag, 1995. [11] Bergman, L. D., Richardson, J. S., Richardson, D. C., and Brooks, F. P. VIEW - An Exploratory Molecular Visualization System with User-Definable Interac- tion Sequences. Proc. SIGGRAPH 1993 (1993), 117­126. [12] Bern, M., Eppstein, D., Agarwal, P. K., Amenta, N., Chew, P., Dey, T., Dobkin, D. P., Edelsbrunner, H., Grimm, C., Leonidas, Guibas, J., Harer, J., Hass, J., Hicks, A., Johnson, C. K., Lerman, G., Letscher, D., Plass- mann, P., Sedgwick, E., Snoeyink, J., Weeks, J., Yap, C., and Zorin, D. Emerging Challenges in Computational Topology. NSF Workshop on Computational Topology, URL: http://www.cis.ohio-state.edu/ tamaldey/paper/topowork/paper.ps.gz, 1999. [13] Bern, M., Eppstein, D., and Gilbert, J. R. Provably Good Mesh Generation. 31st Annual Symposium on Foundations of Computer Science, IEEE (1990), 231­241. [14] Bezdek, K., and Connelly, R. Pushing Disks Apart - The Kneser-Poulsen Conjec- ture in the Plane. http://www.math.cornell.edu/ connelly/kneser.html. [15] Blaney, J., and J.S.Dixon. A Good Ligand is Hard to Find: Automated Docking Methods. Perspectives in Drug Discovery and Design 1 (1993), 301­319. [16] Blinn, J. A Generalization of Algebraic Surface Drawing. ACM Transactions on Graphics 3 (1982), 235­256. 156 [17] Bloomenthal, J., Bajaj, C., Blinn, J., Marie-Paule, Cani-Gascuel, Rock- wood, A., Wyvill, B., and Wyvill, G. Introduction to Implicit Surfaces. Morgan Kaufmann Publishers, Inc, 1997. [18] Bowman, F. Introduction to Elliptic Functions with Applications. Dover, 1961. [19] Cannon, J. W., Floyd, W. J., Kenyon, R., and Parry, W. R. Hyperbolic Geometry. Flavors of Geometry (1997), 59­115. [20] Chen, H., and Bishop, J. Delaunay Triangulation for Curved Surfaces. Proceedings, 6th International Meshing Roundtable (1997), 115­127. [21] Cheng, H. The Skin Software. URL: http://www.cs.duke.edu/ hcheng/Research/ . [22] Cheng, H. The Construction and Display of a Smooth Surface Defined by Weighted Points in Three Dimension. Master's thesis, USTHK, 1996. [23] Cheng, H. Quadratic Patches in NURBS Form. manuscript, URL:http://www.cs.duke.edu/ hcheng/NURBS (1997). [24] Cheng, H., Dey, T. K., Edelsbrunner, H., and Sullivan, J. Dynamic Skin Triangulation. Discrete Comput. Geom. 25 (2001), 525­568. [25] Cheng, H., Edelsbrunner, H., and Fu, P. Shape Space from Deformation. Proc. 6th Pacific Conf.. Comput. Graphics Appl. (1998), 104­113. [26] Cheng, S., Edelsbrunner, H., Fu, P., and Lam, P. Design and Analysis of Planar Shape Deformation. Comput. Geom. Theory Appl. 19 (2001), 205­218. [27] Chew, L. P. Guaranteed-Quality Mesh Generation for Curved Surfaces. Proceedings of the Ninth Symposium on Computational Geometry (1993), 274­280. [28] Connolly, M. Analytical Molecular Surface Calculation. J. Appl. Cryst. 16 (1983), 548­558. 157 [29] Creighton, T. Proteins: Structures and Molecular Principles. Freeman, New York, 1984. [30] Curless, B., and Levoy, M. A Volumetric Method for Building Complex Models from Range Images. Proc. SIGGRAPH 1996 (1996), 303­312. [31] de Berg, M., van Kreveld, M., Overmars, M., and Schwarzkopf, O. Com- putational Geometry: Algorithms and Applications. Springer-Verlag, 1997. [32] de Cougny, H. S., and Shephard, M. S. Surface Meshing Using Vertex Insertion. Proceedings, 5th International Meshing Roundtable (1996), 243­256. [33] Devillers, O. On Deletion in Delaunay Triangulation. Proc. 15th Ann. Sympos. Comput. Geom. (1999), 181­188. [34] Dey, T., Edelsbrunner, H., Guha, S., and Nekhayev, D. V. Topology Pre- serving Edge Contraction. Publications De L'institut Mathematique 66 (1999), 23­45. [35] Eck, M., and Hoppe, H. Automatic Reconstruction of B-Spline Surfaces of Arbitrary Topology Type. Proc. SIGGRAPH 1996 (1996), 325­334. [36] Edelsbrunner, H. Algorithms in Combinatorial Geometry. Springer-Verlag, 1987. [37] Edelsbrunner, H. Smooth Surfaces for Multiscale Shape Representation. Proc. Sympos. Found. Software Techn. Theoret. Comput. Sci. (1995), 391­412. [38] Edelsbrunner, H. The Union of Balls and its Dual Shape. Discrete Comput. Geom. 13 (1995), 415­440. [39] Edelsbrunner, H. Shape Reconstruction with Delaunay Complex. LNCS 1380, LATIN'98: Theoretical Informatics (1998), 119­132. [40] Edelsbrunner, H. Deformable Smooth Surface Design. Discrete Comput. Geom 21 (1999), 87­115. 158 [41] Edelsbrunner, H., Facello, M., and Liang, J. On the Definition and the Con- struction of Pockets in Macromolecules. Discrete Appl. Math. 88 (1998), 83­192. [42] Edelsbrunner, H., and Fu, P. Measuring Space Filling Diagrams and Voids. Tech. Rep. UIUC-BI-MB-94-01, Beckman Inst., Molecular Biophysics Group, Univ. Illinois at Urbana-Champaign, (1994). [43] Edelsbrunner, H., Liang, J., and Woodward, C. Anatomy of Protein Pock- ets and Cavities: Measurement of Binding Site Geometry and Implication for Ligand Design. Protein Science 7 (1998), 1884­1897. [44] Edelsbrunner, H., and M¨ucke, E. Three-dimensional Alpha Shapes. ACM Trans. Graphics 13 (1994), 43­72. [45] Edelsbrunner, H., and M¨ucke, E. P. Simulation of Simplicity: a Technique to Cope with Degeneracies in Geometric Algorithm. ACM Trans. Graphics 9 (1996), 67­104. [46] Edelsbrunner, H., and Seidel, R. Voronoi Diagrams and Arrangement. Discrete Comput. Geom. 1 (1986), 25­44. [47] Edelsbrunner, H., and Shah, N. Triangulating Topological Spaces. Internat. J. Comput. Geom. Appl. 7 (1997), 365­378. [48] Eisenberg, D., and McLachlan, A. Solvation Energy in Protein Folding and Binding. Science 319 (1986), 199­203. [49] Eschenauer, H., Kobelev, V., and Schumacher, A. Bubble Method of Topology and Shape Optimization of Structures. Structural Optimization 8 (1994), 42­51. [50] Foley, J. D., van Dam, A., Feiner, S. K., and Hughes, J. F. Computer Graphics: Principles and Practice in C. Addison Wesley, 1996. 159 [51] Giblin, P. Graphs, Surfaces, and Homology. Chapman and Hall, London, 1981. [52] Gregory, Arthur, State, A., Lin, M. C., Manocha, D., and Livingston, M. A. Feature-based Surface Decomposition for Correspondence and Morphing between Polyhedra. Proceedings of Computer Animation (1998), 64­71. [53] Haber, R., Bendsøe, M., and Jog, C. A New Approach to Variable-topology Shape Design Using a Constraint on the Perimeter. Structural Optimization 11 (1996), 1­12. [54] Hart, J. Computational Topology for Shape Modeling. Proc. Shape Modeling Inter- national '99 (1999), 36­45. [55] Hart, J., Durr, A., and Harsh, D. Critical Points of Polynomial Metaballs. Proc. Implicit Surfaces 98 (1998), 69­76. [56] Hershberger, J., Rauch, M., and Suri, S. Fully Dynamic 2-edge-connectivity in Planar Graphs. Proc. 3rd Scandinavian Worksh. Alg. Theory (1992), 233­244. [57] Hoppe, H., DeRose, T., Duchamp, T., McDonald, J., and Stuetzle, W. Surface Reconstruction from Unorganized Points. Proc. SIGGRAPH 1992 (1992), 71­ 78. [58] Hsu, W. M. Direct Manipulation of Free-Form Deformations. Proc. SIGGRAPH 1992 (1992), 177­184. [59] J.Hershberger, and Suri, S. Kinetic Connectivity of Rectangles . Proc. 15th Annu. ACM Sympos. Comput. Geom (1999), 237­246. [60] Kent, J. R., Carlson, W. E., and Parent, R. E. Shape Transformations for Polyhedral Objects. Proc. SIGGRAPH 1992 (1992). 160 [61] Keren, D., and Gotsman, C. Tight Fitting of Convex Polyhedral Shapes. Interna- tional Journal of Shape Modeling (1998), 111­126. [62] Kim, E., Baker, C., Dwyer, M., Murcho, M., Rao, B., Tung, R., and Navia, M. Crystal Structure of HIV-1 Protease in Complex with VX-478, a Potent and Orally Bioavailable Inhibitor of the Enzyme. J. Am. Chem. Soc. 117 (1995), 117. [63] Kuntz, L., Blaney, J., Oatley, S., Langridge, R., and Ferrin, T. A Geomet- ric Approach to Macromolecule-ligand Interactions. J. Mol. Biol. 161 (1982), 269­288. [64] Laug, P., and Borouchaki, H. Surface Grid Generation on Intersecting Spheres. Numerical Grid Generation in Computational Field Simulations (1998), 717­726. [65] Laug, P., and Borouchaki, H. Automatic Generation of Finite Element Meshes for Molecular Surfaces. European Congress on Computational Methods in Applied Sciences and Engineering, URL:http://www.imamod.ru/jour/conf/ECCOMAS 2000/pdf/70.pdf (2000). [66] Leach, A. R. Molecular Modeling Principles and Applications. Longman, 1996. [67] Lee, B., and Richards, F. M. The Interpretation of Protein Structures: Estimation of Static Accessibility. Journal of Molecular Biology 55 (1971), 379­400. [68] Lee, S., Chwa, K., Hahn, J., and Shin, S. Image Morphing Using Deformation Techniques . J. Visualization Comput. Animation 7 (1996), 3­23. [69] Lerios, Apostolos, Garfinkle, C., and Levoy, M. Feature-Based Volume Meta- morphosis. Proc. SIGGRAPH 1995 (1995), 449­456. [70] Lerios, A., Garfinkle, C. D., and Levoy, M. Feature-Based Volume Metamor- phosis. Proc. SIGGRAPH 1995 (1995), 449­456. 161 [71] Loop, C. Smooth Spline Surfaces over Irregular Meshes. Proc. SIGGRAPH 1994 (1994), 303­310. [72] Markosian, L., Cohen, J. M., Crulli, T., and Hughes, J. Skin: A Constructive Approach to Modeling Free-form Shapes. Proc. SIGGRAPH 1999 (1999), 393­400. [73] Milnor, J. Morse Theory. Princeton U. Press, 1963. [74] Moreton, H. P., and Se´quin, C. H. Functional Optimization for Fair Surface Design. Proc. SIGGRAPH 1992 (1992), 167­176. [75] Munkres, J. Elements of Algebraic Topology. Addison-Wesley, 1984. [76] Piegl, L., and Tiller, W. The NURBS Book. Springer-Verlag, 1995. [77] Provatas, N., Goldenfeld, N., and Dantzig, J. A. Adaptive Grid Methods in Solidification Microstructure Modeling. Phys. Rev. Lett. 80 (1998), 3308­3311. [78] Richards, F. M. Areas, Volumes, Packing and Protein Structure. Annual Review in Biophysics and Bioengineering 6 (1977), 151­176. [79] Rogers, D. F., and Adams, J. A. Mathematical Elements for Computer Graphics. McGraw-Hill, 1990. [80] Rossignac, J., and Borrel, P. Multi-resolution 3D Approximation for Rendering. Modeling in Computer Graphics (1993), 455­465. [81] Rozvany, G. Layout Theory for Grid-type Structures. Topology Design of Structures (1993), 251­272. [82] Santosa, F. A Level-set Approach for Inverse Problems Involving Obstacles. ESAIM: Control Optimization and Calculus of Variations 1 (1996), 17­33. [83] Seitz, S., and Dyer, C. View Morphing . Computer Graphics 507 (1996), 21­30. 162 [84] Sigmund, O., and Petersson, J. Numerical Instabilities in Topology Optimization: A Survey on Procedures Dealing with Checkerboards, Mesh-dependencies and Local Minima . Structural Optimization 16 (1998), 68­75. [85] Singh, K., and Parent, R. Joining Polyhedral Objects Using Implicitly Defined Surfaces. Visual Computing 7 (2001), 417­428. [86] Sokolowski, J., and Zochowski, A. On the Topological Derivative in Shape Op- timization. SIAM Journal on Control and Optimization 37 (1999), 1251­1272. [87] Stofer, E., and Lavery, R. Measuring the Geometry of DNA Grooves. Biopolymers 34 (1994), 337­346. [88] Tam, R., and Fournier, A. Image Interpolation Using Unions of Spheres . The Visual Computer 14 (1998), 401­413. [89] Taubin, G. An Improved Algorithm for Algebraic Curve and Surface Fitting. Pro- ceedings Fourth International Conference on Computer Vision (1993), 658­665. [90] Turk, G. Re-tiling Polygonal Surfaces. Proc. SIGGRAPH 1992 (1992), 55­64. [91] Turk, G., Dinh, H. Q., O'Brien, J., and Yngve, G. Implicit Surfaces That Interpolate. Proc. Shape Modeling International 2001 (2001), 62­71. [92] Turk, G., and O'Brien, J. Shape Transformation Using Variational Implicit Func- tions. Proc. SIGGRAPH 1999 (1999), 335­342. [93] Watt, A., and Watt, M. Advanced Animation and Rendering Techniques. Addison- Wesley, 1992. [94] Wolberg, G. Recent Advances in Image Morphing . Proc. Computer Graphics Internat. (1996), 64­71. 163 [95] Woo, M., Neider, J., and Davis, T. OpenGL Programming Guide. Addison Wesley, 1996. 164 Vita Holun Cheng was born in Hong Kong on July 30, 1972. He completed his Bachelor's degree in computer science at the Chinese University of Hong Kong and M.Phil. in computer science at the Hong Kong University of Science and Technology. He then began his doctoral work at the University of Illinois at Urbana-Champaign, completing a Ph.D. in computer science in May 2002. Holun's interest include computational geometry, parametric surfaces, computer graphics, visualization. He was a teaching assistant during his master degree. He received the Best Teaching Assistant Award in 1996 at HKUST. 165