A U N I O N OF S P H E R E S R E P R E S E N T A T I O N F O R 3D O B J E C T S By Vishwa Ranjan B. Tech. (Elec. Engg.) Indian Institute of Technology, Kanpur, India. M . Sc. (Comp. Sc.) University of Calgary, Calgary, Canada. A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S C O M P U T E R S C I E N C E We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A October 1996 © Vishwa Rarijan, 1996 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Computer Science The University of British Columbia 2366 Main Mall Vancouver, Canada V6T 1Z4 Abstract Reconstruction of an object from a set of points sampled from its boundary is an im-portant problem in graphics and vision. Several methods exist to compute and display surface (e.g., polygonal) and volumetric (e.g., polyhedral) models of objects from the boundary points. In order to display, transform, and compare objects, it is often convenient or necessary to use different representations of objects. Basic desired properties of representations of objects are efficiency of computation, storage, and display. Other important properties include stability (small changes in the data, such as noise or small distortions, cause small changes in the model), the ability to determine the similarities between two data sets, and the computation of simplified models. A survey of common representations of objects (e.g., surface, octrees, etc.) shows that some important properties are lacking in these representations. In this thesis we study a union of spheres representation (UoS) for the volume enclosed by an object's boundary. We present algorithms to generate stable union of spheres mod-els of objects from various sources of data, such as volumetric data (e.g., data from CT or MRI scanners), range data, and other existing models. The spheres can be simplified to obtain multi-scale models. We present a method to establish correspondence between two objects using their union of spheres models and use this to calculate distances be-tween objects, to register objects, and to interpolate between objects. This establishes a measure to study and compare such models. Examples with simple and complex ob-jects show how this measure corresponds closely to the intuitive human understanding of shape. ii Table of Contents Abstract ' ii List of Tables vii List of Figures viii Acknowledgement xv 1 Introduction 1 1.1 Sources of 3D data ' . 4 1.2 Terminology 5 1.2.1 An object 5 1.2.2 Representation versus model 5 1.2.3 Shape of an object . . . 6 1.2.4 Transformations 6 1.2.5 Sampling density or boundary point density 7 1.2.6 Error in a model 8 1.3 Desirable properties in representations 9 1.4 Union of Spheres (UoS) representation 10 1.4.1 Motivation 10 1.4.2 The UoS processing pipeline . 12 1.5 Contributions of the research 12 1.6 Organization 14 iii 2 Survey of 3D Object Representations 15 ' 2.1 Representations 15 2.1.1 Surface-based methods 15 2.1.2 Volume-based methods 18 2.2 Work related to UoS representation 22 2.2.1 Generation of union of spheres models . 23 2.2.2 Simplification of union of spheres models . 27 2.2.3 Computation of mass properties 27 2.2.4 Display of union of spheres models 30 2.3 Estimation of properties for various representations . 32 2.4 Summary 33 3 Generation of a Union of Spheres Model for Objects 34 3.1 General approach 34 3.2 UoS model from volumetric data 37 3.2.1 Algorithm for UoS model from volumetric data 38 3.2.2 Some properties of the generated models 39 3.3 UoS model from range data 42 3.3.1 Algorithm for UoS model from range data 43 3.4 UoS model from existing computer models , 46 t 3.5 UoS model from a polygon-mesh . 47 ' 3.6 Implementation and results . 47 3.7 Summary 57 4 Simplification of Union of Spheres Model 58 4.1 Triangles versus circles 59 4.2 Clustering of spheres 60 iv 4.2.1 Sphericity of a set of spheres 60 4.2.2 Algorithm description 62 4.2.3 Verification of spheres in the simplified model 64 4.2.4 Selection of sphericity 65 4.3 Implementation and results 65 4.4 Summary . 74 5 Approximation properties of the Union of Spheres Representation 75 5.1 General strategy ; 75 5.2 Relevant work . . 76 5.3 Approximating curves with arcs of circles 78 5.4 Approximating 2D objects with union of circles 81 5.5 Approximating 3D objects with union of spheres 86 5.6 Summary 90 6 Matching Objects Using Union of Spheres Models 91 6.1 Matching problem: a perspective 94 6.1.1 The problem . 94 6.2 General strategies 95 6.2.1 Image registration 96 6.2.2 Shape interpolation 97 6.3 Matching with UoS Models 99 6.3.1 Step 3: Model alignment 100 6.3.2 Step 4: Segmentation 101 6.3.3 Step 5: Match spheres in matched parts 104 6.4 Applications of matching 114 6.4.1 Object registration 115 v 6.4.2 Measuring distances between objects 121 6.4.3 Showing the stability of the UoS representation 126 6.4.4 Shape interpolation 132 6.5 Summary 145 7 Conclusion and Future Work 148 7.1 Summary of contributions 149 7.2 Future work 151 Bibliography 154 Notation 166 Glossary 168 vi List of Tables 2.1 Estimation of properties for various representations (1-poor, 5-very good). 33 3.2 Results of our algorithm on 3 D data 53 3.3 Results of our reconstruction algorithm on rabbit, knot, and hand data. . 53 4.4 Various definitions of sphericity and relationships. 61 4.5 Summary of results of the simplification algorithm on 2 D data 65 4.6 Summary of results of the simplification algorithm on 3 D data 66 6.7 Normalized distances between various 2 D objects. . . 124 6.8 Normalized distances between hands. . . 125 vii List of Figures 1.1 Modeling 2D objects using circles, (a) MRI slice (256 x 256) through human head. (b,c) Brain segmented from (a) modeled as union of 228 circles, (d) Boundary points for giraffe scanned from a book. (e,f) Giraffe modeled as a union of 370 circles 2 1.2 Modeling objects using union of spheres, (a) Boundary points for a rabbit obtained from a laser range-finder, (b) Rabbit modeled as a union of 1046 spheres, (c) Boundary points for a polygonal model of a hand, (d) Hand modeled as a union of 2759 spheres. 3 1.3 Example of boundary point data, (a) Points on the surface of a skull in computed tomography head data, and (b) the corresponding iso-surface. 4 1.4 Sampling density 5 7 1.5 Error e in the model R of an object 0 8 1.6 The UoS processing pipeline. . : 13 2.7 A n example of (a) blending function used for blending the sphere model and (b) how blending works. . 24 2.8 (a) Set of nine circles, (b) The associated power diagram, (c) A decom-position using the power diagram, (d) The union of circles 29 3.9 Empty circle property of Delaunay triangulation 35 viii 3.10 Steps of the algorithm to generate the UoC model, (a) A (32 x 32) size C T slice of human head (range 0-255). (b) Boundary points for a threshold of 70.3 (c) Delaunay triangulation of the boundary points, (d) UoT model, (e) UoC model, (f) UoC and UoT models superimposed, and the centers of the circles 40 3.11 (a) A circular cap of volume V\/2 on one of the boundary faces, (b) A configuration of boundary points so that the circles touch each other but the corresponding triangles do not 41 3.12 Classification of spheres as inside or outside, (a) Usual inside and outside cases, (b) A pathological case 44 3.13 Result on MRI brain data, (a) A (256 x 256) size M R I slice of human brain (range 0-255). (b) Boundary points for a threshold of 100.3 (1696 points). r (c) Delaunay triangulation of the boundary points, (d) UoT model (1703 triangles), (e) Union of triangles in the UoT model, (f) Union of circles in the UoC model (1703 circles) 49 3.14 Result on cow data, (a) Boundary points for cow (512 x 512 pixels), (b) Union of circles in the UoC model (2054 circles), (c) Histogram of radius distribution, (d-f) Use of histogram in incremental display of details in the object 50 3.15 Results of the generation algorithm on a (64 x 64 x 128) C T data for femur bone, (a) Boundary points (17384 points), (b) Surface model using marching cubes (34800 triangles), (c) UoS model (49902 spheres), (d) Histogram showing the radius distribution, (e-g) Use of histogram to control display of object (e) 8685 spheres, (f) 17531 spheres, and (g) 29357 spheres 51 ix 3.16 Results of the generation algorithm on a (64 x 64 x 64) P E T brain data. (a) Boundary points for a threshold of 140.3 (range 0-255, 17598 points). (b) Surface model using marching cubes (35216 triangles), (c) UoS model (58194 spheres), (d) UoS model with spheres of radius less than the inter-voxel distance removed (27917 spheres) 52 3.17 Results of the generation algorithm on 3D rabbit data, (a) Point sampling (22652 boundary points), (b) Reconstructed polygonal model [Turk and Levoy 94] (45300 triangles), (c) UoS model (57931 spheres) 54 3.18 Results of the generation algorithm on 3D knot data, (a) Polygonal model (3200 triangles), (b) Point sampling (13639 boundary points), (c) UoS model (38993 spheres) 55 3.19 Results of the generation algorithm on 3D hand data, (a) Polygonal model (1287 triangles), (b) Point sampling (15227 boundary points), (c) UoS model (31883 spheres) • 56 4.20 Figure showing the stability of the U o C representation 59 4.21 Definition of sphericity (a), (a) Large sphericity, (b) Small sphericity. . . 60 4.22 Generation of redundant circles during simplification, (a) Configuration of circles to be simplified, (b) A is simplified to A', C to C, and B to B'\ B' is contained in the union of A' and C. . 64 4.23 Results of the simplification algorithm on 2D brain data, (a) A n M R I slice through a human brain, (b) U o C model of the brain (1703 circles), (c-f) Simplified U o C models (dark) overlapped by the original U o C model (light) for a values of 0.97 (627 circles), 0.95 (470 circles), 0.90 (228 circles), and 0.85 (142 circles) respectively. 67 \ x 4.24 Results of the simplification algorithm on 2D "giraffe" object scanned from a book (512 x 512 pixels), (a) Boundary points (2112 points), (b) UoC model (2110 circles), (c-f) Simplified UoC models (dark) overlapped by the original UoC model (light) for an absolute error of 1.0 pixel (267 circles), 2.0 pixels (157 circles), 3.0 pixels (112 circles), and 4.0 pixels (89 circles) respectively 68 4.25 Results of the simplification algorithm on 3D rabbit data., (a) Polygonal model, (b) The UoS model (57931 spheres), (c-e) Simplified UoS models for a values of 0.97 (7766 spheres), 0.95 (5048 spheres), and 0.93 (3283 spheres), (f) Blended spheres display for the simplified model in (e). . . 69 4.26 Results of the simplification algorithm on 3D knot data, (a) Polygo-nal model, (b) The UoS model (38993 spheres), (c-e) Simplified UoS models for a values of 0.97 (2838 spheres), 0.95 (1061 spheres), and 0.93 (763 spheres), (f) Only, large spheres displayed using the.histogram (8063 spheres). • 70 4.27 Results of the simplification algorithm on 3D hand data, (a) Polygonal model, (b) The UoS model (31883 spheres), (c-e) Simplified UoS models for o~ values of 0.97 (3391 spheres), 0.95 (2319 spheres), and 0.93 (1491 . spheres), (f) Simplified UoS model for absolute tolerance of 3% of largest sphere radius of 1.70 units (1092 spheres) 71 4.28 Sphericity vs. number of spheres for rabbit data 72 4.29 Histogram of spheres for rabbit data before simplification. . . . . . . . . 72 4.30 Histogram of spheres for rabbit data after simplification 73 5.31 A n r-regular object. A n open circle of radius r can freely cover the inside and outside of the object without ever crossing the boundary. . . . . . . 77 xi 5.32 The error-bounding ellipse 79 5.33 Problems can arise in the proposed reconstruction method: (a) if the sam-pling constraint is violated, and (b) if the object is not r-regular 82 5.34 Modeling planar objects as union of triangles 83 5.35 Relation between UoT and UoC models 84 5.36 Figure to assist in the proof of Theorem 2. (a) q is outside the boundary of the object, (b) q is inside the boundary of the object 85 5.37 Error in a cluster 89 6.38 Simplified UoC models of some objects used in this chapter, (a-b) A cow (239 circles), (c-d) A giraffe (370 circles), (e-f) Fish (1) (45 circles), (g-h) Fish (2) (43 circles) 92 6.39 Simplified UoS models of some objects used in this chapter, (a) A rabbit (1046 spheres), (b) A head (262 spheres), (c) Hand (1) (1086 spheres). (d) Hand (2) (1087 spheres) 93 6.40 Segmentation of objects using UoC models. (a,c,e,g) Regular triangulation of circles in the UoC model of objects. (b,d,f,h) "Skeleton" extracted from the regular triangulation along with the circles 102 6.41 (a) Cylindrical portion of an object, (b) A corner, (c) A branching struc-ture. Solid lines show the boundaries, dotted circles the UoC model, and the dashed lines the "ridges". . . 106 6.42 Gradient in direction of neighbor of circle 108 6.43 Feature of a circle . 109 6.44 Feature distances I l l 6.45 (a) Bipartite graph after distances are assigned, (b) Bipartite graph after matching. . ". . . . . . . 113 xii 6.46 Registration of two brain slices using UoC. (a) A 256 x 256 M R I slice through human head, (b) Same image, rotated, translated, scaled, and subsampled at 64 x 64 resolution, (c) UoC model of brain in (a), (d) UoC model of brain in (b). (e) The two UoC models together showing the extent of mismatch, (f) The two UoC models after registration. 118 6.47 Example of registration in 3D. (a) Brain at 64 x 64 x 64 resolution (27917 spheres), (b) Transformed brain at 16 x 16 x 16 resolution (2774 spheres), (c) Brain in (a) after simplification (746 spheres), (d) Brain in (b) after simplification (734 spheres), (e) Data sets (a) and (b) before registration. (f) Data sets (a) and (b) after registration. - 120 6.48 Hand (3) and hand (4): models to be classified 125 6.49 Stability of unsimplified. UoC representation to boundary noise for 2D giraffe data 128 6.50 Stability of simplified UoC representation to boundary noise for 2D giraffe data 128 6.51 Stability of simplified UoC representation to shape-preserving transforma-tions for 2D giraffe data 129 6.52 Stability of simplified UoC representation to some low frequency distor-tions for 2D giraffe data 130 6.53 Stability of UoS representation to boundary noise for 3D hand data. . . 131 6.54 Interpolation from a rectangle to an MRI slice of brain for ws = 0.2, wp = 1.0, and Wf = 0.1 134 6.55 Interpolation from fish (1) to fish (2). The following weights were used in the distance measure: ws = 0.4, wp = 1.0, and Wf = 0.2 136 xii i ( 6.56 Matching, circles in the two fishes using features alone, (a-b), (c-d), (e'-f), and (g-h) show four pairs of matched circles. Notice a "poor" match in (g-h) pair. . 137 6.57 Interpolation from a calf to a cow (without segmentation). Weights used: ws = 0.2, wp = 1.0, and Wf = 0.3 ." 138 . 6.58 Interpolation from a calf to a cow (using automatic segmentation). Dif-ferent weights were used to match different components. . ' 139 6.59 Interpolation from a cow to a giraffe (using automatic segmentation). Dif-ferent weights were used to match different components 140 6.60 Interpolation from hand (1) to hand (3) 142 6.61 Interpolation from hand (1) to hand (2) 143 6.62 Example of matching in 3D. Four pairs (a-b, c-d, e-f, g-h) of matched spheres along with their features for the hand data 144 6.63 Interpolation from a rabbit to a human head 146 xiv Acknowledgement Thanks to Alain; this thesis belongs to him. Thanks to the committee members for their help, expertise, and encouragement. Thanks to other faculty members for the courses I enjoyed. Thanks to the Department of Computer Science, U B C , the Government of British Columbia, and the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support. Thanks to Imager people, past and present, for innumerable things. Thanks to David Watson, Bernd Gaertner, Robert Kennedy, The Geometry Center, and Edelsbrtinner's group at the University of Illinois at Urbana-Champaign, for sharing their code. Thanks to all colleagues and friends who made my stay in Vancouver pleasant. Thanks to Mom, Dad, Lok, Desh, Ritu, Apratim. - • Thanks to my "family" members in Vancouver. Thanks to Xiaoli for sometimes being a good friend in need. Thanks to everyone else who deserves to be thanked! To all of thee, I am grateful, forever. xv Chapter 1 Introduction A concise statement of the problem addressed in this thesis is Given a set of points V known to be on the boundary of an object, what is a good model for the object? We will specifically study the modeling of objects as a set union of overlapping spheres. The union of spheres approximates the volume of the object and the surfaces of spheres not inside any other sphere approximate the surface of the object. The algorithms are described and implemented for 2D objects as well. Therefore, it is also a study of the union of circles model for planar objects. As an illustration of the technique in two-dimensions, Figure. 1.1 shows two examples of objects modeled as a union of eircles. Figure 1.2 shows two examples of three-dimensional objects modeled as a union of spheres. If no information is available on how the boundary points were obtained, the problem is to model an arbitrary cloud of points [C.3.1-Boissonnat 84, C.3.2-Edelsbrunner and Mucke 92, C.1.3-Hoppe et al. 92]. The problem we will consider is a restricted form of this problem because we assume that we know something about how the boundary points were obtained. We use this information to generate a union of spheres model of objects from several types of boundary point data by deriving a test for inclusion or exclusion of a sphere in the defined object. Sometimes the information about an inclusion-exclusion test is available explicitly, as in the case of volumetric data where it is common to use a threshold value to define points on the boundary of the object, and the corresponding iso-surface (Figures 1.3a and 1.3b) [C.1.6-Lorensen and Cline 87]. Sometimes the required 1 Chapter 1. Introduction 2 Figure 1.1: Modeling 2D objects using circles, (a) MRI slice (256 x 256) through human head. (b,c) Brain segmented from (a) modeled as union of 228 circles, (d) Boundary points for giraffe scanned from a book. (e,f) Giraffe modeled as a union of 370 circles. Chapter 1. Introduction 3 Figure 1.2: Modeling objects using union of spheres, (a) Boundary points for a rabbit obtained from a laser range-finder, (b) Rabbit modeled as a union of 1046 spheres, (c) Boundary points for a polygonal model of a hand, (d) Hand modeled as a union of 2759 spheres. Chapter 1. Introduction 4 Figure 1.3: Example of boundary point data, (a) Points on the surface of a skull in computed tomography head data, and (b) the corresponding iso-surface. information is available implicitly, as in the case of data provided by 3D laser range-finders and ray-tracers where an inclusion-exclusion test can be (as we shall see later) derived from the imaging geometry. 1.1 Sources of 3 D data Data in the form of discrete sample points is now common due to the advancement of technology in 3D imaging. A number of techniques exist for acquiring 3D bio-medical images such as computed tomography [B.l-Bates et al. 83], magnetic resonance imaging [B.6-Riederer 89], ultrasound imaging [B.7-Wells 88], confocal microscopy [B.5-Paddock 91], and positron and photon emission tomography [B.2-Jones 90, B.3-Knoll 83]. Differ-ent techniques are good for imaging different types of tissues and they have different reso-lutions. A good survey of these techniques is given by Stytzi et al. [A.9-Stytzi et al. 91]. These 3D imaging techniques generally yield scalar data on a regular grid. Typically, Chapter 1. Introduction - 5 a number of discrete 2D slices (or cross-sectional images) of objects are available. This slice-by-slice digital data or the "3D image" can be conveniently represented as a 3D array of values. Although bio-medical images are the main source of 3D volumetric data, the tech-niques presented here apply to astrophysical, atmospherical, computational fluid dynam-ics (CFD), geological, and oceanographical data [A.6-McCormick et al. 87]. Data from 3D surface scanners such as laser range-finders is also amenable to the techniques presented. Boundary point data similar to range, data can be generated by point-sampling the surface of an object, for instance by ray-casting any computer model of the object from different directions; therefore, the techniques presented here apply to synthetic data as well. 1.2 Terminology 1.2.1 A n object An object O is a closed, bounded subset of three-dimensional Euclidean space E3 (subset of E2 in two-dimensional space). The boundary of an object B(0) is a closed surface. Given any point in space, it is either outside the object, on the boundary of the object, or inside the object. 1.2.2 Representation versus model A representation IZ means a computational method and a data-structure to describe an object. A model is a specific instance of a representation similar to the difference between a class (representation) versus a particular instance (model). Examples of representations are boundary point representations, surface-mesh representations, union of tetrahedra representations, and union of spheres representations. Chapter 1. Introduction 6 1.2.3 Shape of an object While the notion of shape is intuitive, its definition is elusive. Traditionally, the shape of an object is defined in terms of its outward appearance or "look." We are concerned more with the structure of the object, both inside and outside. We borrow the defini-tion of shape from [H.6-Lord and Wilson 84] as follows: "It is the characteristic way in which an object occupies space." The shape of an object in this sense is independent of transformations such as translation, rotation,1 and uniform scaling. Although we do not deal directly with the shape of an object in this thesis, we discuss shape comparison, registration, and interpolation. 1.2.4 Transformations In this thesis we will distinguish between shape-preserving and non-shape-preserving transformations. Shape-preserving (Tsp): These include rigid body transformations (rotation and translation), and uniform scaling. These transformations are a subset of affine transfor-mations and are angle-preserving. Non-shape-preserving (TNSP)'- These include other affine transformations such as non-uniform scaling and shearing. Affine transformations are characterized by a constant Jacobian and they preserve lines and their parallelism. Other commonly used non-shape-preserving transformations in graphics include perspective (which preserves lines but not parallelism) and instances of differentiable homeomorphism (meaning that the Jacobian of the transformation is well defined) such as "free-form deformations" controlled by parametric polynomial meshes [G.15-Sederberg and Parry 86]. In general, these are non-linear in the sense that lines can be transformed to curves. 1 T h e human perception of shape is in fact often sensitive to rotation. Chapter 1. Introduction 7 1.2.5 Sampling density or boundary point density The reconstruction of an object from boundary points is sensitive to the resolution at which the surface was sampled. We measure the sampling density in terms of the largest distance 8 between any point q G B(Q) and the sample point p G V nearest to it (Figure 1.4). The sampling density 8 is defined as the radius of the smallest circle (sphere in object boundary * sampled boundary points Figure 1.4: Sampling density 8. 3D) with the property that a circle of radius 8 centered anywhere on the boundary of the object is guaranteed to contain at least one sampled boundary point. Note that according to this definition of 8 the actual spatial density of sample points is at least ^ in 2D (at least ^ in 3D). Other researchers have used a similar definition to express the density of sample points [C.1.3-Hoppe et al. 92, D.5.2-Brandt 92, D.5.1-Attali and Montanvert 94]. Chapter 1. Introduction 8 1.2.6 Error in a model To define error between an object and its model we choose the Hausdorff distance which is the maximum distance between any point on the surface of the object and the nearest point on the surface of its model and vice-versa [H.ll-Preparata and Shamos 85]. Gen-erally, this definition of error is preferred because it implies a bound on some other error measures as well (e.g., difference in volumes). B(R) B(O) P (^3) ° b J e c t 0 C - - . - : M o d e I R Figure 1.5: Error e in the model R of ah object O. Let the boundary of an object O be noted as B(0). Let 7Z be a model of O and B(7Z) be the boundary of 71. Assume that the object and its model reside in the same 3D Euclidean space. Let dE(p,q) denote the Euclidean distance between two points in this space. Then, an error bound of e between O and 7t means that (Figure 1.5): (A) for any point p G B(0), 3 q £ B(7l) s.t. dE(p, q) < e, and (B) for any point q <E B(7l), 3 p e B(O) s.t. dE(p,q) < e. Chapter 1. Introduction 9 We will refer to these as Conditions A and B respectively. The smallest value" of e that satisfies both these conditions is the Hausdorff distance between the boundaries of O and TZ and is defined to be the error in the model. Conditions A and B are necessary and sufficient to show an error upper bound of e between an object and its model. 1.3 Desirable properties in representations We want a model that possesses some of the intrinsic properties of the object modeled, in addition to basic ones such as ease in generation, manipulation, storage and display. "Visualization" depends on one of these intrinsic properties, namely, the look. 1 In gen-eral, intrinsic properties are properties independent of the resolution of the data, of the imaging geometry, or even of the imaging method itself. Such properties would help in comparing objects, establishing whether two data sets are likely to be from the same object and determining the possible transformations between them if they are. Another property that is useful to us is the ability to simplify the model while keeping, as much as possible, the properties of the object modeled. This helps speed up display, interactive manipulation, and the extraction of models independent of the original level of details available in the data. To achieve these goals, one important property that a representation should have is stability: non-intrinsic changes in the data should not result in significant changes in the model. Non-intrinsic changes include noise and small distortions in the data, changes in the resolution in the case of discrete data, and basic geometric transformations such as uniform scaling, translation and rotation. If we consider two models Ri and R2 of the same object defined by two point sets P i and P 2 respectively, then we call the representation stable if there is a finite positive constant K, such that: d(RuR2) < Kd(PuP2) Chapter 1. Introduction 10 where, d(Ri,R2) measures the distance between two models of the object and d(P\,P2) measures the distance between the corresponding two point sets. In other words, the changes in distance between the primitives of the model are commensurate with the changes in the input data. Intuitively, this means that the primitives in the model of an object (derived from different point sets) depend on the shape of the object and not on the sampling artifacts. Note that we assume that we can compute distance between the models. 1.4 Union of Spheres (UoS) representation 1.4.1 Motivation To give an intuitive example of the kind of model we are looking for, assume that we have data sampled from a human head. The head can in a first approximation be modeled by a sphere. The surface area and the volume of the sphere wi l l provide us rough but useful estimates of the surface area and the volume of the head. At the same time the position and radius of the sphere will give us an idea of the translation and scaling to apply to transform the head from its canonical position. If we fit an ellipsoid of revolution instead, then the two additional degrees of freedom available in the ellipsoid might allow us to obtain the parameters of the rotations to orient the head in addition to position it. Under the right circumstances, many properties such as surface area, volume, even topology or connectivity can be obtained from primitives that are not intrinsically related to the object they model. For instance, surface area can be computed from a triangulation of a surface even if no information is available on the adjacencies among the triangles. However, if we would like the model to be: • complete, in the sense that both the boundary and the inside of the object are modeled, Chapter 1. Introduction 11 • accurate, that is, the measurable properties and the "look" of the model should be close to that of the object, • stable, as defined above, • multi-scale, so that it is possible to work with it at a variable level of detail, • displayable, so that it can be rendered at a reasonable speed on machines with and without special graphics hardware, • transformable, that is, the object can be geometrically transformed (e.g, rotated, scaled, translated) by transforming the model, • manipulable, that is, it permits control over the object by providing control over its primitives, e.g., altering shape of the object, Boolean combinations of parts, etc., then the choices of primitives to use are relatively few. We need volume primitives to satisfy completeness, simple volume primitives to reduce the search space when fitting to the data and facilitate visualization, and more than one primitive to hope for accuracy. In this thesis we study a representation of objects using one of the simplest volume primitives, namely, spheres. The idea is to explore a union of spheres representation (UoS) in detail, examine its strengths and weaknesses, and consider its possible applica-tions. In the UoS representation, the spheres are derived from a Delaunay tetrahedral-ization of boundary points; this choice of spheres allows us to do efficient computation and simplification of union of spheres models of objects.' We describe algorithms to ex-tract UoS models of objects from the given boundary points, to simplify and display these models, to compare and interpolate objects using these models, and to retrieve transformations between related data sets that use these models. Chapter 1. Introduction 12 1.4.2 The UoS processing pipeline Figure 1.6 shows the sequence of processing steps in our approach. In the case of a single object these are: 1. Obtain or generate the boundary points V on the surface of the object. 2. Generate the UoS model from these boundary points (TI). 3. If needed, simplify the UoS model (new TV). . 4. Use 1Z for visualization and other applications. 1.5 Contributions of the research The results presented in the body of the thesis are: • We show how we can generate a UoS model for a 3D object from a set of points on its boundary. • We provide an algorithm for a stable simplification of the UoS model. • We determine the approximation properties of the UoS representation. • We demonstrate how we can compare two objects by matching the spheres in their UoS models. • We present a method to extract a transformation between two objects, again, by matching the spheres in the UoS models of two objects. We then apply this tech-nique to object registration and shape interpolation. • We discuss various operations (e.g., display, computing mass properties) and how easily or not they can be computed with the UoS model. Chapter 1. Introduction 13 REAL OBJECTS Surface scan e.g., range-finder data * * * \ 1 » \ N < > « CCG MODELS A Ray casting ) Sampled boundary points Calculate the UoS model Simplify REAL OBJECTS Volumetric scan e.g., CT or MRI (isovalue) Simplified UoS model Figure 1.6: The UoS processing pipeline. / Chapter 1. Introduction 14 1.6 Organization The thesis is organized as follows. Chapter 2 presents a survey of different representations for three-dimensional objects. It also describes previous work related to union of circles (in 2D) and union of spheres (in 3D) representations of objects. Chapter 3 describes how we can generate the UoS models of objects from available 3D data, such as C T scans, M R I scans, range data, and other computer graphics models. Chapter 4 presents an algorithm to simplify the generated UoS model. Chapter 5 discusses the approximation properties of the UoC and UoS representations. In particular, it points out the conditions under which the error in the UoS model obtained from the described reconstruction method is bounded. Chapter 6 explains how we can match and compare two objects using their UoS models. It illustrates the usefulness of matching by presenting algorithms (based on matching) to compute distances between two objects, to register two objects, and to interpolate between two objects. Finally, Chapter 7 summarizes the contributions and outlines directions for future work. This order of presentation tries to be consistent with the usual discussion of the classic problem of reconstruction of a signal from a set of samples (namely, discrete samples —r continuous version of the signal using a particular reconstruction method —>• study the approximation properties of the proposed reconstruction method). Readers who are curious about the conditions under which the presented UoS reconstruction method results in a bounded error, could read Chapter 5 before reading Chapters 3 and 4. Chapter 2 Survey of 3D Object Representations In this chapter we present a survey of different representations for 3D objects. The rep-resentations used can be divided into three categories: representations that construct the surface explicitly, representations that use volume primitives whose surfaces (or a subset of the surface) define the object's surface, and representations that do not construct the surface explicitly, but rely on direct volume rendering and (human) visualization to "see" the boundary. The last category of representations will not be discussed here, because an explicit computation of the surface is required to meet many of the criteria given in Chapter 1. Of course if visualization is the only goal, these representations can be quite effective and powerful [G.6-Drebin et al. 88, G.14-Sabella 88]. We focus instead on surface-based and volume-based representations. 2.1 Representations 2.1.1 Surface-based methods Surface methods model the boundary of an object by fitting surface primitives to point, data. Usually the primary concern is fast surface display from different viewpoints. The complexity of the method depends on the surface primitives used. One advantage of a surface model is that it is easy to display, especially if it can be converted into polygons, because many graphics workstations support polygon display in hardware. The disadvantages are that surface models are usually piecewise, and lack information about 15 Chapter 2. Survey of 3D Object Representations 16 the global properties of the object such as total volume, connectivity, and inclusion, (a) Marching cubes In medical imaging, the most commonly used method of this type is the so called march-ing cubes technique, which fits triangles to the boundary points in every boundary cell [C.1.6-Lorensen and Cline 87] of the given volumetric data. The algorithm is sim-ilar to the one presented by Wyvil l et al. in the context of modeling "soft" objects [D.4.3-Wyvill et al. 86]. The algorithm assumes that the voxels have been classified as inside or outside the object and that the surface of the object passes through an edge of a cell if the voxels defining that edge are of different types. Wyvil l et al. noted that the 256 different configurations that are possible for a cell (each of the 8 defining voxels can be inside or outside) can be reduced to 14 basic configurations—others are permutations of these. The idea then is to first classify each cell into one of the fourteen different configurations possible (based on the inside-outside configuration of the defining voxels), and then to use a table look-up to fit triangles in every cell to boundary points obtained by interpolating between all neighboring inside-outside pairs of voxels. One advantage of the algorithm is that the model is quick to compute. The method's success is due to the fact that interactivity (that is, the ability to compute a new image at or near video rates) plays a crucial role in visualization. One problem with the algorithm is that some configurations of cells are ambiguous, so the model that is created may not be a manifold. There may be a large number of triangles in the model and their number and distribution are highly resolution and sampling dependent. A l l of these mean that the triangles cannot be used to characterize the object or compare it to other objects. Chapter 2. Survey of 3D Object Representations 17 (b) Hoppe's method A similar approach, based on different input data, fits a surface to the set of bound-ary points considered as an arbitrary cloud of points ignoring whatever neighborhood information could be kept from the original volumetric data. A n effective algorithm to compute such surfaces has been developed by Hoppe et al [C.1.3-Hoppe et al. 92]. The fact that they do not assume any structure in the set of points is both a strength • ( and a weakness of the method for our purposes, since useful information about in-side/outside classification may be lost. The method also shares with the marching cubes the problem of possibly generating a large number of triangles—in fact they use marching cubes as one step in the overall method. Methods to simplify this kind of mesh exist ( [C. l . l l -Turk 92, C.1.9-Schroeder et al. 92, C.1.4-Hoppe et al. 93]), but they can add problems of their own (in particular, stability). A recent paper by Eck et al. shows how to conveniently perform a multiresolution analysis on arbitrary meshes of this type [C. l . l -Eck et al. 95]. (c) Surfaces from contours One other approach for generating the surface-mesh of an object from volumetric data is to derive a surface from contours [C.1.2-Fuchs et al. 77, C.l.lO-Meyers et al. 91]. In this approach, contours on different 2D slices of the volumetric data are first gen-erated and then the contours on the neighboring slices are connected. This reduces the 3D surface reconstruction problem to a 2D contour-tracing problem followed by a "stitching together" phase. Higher order surfaces have also been fitted using this ap-proach [C.1.5-Lin et al. 89]. The method works well when slices have a single contour each, but difficulties arise when multiple contours appear in a slice (for example, due to branching in the object). Several researchers have addressed this problem in detail "A Chapter 2. Survey of 3D Object Representations 18 [C.l.lO-Meyers et al. 91]. (d) Surface models from range data A large number of techniques exist in the graphics and vision literature to reconstruct surface models of objects from range data. Surface methods first deduce a surface model (generally polygonal) for each of the views individually, and then merge the informa-tion from all the views by identifying and taking care of the overlapping portions of the surfaces in the different scans. Some of the successful methods can be found in [F.ll-Schmitt 91, E.6-Soucy and Laurendeau 92, E.7-Turk and Levoy 94]. Another re-cent technique for creation of detailed surface models from range data uses volumetric data as an intermediary step for reconstruction [E.3-Curless and Levoy 96]. 2.1.2 Volume-based methods I. Voxel-based methods Voxel methods store directly or indirectly some or all of the voxels (i.e., their position coordinates and values) in discrete volumetric data. Three methods will be discussed here. (a) Octrees Octrees are hierarchical data-structures used to model volumetric objects [C.2.1-Meagher 82], assuming a binary classification of each voxel as being inside or outside the object. A bounding box around the object forms the root node for the tree. Beginning from the root node, each node is recursively subdivided into eight sub-nodes (cubes or octants) if it is mixed (not entirely inside or entirely outside the object). The subdivision process stops when either there are no more nodes of mixed type in the tree or when a predefined Chapter 2. Survey of 3D Object Representations !9 finest level of subdivision has been reached. Octree models are generally more compact than the voxel model, are hierarchical, and there exist elegant algorithms to generate them and apply some operations (e.g. Boolean operations) on them. For display, voxels can be displayed directly (using volume rendering), or triangles can be generated within all boundary cells (using a method similar to marching cubes). One major disadvantage with the octree representation is that'it is highly dependent on the coordinate system. A small change in value, especially near the boundary between nodes deep in the hierarchy, will result in very different octrees, and this hampers the use of octrees in problems such as characterization of shapes and comparison of data sets. In 2D, a similar technique using box subdivision provides a quadtree model for the area of a 2D object. (b) Shells An interesting model, with characteristics of voxel fuzzy classification and boundary representation, uses shells [C.2.5-Udupa and Odhner 93]. A shell is a set of voxels in the vicinity of the boundary, with "fuzzy" membership values associated with them. Shells can be used directly for rendering, but Udupa and Odhner show that they can also be used for various operations such as volume computations and measurements on surfaces. As useful as they seem to be, especially because of the speed of rendering obtained, shells do not meet some of our criteria, mostly because they are the result of a classification scheme of voxels that is dependent on the resolution and geometry of the original imaging system. (c) Wavelets Wavelets have been developed in the last few years as powerful analytic and numeric tools for a variety of tasks [H.2-Chui 92]. They allow the modeling of signals in a hierarchical, Chapter 2. Survey of 3D Object Representations 20 multiscale basis, sometimes with a high level of compression. They have been used in the context of volume data representation by Muraki [C.2.3-Muraki 93]. Even though the intrinsically hierarchical nature of wavelet transforms allows several levels of modeling, it is not clear that each level is in fact useful either visually or for the purpose of the analysis of global properties. In particular, the fact that the wavelet transforms used are not invariant under translation or rotation, makes their use for the comparisons of objects and the extraction of transformations problematic. II. Fitting volume primitives The problem of deriving a geometric model from discrete boundary data of an object for classification or comparison is a well-studied problem in computer vision. A common technique in vision is to fit a single volume primitive (such as a deformable superquadric [C.3.7-Terzopoulos and Metaxas 91]) to the boundary points by posing the problem as an optimization problem. A fairly complex shape can be defined with relatively few parameters; the hope is that the method is stable, making it useful for shape comparison. Because the inside-outside property is well-defined for these primitives, they can be used in ray-casting or ray-tracing algorithms for display. This approach is attractive for use in visualization, but its main drawback in our context is that a single primitive cannot model effectively most objects of interest. Some a priori decision has to be made about the expected complexity of the object(s) to be modeled. We will see that our approach is a variation with simpler primitives, but using a number of primitives that is not pre-determined. III. Assemblage of volume primitives Various methods have been used to model an object with a collection of volume primitives. We will mention three different ones. Chapter 2. Survey of 3D Object Representations 21 (a) BSP trees Binary Space Partitioning (BSP) trees are structures made of the union and intersection of half-spaces [C.3.5-Naylor 92b]. They are widely used not only as data structures for efficient spatial query and navigation, butalso as modeling tools [C.3.4-Naylor 92a]. For boundary point data, the strategy is to generate the BSP trees by fitting planes to the boundary points. This needs good heuristics to be effective. Once this is done, the structure produced has the advantages of BSP trees (quick Boolean operations, good visualization). In addition, the simplification of the model is relatively easy. Stability is not guaranteed, being mainly dependent on the methods used to produce the original planes. (b) Delaunay tetrahedra Some methods such as that of Boissonnat [C.3.1-Boissonnat 84] use the Delaunay tetra-hedralization of points followed by an inside-outside classification of the tetrahedra. The faces where the inside and outside tetrahedra meet define the boundary of the object. Of-ten, it is not clear which tetrahedra should be kept and which ones should be discarded; heuristics are generally used. From the tetrahedral model it is possible to obtain the polygonal mesh for the surface of the object which can then be simplified by any of the mesh-simplification methods mentioned previously. In general, the number of tetrahedra can be as large as 0(n 2 ) , where n is the number of boundary points, hence the model, unless simplified, may not be compact enough to store, or easy to manipulate. Delau-nay tetrahedralization is not stable to boundary noise, therefore it may not be possible to use the tetrahedra in shape comparison applications. Algorithms that simplify the tetrahedral model exist but have found limited use [D.5.1-Attali and Montanvert 94]. Chapter 2. Survey of 3D Object Representations 22 (c) Alpha-shapes One innovative computational geometry technique to determine "shapes" from an ar-bitrary cloud of points is alpha-shapes [C.3.2-Edelsbrunner et al: 94]. The idea is to assume a spherical eraser of radius a and move it everywhere in 3-space where it can go without including any of the original points. The part of 3-space that has not been erased is the alpha-shape, and models the object at the scale a. At one extreme, for an infinite a, the result is the convex hull of all the points. At the other extreme, for an a of 0, the re-sult is the set of points itself. In-between, the shape is normally bounded by arcs of circles and spherical segments, but they are normally replaced by straight line and plane seg-ments. In the algorithms designed and implemented by [C.3.2-Edelsbrunner et al. 94], the model is derived efficiently using Delaunay tetrahedralization. The model obtained has the great advantage of being at a controllable scale. One disadvantage is that the same scale applies to all features. This shortcoming was removed by the definition of weighted alpha-shapes [C,3.3-Edelsbrunner 92], but the problem of automatically select-ing different weights for different features still remains. Moreover, there can be "hanging" faces and edges in the model (i.e., the surface may not be a manifold). The most serious drawback for our purpose is the fact that the intuitive notion of shape is not captured by the model, and the volume primitives used, tetrahedra, are not sufficiently stable to be directly usable. 2.2 Work related to UoS representation This section deals with previous work related to the UoS representation. We will discuss methods of generating union of spheres models of objects, computation of mass properties from the union of spheres model, and the display of objects that are modeled as a union of spheres. Chapter 2. Survey of 3D Object Representations 23 2.2.1 Generation of union of spheres models (a) Union of spheres Union of spheres models have been used extensively in visualizing molecules. Each atom can be modeled by an individual sphere; atoms of different elements may have different sizes corresponding to their Van der Waal's radii or other physical properties. Some systems use cylinders to denote the bonds between atoms whereas other systems are based on overlapping spheres only [D.2.8-Knowlton and Cherry 77, D.2.9-Max 79, D.2.2-Connolly 83a]. The use of a union of spheres for modeling objects for com-puter graphics and vision applicafions was suggested by Badler and Bajcsy in 1978 [A.l-Badler and Bajcsy 78]. In 1979, O'Rourke and Badler presented an algorithm to de-compose an object into spheres from discrete contour data [D.1.5-0'Rourke and Badler 79] Their algorithm is unstable and there is no guarantee that all of the object volume will be n. included in the model. In another paper, Badler et al. proposed using a spherical decom-position of a human body for animation and visualizing motion [D.l . l-Badler et al. 79]. Mohr described a technique based on the distance transform to obtain a model using non-overlapping, tangential spheres from which he derived a network (graph) describ-ing the object for shape comparison [D.1.4-Mohr and Bajcsy 83]. The union of spheres and ellipsoids representations have also been suggested for modeling objects in com-puter graphics by Herbison-Evans [G.12-Herbison-Evans 82]. Hubbard derived a union of spheres model ("sphere-tree") from an octree model for the purposes of efficient colli-sion detection between objects [D.1.3-Hubbard 95]. (b) Union of blended spheres or "blobbies" Blinn showed how blended spheres can be used to model trees and human beings [D.4.1-Blinn 81] and Wyvil l et al. used Blinn's technique to model "soft" objects [D.4.3-Wyvill Chapter 2. Survey of 3D Object Representations 24 et al. 86]. Blended spheres are sometimes called "blobbies". Blobbies model a 3D surface as an isosurface of a scalar field which is produced by a number of field generating primitives or "atoms." In other words, blobbies define implicit surfaces using atoms that Figure 2.7: A n example of (a) blending function used for blending the sphere model and (b) how blending works. have an associated blending function f(r) (a monotonically decreasing field, symmetric in space for spherical atoms) surrounding them (Figure 2.7a). These atoms interact with each other through these blending functions. The sum of the field functions, £ ) i fifr) — T, defines a "smooth" implicit surface. Whereas Blinn used a Gaussian of infinite extent as the blending function, Wyvil l et al. modified the atoms in blobbies to use blending functions of finite support for reasons of efficiency. Figure 2.7b shows how this works. Each atom has five degrees of freedom: namely, three for the center, one for the radius of isolation riso, (i.e., radius of the atom in the absence of other atoms) and one for the radius of influence r,„yj. The use of blobby models in computer graphics is becoming increasingly common; in the past few years blended spheres have been used to model, render, and animate a wide Chapter 2. Survey of 3D Object Representations 25 variety of objects, e.g., smoke, killer whales, etc. Muraki [E.4-Muraki 91] described a method to derive "blobby" models from range data. In some sense, his method tries to fit n overlapping spheres which when blended together best approximate the surface of the 3D object. One difference with our approach is that whereas Muraki's reconstruction technique works top-down (it starts with a single primitive and recursively splits it to fit the data), our approach is bottom-up (it starts with numerous spheres, and clusters them into larger approximating spheres). (c) Skeletons Some union of spheres models are based on the notion of skeleton or medial-surfaces (axis in 2D) of an object. These are the locus of the centers of all optimal spheres (circles in 2D) inside an object. An optimal sphere is a sphere inside the object that cannot be contained within any other sphere inside the object; it is defined by at least three points on the boundary of the object. The skeleton of an object can be described in a space that is a dimension less than that of the object itself, thereby hopefully facilitating comparison of objects. Several methods exist to compute discrete (hence approximate) version of skeletons from voxel data or from boundary points. Some methods are based on thinning of the binary volumetric voxel model of the object [D.5.7-Tsao and Fu 84], other methods are based on the computation of the Voronoi diagram of the boundary points [D.5.4-Goldak et al. 91, D.5.2-Brandt 92], while some are based on ridge-tracking in the Euclidean distance field for the object [D.1.4-Mohr and Bajcsy 83]. Attali studied in detail the skeletal properties of the union of spheres model of objects [D.5.1-Attali and Montanvert 94]. The starting point of their method is similar to ours. One problem with skeletal representations for our purposes is that they are unstable with respect to boundary noise [H.7-Marr 82]. Several researchers address the issue of how to compute stable skeletons (mostly in 2D) for objects. These methods can be broadly classified into Chapter 2. Survey of 3D Object Representations 26 two categories: (i) methods that calculate the (unstable) skeleton from the data and then apply a simplification procedure which stabilizes it, and (ii) methods which pre-process the data using a blurring operator which hopefully takes care of the boundary noise, after which they compute the skeleton. One distinction between the skeletal model and the UoS model studied in this thesis is that we are not concerned with the connectivity or other topological properties of the skeleton. The union of spheres is treated as a union of simple volume primitives and not as a graph denoting the connectivity between the spheres. As we will see in later chapters, this actually reduces the complexity of the generation and simplification of the UoS model and facilitates comparison of objects using the model. (d) Cores The core representation generalizes the notion of medial surfaces (axes) for gray-scale images [D.5.6-Morse et al. 94]. The model is obtained by first blurring the voxel data followed by application of "boundariness" and "medialness" operators to the resulting data. The medial points obtained by ridge-tracking form the "core" of the object at that level of blurring. For every core point, both the location of the core point and its distance from the estimated boundary at that level are stored. This process is repeated for different (but exponentially related) blurring factors. Some of the cores at lower levels are reinforced by cores at a higher level while some new ones are created. The cores from ' different levels are finally merged to form a hierarchical model of the object using cores. One advantage of this method is that the object does not have to be explicitly defined in the volumetric data. Although the representation seems promising for data where it is not easy to define the object boundary, and Morse et al. have presented preliminary results of its application, the representation is still in experimental stages. The merging process in the computation of cores seems difficult. When the object boundary can be Chapter 2. Survey of 3D Object Representations 27 denned, there are easier ways to derive the "core" or "skeletal" model. 2.2.2 Simplification of union of spheres models Recently, Hubbard proposed a method to simplify the union of spheres model obtained from the Voronoi approach to compute approximate skeletons [D.1.2-Hubbard 96]. This is in the context of building compact bounding-sphere hierarchies for objects for real-time collision detection. The simplification method merges two neighboring spheres in a Voronoi skeleton, if the merger does not increase a cost function by some pre-specified threshold. Two spheres are merged by replacing them with the smallest sphere enclosing the points defining the two spheres. The replacing sphere inherits all these points for possible successive mergers. One limitation of this simplification method for our purposes is that it may not result in a stable model, mainly because the volume enclosed by the smallest sphere enclosing the points defining two spheres a and b can be quite different from the volume enclosed by the union of a and b. Our simplification method is significantly different from Hubbard's method in at least two respects: (i) the criteria that decide which spheres should be grouped together, and (ii) the computation of the sphere that replaces a group of spheres. Our algorithm produces a stable model of an object and can merge any number of spheres as long as those spheres together satisfy our criteria. 2.2.3 Computation of mass properties It is possible to approximate the mass properties of objects by voxelizing the space containing them, that is, overlaying them with- a grid and checking which grid points are inside the object model and which ones are outside [D.3.6-Lee and Requicha 82]. Connolly's method measures analytically the surface area of the molecular surface de-fined by union of spheres [D.3.4-Connolly 83b]. It can also approximate other mass Chapter 2. Survey of 3D Object Representations 28 properties such as volume and center of mass. The topological and mass properties for a union of balls can be elegantly computed using computational geometry techniques [D.3.3-Avis et al. 88, D.3.5-Edelsbrunner 93]. These are based on the concept of a power diagram [D.3.1-Aurenhammer 87]. (a) Power diagram The power diagram or Laguerre-Voronoi diagram is a particular generalization (for spheres) of the Voronoi diagram for points [D.3.1-Aurenhammer 87]. Let S = {si, s 2 , . . . , s„} de-note a set of n spheres in a d-dimensional Euclidean space Ed. For this set of spheres, the power diagram PD(S) partitions Ed into convex polyhedra (i.e. a cell complex in Ed) defined as follows. Each s E S is associated with its power cell C(s), C(s) = {x e Ed | pow(x, s) < pow(x, t), Vt e S} where, pow (x, s) = of|;(x, c) — r2, c is the center of the sphere s, r is its radius, and dg(,) is the Euclidean distance. Figure 2.8 shows the power diagram for nine circles. The power diagram and Voronoi diagram are identical if the spheres in S have identical radii [D.3.1-Aurenhammer 87]. Let U(S) denote the union of spheres in S. The usefulness of power diagrams for visualizing unions of spheres lies in a theorem which converts the union to a summation [D.3.2-Aurenhammer 88]. It states that the union of spheres can be computed by sum-ming up the intersections of the spheres with their corresponding power cells. In other words, W(5) = s i U s 2 U . . . U s n = 53" = 1s inC7(s i) Note that some of these intersections may be empty (e.g., Figure 2.8, circle number 9). Chapter 2. Survey of 3D Object Representations 29 (c) (d) Figure 2.8: (a) Set of nine circles, (b) The associated power diagram, (c) A decomposition using the power diagram, (d) The union of circles. Chapter 2. Survey of 3D Object Representations 30 (b) Some observations , For a set of n spheres in 3D space, the following hold [D.3.2-Aurenhammer 88]: • Their power diagram can be computed in 0(n2) time and 0(n2) space. • Intersection of a sphere s with its power cell C(s) can be computed in 0{k) time, where k is the number of vertices in C(s). • The surface and volume of the union (or intersection) of spheres can be found in 0(n2) time. • The point location problem (i.e., determining whether a point belongs to the union) can be solved in 0(log 2n) time. • Various topological properties can be determined in 0{n2) time [D.3.5-Edelsbrunner 93]. 2.2.4 Disp lay of union of spheres models There are two issues concerning the display of a union of spheres: (i) speed, and (ii) smooth display. Both of these issues have been addressed before in the context of the union of spheres model for displaying molecules. (a) U n i o n of spheres Fast display techniques for a large number of overlapping or non-overlapping spheres by sorting the spheres in depth followed by their projection on the screen have been dis-cussed in literature [D.2.8-Knowlton and Cherry 77, D.2.9-Max 79, D.2.7-Knowlton 81, D.2.4-Franklin 81]. For overlapping spheres, one method for fast display is to use con-tours at regular intervals. Another method is to obtain a surface tessellation from the Chapter. 2. Survey of 3D Object Representations 31 union of spheres—this can take advantage of the current computer architectures that sup-port polygon display in hardware. If the spheres are polygonized, the polygons that are completely inside the union can be discarded. Connolly's method [D.3.4-Connolly 83b] and Edelsbrunner's method [D.3.5-Edelsbrunner 93] can compute the surface of the union of spheres. They can also give an approximate tessellation using triangular facets (the spherical triangles are approximated by planar triangles) [D.2.2-Connolly 85, D.3.5-Edelsbrunner 93]. Fournier and Ranjan developed a new method of obtaining an exact tessellation of union of spheres using generalized spherical triangles [I.2-Fournier and Ranjan 96]. Some special hardware architectures such as PixelPlanes™ can display spheres faster than polygons [D.2.5-Fuchs et al. 85], but these architectures are not commercially avail-able. (b) Blended spheres As discussed before, spheres can be blended together for smooth display using "blobbies" [D.4.1-Blinn 81, D.4.3-Wyvill et al. 86]. The radius of influence of atoms-provides a control on the smoothness of the resulting surface. If the union of spheres model is obtained from volumetric data, for all atoms we can set rinfl = riso + 0.5^, where v is the inter-voxel distance. This ensures that sufficiently close spheres blend with one another [I.6-Ranjan and Fournier 94]. Connolly studied the smooth display of union of spheres in terms of solvent-accessible surfaces of molecules [D.2.2-Connolly 83a]. The idea is to consider a molecule as a union of spherical atoms of different sizes, and find where a probe atom of fixed size (the solvent) can go without intersecting any atom. The boundary of the portion of the 3-space which is untouched by the probe defines the Connolly surface. Several efficient implementations to compute the solvent accessible surfaces of molecules exist [D.3.4-Connolly 83b, D.2.11-Sanner et.al. 95]. Chapter 2. Survey of 3D Object Representations 32 A number of other methods exist for polygonizing implicit surfaces defined by the union of blended or non-blended spheres. One way to do this is to voxelize the space and then use an algorithm similar to marching cubes [D.4.3-Wyvill et al. 86]. Since this tes-sellation depends on the voxelization artifacts and not on the complexity of the surface, generally, a large number of triangular facets are required. This can be improved by using a better voxelization data structure such as an octree [D.2.1-Bloomenthal 88]. A n adap-tive algorithm that depends on the complexity of the surface (curvatures of arcs on the surface) has been presented by Overweld and Wyvil l [D.2.10-Overweld and Wyvi l l 93]. Recently, Wyvi l l et al. have extended their polygonizer to combine CSG (Constructive Solid Geometry) with implicit surfaces. 2.3 Estimation of properties for various representations Table 2.1 shows our subjective evaluation of the representations discussed in this chapter with respect to our criteria. A '1' for a particular property for a representation means that the representation is "poor" (inconvenient or inefficient) for that property; a '5' means that the representation is "good" for that property. A '?' means that we do not have enough evidence to judge the representation for that property. For instance, in the category "stability," octrees are a '1' because a small change in the position can lead to a big change in the model; surface-mesh is a T in this category because a small change in the data can change significantly the primitives in the model. Surface-mesh is a '5' for "generation" because it can be generated quickly (for example, using marching cubes on volumetric data). Octrees and wavelets, being inherently hierarchical, are.a '5' for "scalability." Chapter 2. "Survey of 3D Object Representations 33 P r o p e r t y R e p r e s e n t a t i o n Generation Storage Display Manipulation Completeness Accuracy Stability Scalability Interactivity Transformability Texture mapping Surface-mesh 5 2 5 1 1 4 1 •3 5 5 5 Octree 4 '5 2 1 5 2 1 5 3 2 1 Shell 2 4 2 1 4 2 ? 5 2 2 1 Wavelet 2 5 2 1 4 4 ? 5 2 2 1 Volume primitives 1 5 3 5 5 3 3 1 1 5 ? BSP trees A 5 4 3 5 4 2 4 4 5 4 Tetrahedra 4 4 4 1 5 4 2 2 3 5 4 Alpha-shape 3 4 4 1 1 2 2 5 1 5 1 Union of spheres 4 4 2 4 5 3 4 4 1 5 ? Table 2.1: Estimation of properties for various representations (1-poor, 5-very good). 2.4 Summary In this chapter we discussed some of the existing representations of 3D objects. It should be clear that different representations are good for different applications; no representa-tion is superior to others in all respects. The choice of a representation depends on the application at hand. However, most representations lack the property of stability; this prohibits their use in applications involving object comparison. We will see later how we can obtain stable UoS models of objects for such applications. Chapter 3 Generation of a Union of Spheres Model for Objects We generate stable1 UoS models of objects from boundary points in two steps. In the first step, we generate a UoS model which is faithful to the input data in the sense that all the boundary points lie on its surface. In this chapter, we present algorithms to derive such a model for 3D objects from volumetric data (such as CT and MRI scans) and from the depth data available from laser range-finders and ray-tracers. Essentially, this means that we can obtain the UoS model for objects using several available sensors and from virtually any other existing computer model of the object. In the second step, we simplify the "dense" UoS model generated by the first step to stabilize it and to reduce the number of primitives. This is done by clustering similar or near-by spheres and eliminating redundant or non-significant spheres (in terms of contribution to the volume of the object). Together, the two algorithms yield stable UoS models of objects from their boundary points. The simplification algorithm will be described in the next chapter. \ 3.1 General approach Our algorithm is based on the Delaunay tetrahedralization (triangulation in 2D) of the boundary points. We use Delaunay tetrahedralization because of its empty sphere prop-erty, which guarantees that we will fit the largest spheres possible defined by four bound-ary points and not including other boundary points. 1 We will show the stability property of the UoS representation in Chapter 6. 34 Chapter 3. Generation of a Union of Spheres Model for, Objects 35 Empty Circle (Sphere) Property The Delaunay triangulation of a set of sites in a plane has the property that the cir-cles circumscribing the Delaunay triangles are defined by three sites (assuming points in general position) and cannot contain any other site (Figure 3.9). Similarly, in 3D, Figure 3.9: Empty circle property of Delaunay triangulation. the spheres circumscribing the Delaunay tetrahedra are defined by four sites and do not contain any other site. The steps of the algorithm are: 1. Generate the boundary points (if not available). 2. Compute the Delaunay tetrahedralization of boundary points. 3. Compute the spheres circumscribing the Delaunay tetrahedra. 4. Discard the spheres "outside" the object and keep the ones that are "inside". This requires a test for determining the outside and inside spheres. Chapter 3. Generation of a Union of Spheres Model for Objects 36 The steps discussed here are 1 and 4. In the following sections, we will describe these steps for several types of data: (i) volumetric data, (ii) range data, (iii) ray-casting data, and (iv) polygonal-mesh data. The complexity of the algorithm is 0(n?) in the worst case, where n is the number of boundary points, because Delaunay tetrahedralization has the complexity 0(n2) and all other steps are linear as a function of the number of spheres produced. In our experiments, the maximum number of spheres obtained so far is about 3n. The technique of using the Delaunay tetrahedralization of boundary points, followed by an "inside-outside" classification of .the spheres circumscribing the tetrahedra has been used before to obtain skeletons of 3D objects (for example, [D.5.4-Goldak et al. 91, D.5.2-Brandt 92]). Most of these methods assume that a polyhedral model of the object is given and that a "point-in-polyhedron" test can be used to classify the tetrahedra (and the corresponding spheres) as inside or outside. For example, we can check if the centers of the spheres circumscribing the Delaunay tetrahedra are inside the polyhedral model or outside, and we keep only the spheres for which the centers are inside. However, there are three problems with this approach: (a) How do we classify the spheres if a polyhedral model of the object is not available? More specifically, what do we do for other types of data, e.g., volumetric data or range-finder data? (b) How do we ensure that the error in the reconstruction is bounded? In particular, how do we know that some tetrahedra (spheres) are not partially inside and partially outside the object? (c) Does the reconstruction converge to the actual object in the limit as the density of boundary points increases? Chapter 3. Generation of a Union of Spheres Model for Objects 37 We will address these issues in detail in this chapter and in Chapter 5. The next four sections provide answers to (a). One common approach to answer (b) and (c) is to assume that the polygonal boundary of the object is .contained in the faces of Delaunay tetrahedralization [C.3.1-Boissonnat 84]. Then the tetrahedra can be classified as completely inside or outside the object and partially inside, partially outside tetrahedra are avoided. This is generally a reasonable approximation (there is no guarantee however) if the boundary point density 8 is small compared to the size of object features. Under the limit as the density of boundary points increases indefinitely, it can be shown that the centers of Delaunay spheres approach the skeleton of an object and that their union tends to the actual object [D.5.4-Goldak et al. 91]. 3.2 UoS model from volumetric data To generate UoS models from volumetric data (such as C T or M R I data), we noticed that if we derive the. boundary points by interpolating between all pairs of inside-outside voxels, the empty-sphere property of Delaunay spheres can be used to classify the tetra-hedra (spheres). This provides a fast and easy method of classification of tetrahedra (spheres) as inside or outside and shows that for volumetric data partially inside and partially outside tetrahedra are avoided. The input to our algorithm is the voxel data for the object and a threshold (or isovalue) or a pair of thresholds (in general, any binary classifier) which defines the object of interest in the volumetric data. From this binary voxel data we generate a set of points on the surface of the object to be reconstructed. We assume that the surface of the object has been reasonably sampled by this procedure, i.e., 8 is small compared to the size of the features in the object. This is because the quality of reconstruction will be determined by the errors induced by sampling, namely, we will not reconstruct features that are small with respect to 5. Chapter 3. Generation of a Union of Spheres Model for Objects 38 3 . 2 . i Algorithm for UoS model from volumetric data The generation algorithm creates a UoS model which is faithful to the input volumetric data in the following sense: (a) all inside voxels are included in at least one sphere, (b) no outside voxel is within any sphere, (c) all boundary points are on the surface of at least one sphere, (d) no boundary point is inside any sphere. For sake of clarity, the algorithm will be described in two dimensions. For 3D, substitute "tetrahedron" for "triangle," and "sphere" for "circle." We will explain the steps with the help of the object shown in Figure 3.10a, a computed tomography (CT) slice through a human head at a 32 x 32 resolution. The steps in the model generation algorithm are: 1. Compute the boundary points: A l l points at the value of the threshold are defined to be on the boundary of the object. In this example, we calculate a discrete set of boundary points using linear interpolation between all pairs of neighboring (in the 4-connected sense) inside and outside voxels. In Figure 3.10b, the selected regions correspond to a vertebra and the lower jaw. 2. Compute the Delaunay triangulation of the set of boundary points (Figure 3.10c) generated above. We assume here that the points are in general position (see glossary). This is not a limitation because there are reliable ways to ensure this [C.3.2-Edelsbrunner et al. 94]. 3. Compute the circumscribing circle of each triangle. 4. Test for inside and outside circles: By the empty circle property of the De-launay triangulation, a circumscribing circle does not contain any other boundary point. This also means that these circles contain either only outside voxels or only inside voxels. This is true because if any circle contains both inside and outside voxels, it will contain a boundary point as well and this is forbidden by the empty Chapter 3. Generation of a Union of Spheres Model for Objects 39 circle property of the Delaunay triangulation. Therefore, i f any voxel wi th in a circle is known to be inside (or outside), the whole circle can be defined to be inside (or outside). Use the original volumetric data to verify which circles are inside. This is done by checking the value of the voxel if any voxel is inside the circle. Instead of working wi th voxels, we can check if the median points of triangles, or the centers of circumscribing circles are inside or outside. B y using values at voxels we avoid unnecessary interpolation. If there is no voxel wi thin the circle, we use interpolation to see if the center of the circle is inside (or outside). For consistency, we use the same interpolation scheme for the inside/outside test of the non-voxel point as was used to derive the boundary points. It is important to note that inside and outside is here a classification of the whole primitive (triangle or circle), but not of a l l the points inside. In particular an outside circle can intersect an inside one (but it cannot happen with the triangles). The inside circles define our union of circles (UoC) model (Figures 3.10d and 3.10e). The model based on inside triangles wi l l be called union of triangles (UoT) subsequently. S im-ilarly, from 3D data, we obtain union of spheres (UoS) and union of tetrahedra (UoTet) models. A s is shown in the next section, the differences in the properties of the U o T and U o C models depend on the boundary point density 8 on the surface of the object. If 8 is small compared to the "size" of the object features (we wi l l make this more precise in Chapter 5), the two models are close approximations of each other. 3.2.2 Some properties of the generated models This section presents some observations (and their brief justifications) relevant to the models described above. The results about error w i l l be discussed in Chapter 5. Chapter 3. Generation of a Union of Spheres Model for Objects 40 Figure 3.10: Steps of the algorithm to generate the UoC model, (a) A (32 x 32) size C T slice of human head (range 0-255). (b) Boundary points for a threshold of 70.3 (c) Delaunay triangulation of the boundary points, (d) UoT model, (e) UoC model, (f) UoC and UoT models superimposed, and the centers of the circles. Chapter 3. Generation of a Union of Spheres Model for Objects 41 Observation 1 Let VT denote the volume of the UoTet model and Vs the volume of the UoS model' Let A be the maximum size of the edge in any boundary face. Let V\ denote the volume of a sphere with diameter A. Then, for an object with one component and no holes, VT < Vs < Vr + (n — 2)V\, where n is the number of boundary ^points. Since we have a circumscribing sphere to every tetrahedron, obviously, Vs > VT- The second inequality follows from the facts that for an object with one component and no holes, tetrahedralization of n boundary points has (2n — 4) boundary faces, and for every boundary face, there is one spherical cap of volume at most Vx/2 (Figure 3.11a). This bound is quite pessimistic. If the diameter of a sphere passing through the three vertices of a face is greater than A (as is mostly the case), the volume of the spherical cap is less than Vx/2. (a) (b) Figure 3.11: (a) A circular cap of volume Vx/2 on one of the boundary faces, (b) A con-figuration of boundary points so that the circles touch each other but the corresponding triangles do not. Chapter 3. Generation of a Union of Spheres Model for Objects 42 Observation 2 The topology of the UoTet model and the topology of the UoS model may differ. The topology of the UoTet model and that of the UoS model differ when two spheres touch each other but the corresponding tetrahedra do^not (Figure 3.11b). When this happens, it is indicative of the fact that sharp features exist in the object. In these places, either topology (UoS or UoTet) can be correct; in fact, this case is ambiguous for other algorithms (such as marching cubes) as well. One related consequence of this is that two components or branches of the same component separated by more than A cannot join. This is because for a maximum edge size of A, the maximum a sphere can protrude is A/2 (see Figures 3.11a and 3.11b). 3.3 UoS model from range data In this section, we present a method to obtain a UoS model of an object from range data. A number of attempts have been made to model 3D objects by scanning their visible surface from different directions using a range-finder and then integrating the depth information (generally available as points on the surface of the object) provided by the different views. Some surface-based and volume-based methods for reconstruction of models from range data were discussed in Chapter 2. Model acquisition by this tech-nique can be broken down into three sub-problems: (i) data acquisition, (ii) registration of data captured from different views, and (iii) integration of information provided by the different views. This section deals with (iii) only. Registration of different images is crucial (although sometimes it can be avoided), but there are fairly reliable ways to en-sure this [E.l-Besl and McKay 92, E.2-Chen and Medioni 92, E.7-Turk and Levoy 94]. Although we do not address the issue of registration, our method is relatively insensitive to small registration errors; this can be shown experimentally using a method similar Chapter 3. Generation of a Union of Spheres Model for Objects 43 to the one described in Section 6.4.3 to study the stability property of the UoS repre-sentation (namely, slightly misaligning the different scans, computing a new UoS model using the misaligned point set, and comparing the new UoS model with the original UoS model). 3.3.1 Algorithm for UoS model from range data As in the case of volumetric data, we use the spheres circumscribing the tetrahedra computed from a Delaunay tetrahedralization of the boundary points (which are available as the input). We need a test to determine if a sphere is inside the object or outside. Test for inside and outside spheres: In volumetric data, classification of spheres is possible using the voxel values (as described in the previous section). In the present case, we make use of the visibility information to . classify the spheres as follows: • Find the points that define a given sphere (namely, the vertices Vj, % = 1..4 of the Delaunay tetrahedron). We assume again that the points are in general position. • For each of the four points, do the following: - Find the direction from which this point is visible (i.e., direction of scanner —* for the range-finder). Call that direction d^ . - Define another vector n; as the direction from the center of the sphere to the point. This is the normal to the surface. —* - Take the dot-product n^ • dj of the two vectors; a positive value means that we see the outside of the sphere whereas a negative value indicates that we see the inside (Figure 3.12a). Chapter 3. Generation of a Union of Spheres Model for Objects 44 Inside circle third scan A Direction of third scan (b) Figure 3.12: Classification of spheres as inside or outside, (a) Usual inside and outside cases, (b) A pathological case. Chapter 3. Generation of a Union of Spheres Model for Objects 45 • If all four dot-products (corresponding to four points) are positive, the sphere, is temporarily classified as inside. Otherwise, if any of the dot-products is negative, the sphere is classified as outside (Figure 3.12a). This is because from the view-point we are supposed to see the outside of the object. If we see the inside of the sphere, there cannot be an element of surface in front of it, and therefore the sphere is outside (note that this test cannot fail except for reasons of numerical inaccuracies). Unfortunately, in certain cases it is possible to see the outside of an outside sphere; see an example of a possible problem in Figure 3.12b. Therefore, this local test is not enough to correctly classify all the outside spheres. To resolve "these cases, if necessary,2 we perform an additional verification step for the spheres temporarily classified as inside. Veri f icat ion of inside spheres The "inside" spheres can be verified by checking the condition that any point on the visible surface of the model should be within a distance 5 of a sampled boundary point. Any sphere that violates this condition is eliminated. The power diagram of the spheres in the model can be used to compute the contributions of the spheres to the surface of the model. However, this approach may not be very efficient as recomputation of (parts of) power diagram may be necessary after each elimination. To test which spheres should be kept and which ones should be discarded we perform a ray-cast of the spheres in the generated model from all the directions from which it was scanned. One underlying assumption for this test is that 5 is known; the density of rays that we use depends on the resolution of the laser range-finder (which in some sense determines 5). For each ray, we find the first intersection point p between the ray and the spheres and check if there is a sampled boundary point within a distance 5 of p. If for 2Results on practical data for several objects have few large undesired spheres. Chapter 3. Generation of a Union of Spheres Model for Objects 46 a ray, no such boundary point is found, we discard the intersected sphere and continue with the next intersection of the ray until either a valid intersection is found, or the ray crosses the bounding box. When 5 is unknown^ we traverse the rays beginning from the known boundary points in the direction of their respective viewpoints. A l l spheres intersected by these "visible" portions of rays (the part of rays from their origin to the visible points) are removed; this may remove some "inside" spheres that are grazed by these rays but this seems to be the most reasonable thing to do. . The verification algorithm can be applied directly to the model containing all inside and outside spheres, but the combined approach is much faster in practice (the first test based on 4 dot products eliminates most of the undesired spheres). Note that this method may not work on the tetrahedra because Delaunay tetrahedra can be arbitrarily thin and rays can miss them. The time complexity of the verification step using the ray-casting approach depends on the number of rays that we use. 3.4 UoS model from existing computer models e This section deals with the problem of deriving UoS models of objects from their existing computer models. Generating a union of spheres model of an object can be posed as an optimization problem as follows: Given n, find n (possibly overlapping) spheres inside the object such that they optimize a property, e.g., cover the maximum volume of the object. This is a complex optimization problem and an efficient solution seems unlikely. Techniques such as simulated annealing can be tried but they will be inefficient because they will require computation of the volume of the union of spheres (which has a time complexity of 0(n 2 ) , n being the number of spheres) several times in every iteration. As discussed in the previous chapter, several other techniques exist which yield a union of spheres model of an object from the boundary points in a reasonable time, therefore, the Chapter 3. Generation of a Union of Spheres Model for Objects 47 problem is rarely (if at all) posed as an optimization problem. We/generate a UoS model of an object from its existing computer model using the general approach described in Section 3.1. We first generate a set of points on the boundary of the object by ray-casting it from different directions (the number of views to use, the choice of viewpoints, and the choice of resolutions depend on the object). This provides data similar to the range-finder data—a set of points on the boundary of the object, and the corresponding directions from which they are visible. Given this data, the generation algorithm is exactly the same as that for range data. Sometimes, the original model may aid in the verification of inside spheres. 3.5 U o S model from a polygon-mesh The input in this case is a set of points on the boundary of an object, a list of polygons defined by those points, and the "outward" normals for the polygons. One way to generate the UoS model from this data is to ray-cast this polygonal model and then use an algorithm similar to the one discussed in the previous section. However, this will (in general) create points on the polygonal model that is itself an approximation of the original object. We can use directly the points provided as vertices of the mesh. We use a test similar to the test for range data except that instead of using a visibility direction of a point, we use the "outward" pointing normal at every vertex. 3.6 Implementat ion and results The above-described algorithms (except UoS model from polygon-mesh) have been im-plemented in the C language using the G L graphics library and compiled for Silicon Graphics and I B M RS/6000 platforms. We will first present the results for 2D objects and then show the results of reconstruction in 3D. Chapter 3. Generation of a Union of Spheres Model for Objects 48 Figure 3.13 shows the results of applying the volumetric data algorithm on an M R I slice through a human head (256 x 256 pixels with physical resolution in x and y of about 0.1 mm). The region selected for a threshold of 100.3 corresponds to a section through the brain. There are 1696 boundary points. The inside-outside classification scheme resulted in 1703 inside triangles and circles in the UoT and UoC models. The computation time was 0.3 seconds on an IRIS Indigo 2 machine. Note how close the UoT and UoC models actually are (this will be formalized in the next chapter). Figure 3.14 shows the results of the same reconstruction technique on a 2D "cow" object scanned at a resolution of 512.x 512 pixels. Figure 3.14a shows the boundary points and Figure 3.14b shows the union of circles in the. UoC model (2054 circles). Figure 3.14c shows the histogram of distribution of circles in the model in terms of their size (measured by radius). One use of the histogram can be in displaying the structure of the object. First, circles with relatively larger radii can be displayed to show the structure of the object at a "coarse" level. Then circles with relatively smaller radii_can be incrementally added to show the "fine" level details. Figures 3.14d-f show the results of this process. Figures 3.15 and 3.16 show the results of the algorithm on two 3D objects: (i) a C T scan of the femur bone3 (courtesy Insituti Ortopedici Rizzoli, C I N E C A , Cray Research, and University of Bologna), (ii) P E T brain data (courtesy, U B C hospital). Table 3.2 summarizes the results for these objects. Delaunay tetrahedralization being an 0(n2) process is the bottleneck for our algorithm in terms of speed: 1 For construction of UoS models from range data and other existing models, we show the results on three data sets: (i) the range-finder data for a clay rabbit (courtesy Cy-berware and Stanford University [E.7-Turk and Levoy 94]), (ii) for the polygonal model of a knot (an algorithmically generated object, courtesy Rob Scharein), and (iii) for the 3 U R L address: h t t p : / / s i r i o . c i n e c a . i t / p r o m e t e o / d a t a s e t s . h t m l Chapter 3. Generation of a Union of Spheres Model for Objects 49 Figure 3.13: Result on MRI brain data, (a) A (256 x 256) size MRI slice of human brain (range 0-255). (b) Boundary points for a threshold of 100.3 (1696 points), (c) Delaunay triangulation of the boundary points, (d) UoT model (1703 triangles), (e) Union of triangles in the UoT model, (f) Union of circles in the UoC model (1703 circles). Figure 3.14: Result on cow data, (a) Boundary points for cow (512 x 512 pixels), (b) Union of circles in the UoC model (2054 circles), (c) Histogram of radius distribution, (d-f) Use of histogram in incremental display of details in the object. Chapter 3. Generation of a Union of Spheres Model for Objects 5 1 (S) <D s-, <u ft 0 (HI 1 3 2 4 6 8 10 12 14 16 18 20 radius (voxel units) (d) (g) Figure 3.15: Results of the generation algorithm on a (64 x 64 x 128) C T data for femur bone, (a) Boundary points (17384 points), (b) Surface model using marching cubes (34800 triangles), (c) UoS model (49902 spheres), (d) Histogram showing the radius distribution, (e-g) Use of histogram to control display of object (e) 8685 spheres, (f) 17531 spheres, and (g) 29357 spheres. Chapter 3. Generation of a Union of Spheres Model for Objects 52 (d) Figure 3.16: Results of the generation algorithm on a (64 x 64 x 64) P E T brain data, (a) Boundary points for a threshold of 140.3 (range 0-255, 17598 points), (b) Surface model using marching cubes (35216 triangles), (c) UoS model (58194 spheres), (d) UoS model with spheres of radius less than the inter-voxel distance removed (27917 spheres). Chapter 3. Generation of a Union of Spheres Model for Objects 53 Femur Bone B r a i n Voxel resolution 64 x 64 x 128 64 x 64 x 64 Number of surface points 17384 17598 Time to tetrahedralize 8 minutes 8 minutes (Iris Indigo 2) Inside-outside check 25 seconds 27 seconds Number of inside spheres 49902 58194 Table 3.2: Results of our algorithm on 3D data. polygonal model of a hand (a modeled object, courtesy Sidi Yu). The boundary points P rope r ty R a b b i t K n o t H a n d Directions of scans 10 14 6 Number of surface points 22652 13639 15227 Time to tetrahedralize > 10 minutes 3 minutes 4 minutes (on Iris Indigo-2) Inside-outside check 35 seconds 20 seconds 22 seconds Number of inside spheres 57931 38993 42258 Table 3.3: Results of our reconstruction algorithm on rabbit, knot, and hand data. for (ii) and (iii) were obtained from the ray-tracer optik [G.l-Amanatides 87]. The in-formation about the three objects is included in Table 3.3 and the rendered images are shown in Figures 3.17, 3.18, and 3.19. In Figures 3.17c and 3.18c we can see some of the "misclassified" (relative to the polygonal model) spheres for the rabbit and knot data. Of course, some of these spheres may be as faithful to the input data as are the "properly" classified spheres. Generally, when the boundary of the object is densely sampled, there are a large number of small spheres. Our simplification algorithm eliminates most of the small spheres (e.g., see the results of simplification on the rabbit data in Figure 4.25). If we know that the underlying object has no features of size smaller than A, it makes sense to remove all spheres with size less than A; they are merely the artifacts of the reconstruction Chapter 3. Generation of a Union of Spheres Model for Objects 54 Figure 3.17: Results of the generation algorithm on 3D rabbit data, (a) Point sampling (22652 boundary points), (b) Reconstructed polygonal model [Turk and Levoy 94] (45300 triangles), (c) UoS model (57931 spheres). Chapter 3. Generation of a Union of Spheres Model for Objects 55 Figure 3.18: Results of the generation algorithm on 3D knot data, (a) Polygonal model (3200 triangles), (b) Point sampling (13639 boundary points), (c) UoS model (38993 spheres). Chapter 3. Generation of a Union of Spheres Model for Objects 56 (c) Figure 3.19: Results of the generation algorithm on 3D hand data, (a) Polygonal model (1287 triangles), (b) Point sampling (15227 boundary points), (c) UoS model (31883 spheres). Chapter 3. Generation of a Union of Spheres Model for Objects 57 process.4 In Figure 3.19c we show the UoS model for the hand after small spheres were removed using the histogram. In the original model there were 42258 spheres and in the shown model there are 31883. It is interesting to note that the UoS model looks fairly smooth when the artifacts of reconstruction are removed. 3.7 Summary In this chapter we described algorithms to generate UoS models of objects from volu-metric data such as obtained from MRI and C T scanners, and surface data as obtained from laser range-finders. A similar method can also be used to generate the UoS model of an object from virtually any other existing model. The approximation properties of the UoS representation will be presented in Chapter 5. 4 Even results on smooth objects have small spheres "hanging" out. Some of these spurious spheres can be. explained by the perturbation of boundary points (to maintain the generality condition) in the Delau-nay tetrahedralization code that we use. These will not appear if some other techniques for maintaining the generality condition of points (e.g., Simulation of Simplicity (SoS) [G.7-Edelsbrunner and Miicke 90]) are used. Chapter 4 Simplification of Union of Spheres Model The main motivations for simplification of UoS models are to remove nearly identical or redundant spheres that do not contribute much to the surface or volume of the object, to speed up manipulation and display, and to facilitate comparison of objects by reducing the models to almost the same number of primitives. A n ideal algorithm to do this should be efficient, should give control over errors and the number of primitives, should yield a stable model, and should preserve topology. These properties are difficult to achieve simultaneously. Clustering n points into c clusters so that they minimize some function (e.g., sum of squared distances of points from the centers of their respective clusters), is a computation-ally expensive problem (the time complexity is exponential in n [H.3-Duda and Hart 73]). A n efficient deterministic algorithm for the clustering problem is therefore not possible; almost always (unless n is small), heuristics are used. In this section we present an algo-r i thm which yields a stable model, gives some control over the simplification error and the. number of primitives used. Unlike many point-clustering methods that can be used (e.g. median cut [G . l l -Heckbe r t 82]), the algorithm makes use of the fact that we are clustering spheres and not points. It is fairly efficient but does not guarantee preservation of topology. 58 Chapter 4. Simplification of Union of Spheres Model 59 (a) (b) Figure 4.20: Figure showing the stability of the UoC representation. 4.1 Triangles versus circles Figure 4.20 shows an example where the triangles in the Delaunay triangulation are un-stable but the circumscribing circles are stable. Figure 4.20a shows Delaunay triangles and the corresponding circumcircles for a set of eight nearly co-circular points (6 triangles and 6 circles). Figure 4.20b shows the Delaunay triangles and the corresponding circum-circles when the points were perturbed by small amounts. This example demonstrates that: 1. There can be cases when the UoT representation is unstable but the UoC repre-sentation is stable. Note that if the UoT representation is stable, so is the UoC representation, but not vice-versa. 2. There can be redundancy in the UoC model in the case of nearly co-circular points, and the circles generated from these points can be replaced by one circle, with little added error. This is somewhat analogous to replacing a number of adjacent nearly Chapter 4. Simplification of Union of Spheres Model 60 co-planar polygons by one polygon in the simplification of a polygonal model. The second point motivated our approach to simplification. 4.2 Clustering of spheres 4.2.1 Sphericity of a set of spheres The simplification algorithm is based on the concept of sphericity, which is a measure of how well a set of spheres can be modeled by a single sphere. Sphericity is defined as the ratio of a dimension of the union of spheres to a dimension of a sphere including them. The obvious choices for the dimension are linear (the radius or diameter), quadratic (the (a) (b) sphericity a = r / R Figure 4.21: Definition of sphericity (a), (a) Large sphericity, (b) Small sphericity. surface area) or cubic (the volume). For computational convenience we have chosen a linear measure, specifically the ratio of the radius of the largest sphere, in the set to the radius of the smallest enclosing sphere (Figure 4.21). According to this definition, sphericity is a number between 0 (very non-spherical) and 1 (perfectly spherical). This Chapter 4. Simplification of Union of Spheres Model 61 measure is related to the other ones, by noting that the surface area of the largest sphere is less than or equal to the surface area of the union 1, and similarly the volume, of the largest sphere is less than or equal to the volume of the union. We summarize these relationships in Table 4.4. Note that the notion of sphericity is scale-independent (as it is a ratio); this helps in comparison of objects at different scales after appropriate simplification. The linear measure of sphericity can be related to how close the boundary points are to being co-spherical. Assume, for instance, that we have a set of points that are on the surface of a sphere of radius r with a deviation such that they are all at distance between r ( l ± k) from the center of the sphere. The radius of the minimum enclosing sphere is r ( l + k). By the property of the Delaunay triangulation (it contains the largest empty sphere), there is a sphere inside of radius at least r ( l — k). Therefore the sphericity of Order Measure of UoS Measure of Largest Sphere Measure of Enclosing Sphere Sphericity Approx. Relation Linear rc oL = ^ rr, • aL = 1 - e Quadratic A ATTTL2 < A Aurc2 OA ~ 4 n r p 2 > ° L aA>l-2e Cubic V 47rr L 3/3 < V 4.irrc3/3 °V ~ 47rrc3/3 ^ °" L oy > 1 — 3e Table 4.4: Various definitions of sphericity and relationships. the UoS is at least (1 — k)/(l + k) which is approximately 1 — 2k if k <C 1. This means that if we chose a sphericity of 1 — e, we are sure to be able to enclose a set of boundary points which deviate by less than 1 ± 0.5e from the true sphere.- Stability is achieved in that case because any initial Delaunay tetrahedralization of the set of nearly co-spherical points with noise of that magnitude will eventually give the same (within comparable variations in the position of center and radius) clustered sphere. : Note that the surface area of the union can be more than the surface area of the enclosing sphere. Chapter 4. Simplification of Union of Spheres Model 62 4.2.2 A l g o r i t h m descript ion The simplification algorithm processes the spheres from the largest to the smallest. In every iteration, using the largest remaining sphere a, the algorithm calculates a new sphere c such that c encloses a, the set of spheres it represents has a sphericity > o (a user defined threshold, generally close to 1), and c tries to enclose as many other spheres as possible. A l l the spheres within c are then deleted and the iterations continue until all the spheres in the original model have been taken care of. Let S denote the set of spheres to be simplified. The steps in the simplification algorithm are: 1. Mark all the spheres in S as "unclustered." 2. Select the largest "unclustered" sphere a G S. Let r a be its radius. Calculate a sphere b of radius ra(2/a — 1) concentric with a. The cluster c of radius r c we want to obtain has to be within b to satisfy the constraint that ra/rc > a. 3. Mark all "unclustered" spheres in S that are inside b as possible candidates for . inclusion in c. Denote this set as G. Find the smallest bounding sphere c for spheres in G [G.16-Welzl 91]. 4. Compute the sphericity of G, a a = ra/rc. If O~Q > cr, goto step 5 else goto step 6. 5. Keep c as the new cluster. Mark all the "unclustered" spheres in G as "clustered." If there are "unclustered" spheres in S, goto step 2, else quit. 6. Sphere c is defined by 2, 3, or 4 spheres [G.16-Welzl 91]. Out of these defining spheres, eliminate the sphere, different from a, from G whose removal increases the sphericity of G the most. Goto step 4. Chapter 4.. Simplification of Union of Spheres Model ^ 63 The loop in Step 6 wi l l terminate because when a sphere is removed from G, the sphericity either increases or stays the same and when only a is left in G, OQ = 1. Note that it uses a greedy approach, which eventually leads to a sphere c which satisfies ra/rc > a but may not be the op t imal 2 one globally. The worst case time complexity of the simplification algorithm is 0(n3), n being the number of spheres in the model to.be simplified. The actual implementation uses a grid data-structure to aid in searching the spheres contained in a given sphere [G. l -Amanat ides 87]. O n the tested objects, the running time of the simplification algorithm is about 3-20 times faster than the time for tetrahedralization. We have used the algorithm to simplify a variety of objects. It can be applied wi th different values of sphericity threshold a and the simplified model always contains the original model. Since the algorithm ensures that a l l the sets of clustered spheres have a sphericity > a, it puts a A bound on the simplification error. The algorithm does not guarantee preservation of topology (similar to most clustering algorithms) but it w i l l not disconnect connected components, or create cavities in the object, because the enclosing sphere contains entirely the spheres of the cluster. O n the other hand it may connect components previously disconnected, or make existing holes and cavities disappear. The enclosing sphere can extend at most 2 r ( l — a)jo beyond the enclosed sphere, r being the radius of the largest enclosed sphere. It is easy then to see that for an intersection to be created where there were none, the gap has to be less than the sum of these terms for each cluster. In particular, i f we know the radius of the largest sphere, we know that nothing at a distance more than 4 r ( l — a)/a w i l l become connected i f it is not already. For example if o = 0.98, this distance is about 8% of the largest sphere, and can be controlled. One other "flaw" wi th the simplification method is that the volume is not preserved during simplification; it increases as the sphericity threshold o decreases. This 2 Optimal here means the smallest sphere that includes the largest cluster above the sphericity threshold. Chapter 4. Simplification of Union of Spheres Model 64 can be taken care of by shrinking the simplified spheres such that the original volume is preserved, but it involves computing that volume, and can create more disconnected components. 4.2.3 Verification of spheres in the simplified model One problem with the simplification method is that it can create spheres that are inside the union of other generated spheres. See Figure 4.22 for an example. To take care Figure 4.22: Generation of redundant circles during simplification, (a) Configuration of circles to be simplified, (b) A is simplified to A ' , C to C, and B to B'. B' is contained in the union of A' and C. of these redundant spheres we can use the power diagram [D.3.1-Aurenhammer 87]. According to a property of power diagrams, for spheres that do not contribute to the volume of the union of spheres, their intersection with their power cell is empty. Therefore spheres which have an empty intersection with their power cell can be safely eliminated. The calculation of the power diagram and its use to compute the volume of union of spheres (area of union of circles in 2D) is publicly available [D.3.5-Edelsbrunner 93] and we used it in our implementation. Chapter 4. Simplification of Union of Spheres Model 65 4.2.4 Selection of spherici ty The choice of sphericity is generally governed by two factors: (i) the error tolerance, and (ii) the number of primitives. For (i) we desire that any point on the surface of the new model be within a distance e from the closest point on the surface of the original model. To achieve this, the sphericity threshold a is not kept constant but is allowed to change; it is different in every iteration and is denned by the radius of the largest remaining sphere and the user-specified tolerance e (see the discussion on sphericity above). For (ii) we find the appropriate value of a by the trial and error approach. The number of primitives in the simplified model is monotonically non-increasing with a. In the extreme cases, when a = 0, we get a single sphere as the output, and when a = 1, we get the original UoS model itself as the output. 4.3 Implementat ion and results The algorithm has been implemented in 2D and 3D in the C programming language using the G L graphics library and compiled for SGI and I B M workstations. The results of the simplification algorithm are first shown in 2D for clarity of illustration. Table 4.5 provides some statistics regarding simplification of the 2D example objects. Brain Giraffe (n = 1703 circles) (n = 2110 circles) a circles ratio time abs. error circles ratio time s n/s (sec.) (pixels) s n/s (sec.) 0.97 627 2.72 0.7 1.0 267 7.90 2.2 0.95 470 3.62 0.3 2.0 157 13.44 2.1 0.90 228 7.46 0.1 3.0 112 18.84 1.2 0.85 142 11.99 0.2 4.0 •89 23.71 1.0 Table 4.5: Summary of results of the simplification algorithm on 2D data. Figure 4.23 shows the result for the UoC model of a slice through the human brain. Chapter 4. Simplification of Union of Spheres Model 66 It is simplified for different values of sphericity threshold and the discrepancy between the original model and the simplified model is also shown (Figure 4.23). Figure 4.24 shows the results for the UoC model of a giraffe scanned from a book. It is simplified for different values of absolute error and the discrepancy between the original model and the simplified model is also shown (Figure 4.24). Simplification results for 3D data are presented in Figures 4.25, 4.26, and 4.27. We show the results on objects obtained from range-finder and ray-tracer data—this is be-cause these objects are more familiar to human beings and we can see that the simplifi-cation algorithm tries to preserve shape. Also, they present a varied spectrum of shapes (flat regions, sharp edges, corners, conical and cylindrical regions) to be simplified. Table 4.6 provides some statistics regarding simplification of the 3D example objects. As ex-Rabbit Knot Hand (n = 57931 spheres) (n = 38993 spheres) (n = 31883 spheres) o spheres ratio time spheres ratio time spheres ratio time s n/s (min.) s n/s (min.) s n/s (min.) 0.97 7766 7.45 5.7 2838 13.74 2.2 3391 9.40 1.1 0.95 5048 11.47 6.3 1061 36.75 1.2 2319 13.75 0.75 0.93 3283 17.64 12.1 763 51.10 1.0 1491 21.38 0.51 Table 4.6: Summary of results of the simplification algorithm on 3D data. pected, the algorithm results in fewer spheres (after simplification) in the regions where the object is "fat" (or the curvature is low), e.g., the rabbit's body, and more spheres where the features are "thin", e.g., the rabbit's ears. Figure 4.25f shows the result of rendering the simplified model shown in figure 4.25e as blended spheres (or "blobbies"). To render an object with about 1000 spheres, as blobbies, our ray-tracer takes about 4 minutes on an IRIS Indigo-2 at 512x512 resolution. Figure 4.26f shows the results of displaying only relatively larger spheres using the histogram. As can be seen, these spheres do a good job of modeling the knot in this case. Chapter 4. Simplification of Union of Spheres Model 67 (e) (f) Figure 4.23: Results of the simplification algorithm on 2D brain data, (a) An MRI slice through a human brain, (b) UoC model of the brain (1703 circles), (c-f) Simplified UoC models (dark) overlapped by the original UoC model (light) for a values of 0.97 (627 circles), 0.95 (470 circles), 0.90 (228 circles), and 0.85 (142 circles) respectively. Chapter 4. Simplification of Union of Spheres Model 68 Figure 4.24: Results of the simplification algorithm on 2D "giraffe" object scanned from a book (512x512 pixels), (a) Boundary points (2112 points), (b) UoC model (2110 circles), (c-f) Simplified UoC models (dark) overlapped by the original UoC model (light) for an absolute error of 1.0 pixel (267 circles), 2.0 pixels (157 circles), 3.0 pixels (112 circles), and 4.0 pixels (89 circles) respectively. Chapter 4. Simplification of Union of Spheres Model 69 ( e ) ( D Figure 4.25: Results of the simplification algorithm on 3D rabbit data, (a) Polygonal model, (b) The UoS model (57931 spheres), (c-e) Simplified UoS models for o values of 0.97 (7766 spheres), 0.95 (5048 spheres), and 0.93 (3283 spheres), (f) Blended spheres display for the simplified model in (e). Chapter 4. Simplification of Union of Spheres Model 70 Figure 4.26: Results of the simplification algorithm on 3D knot data, (a) Polygonal model, (b) The UoS model (38993 spheres), (c-e) Simplified UoS models for a values of 0.97 (2838 spheres), 0.95 (1061 spheres), and 0.93 (763 spheres), (f) Only large spheres displayed using the histogram (8063 spheres). Chapter 4. Simplification of Union of Spheres Model 7 1 Figure 4.27: Results of the simplification algorithm on 3D hand data, (a) Polygonal model, (b) The UoS model (31883 spheres), (c-e) Simplified UoS models for a values of 0.97 (3391 spheres), 0.95 (2319 spheres), and 0.93 (1491 spheres), (f) Simplified UoS model for absolute tolerance of 3% of largest sphere radius of 1.70 units (1092 spheres). Chapter 4. Simplification of Union of Spheres Model 72 60000 50000 a 4 0 0 0 0 30000 B 3 Z 20000 10000 1 1 r 0 I 1 1 1 1 1 1 1 1.00 0.99 0.98 0.97 0.96 0.95 0.94 0.93 sphericity Figure 4.28: Sphericity vs. number of spheres for rabbit data. 7000 6000 -C 5000 4000 B 3 Z 3000 2000 1000 0 Ita: 0.00 0.50 1.00 1.50 2.00 radius 2.50 3.00 3.50 4.00 Figure 4.29: Histogram of spheres for rabbit data before simplification. Chapter 4. Simplification of Union of Spheres Model 73 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 radius Figure 4.30: Histogram of spheres for rabbit data after simplification. Figure 4.27f shows the simplification result for hand data for an absolute tolerance of 0.05 units which is about 3% of the largest sphere radius (1.70 units). This results in a lower sphericity of resulting clusters in the fingers and a higher sphericity (i.e., relatively more accuracy) in the palm. Figure 4.28 shows the number of spheres for different values of sphericity for the rabbit data. The model has too much redundancy in the beginning (indicated by the steep portion of the curve) but starts to stabilize for sphericity close to 0.95 (the curve begins to flatten). One interpretation of this is that while in the beginning the number of spheres depends on the number of boundary points, after proper simplification it depends on the shape of the object. Figures 4.29 and 4.30 show the histograms of radius distribution for the UoS model for the rabbit data before and after simplification. Before simplification there are large numbers of spheres of all sizes; this is because the surface of the object is densely sampled in all places. After simplification, there are a relatively large number of small spheres Chapter 4. Simplification of Union of Spheres Model 74 (due to ears and legs) and small number of large spheres (corresponding to spheres in the body and head). 4.4 Summary In this chapter we showed how to simplify the UoS models generated using the algorithms presented in Chapter 3. The algorithm uses a new clustering method based on the concept of sphericity. The algorithm stabilizes the UoS model, removes redundancies from the UoS model, and provides a control over the desired error. We will discuss some applications of the stable UoS models in Chapter 6. Chapter 5 Approximation properties of the Union of Spheres Representation In this chapter we study the approximation properties of the proposed UoS representa-tion. This study is necessary as a guide to when the use of the UoS representation will be appropriate. We will not discuss the approximation properties of the UoS representa-tion in the most general context, but rather in a manner which will show the connection between the boundary point density 5 and the error e in the UoS model obtained from those points. The main result of this chapter is that under certain restrictions on the boundary of the underlying object, the error e in the reconstructed UoS model is 0(5). In other words, the error in the reconstruction depends on how densely the boundary of the object was sampled. For data sets that do not satisfy the imposed constraints on the boundary of objects or the sampling criterion, the error in the reconstructed UoS model may not be bounded. 5.1 General strategy We first discuss the approximation properties of the union of arc representation of curves and relate it to the approximation properties of the polyline representation of a curve. Next, we consider the union of circles (UoC) representation for 2D objects—this study makes use of the union of triangles (UoT) representation of a 2D object. Finally, we derive the error bounds for the union of spheres (UoS) representation of 3D objects— here, we make use of the union of tetrahedra (UoTet) representation of a 3D object. This approach has been taken for the following reasons: 75 Chapter 5. Approximation properties of the Union of Spheres Representation 76 • The approximation properties of the polyline (polygonal in 2D, polyhedral in 3D) representations have been widely studied and we can compare the two representa-tions in terms of the approximation properties. In particular, we will show that in reconstruction from boundary point data, both representations incur the same order of error. • The UoC (UoS) model described in Chapter 3 is actually derived from the UoT (UoTet) model of an object. Therefore, it will be easy to relate the results of this chapter to the results of the previous chapter. In both 2D and 3D our approach is to obtain the error bounds as follows. We first put some constraints on the complexity of the boundary of objects (this is somewhat analogous to the concept of bandlimitedness of a signal in sampling theory). We then define a proper sampling of the boundary of the object and assume that the boundary has been properly sampled. Based on this definition we derive error bounds for the polygonal (polyhedral) description of the object. Finally, we present a construction of the proposed model which shows that the error between the polygonal (polyhedral) model and the corresponding proposed model.is closely related, which provides the appropriate error bounds. To show an error bound of e we will show that Conditions A and B (defined in Section 1.2.6) are satisfied for that value of e. 5.2 Relevant work There is some overlap in the work presented in this chapter and the work of Brandt and Algazi who studied the modeling of 2D objects by their approximate skeletons derived by computing the Voronoi diagram of the boundary points [D.5.3-Brandt and Algazi 92]. In particular, we make use of the concept of r-regular object, which is taken from the field of mathematical morphology [H.13-Serra 82]. In simple terms, an object is r-regular if Chapter 5. Approximation properties of the Union of Spheres Representation 77 it is possible to cover the entire space occupied by the object's interior and exterior by a union of open circles (spheres in 3D) of radius r without touching the object's boundary. See Figure 5.31 for an example. This means that: (i) the curvature of the boundary Figure 5.31: An r-regular object. A n open circle of radius r can freely cover the inside and. outside of the object without ever crossing the boundary. of the object cannot be greater than 1/r anywhere, and (ii) the object cannot have an isthmus or a canal narrower than 2r. Another consequence of this definition is that a circle defined by any three points on the boundary of an r-regular object has a radius of at least r [D.5.3-Brandt and Algazi 92]. Brandt and Algazi showed that for such objects if the sampling density 8 < r, the maximum error in the UoC model obtained from the Delaunay triangulation of boundary points is 0(5) [D.5.3-Brandt and Algazi 92]. Kirkpatrick was the first to note the direct relationship between the skeleton of a poly-gon and the Voronoi diagram of its vertices [D.5.5-Kirkpatrick 79]. Boissonnat noted that if the boundary of a polyhedral object is contained in the Delaunay tetrahedralization of its boundary points (called the "contour containment criterion") then the tetrahedra can be classified into two types: inside or outside [C.3.1-Boissonnat 84]. He also men-tioned, in passing, that the spheres circumscribing the inside tetrahedra can be used to Chapter 5. Approximation properties of the Union of Spheres Representation 78 approximate the skeleton of an object. Goldak et al. used this approach to derive a skeletal model of objects and showed that as the density of boundary points increases indefinitely, the centers of the spheres define exactly the medial surface (or skeleton) of the object [D.5.4-Goldak et al. 91]. Several papers touch on the issue of distribution and placement of points on the surface of the object but do not go into details [D.5.4-Goldak et al. 91, D.5.2-Brandt 92]. Turk discusses this issue with the objective of deriving "good" surface triangulations of objects [C . l . l l -Tu rk 92]. We will not deal with this problem in this thesis but will require that the boundary points satisfy our sampling criterion. 5.3 A p p r o x i m a t i n g curves w i t h arcs of circles This section deals with the approximation of a planar curve from a set of points (samples) on the curve. We begin this discussion by assuming that the curve has been properly sampled. This means that the sampling density is 8, i.e., for any point q on the curve, there is a sample point p such that d£(p ,q ) < 8. One consequence of this is that for any two neighboring sample points p ;_i and pj, d#(Pi-i) p;) < 28. The results for planar curves are based on the following theorem. Theorem 1 If the curve segment Ci, whose end points are two samples pj_i and'Pi, has a length 21 and the distance between pj_i and Pi is 2d, then Ci stays on or within the ellipse x2 y2 aligned with line-segment p7~Tp7 and having o as center, where o is the mid-point between P i _ i and Pi-P r o o f The ellipse is shown in Figure 5.32. Because I2 is greater than (I2 — d2) the Chapter 5. Approximation properties of the Union of Spheres Representation 79 (-d, 0), (-1, 0) 1 \2-d2 ) ^ < 0 0) \(l,0) i-1 (0,0) i-1 (0.-Ji2-d 2 ) Figure 5.32: The error-bounding ellipse. principal foci are on the x-axis as shown. In fact, pj_i and p; are the foci of the ellipse. The length of the semimajor axis of the ellipse is 21. This means that the sum of the distances of any point on the ellipse to the two foci is equal to 21. We will use this fact to prove the theorem by contradiction. Assume that the curve segment c"i crosses the ellipse at some point q. Let q' be one point on c~i that is outside the ellipse. Note that d#(q, q') > 0 because q ^ q'. Then, 21 - dc(pi-i,Pi) = dc(Pi-i,q) + dc(q,^)+dc(qL',Pi) > dE(pi_1,q)-rdE(ci,qL') + dE(q.',pi) > rfjg(pi_i,q) + d E (q ,p j ) (by triangle inequality) > 21 (by definition of the ellipse) which cannot be true. Hence the curve segment cannot go out of the ellipse. Q.E.D. Chapter 5. Approximation properties of the Union of Spheres Representation 80 If we can generate an arc which stays within this ellipse, then Conditions A and B will be satisfied and the error will be bounded by the width of the ellipse (i.e., 2y/l2 — d2). Because d < 8, the error essentially depends on I. If there is no restriction on the length of the curve between any two consecutive samples, by Theorem 1 the error will not be bounded. We will say that a curve is well-behaved and is properly sampled if: • the sampling density, as before, is 8. • for every point on the curve, there is a path of length < 8 along the curve to the nearest sampled point (i.e., / < 8; in general, we require that / be bounded). • as 8 tends to zero, the ratio l/d tends to one. This is necessary to ensure that the curve is modeled exactly in the limit. These conditions imply that for any two consecutive samples Pi-i arid p;, 2d = dE(Pi-i, Pi) < dc(Pi-i, Pi) = 21 < 28. (a) Polyline representation of a curve: One arc that stays within the ellipse is the straight line joining the sample points (arc of circle of infinite radius). Generating such an arc for every segment yields a polyline model of a curve. For such a model, Conditions A and B are satisfied for e = 8. When l/d is close to one, Theorem 1 lets us provide tighter error bounds on curves (e.g., for l/d < 1.1, e =.0.4.8) without putting any extra restriction on the curve. The maximum error in the segment c\ according to Theorem 1 can be yf{l2 - d2). The maximum error for the whole curve is the maximum of these error maxima for different segments of the curve. (b) Union of arc representation of a curve: Another arc within the ellipse is an arc of length exactly equal to 21. This arc is guaranteed to stay within the ellipse accord-ing to Theorem 1. Generating such an arc for every segment of the curve yields a union of Chapter 5. Approximation properties of the Union of Spheres Representation 81 arc model of a curve. It also has an advantage that if the curvature of the curve segment Ci is a constant, the found arc models the curve exactly [I.3-Fournier and Ranjan 95a]. According to Theorem 1 the maximum error in the segment Cj is 0(6). Note that the error in the union of arc model in the worst case cannot be worse than the error in the polyline model because the latter is a special case of the former. Although it is possible to obtain slightly tighter error bounds by putting restrictions on the curve such as—(i) the curve is allowed to be on one side of the line joining the two neighboring samples, or (ii) the curvature k in the curve segment is bounded between kmin and kmax—we do not present the analysis here because the error bound does not change significantly (in particular, the order of the error does not change [I.3-Fournier and Ranjan 95a]). 5.4 Approximating 2D objects with union of circles The boundary of a 2D object consists of one or more closed curves. If the goal is to model only the boundary of the 2D object, the error bounds of the previous section are valid whenever the boundary of the object is well-behaved and is properly sampled. If, however, we are interested in modeling only the interior of the object; (i) we cannot use arcs of low curvature (i.e., the radius of the circle is limited), and (ii) even the concave portions of the object's boundary have to be approximated by arcs of "inside" circles. Figure 5.33 shows the two basic problems that lead to large reconstruction error using our approach. In Figure 5.33a, the boundary of the object is not sampled well enough. This problem is avoided by demanding a sampling density of 6, and that the length of the boundary curve be bounded between any two successive samples. In Figure 5.33b, the underlying object is not r-regular; the features in the objects are too sharp with respect to the sampling. If the dotted circles in Figure 5.33b are classified as inside, the error can be large; if they are classified as outside, we can break the model. This situation Chapter 5. Approximation properties of the Union of Spheres Representation 82 (a) (b) Figure 5.33: Problems can arise in the proposed reconstruction method: (a) if the sam-pling constraint is violated, and (b) if the object is not r-regular. can also occur near the sharp edges and corners of an object. To take care of this, we assume that the objects are r-regular. Note that r can be small, but then so should be the corresponding 6. We will present the error analysis for the case when the underlying object is r-regular, where r > 25. (a) U n i o n of triangles ( U o T ) representation: Consider a union of triangles model of the object derived from the Delaunay triangulation of the boundary points by deleting all triangles for which the center of the circumscribing circle is outside the object1 (Figure 5.34a). For such a model, all the vertices are on the boundary of the object and the unshared edges of triangles form the boundary. The length of each boundary edge is less than 26. Together, they give a polyline model of the boundary curve and the error bounds for such a representation were derived in the last section. Figure 5.34b shows the error ellipse2 associated with one segment of the boundary. (b) U n i o n of circles ( U o C ) representation: Now consider a union of circles 1 Recall that from our assumptions this inside-outside test is available to us. 2 Brandt and Algazi showed that for an r-regular object, the boundary curve between any two suc-cessive samples must be contained in the lune formed by the intersection of two circles of radius r intersecting at the sample points [D.5.3-Brandt and Algazi 92]. Chapter 5. Approximation properties of the Union of Spheres Representation 83 Figure 5.34: Modeling planar objects as union of triangles. model obtained by taking the set union of circles circumscribing all the "inside" Delaunay triangles. By the empty circle property, the circles do not contain any other sampled boundary point. We will show that under the imposed restrictions on the object and its sampling, the error incurred by such a construction is 0(5). In 2D, notice that all the boundary edges in the UoT model are less than 25, and that the boundary of the UoC model is denned by arcs of circles "poking out" from such edges. When the radius of the circumscribing circle is 25, the error between the UoT model and the UoC model is 0.26 5 (figure 5.35). As this circle radius increases (or alternatively, 6 r . decreases), the two models get closer and closer. Figure 5.35 shows what happens at a concave portion of the object as the sampling density 5 is increased. Since the error in the UoT model is 0(5), and the error between the UoT and UoC models is also 0(5), the overall error in the UoC model is also 0(5). Now we will present a more rigorous proof that the error in the reconstructed UoC model is bounded. Chapter 5. Approximation properties of the Union of Spheres Representation 84 Figure 5.35: Relation between UoT and UoC models. Theorem 2 / / the underlying object being reconstructed is r-regular, and the boundary point density condition of 8 < r holds, the error in the UoC model obtained from the Delaunay triangulation of the boundary points is less than or equal to 8. Proof: We need to show that Conditions A and B are satisfied for e = 5. C o n d i t i o n A : Since all sampled points lie on the boundary of the UoC model, Condition A is satisfied for e = S by the definition of the sampling density. C o n d i t i o n B : To show that for any point q on the boundary of the UoC model, there is a point p belonging to the boundary of the object and within a distance e from q, we recall that the boundary of the UoC model is formed by the arcs of the circumscribing circles that are not contained in the UoT model. Consider a point q on the boundary of the UoC model. There are three cases: (i) q is on the boundary of the object, (ii) q is outside the boundary of the object, and (iii) q is inside the boundary of the object. For (i) Condition B is obviously satisfied. For (ii), consider the line joining q to the center of Chapter 5. Approximation properties of the Union of Spheres Representation 85 (a) (b) Figure 5.36: Figure to assist in the proof of Theorem 2. (a) q is outside the boundary of the object, (b) q is inside the boundary of the object. the circle (Figure 5.36a). If it intersects the boundary within 6 from q, we have found a point p such that d(p, q) < 5. If it intersects at a distance greater than 6, we have found a point p for which the sampling condition fails (because the circle cannot contain any sampled boundary point). For (iii), consider a circle of radius 8 centered at q (Figure 5.36b). If it intersects the boundary of the object, we have found a point p satisfying the Condition B. If it is totally contained in the object, the sampling condition is violated between the two boundary points defining the arc on which q lies. Q.E.D. Chapter 5. Approximation properties of the Union of Spheres Representation 86 5.5 Approximating 3 D objects with union of spheres Analogous to.the 2D case, unless we restrict the objects to be r-regular, and the sampling density, there can be arbitrarily large error in the UoS model constructed from the Delaunay tetrahedralization. In our discussion, we will assume both of these conditions, and that 5 < r. Note that r-regularity implies that any sphere defined by four boundary points will have a radius of at least r. (a) Union of tetrahedra (UoTet) representation: Consider a union of tetrahe-dra model of a 3D object derived by a Delaunay tetrahedralization of boundary points followed by an inside-outside classification of the tetrahedra based on the centers of the circumscribing spheres. For r-regular objects, it leads to a UoTet model such that all the vertices are on the boundary of the UoTet model and the unshared faces of tetrahe-dra form the boundary. If all the boundary faces have edges of length 0(5), then it is straight-forward to extend the arguments for the 2D case to show that the error bound e is 0(5). However, unlike the 2D case, the sampling condition does not guarantee that all boundary edges < 25. Although several papers claim that the error in the UoTet model is bounded provided 5 is small enough [C.3.1-Boissonnat 84], it is still an open question to show that rigorously. The empty sphere property of Delaunay tetrahedralization may be useful in this regard. (b) Union of spheres (UoS) representation: Consider a union of spheres model obtained from the spheres circumscribing all the tetrahedra in the UoTet model. The error incurred by such a construction is 0(5) if the error incurred by the UoTet repre-sentation is 0(5). The argument for the 2D case holds exactly by replacing "circle" with "sphere", and "triangle" with "tetrahedron." For the UoS model derived from volumet-ric data, we can show that the error is bounded without putting any constraint on the length of the edges of the boundary faces in the UoTet model. Chapter 5. Approximation properties of the Union of Spheres Representation E r r o r bound for U o S model from volumetr ic data 87 Theorem 3 Let there be a volumetric sampling of an r-regular object O on a grid of cell size (v x v x v), where each voxel has been classified as "inside" or "outside" the object. Let the boundary of the object be sampled by computing boundary points between all neighboring (in the 6-connected sense) inside and outside pairs of voxels. If r > the error in the UoS model derived by our construction is less than or equal to \^3u. Proof: To show the error bound we first note that if r > for our sampling procedure, the sampling density 8 < y/3v. This is due to the facts that every sphere of radius > r necessarily contains at least one voxel, and that for an r-regular object it is possible o to put a sphere of radius r completely inside (or completely outside) the object and touching the boundary at any boundary point. To show that 8 < \J?>v consider any point q G B(0) and the cell containing q. There are three possible voxel configurations for the eight voxels defining this cell: (i) some voxels are inside, some outside (i.e., cell is of mixed-type), (ii) all voxels are inside the object, and (iii) all voxels are outside the object. For case (i), since we compute a boundary point between all neighboring inside-outside pairs of voxels, there is definitely a sampled boundary point within \[Zv of q. For case (ii), consider a sphere of radius r completely outside the object and touching the object's boundary at q. This sphere contains at least one voxel defining the neighboring cell, and because the sphere is completely outside the object, this voxel has to be outside the object. Therefore the neighboring cell is of mixed type and it must define a boundary point and this boundary point cannot be further than \[Zv from q. Case (iii) is similar to case (ii)—just the role of inside and outside are reversed. Now we will show that Conditions A and B are satisfied for e = y/3v. Chapter 5. Approximation properties of the Union of Spheres Representation 88 C o n d i t i o n A : This is satisfied in our model for e = \/Zv > 8 because the sampled boundary points lie on the surface of the UoS model also. C o n d i t i o n B : To show that Condition B is satisfied for e = \/oV, we use an argument similar to the one presented above to show that 8 < \[Zv. Consider a point q G B(7l) and the cell that contains it. There are three possibile voxel configurations for the eight voxels defining this cell: (i) some voxels are inside, some outside, (ii) all voxels are inside the object, and (iii) all voxels are outside the object. For case (i) there is necessarily a sampled boundary point within y/Zv from q. If (ii) is true, q cannot be a boundary point (contradiction). For (iii), we note that the sphere on which q lies has a radius greater than or equal to r (all spheres generated from samples on the boundary of an r-regular object have a radius greater than r), therefore it contains at least one of the voxels in the neighboring cell. This voxel has to be inside the object because each sphere contains either only inside voxels or only outside voxels and this sphere was defined to be inside the object based on the type (inside or outside) of voxels it contains. Again, in ' this situation there must be a sampled boundary point within y/3i/ of q. From the above arguments we conclude that the error in the UoS model is less than or equal to y/Zu. Q.E.D. One practical use of this theorem is that if we know that an object is r-regular, for the error in the reconstruction of its UoS model to be bounded, the object should be sampled with a linear grid resolution of v < E r r o r bound for simplified U o S model Finally, we provide an upper bound for the error in the simplified UoS model. Chapter 5. Approximation properties of the Union of Spheres Representation 89 Observat ion 3 Let S denote the set of spheres in the original model for object O. Let the error in U(S) be bounded by e. Let S' denote the set of spheres in the simplified model. If S is simplified for a sphericity threshold of a to obtain S', then the error bound forU(S') < e + 2r^(l — <j)/a, where is the radius of the largest sphere in S. Consider a cluster c G S'. Let G be the set of spheres from S contained in c. We first show that c introduces an error of at most 2rc( l — a)/a, where ra is the radius of the largest sphere in G. The worst case scenario for c in terms of error is shown in Figure Figure 5.37: Error in a cluster. 5.37. In this situation, the error in c is 2{rc — rc). This is maximum for the largest value of rc, which is r ^ / a , so maximum error in c (emax in Figure 5.37) is 2 rc ( l — o)ja. Therefore, in the simplification of the whole set S, the Hausdorff distance between U(S) and U(S') is bounded by 2TL(1 — o~)/a, where TL is the largest sphere in S. Since Conditions A and B are satisfied between O and U(S) for an error bound of e (given), and between U(S) and U(S') for an error bound of 2TT,(\ — o)jo (shown above), we conclude that for any e'"> e + 2rx,(l — cr)/a, the Conditions A and B will be satisfied between O and U(S'). Note that the actual Hausdorff distance between O and U(S') Chapter 5. Approximation properties of the Union of Spheres Representation 90 may be much less than e'. 5.6 Summary In this chapter we looked at the approximation properties of the union of arc representa-tion of curves,, union of circles (UoC) representation of 2D objects and union of spheres (UoS) representation of 3D objects. If the underlying object is r-regular, the error in the reconstruction is 0(5), where 5 is the sampling density, provided 5 < r. UoS may not be suited for compactly modeling flat portions of an object or sharp edges and corners with high accuracy. However, for smooth curved objects with not many sharp features, the UoS model yields a compact model for a reasonable error tolerance. In addition, we will see in the next chapter that it is stable enough to be used in problems involving shape comparison. An interesting question is how local bounds on sampling density and curvature affect the error bound in the reconstruction. This could be explored in the future. Chapter 6 Matching Objects Using Union of Spheres Models The notion of shape is both intuitive and elusive. We all think we know what shape is, but when asked to define the notion precisely, we have great trouble avoiding circular definitions. Similarly, it has been difficult to capture the concept of shape in programs that would recognize or compare objects from either images or three dimensional models. In this chapter we show the potential for using the UoS (UoC in 2D) model of objects for this purpose. Given two UoS models A and B, we first deal with the problem of matching spheres from one UoS model to the other. We use the matches to show experimentally that the UoS representation is stable with respect to some transformations and also to boundary noise. We use the UoS model to compute distances between objects, and to extract transformations between objects for object registration and later for the problem of object shape interpolations. The results obtained demonstrate that the UoS representation is indeed stable and that it captures, to some extent, the intuitive notion of shape. In all of the discussion, we assume that the UoS models of the objects are available to us. Again, for clarity of illustration, we will describe some algorithms in 2D but show the results both for 2D and 3D. Whenever the extension of an algorithm to three dimensions is non-trivial, the difficulties will be pointed out along with the proposed solution. Figures 6.38 and 6.39 show some of the objects that will be used in our examples. 91 Figure 6.38: Simplified UoC (239 circles), (c-d) A giraffe circles). models of some objects used in this chapter, (a-b) A cow (370 circles), (e-f) Fish (1) (45 circles), (g-h) Fish (2) (43 Figure 6.39: Simplified UoS models of some objects used in this chapter, (a) A rabbit (1046 spheres), (b) A head (262 spheres), (c) Hand (1) (1086 spheres), (d) Hand (2) (1087 spheres). Chapter 6. Matching Objects Using Union of Spheres Models 94 6.1 Matching problem: a perspective 6.1.1 The problem The problems of object registration, shape interpolation, object classification, or object recognition are all closely related. They all involve matching two models of the same or two different objects. 1. When A and B model the same object with a transformation between the two images, the aligning of the two models is referred to as registration. The transfor-mation can be affine or non-affine. This is a common situation in medical imaging, with images produced at different times, and/or through different imaging geome-tries. i 2. When A and B model different objects, the aligning and interpolation of one into another is a problem of shape interpolation. This can again be useful in the bio-medical field, for example, in visualizing the growth of an organ, or of a tumor. 3. When A is a model of an object that needs to be classified as one from a set (B\, B2, • •. Bn) of n objects, the problem is that of object classification. 4. When A is a model of an unknown object, and B is a template of a known object, the problem is that of object recognition. This is similar to 3 except that a threshold is used to decide whether A and B model the same object" or not. We discuss all of these problems in this chapter, but will focus mainly on the first two. The problem that we mostly want to address in this chapter is how to obtain a stable matching, and how to obtain as automatic a process as possible or required. Chapter 6. Matching Objects Using Union of Spheres Models 95 6.2 General strategies There are a number of ways of comparing shapes through different representations of objects [H.3-Duda and Hart 73, A.5-Hebert et al. 94]. We will not attempt to present an exhaustive survey of these techniques in this thesis, but merely summarize. We can identify many ways of matching objects using their computer models: 1. By comparing global properties. This might involve, for example, computing the surface area, volume, and/or moments or fitting a single volume primitive to the model. 2. By establishing local correspondence and comparing properties in that local region. This might involve feature detection and matching between models. Features can be points, edges, corners, etc. 3. If the object is too complex, by segmenting the object into parts and using ap-proaches 1 or 2 applied to the parts. 4. Hierarchical descriptions of shapes or decomposition of shapes using some "ba-sis" functions (analogous to Fourier decomposition of a signal) could be used. See for example the work of Hughes, Ghosh et al. and Muraki ([F.5-Hughes 92], [G.lO-Ghosh and Jain 93], [C.2.3-Muraki 93]). , -5. Other methods include reducing the dimension by using skeletons [H.3-Duda and Hart 73], [H.l-Bookstein 91]. Generally, the first method is rather easy to apply, but has limited flexibility and gener-ality. It can, however, be used as a first step for more specific matching, as we will indeed do in our case. The third approach, in a way, circumvents the problem, since an effective segmentation implies that some matching has already taken place. In particular, if the Chapter 6. Matching Objects Using Union of Spheres Models 96 segmentation requires human intervention, that effectively means that the human does some or most of the shape recognition, and we know that we (humans) can do that. 6.2.1 Image registrat ion The problem of registration can be defined as: given two data sets of the same object, and a class of transformations T , find the transformation from one to the other. The amount of research on registration is also considerable and an exhaustive survey of techniques is beyond the scope of this thesis. It is hard to even distinguish between the image warping techniques and image registration; there is a significant overlap in the objectives and therefore the techniques. The reader is referred to [A.3-Brown 92] for a good survey of registration methods. There is certainly no unique solution to all registration problems. This is mainly because there are many known (some of which can be modeled, some cannot) and unknown reasons for mis-registration to begin with. To find a provably "optimal" transformation, which minimizes some difference measure between the two models after T has been applied to one is hard (at least NP-complete [A.3-Brown 92]). Therefore, almost always, the problem is not formulated as a global minimization problem but rather as a problem of finding a good match. Most registration algorithms have three parts [A.3-Brown 92]: • segmentation: where the features to be used for matching the two data sets are detected. A wide variety of features have been used such as points (voxels), lines (edges), curves (contours), etc. • matching: where the features are matched and the transformation parameters are estimated using the matches. Again, there are a number of ways of matching the detected features (from manual to semi-automatic to totally automatic). Fur-thermore, several ways of finding the transformation from the matches have been Chapter 6. Matching Objects Using Union of Spheres Models 97 . employed, and several classes of transformations have been used. • verification: where the quality of the registration is assessed. This implies a. measure for the difference between the two models, and this is difficult to define in a way not too influenced by the matching process. 6.2.2 Shape interpolat ion Shape interpolation is the process of continuously transforming one object model into another. In addition to its aesthetic appeal (e.g., in "morphing" or metamorphosis be-tween objects [F.l-Beier and Neely 92]), it can be used in visualizing the evolution of objects (e.g., growth of a tumor), or in automatic in-betweening in key-frame anima-tion [F.lO-Reeves 81]. It is difficult to assess objectively the quality of interpolation algorithms in terms of the intermediate shapes that are generated, because it is already difficult to compare shapes themselves. Generally, the criteria are application-dependent. For example, in animation it may be the aesthetic appeal of the interpolation and the consistency of the intermediary shapes (to read an early discussion of these, see the pa-per by Catmull [G.5-Catmull 78]); in visualizing the growth of a tumor, it might be the plausibility of the intermediate shapes and of the evolution of the phenomenon. However, it may be. useful to list a few desirable properties (not necessarily all required at the same time): • The boundary of the interpolated models should remain smooth and continuous, without tears or self-intersections. • There should not be unnecessary distortions. This means, for example, that a large volume of one object should not transform into a small volume of the other. Chapter 6. Matching Objects Using Union of Spheres Models 98 • The algorithm should not generate too many intermediate primitives (only of the order of the number of primitives in the models). • It should allow some (intuitive) control over the interpolation. The user should be able to select or modify some aspects of the interpolation. • It should preserve as much as possible of what is common between the two models. For instance if A and B have the same topology, the intermediate shapes should have it too; if A and B are convex, so should be all the in-between models. • Neighborhood and context should be preserved. If two matched pairs of primitives are adjacent, they should remain so in the interpolation. • It should allow other attributes such as color or texture to be interpolated whenever possible (sometimes different shape or topology make it impossible). • It should be efficient and easy to implement. Essentially, there are three classes of algorithms to perform shape interpolation. Some algorithms first establish a correspondence between the primitives of the two models (often after bringing them to some canonical position) and then interpolate between the matched primitives [F.4-Hong et al. 88, F.7-Kent et al.]. The second class of algorithms does not match explicitly, but uses a global transformation acting on all the primitives at once [F.5-Hughes 92, F.9-Payne and Toga 92]. The third class of algorithms uses specialized matching primitives. This is probably the oldest and most popular way;, the primitives can be strokes, skeletons, members of an articulated model, or segments drawn by the user who determines the matching. This has been used in key-frame interpolation [F.2-Burtnyk and Wein 76], in-betweening [F.lO-Reeves 81] and morphing [F.l-Beier and Neely 92]. The main difficulty with the first class is that the matching Chapter 6. Matching Objects Using Union of Spheres Models 99 process is difficult and can be unstable (a small change in the input data triggers a large change in the primitives and their matches). The main difficulty with the second is that there is little control over the intermediary shapes. None of these.approaches deal with "shape" primitives, and as a result there is little control over the result, except in the last category, where the control is obtained by direct human intervention in defining and/or matching the primitives. 6.3 Matching with UoS Models In this section, we present a method to establish the matches between the spheres of A and spheres of B. Let A be modeled by a set of spheres, {ai , a 2 , . . . , a n } and B by {bi , b 2 , . . . , bfc}. The number of spheres n and k, and the topologies of A and B can be different. Without loss of generality, we assume n < k. Our approach can be summarized by the' following steps: 1. Compute a UoS model for the two objects. 2. Simplify the model, to stabilize it and remove redundancies. 3. Align the models globally if needed (optional). 4. (a) Segment the objects into parts if necessary (optional). (b) Match parts between objects if they have been segmented. 5. Match spheres in matched parts.-(a) Compute distances between spheres in the two models. (b) Match spheres based on the distances. 6. Use the matches in a specific application. Chapter 6. Matching Objects Using Union of Spheres Models 100 Steps 1 and 2 have been described in Chapters 3 and 4 respectively. Steps 3, 4, and 5 will be described in this section. Specific applications of matching (step 6) will be described in the next section. 6.3.1 Step 3 : Model alignment It is useful to bring the objects to some canonical size and position before matching. This is the purpose of step 3, and can be manual or automatic. In the manual case, the user is allowed to uniformly scale, rotate, and translate the two models in the world space. In the absence of user intervention, we use the method of central moments to align the two models [C.3.6-Solina and Bajcsy 90]. The models are translated such that their centers of masses match, they are scaled such that their volumes1 match and they are rotated such that their local axes established from second order moments match. To estimate the moments and volumes we lay a 3D grid over the two models keeping track of grid points that are inside each model. The center of masses of the models are estimated by taking the average of all the corresponding m interior grid points. j m m m x = —Y,xi, y=—J2vi, z ^ — ^ Z i m~1 m~\ m~1. i=i i=i i=i The eigenvectors of the following "covariance" matrix M are used to attach a local coordinate system to the models. -i m mi=i (Vi - V)2 + (zi - z)2 (yi - y) (x{ - x) (zi - z) (x{ - x) (xi-x)(yi-y) (x{ - x)2 + (z{ - z)2 (zi-z)(yi-y) (xi - x)(zi - z) ' (yt - y)(z{ - z) (x{ - x)2 + (y, - y)2 Because the matrix M is symmetric, its eigenvectors are orthogonal. The eigenvalues of the matrix are used to put an ordering on the three eigenvectors. The z-axis is aligned 1 Sometimes it is more convenient or useful to scale models by their linear dimension instead of volume. One such approach is discussed in Section 6.4.2. Chapter 6. Matching Objects Using Union of Spheres Models 101 with the eigenvector corresponding to the smallest eigenvalue (axis of least inertia), y-axis to the next, and rr-axis to the axis of largest inertia. This approach works well for objects that have similar shapes, because they have similar moments. However, when the two objects have different shape, there are instances where the moments are not sufficient to obtain a good preliminary alignment of the models (one that looks reasonable to a human observer). In this case the user can scale, rotate, and translate the models manually before the matching step. 6.3.2 Step 4: Segmentation The U o C model allows us to compute semi-automatically an effective segmentation, al-beit wi th some heuristics. The first step is to extract skeletons from the U o C model. This problem has been studied before, in particular by A t t a l i [D.5 .1-At ta l i and Montanvert 94], following earlier work by Brandt [D.5.2-Brandt 92], Kirkpat r ick [D.5.5-Kirkpatr ick 79] and Boissonnat [C.3.1-Boissonnat 84]. Our approach is to compute first the power d i -agram [D.3.1-Aurenhammer 87] for a simplified version of the U o C model. The dual of the power diagram is a triangulation (the regular triangulation), where an edge ex-ists between two circles that are neighbors in the power diagram. The next goal is to eliminate edges in this triangulation unti l most circles have only two neighbors, wi th the ones at remaining branches of the skeletons having three neighbors, and the circles at the extremities (leaves) of the skeleton having only one neighbor. Since we do not want edges between non-intersecting circles, the obvious candidates for elimination, are the edges between circles that do not intersect, even though they are neighbors in the power diagram. The next candidates are the skinny triangles (i.e., triangles wi th one angle close to 180 degrees), which are likely to link circles that are already linked more directly. In this case, the largest edge of the skinny triangle is eliminated. These are fairly simple steps, and do not by themselves guarantee that the output w i l l be a tree but do produce Chapter 6. Matching Objects Using Union of Spheres Models 102 Figure 6.40: Segmentation of objects using UoC models. (a,c,e,g) Regular triangulation of circles in the UoC model of objects. (b,d,f,h) "Skeleton" extracted from the regular triangulation along with the circles. Chapter 6. Matching Objects Using Union of Spheres Models 103 a reasonable segmentation most of the time. Figure 6.40 shows the result with the four objects already introduced. Once the skeleton is determined, one can form the segments by grouping the circles with only two neighbors, and stopping when either a dead-end or a three-way branch is found. The circles found along the way are marked, and the process is continued until there are no unmarked circles left. The last step is then to match the segments between the two objects. This can be done by hand and that is what we have done in our examples. A n automatic process could be used, with the same approach that we will present for the matching of the primitives (see below), given a measure of the distance between two segments. In 3D, the skeleton of an object can have surfaces as well as curves. So in the continuous version of the skeletal surface, each skeletal point has an infinite neighborhood of skeletal points (unlike in 2D, where the skeleton is made of curves and each skeletal point has a finite number of neighbors). It is not clear how many neighbors should be used for each skeletal point in the discrete version of a skeletal surface. This makes the automatic segmentation procedure somewhat difficult, although the direct extension of the 2D approach works effectively in the cylindrical and conical sections of the object. It should be stressed that segmentation was introduced here in order to determine if it really improved the quality of the results, not to present a new, improved segmentation technique. x The next step, matching spheres, is central to this chapter, and will be described in detail' in the following sections. Chapter 6. Matching Objects Using Union of Spheres Models 104 6.3.3 Step 5: Match spheres in matched parts Step 5a: Distance computation After computation of UoS model and alignment, we define and calculate the distances d(a, b) between every sphere a G A and b G B. The distance between two spheres is a function of their location (coordinate of centers after alignment, (xa,ya,za) and (xb, yb, Zb)), size (radius after alignment, ra and rb) and the difference between their characteristic features. Specifically: d(a, b) = wpdp(a, b) 4- wsds(a., b) + w/df(a, b) where, dp(a, b) = ( ( x a - xb)2 + (j/0 - yb)2 + (za - zb)2) di(a,b) = ( r a - r b ) e d/(a, b) will be explained in the next sub-section. wp, ws, and Wf are adjustable weights, e is an exponent (2 if the areas of circles are to be emphasized, 1 if we want to consider linear scales). We chose e = 2 to ensure that the distance rf(a, b) is dimensionalfy consistent. By varying wp, ws, and wj we can assign different weights to the size, location and feature of the primitives. For example, when ws = 0 and Wf = 0 it uses position alone, and when wp = 0 and Wf = 0 it uses size alone. This function strictly depends on the geometry of the circles. Note that d p(,), ds(,), and df(,) will be normalized so that by default we can set wp = ws = Wf = 1, and have all three components have equivalent weights (see Section 6.4.2). The units of our distance measure are units of distance squared (e.g., square centime-ters, or square inches). This measure emphasizes differences .for both 2D and 3D objects. If we treat the three components of our distance measure as three components of an Chapter 6. Matching Objects Using Union of Spheres Models 105 orthogonal 3-dimensional space, the formula for total distance is analogous to computing the Euclidean distance2 squared between two points in this space. Using squared units makes the mathematics a bit easier and implementation a bit faster. Features and distances between features We now describe the computation of dj(a, b). We will explain it in four parts: (a) the concept of radius field, (b) how to estimate the gradient in the radius field from the UoS model, (c) how to calculate features from this gradient, and (d) how to calculate distance between features. (a) Radius field For every closed object, it is possible to assign to every point in space a value which is its signed distance (positive if the point is inside the object, negative if outside) to the nearest boundary point. This is the well-known concept of the Euclidean distance field—we prefer to call it the radius field because the value at any point is the radius of the largest sphere centered at that point and not containing any boundary point. The gradient of the scalar radius field dr/dx at any point captures how the radius of the largest such sphere varies in space, where dx is the distance in a specific direction, and dr is the variation in the value of the radius field in that direction. The radius field and its gradient have been studied in computational vision in the context of their ability to characterize how the shape of an object is changing at any point. The concept is closely related to the concept of skeleton or medial axis (surface) of an object. Several papers describe how ridge-tracking in the radius field (generally discrete) can be used to approximate the skeleton of an object. The direction of the ridge at any ridge-point 2Technically speaking, our. distance measure and its components are not distances; they are pseudo-distances because they do not obey the triangle inequality, but for convenience we will continue using the term "distance" for them. Chapter 6. Matching Objects Using Union of Spheres Models 106 (a) (b) (c) Figure 6.41: (a) Cylindrical portion of an object, (b) A corner, (c) A branching structure. Solid lines show the boundaries, dotted circles the UoC model, and the dashed lines the "ridges". is represented by the directions of maxima in the signed gradient of radius field at that point. The radius gradient along a ridge in the radius field is critical in identifying the features or immediate context (defined more precisely in the next sub-section) of optimal spheres in an object. For example, Figure 6.41a shows the direction of the ridge for a cylindrical section of an object (dr/dx is 0 along this ridge which is shown as a dashed line). Figure 6.41b shows a corner of an object—dr/dx along this ridge actually can estimate the corner angle or how sharp the corner is (to be precise, a = &icsm(dr / dx)). Finally, Figure 6.41c shows a branching structure. At the junction, three ridges meet, so if dr/dx is considered in every direction it will have three local maxima corresponding to the three ridge directions. Chapter 6. Matching Objects Using Union of Spheres Models 107 (b) Radius gradient from UoC and UoS models We will first describe the computation of the radius gradient for 2D objects. As mentioned before, the features that we use are derived from the gradient in the radius field. There are several possible definitions of radius gradient in the case of a UoC model. In general, we want a measure related to the distance field gradient, where dr/dx is defined for every direction. In that case, because the centers of the Delaunay circles are maxima in that field, the values of the gradients so defined would always be negative. The gradient will be maximum (in signed value) when moving towards the neighbors of the circle considered, and minimum when moving towards the boundary points which define the circle. In the 2D Delaunay case, we have a maximum of three neighbors, and in the 3D case a maximum of four neighbors. For our purpose, it is better to have a small number of discrete values, so we will consider only the values in the direction of the neighboring circles. One possible definition for the gradient is to consider the variation of the radius of the circle that is totally inside the model as we move from the center of the circle considered to the center of the next circle (Figure 6.42a). The radius of this circle is smaller than the radius of either circle as their intersection is smaller than any of them. The advantage of such a definition is that it would more accurately represent the case where the size of intersection (defined by the line segment joining the two intersection points of the two circles in 2D, and by the diameter of the intersecting circle in 3D) is quite a bit smaller than any of the two circles. The main drawback is that the gradient is always negative, and that it is more local a characteristic than we would like: consider for instance a conical object, this would not adequately capture the fact that the radii of the circles increase regularly in one direction and decrease at the same rate in the other direction. We define the gradient as the ratio of the difference in radii of neighboring circles to Chapter 6. Matching Objects Using Union of Spheres Models 108 Cirdie 1 (a) (b) Gradient computed for circle 1 in direction of circle 2 Gradient computed for circle 1 in direction of circle 2 dr/dx = (r — r})/x1 dr/dx = ( r i - r 1 y x Figure 6.42: Gradient in direction of neighbor of circle. the distance between their centers (Figure 6.42b). In general, the gradient according to this definition is in the range (—oo, oo) but in our case because a circle does not contain any of its neighbors, the value of dr/dx is limited to the range (—1,1). This follows from the relation r i < r2 + x, which is true if circle 1 does not contain circle 2 (see Figure 6.42b). Our definition does differentiate between the cases of increasing or decreasing radii between the neighbors but at the expense of misrepresenting cases where there is a significant narrowing on the model between two adjacent circles. The justification here is that, in general, the UoC model will have one or more circles of a diameter commensurate with the narrowing. (c) Features from radius gradient Since each circle in the original model is circumscribing a Delaunay triangle, this can be used effectively to characterize the immediate context of the circle. If we consider from the center of the circle in the direction of its (at most) three neighboring circles, the gradient dr/dx, where dr is the signed difference between the radius of the neighboring Chapter 6. Matching Objects Using Union of Spheres Models 109 circle and the circle we consider, and dx is the distance (unsigned) between the centers of the circles (Figure 6.42), we obtain a feature characteristic of the circle and its neighbors, formed of three values in three directions (see Figure 6.43). If the circle has less than O N e g a t i v e Figure 6.43: Feature of a circle. three neighbors (for circles whose triangles have one or more sides on the boundary) the value for missing neighbors is set at the minimum of —1, since it means that in this direction the neighboring circle shrinks to 0 for any distance moved. Because we want to match the circles depending on the direction which are not blocked (that is where the neighbors are not shrinking too fast), we map the dr/dx value to a (0, 2) range, where —1 is mapped to 0, 0 to 1 and +1 is mapped to 2. The asymmetry introduced in the range is legitimate because bigger circles are more important (are less sensitive to noise and errors). The range (—1,1) can be mapped to the range (0,2) by a direct shift (i.e., adding 1) but we used the mapping (1 + sin(arctan(dr/cfo))) which maps the range (-co, +oo) to the range (0,2) so that if we allow overlapping circles later the mapping will still work. The range (—1,1) is mapped to approximately (0.3,1.7) according to this mapping. This is slightly better than a direct shift also because we do not lose all the information in the Chapter 6. Matching Objects Using Union of Spheres Models 110 directions of gradient—1. In a simplified UoC model, the simplification algorithm used might not preserve the Delaunay properties, so the circles can have more than three neighbors. In that case we use the power diagram of the circles in the UoC model [D.3.1-Aurenhammer 87]. This data structure is to circles what the Voronoi diagram is to points, and once computed gives easily all the neighbors of a circle. We then take the three largest circles of the neighbors (if more than three) and compute the feature as before. The choice of three "most significant" neighbors is rather arbitrary, the justification being that in the original UoC model each circle has exactly three neighbors, and that the matching of features is more manageable. (d) Distance between features Once the feature of each circle is computed, it remains to compute a distance between the features of two circles. To define the distance we use a physical analogy. Imagine that the two'features have a common center, are free to rotate around it, and that between the extremities of each pair of "branches" there is a spring pulling with a force proportional to the distance between the extremities (Hooke's law). The system will be at rest when all the forces sum to 0, and the most stable equilibrium will be when the potential energy is at a minimum. From Hooke's law, the potential energy of a spring = |/ec 2, where K is the force constant of the spring, and x is its extension. So, assuming the same force constants for the springs, the total potential energy of the system will be minimized when tfj is minimum, where the summation is over the matched pairs, and the distances between the extremities of the matched pairs. At that point the residual values of the sum of the distances squared is taken to be the distance between the two features. This has the characteristics of a pseudo-distance: it is non-negative and it is zero if the features are identical, but it does not obey the triangle inequality. Chapter 6. Matching Objects Using Union of Spheres Models 111 Rot = 3.14 rad, Dist = 0.0130 Rot = 5.75 rad, Dist = 0.8022 Rot = 2.09 rad, Dist = 0.0206 Rot = 5.23 rad, Dist = 0.0089 Rot = 5.41 rad, Dist = 0.7506 Rot = 0.00 rad, Dist = 0.0095 Figure 6.44: Feature distances Chapter 6. Matching Objects Using Union of Spheres Models 112 The function d/(a, b) is computed as: (a, b) = minimum over all rotations of SJ=I xij2 Figure 6.44 gives some examples of pairs of features and distances between them. See Figure 6.56 in Section 5.4.4 for examples of circles with matched features. The relative position of the two features at the minimum also gives the rotation that brings the.first feature to the best match to the second. One therefore obtains a translation from the distance between centers of circles pairs, a scaling from the ratio of their radii, and a rotation from their features. This is enough to completely specify a similarity that would transform one circle to another while preserving their local contexts. Step 5 b : M a t c h i n g of spheres One way to match 3 spheres is to pick the nearest neighbor match in the other set for every sphere; this may be fine for certain applications. However, several spheres of A might get matched to a single sphere of B, which should be avoided when there are no redundant spheres in the two models. To avoid this situation, we define the problem as follows: Find a maximum match such that: d(A, B) = J2&eA,beB d(a., b) is minimum. We define a bipartite graph whose nodes correspond to the spheres in A and B respectively, and the weight on the edges are the distances between them (Figure 6.45a). A maximum match in a bipartite graph means that all n nodes of A should get matched to exactly n nodes of B. This problem has been widely studied and can be solved exactly in 0(k3) time by reducing it to a network flow problem [H.14-Tarjan 83], where 3 We use the term "match" to imply pairing of spheres from object A to object B. Although we use maximum matching as one of the steps in this pairing process, the pairing that we finally establish is not a "match" in the graph theoretical sense. Chapter 6. Matching Objects Using Union of Spheres Models 113 Figure 6.45: (a) Bipartite graph after distances are assigned, (b) Bipartite graph after matching. Chapter 6. Matching Objects Using Union of Spheres Models 114 ,2k is the number of nodes in the bipartite graph (in our case, k is the larger number of spheres in the simplified models). The remaining k — n spheres of B can be matched in different ways. In most cases matching them to the sphere in A that matches their nearest neighbor in B has worked out to be a good strategy (in some sense, this tries to preserve local context). After matching, every sphere of A is matched to at least one sphere of B, and every sphere of B is matched to exactly one sphere of A (Figure 6.45b). If the UoS model has been simplified, the spheres in the original model can be matched as a function lof their positions within their representative cluster. The relative weights of position, size, and feature distances in the distance measure can be used to guide the desired matching. For example, a relatively big weight for feature distance will ensure that the spheres with similar local context in the two mod-els will get matched independent of their size, position, and feature orientation (since distance between features is orientation independent). See Figure 6.56 for an example. A large weight for only position distance means that the information about size and features plays little role in determining the matches; only the proximity of spheres is important. Generally speaking, there is no golden rule to determine the weights for a desired matching (the problem indeed is interesting and could be investigated in future). In our experiments, we used equal weights for the three distances initially, and then fine-tuned the weights to get the desired matches. 6.4 Applications of matching In this section we will discuss four applications of the matching procedure described above: (i) object registration, (ii) measuring difference in the shape of objects, (iii) showing the stability of the UoS representation, and (iv) shape interpolation. Chapter 6. Matching Objects Using Union of Spheres Models 115 A l l algorithms described above have been implemented in the C programming lan-guage using the G L graphics library and compiled for SGI machines. The code to compute the power diagrams was obtained from the public domain package by Edelsbrunner et al. [G.8-Edelsbrunner et al. 92]. We used the F O R M S library [G.13-Overmars 92] for the interface. 6.4.1 Object registration The matching process gives us a series of pairwise transformations between circles (spheres) in the UoC (UoS) model of objects. To register objects using these, we have to consider a specific class of transformations. (a) Deriving affine transformations from matched circles A straightforward class to use, and one that can be justified in many applications is the class of affine transformations. In this case (in 2D) we are looking for 6 unknowns (the 2x2 matrix and a translation), and each match gives us 2 equations for the position, 1 for the scaling (the circles are matched with uniform scaling), and 1 for the feature (an angle of rotation). So, in general, this gives an over-determined system of linear equations that can be solved by standard least-squares techniques. Of course, assuming only rigid transformations reduces the number of unknowns to 3, and other restrictions are possible. In 3D, there are 12 unknowns (the 3x3 matrix and a translation) instead of 6; the rest of the procedure is exactly the same. It would be interesting to consider larger classes of transformations, in particular, subsets of differential homeomorphism. Each match between circles gives a Jacobian at the center of the circles. The task would be to find a vector field and a transformation that best fit the transformations given at the centers of circles. Chapter 6. Matching Objects Using Union of Spheres Models 116 (b) Deriving parametric polynomial transformations from matched circles Consider a rectangular region of image; P(X, Y), parametrized in terms of u and v such that X(u, v) = u, Y(u, v) = v, and therefore t u = (Jff^) = (1,0), and t„ = = (0,1). Let the transformed image be denoted by P(X'(u,v),Y'(u,v)) such that: X'(u,v) ud ud~l v 1 Y'(u,v) ud ud~l [M][G y ] [M T ] v ,4-i where [M] is p x p basis matrix corresponding to the transformation model chosen and [G x] and [G y] are p x p matrices of unknowns, p = d 4- 1, d being the degree of the polynomial. So, there are p2 unknowns each for X and Y. We can treat each coordinate separately. Let us see how we can get p2 equations for the X component. Equations for other components can be derived similarly. Let us say that circle a; in model A is matched to circle b; in model B. Let (Xi, YA = (u^ vA denote the coordinates of the center of a*. Since we know the corresponding (X!-,Y-), namely, the center of the transformed circle bj, we get one equation from the above relationship for X'(u, v). We can use center information from p2 pairs of matched circles (if possible) to get p2 equations. However, recall that each match provides us with the information about rotation and scaling at the center of the circles. One way to use the scaling and rotation information is to consider ( ^ - ) u . v . , which is the X-component of t'u, and equate it to the X-component of tu = (1,0) transformed by the scaling and rotation. Similarly, Chapter 6. Matching Objects Using Union of Spheres Models 117 we get another equation by considering the X-component of t'v = ( ^ r ) u . v. and equating it to the X-component of t„ = (0,1) transformed by the scaling and rotation. For example, let us consider the case of biquadratic Bezier function as the transfor-mation model. In this case, if a circle centered at (X;, Yj) = (ui,Vi) is transformed to a circle centered at (X-, Y/ ) , then: Poo Poi P02 u\ u{ 1 [B] PlO Pl l Pl2 [BT ] Vi P20 P21 P22 i 'dXr du 'dXr . dv UitVi Ui,Vi Poo Poi P02 = 2ui 1 0 [B] Pw Pl l Pl2 [B T ] Vi P20 P21 P22 1 Poo ' Poi P02 2vi = U2 Ui 1 [B] PlO Pl l Pl2 [B T ] 1 P20 P21 P22 0 where B is the 3 x 3 (constant) matrix corresponding to Bezier basis functions. There are 9 unknowns, p0o • • - P22- Each relationship above gives one linear equation in the 9 unknowns—so there are 3 equations per circle per coordinate. Since there are 9 unknowns we need at least 3 matched pairs of circles. Of course, if we use more than 3 matched pairs, the system of equations is overdetermined as we have more equations than unknowns; in this case we can use the method of least squares to get the best (in the least squared sense) solution. Examples For registration, the user loads in two images and the corresponding UoC models of the segmented regions to be registered. The program then computes the distances Chapter 6. Matching Objects Using Union of Spheres Models 118 (e) (f) Figure 6.46: Registration of two brain slices using UoC. (a) A 256 x 256 MRI slice through human head, (b) Same image, rotated, translated, scaled, and subsampled at 64 x 64 resolution, (c) UoC model of brain in (a), (d) UoC model of brain in (b). (e) The two UoC models together showing the extent of mismatch, (f) The two UoC models after registration. Chapter 6. Matching Objects Using Union of Spheres Models 119 between primitives, finds the minimum-weight match, and computes the transformation using a user-specified number of good matches for the least-squares fit (we implemented only the case of affine transformations). This transformation can then be applied to the UoC model or directly to the images4. Figure 6.46 shows an example of registration with the UoC method. The first image (a) is an MRI slice of a human brain, at 256x256 resolution. The second image (b) is from the first, subsampled (64x64), rotated by 20 degrees and scaled by a factor of 1.2 (it has also been translated to center the image in the frame, but we will not worry about that). The outline was computed from a range Of threshold (100.3-255). From the original number of circles (1703 and 432) each UoC model was simplified to about 200 circles each (228 and 218) and these were matched following the method described above. Images (c) and (d) show the UoC model outlines, image (e) the superposed outlines without registration and (f) shows the superposed outlines with registration. The scaling and rotation retrieved by computing moments were in this case 1.146 and 19.43 degrees. From the UoC method the scale was 1.166 and the rotation 19.27 to 20.05 degrees (since we did not constrain the matrix to be uniform scaling and rotation, the least-squares solution includes some amount of shearing and non-uniform scaling). Figures 6.47 shows an example of registration in 3D. These results are on P E T brain data. The original data has a resolution of (64 x 64 x 64) and the brain was segmented 4 This is an important point, because in many cases the UoC model is not intended as a substitute for the object itself or some other model, but as a tool to compute registration, matching and/or shape distances. Chapter 6. Matching Objects Using Union of Spheres Models 120 ( e ) ( D Figure 6.47: Example of registration in 3D. (a) Brain at 64 x 64 x 64 resolution (27917 spheres), (b) Transformed brain at 16 x 16 x 16 resolution (2774 spheres), (c) Brain in (a) after simplification (746 spheres), (d) Brain in (b) after simplification (734 spheres), (e) Data sets (a) and (b) before registration, (f) Data sets (a) and (b) after registration. Chapter 6. Matching Objects Using Union of Spheres Models 121 using a threshold of 140.3. This data was transformed by the following matrix: 0.9659 0.2578 -0.0226 0.0000 -0.2588 0.9623 -0.0842 0.0000 0.0000 0.0872 0.9962 0.0000 -5.0000 5.6000 0.2000 1.0000 It was resampled at a resolution of (16 x 16 x 16). There were 27,917 spheres in the UoS model for the original brain data and 2774 spheres in the transformed brain. Both models were simplified to about 750 spheres each. Since there are no prominent features in the two brains, Wf was set to zero while computing distances. The affine matrix obtained from the least squares method using the matches between the spheres in the two models was: 0.9608 0.2332 -0.0843 0.0000 -0.2731 1.0691 -0.0291 0.0000 0.0800 0.0833 1.0062 0.0000 -4.7759 4.8346 0.1792 1.0000 Figure 6.47e shows the two models together before registration and Figure 6,47f shows the two models after transformation of the first UoS model by the above matrix. 6.4.2 Measuring distances between objects The residual distance d(A, B) after registration using the optimal match is defined as the distance between the two objects. It is a measure of how well the objects match or, equivalently, how different the two objects are. When the two models are exactly the same, d(A, B) will be zero for a perfect match (assuming no numerical errors). The residual distance has some advantages over global measures such as moments. For exam-ple, two objects with different shapes (according to the human notion of shape) can have similar moments [H.3-Duda and Hart 73]. So, if we get similar moments, there is an Chapter 6. Matching Objects Using Union of Spheres Models 122 ambiguity about the object shape (see, for example, [H.3-Duda and Hart 73], page 212). Our method does not face such ambiguity because our measure for shape difference uses local comparisons whereas moments are results of global computation. Sometimes, the distribution of individual distances d(a, b) may be helpful in characterizing the extent of difference in shape. To make sure that the distances that we get between different objects are consistent, we normalize the distances. Normalizing the distance measure Let us say that we want to compare objects 0\, 0 2 , . . . , Om using their UoS models. To normalize our distances we do the following: • We scale down all the models by their respective diameters. The diameter of a model is defined as the diameter of the smallest enclosing sphere (circle). This ensures that the scaled model fits in a sphere of diameter 1, and therefore the largest sphere (circle) in a representation has a diameter of at most 1. • We position the models such that their centers of masses are at the same location. • Let (0, r p ) , (0, r s ) , and (0,r/) denote the ranges of the three distance components (namely, position, size, and feature). Since in general these are different for different pairs of objects, for the distance measure to be consistent between objects, we should use the same ranges across all objects. Ideally, rp, rs, and rj should be set to the theoretical maximum of the corresponding distances. Also, if we want equal weight of wp, ws, and Wf to have equal effect on the overall distance, we should also scale the distance components such that rp, rs, and 7 7 are the same (we choose this arbitrarily as 1). In our distance measure we let wp, ws, and uif vary over the Chapter 6. Matching Objects Using Union of Spheres Models 123 range (0,1) as well. What all this means is that the maximum distance between any matched pair of spheres (circles) can at most be 3. It still remains to compute the theoretical maxima of the three distances. Position distance after normalization has a theoretical maximum of 4 (when the center of masses of the models are close to the edge of the bounding circles and the two circles are a diameter away in opposite directions). Size distance after normalization of models by their diameters has a theoretical maximum of 0.25 (one circle has a radius of 0;5, the other 0.0). Feature distance has a theoretical maximum of at least 16 (one set of features has 4 vectors of magnitude 2 each, and the other set has 4 vectors of magnitude 0 each). • Sometimes, average distance between the matched pairs makes more sense than the overall distance between the models. To achieve this, we divide the total distance by the number of matches. If we have performed the normalization as described above, average distance between the matched pairs compared to 3 is a reasonable indicator of difference in shape. Note that there can be several other ways of normalizing the distance measure but we chose this one since it seemed to be simple and without any artifact. In particular, we could have normalized the models by scaling them such that they have the same volume (area in 2D). A careful investigation of the approach shows that it would be useful but one drawback is that unlike our method which yields a theoretical maximum for all three components of the distance measure, the theoretical maximum for the position distance for this approach is infinity. Although our measure does not explicitly take into account the difference in topology of the two objects (e.g., number of holes, cavities, connected components etc.), we believe that these are taken care of implicitly in the UoS model. This is because two objects Chapter 6. Matching Objects Using Union of Spheres Models 124 Table 6.7: Normalized distances between various 2D objects. Fish (1) Fish (2) Calf Cow Giraffe Fish (1) 0.0 0.0542 0.1382 0.2438 0.2495 Fish (2) 0.0542 0.0 0.1593 0.2268 0.2517 Calf 0.1382 0.1593 0.0 0.1019 0.1321 Cow 0.2438 0.2268 0.1019 0.0 0.0823 Giraffe 0.2495 0.2517 0.1321 0.0823 0.0 with "similar" shape—for example, one with a hole, the other without—have significantly different UoS models. Examples Table 6.7 gives the distances for the objects used in our 2D examples (the calf is in-i troduced later in Figure 6.57). These distances were computed as the average of the distances between pairs of circles after matching (with equal weights wp = ws = Wf = 1). These are somewhat dependent on the initial relative position of the two objects. A n obvious procedure is to use the match to register the two objects, and repeat the process until convergence. While we do not have at present any proof that it will converge, it is reasonable to assume that it will for good initial positions. The distances, as they are, indeed capture what one would expect to be the differences between shapes. The closest to fish (1) is fish (2) and the closest to the giraffe is the cow. The closest to the calf is the cow, but the distance is relatively large. This is possible because the two front legs of the calf overlap whereas for the cow they do not. Again, a further exploration of the role of the initial position may be useful to explain this situation. Distances calculated in this fashion can be used in shape classification where the goal is to classify an unknown object as one of m known objects. We will illustrate this with Chapter 6. Matching Objects Using Union of Spheres Models 125 Figure 6.48: Hand (3) and hand (4): models to be classified. the simple example of gesture recognition using the hand models that we have. Two hands (hand (3) and hand (4)) are shown in Figure 6.48. We need to find if they are Table 6.8: Normalized distances between hands. Hand (1) Hand (2) Hand (3) 0.0797 0.1628 Hand (4) 0.3805 0.1270 closer in gesture to hand (1) or hand (2) (Figure 6.39). Table 6.8 shows the distances between hand (3), hand (4) and hand (1), hand (2). These are the normalized distances as described before. For matching, the models were simplified to about 100 spheres each. According to this table, hand (3) is closer to hand (1) and hand (4) is closer to hand (2) which matches with the classification that humans do. Chapter 6. Matching Objects Using Union of Spheres Models 126 6.4.3 Showing the stability of the UoS representation The matching process can be used to demonstrate the stability of the UoS representation. It is important while discussing stability to distinguish between various causes of varia-tions in the data such as noise, affine transformations, or other distortions. Here we will show experimental evidence that the UoS representation is stable at least to boundary noise, shape-preserving transformations, and some other low frequency distortions. This will justify its use in comparing objects under these limited cases of distortions. Recall that a representation is stable if the following relationship holds for some bounded positive constant K: . d(R1,R2) < K,d{Px,P2) where Pi and P2 are two sets of boundary points for the same object, d(Pi,P2) is the distance between the point sets (which can be computed by optimal bipartite matching), and d(Ri, R2) is the distance between the two models. In other words, a representation is stable if the changes in distance between the primitives of the models are commensurate with the changes in the input data. To show the stability of the UoS representation to boundary noise, transformations, and other distortions we conduct the following experiment for a few objects (in 2 D and 3D) : 1. From the given boundary point data P i for an object, compute its model Ri. 2. Apply a perturbation to the boundary points to get P2 and compute the new model Ra-The perturbation can be applied to the boundary points or to the underlying vol-umetric data (if available). For example, it can be a combination of random affine Chapter 6. Matching Objects Using Union of Spheres Models 127 transformation and noise. Further, different number of boundary points can be used by appropriate resampling. 3. Perform matching between Pi and P 2 , and Ri and R2 as described in Section 5.3. Calculate the residual distances d(Pi,P2), and d(Ri,R2) for the best matches. 4. Repeat steps 2 and 3 several times. . 5. Plot d(Ri,R2) against d(Pi,P2). If our method of generating the UoS model is stable, the residual d(Ri,R2) should be bounded. Examples We ran the experiment described above on several data sets. First, we show the results for the 2D giraffe data introduced before under the following four conditions: A . Effect of noise on unsimplified UoC model: We perturbed the boundary points for giraffe data (512 x 512 pixels) randomly by increasing amounts and noted the corresponding average distance between the circles in the UoC models. Figure 6.49 shows the plot of average circle distance versus the average point distance. The plot indicates that the changes in the UoC models are commensurate with the changes in the underlying point sets. B. Effect of noise on simplified UoC model: In this case, the UoC models for the giraffe data (512 x 512 pixels) above were simplified to about 100 circles each. The average point distances and average circle distances were computed as before. Figure 6.50 shows the plot of average circle distance versus the average point distance. Note that the average circle distance is almost constant after simplification. This indicates that the model is determined by the shape of the giraffe and is fairly independent of sampling artifacts. It also reflects the smoothing effect of simplification. Chapter 6. Matching Objects Using Union of Spheres Models 128 o.io T3 1> 0.08 0.06 0.04 0.02 0.00 0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 average point distance Figure 6.49: Stability of unsimplified UoC representation to boundary noise for 2D giraffe data. o.io 0.08 0.06 0.04 0.02 0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 average point distance Figure 6.50: Stability of simplified UoC representation to boundary noise for 2D giraffe data. Chapter 6. Matching Objects Using Union of Spheres Models 129 C. Effect of shape-preserving transformations on simplified UoC model: For a shape-preserving transformation Tsp, UoC(Tsp(P)) ~ Tsp(UoC(P)). This is true if we transform the boundary points by a shape-preserving transformation and recompute the UoC model (this is because Delaunay triangulation is stable under shape-preserving transformations). However, in our case we transformed the underlying image for the giraffe data (256 x 256 pixels) by some random combination of rotation, translation, and uniform scaling. Then we subsampled it at different resolutions. The UoC models were computed for these transformed images and were simplified to about 100 circles each. The point distances and circle distances were computed by optimal bipartite matching. Figure 6.51 shows the plot of average circle distance versus the average point distance. 0.10 0.08 o 0 1 • | 0.06 "u I— & 0.04 ra u i> > 0.02 0.00 0.50 0.55 0.60 0.65 0.70 0.75 0:80 0.85 average point distance (x 100) Figure 6.51: Stability of simplified UoC representation to shape-preserving transforma-tions for 2D giraffe data. Again, note that the average circle distance is relatively constant for the simplified UoC model. 1 1 : 1 1 1 r J I I I I L Chapter 6. Matching Objects Using Union of Spheres Models 130 D . Effect of other distortions on simplified U o C model : In this case, we distorted the underlying giraffe image (256 x 256-pixels) by various distortion tools in PhotoShop™ such as shear, ripple, wave, etc.. Some of these distortions preserve lines while most of them do not and are non-linear. Again, the UoC models were computed for these transformed images and were simplified to about 100 circles each. The point distances and circle distances were computed by optimal bipartite matching. Figure o.io 0.08 •o •p 0.06 0.04 0.02 0.00 0.0000 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 average point distance Figure 6.52: Stability of simplified UoC representation to some low frequency distortions for 2D giraffe data. 6.52 shows the plot of average circle distance versus the average point distance. Again, note that the average circle distance is relatively constant for the simplified UoC model. Although stability depends on the level of simplification, the fact that we can achieve this level of stability justifies the use of this representation for shape comparison applications. To show the results of stability in 3D, we use the hand(l) data introduced before. This data set was obtained by ray-casting a polygonal hand model from 6 different views and using the algorithm described in Chapter 3. There are 15227 boundary points in Chapter 6. Matching Objects Using Union of Spheres Models 131 the original data set. We perturbed the boundary points by about 0.1% of the smallest bounding box side to about 1.0% of the smallest bounding box side. The dense UoS models for these data sets contained about 40000 spheres. These were simplified to have about 100 spheres each, and were then matched to the simplified version of the original model as described in the previous section. Calculating distances for data sets of this size takes about 9 minutes on an Iris Indigo 2. Figure 6.53 shows the result of this a as a. 0.20 0.15 0.10 0.05 -0.00 0.00 5.00 10.00 15.00 20.00 average point distance (x 1000) 25.00 Figure 6.53: Stability of UoS representation to boundary noise for 3D hand data. experiment as a plot of the average sphere distance against the average boundary point distance. Recall that stability demands that the distance in the underlying model be at most proportional to the distance in the boundary samples. The average sphere distance, being nearly constant, indicates that the UoS model can be stabilized to a level where it is almost independent of the perturbation of the boundary points arid hence can be used for shape comparison. We could not conduct experiments in 3D for other types of distortions because com-puting the distance between points involves a bipartite matching on 10,000 points or so. Chapter 6. Matching Objects Using Union of Spheres Models 132 This is not feasible for our implementation given that the process is 0(n2) in space and 0(n 3 ) in time. However, based on the above experiments in 2D and 3D we believe that the UoS representation is fairly stable. This is because the algorithms for generation of the UoS model and its simplification in 3D are based on exactly the same principles as in 2D. 6.4.4 Shape interpolation For interpolation, the goal is to determine a path from every circle in model A to its matched counterpart in model B. Again, in general, the match implies a translation, a scaling, and a rotation. The easiest solution is to interpolate the circles linearly, i.e., the centers of the circles move in a straight line, and the size varies linearly as a function of time. It has been observed from the beginning of in-betweening that there are many cases where such a linear approximation is not appropriate. Two strategies are possible given our approach. One is to apply the effect of rotation through the path between center of circles. Since there is no particular reason to vary the rate of rotation, a curve of constant curvature, namely an arc of circle, is sufficient, and is fully defined given the two end points (the centers of the circles, matched) and the rotation angle (the angle of rotation necessary to match the features of the two cir-cles). The other strategy is to use the skeletons of section 5.3.2. Several papers have considered a similar problem, and the reader is referred to the work 'of Sederberg et al. and Shapira et al. [F.13-Sederberg and Greenwood 92] [F.12-Sederberg et al. 93] [F.14-Shapira and Rappoport 95]. , Examples To achieve shape interpolation with our system, the user loads in the UoC models of two objects. If segmentation is not wanted, the program proceeds to match the UoC models Chapter 6. Matching Objects Using Union of Spheres Models 133 based on the weights wp, ws, and wj given by the user. If segmentation is requested, the program automatically decomposes the objects into segments using the heuristics described earlier. The user then specifies which segments of one object should match to which segments of the other (given a measure for distance between the segments this can be done automatically by using the matching algorithm described above). The user is also allowed to break/merge parts using the skeletons of the objects. The program then matches the circles inside the pairs of matched segments. If not satisfied with the automatic matching, the user is allowed to interfere and govern the matching process by manually aligning the two parts. Finally all the matches are concatenated. The whole process of shape interpolation, beginning from the scanned images of two objects for the 2D examples presented next, took less than 30 minutes each on an IRIS Indigo 2. The longest step is to compute the distances in the bipartite graph, since there are 0(n2) of them. Even though the optimal graph matching is 0(n3), in practice the running time is only a few seconds, as for all the other steps. This includes the time for regular trian-gulation of circles, which is done by a separate program [G.8-Edelsbrunner et al. 92]. Figure 6.54 shows the matching and the subsequent interpolation between a rectangle and a MRI brain slice. This is to show that even polygons can be treated well with UoC model (even though they are not the objects UoC model is best suited for) and that the interpolation process is insensitive to changes in topology. Note also that the common symmetry between the objects is respected. In the absence of any semantics, this is a good solution to this particular problem. Figure 6.55 shows the interpolation between two hand-drawn "fish" objects (their models simplified to 45 and 43 circles, respectively). These objects were selected to show how the matching process behaves on generic cases of branches, taperings and "curling". It should be noted that humans might hesitate to recognize a fish in fish (2) and instead see legs and a curly tail. Nevertheless it is not difficult to correctly match the parts of the Chapter 6. Matching Objects Using Union of Spheres Models (g) ( h ) Figure 6.54: Interpolation from a rectangle to an MRI slice of brain for wa = 0.2, and Wf - 0.1. Chapter 6. Matching Objects Using Union of Spheres Models 135 tails as well as the rest of the body.. The animation shows clearly the potential problems when interpolating in a straight line between positions when the tail curls. Again we should stress that since the rotation is known when the features are matched, we have enough information to follow an arc of circle instead of a straight line. Figure 6.56 shows for the same two objects examples of the circles matched by features alone; this measure is a good one for its purpose. The two main branching features (tail and leg) are correctly matched, and the tapering, "terminating" parts are matched also. Pairs (g) and (h) shows that without position distance, end-parts can be matched "incorrectly". The features indeed match, and in fact so do the sizes, but the positions do not. In this case a non-zero weight for positions will help. Figure 6.57 shows a more complex example, a match and interpolation between a calf and a cow. The images were scanned at a resolution of 512 x 512 pixels. Again this is done from simplified' UoC models, with 154 and 239 circles, respectively. The weights chosen are such that the legs match well even without segmentation. Note that since these are two-dimensional images, the front legs of the calf are overlapping, and therefore there is a many-to-one mapping of "segments" here, and humans will use semantics and knowledge of the three-dimensional shape to decide what is correct. In the interpolation computed, part of the legs of the cows separate during the process. Also, a circle migrates from the front rear leg of the calf to the rear rear leg of the cow. This is because the features and the sizes are similar enough to cancel the effect of the distance in position. Segmentation can help solve these problems as Figure 6.58 illustrates. There are 10 segments in the cow and 9 in the calf. The segments were automatically extracted, but matched manually. One can see that there are no "wandering" circles in this case. Figure 6.59 shows steps in an interpolation between images of a cow and a giraffe. These images were again scanned from a book. The cow has 239. circles in its simplified form, and the giraffe has 370. In this example segmentation was used and the segments Chapter 6. Matching Objects Using Union of Spheres Models 136 (a) (b) Figure 6.56: Matching circles in the two fishes using features alone, (a-b), (c-d), (e-f), and (g-h) show four pairs of matched circles. Notice a "poor" match in (g-h) pair. Chapter 6. Matching Objects Using Union of Spheres Models 138 Chapter 6. Matching Objects Using Union of Spheres Models 139 (a) (b) Figure 6.58: Interpolation from a calf to a cow (using automatic segmentation). Different weights were used to match different components. Chapter 6. Matching Objects Using Union of Spheres Models 140 Figure 6.59: Interpolation from a cow to a giraffe (using automatic segmentation). Dif-ferent weights were used to match different components. Chapter 6. Matching Objects Using Union of Spheres Models 141 were matched manually. There are 10 segments each in the two models. The interpolation is as one would hope, except that the right ear of the giraffe comes from the neck of the cow. This is reasonable given that there is no feature in the cow corresponding to it. Figure 6.60 shows an example of shape interpolation in 3D—the hand (1) model shown before interpolates to hand (3). Notice that the spheres in the fingers and thumb of hand (1) get matched well to the the spheres in the fingers and thumb of hand (3) respectively. The matching and interpolation were totally automatic in this case. The weights that were used in the distance measure were: ws = 0.8, wp = 1.0, and Wf = 0.7. Since the matching and interpolation in our examples is done on the simplified UoS model, the rendered version shown here lacks details. The display of the UoS model can be smoothed by using blending functions (i.e., using blobby display [I.6-Ranjan and Fournier 94]). Figure 6.61 shows an example of shape interpolation between hand (1) to hand (2). In this case, because the fingers and thumbs are bent relatively more, our matching process fails to come up with a 100% "good" matches for the spheres as judged by human standards. There are two spheres which migrate from one finger to another. The matching and interpolation were totally automatic in this case. The. weights that were used in the distance measure were: ws = 0.5, wp = 1.0, and Wf = 0.5. One solution to this problem, as in 2D, is segmentation of objects into parts, matching parts between the two objects, and then matching spheres in the matched parts. For example, in this case, the hands can be segmented into palm, four fingers, and a thumb. Another solution can be to use another hand model in-between hand (1) and hand (2). So for example, in our case, we can first interpolate between hand (1) and hand (3) and then from hand (3) to hand (2). Of course, this is not always possible for other cases. Figure 6.62 shows four pairs of matches between the two hands. Note that the spheres in the body of the finger and at the tip are matched correctly (the features in the tapering portions of object are distinctly defined). It is interesting to see that even in the flat Chapter 6. Matching Objects Using Union of Spheres Models 142 Figure 6.60: Interpolation from hand (1) to hand (3). Figure 6.61: Interpolation from hand (1) to hand (2). Chapter 6. Matching Objects Using Union of Spheres Models 144 Figure 6.62: Example of matching in 3D. Four pairs (a-b, c-d, e-f, g-h) of matched spheres along with their features for the hand data. Chapter 6. Matching Objects Using Union of Spheres Models 145 portions of objects where the features may not be distinctive, the spheres are well matched because a combination of position, size, and feature is used. Finally, Figure 6.63 shows a funny example of shape interpolation between a rabbit (obtained from multiple view range data, 1046 spheres) and a human head (obtained from multiple-view ray-traced data, 262 spheres). The alignment procedure was semi-automatic for this case. First we aligned the models by moments and then rotated one of the models to proper (that is, semantically appropriate) orientation. The estimation of moments for both models took a total of about 1.5 seconds. The distances calculated in this example took into account position and size only (wp = 1.0, ws = 0.25, wj = 0.0). It took about 4 seconds to calculate the distances on an IRIS Indigo 2 machine. Whereas this measure seems to work quite well, other context or feature-based measures will do better for specific purposes, for instance to match the rabbit's ears to the human ears. To match a total of 2000 spheres it took about 1.5 minutes on an IRIS Indigo 2. Figures 6.63b-e show the in-betweens for the rabbit interpolating to a human head. We chose the linear interpolation in radius so that the linear dimensions of the model increase by about equal amounts for each time step. We found this to be more pleasing than an interpolation where the volume of the primitives is increased linearly. The display rate for 1000 spheres (100 polygons per sphere, Gouraud-shaded) is about 2 frames per second on an Indigo 2. 6.5 S u m m a r y In this chapter, we showed how to match spheres between the UoS models of two objects. These matches were then used to define distances between objects. We used the distances to show that the UoS representation is stable to boundary noise, some transformations, and to small distortions in the data. We also showed how we can derive transformations Figure 6.63: Interpolation from a rabbit to a human head. Chapter 6. Matching Objects Using Union of Spheres Models 147 between data sets using the matches and used it to register objects and to do shape interpolation between objects. The shape distance measure deduced from the matches has merit and captures quantitatively the concept of difference in shape. The uses of this are to be investigated further, especially in the area of object recognition. Chapter 7 Conclusion and Future Work This thesis deals with the issue of computer representation of 3D objects. Different representations are good for different applications; no representation will be superior to others in all respects. The choice of a computer representation of an object depends on the application at hand. The major goal of this thesis was to study the union of spheres (UoS) representation of objects, where a union of overlapping sphere primitives is used to model an object. A careful investigation of current existing representations shows that they lack stabil-ity. A stable representation of an object is independent of the sampling artifacts such as the method of generation of boundary points, the resolution and sampling geometry, the imaging methodology, or some noise in the data. Stable representations help character-ize shapes for comparison or recognition; skeletal (or medial axis) and volume primitive representations have been popular in vision for the same reason. Piecewise polyhedral representations, e.g., tetrahedra, and voxel representations, e.g., octrees, generally tend to be unstable. This thesis shows that it is possible to obtain stable UoS models of objects for applications requiring shape comparison, and therefore the UoS representation makes a viable alternative for such applications. Although some operations become difficult to perform because of the overlapping primitives, efficient algorithms exist to perform a wide variety of operations on and using the UoS model. 148 Chapter 7. Conclusion and Future Work 149 7.1 Summary of contributions The following are the contributions of the research work presented in this thesis: • Generation of UoS model for objects: Algorithms to generate UoS models for objects from volumetric data such as C T and MRI scans, and from boundary point data available from range-finders and ray-tracers are developed. What this means is that a wide variety of object models in computer graphics can be converted to the UoS model if deemed useful. The UoS models are derived from a Delaunay tetrahedralization of boundary points followed by an inside-outside classification of the spheres circumscribing the tetrahedra. • Simplification of UoS model: A simplification algorithm based on the concept of sphericity is defined. Sphericity is a measure of how well a set of spheres can be approximated by a single sphere. The simplification method stabilizes the model and reduces the number of primitives while preserving the object shape. It also gives control over the error and the number of primitives, and is efficient and easy to implement. By varying a parameter (the sphericity threshold), different levels of simplification can be achieved. The concept of sphericity is scale-independent and hence objects derived from data sets at different resolutions can be compared using their simplified UoS models. • Error bounds for UoS representation: For an r-regular object, it is shown that if its boundary is sampled properly (i.e., the sampling density 8 < r), the error in the model reconstructed frorh boundary points (using union of arcs representation for curves, and UoC representation for areas) is < 8. This is the same order of error as their corresponding linear counterparts, namely, the polyline model of curves, and the UoT model of 2D objects. The error in the UoS model is 0(8) if the error Chapter 7. Conclusion and Future Work 150 in the UoTet model is 0(6). For volumetric data, the error in the UoS model is <-y/Zv, where v is the linear dimension of the grid used in volumetric sampling. For curves, it is obvious that the union of arcs model cannot do worse than the polyline model (since polylines are a special case of the union of arcs model). For closed objects in 2D and 3D, since there is a limit on the curvature of primitives used (due to finite size of objects, and our demand that the whole primitive be contained in the object), the results summarized in the previous paragraph are not obvious. These results are important because whereas the UoTet representation can be unstable due to boundary noise, the UoS representation can be made relatively immune.to it and hence made useful for comparing objects. In other words, without compromising on the order of error, we can obtain stable UoS models of objects. • M a t c h i n g w i t h U o S models: A n algorithm to match spheres in the UoS models of two objects is presented. Matching is based on a distance measure which captures the difference between the primitives of the two models based on their size, position, and their immediate context in the model; this provides control during the matching process. The matches are obtained by an optimal bipartite graph matching. • Distance between objects and s tabi l i ty of U o S representation: The resid-ual distance from the optimal matching of the UoS models of two objects is used to define distance between them. Distance between objects.can be used, for example, in object classification, determination of growth, and measurement of change in shape. The notion of "shape" distance is then used to show the stability property of the UoS representation. It is shown that the UoS representation can be made relatively insensitive to boundary noise, to shape-preserving transformations, and to other small distortions. Chapter 7. Conclusion and Future Work 151 • Extracting transformations using UoS models: Algorithms to extract trans-formations between data sets from the matched spheres in their UoS models are described. \These are then used to register and interpolate objects. • Software: Programs to compute the proposed volume models (from volumetric data, and from range-finder and ray-tracer data), to simplify them, and to match and interpolate objects using them have been implemented in both two and three dimensions. The research has resulted in extensive tools to study and visualize objects modeled as UoC (in 2D) and UoS (in 3D). 7.2 Future work Ideas for further research fall into two categories: (A) work related directly to UoS model which will enhance its use in graphics, and (B) application of the UoS representation to other interesting problems. Some of them are: • A . l Fast smooth display: Fast and smooth display of objects modeled as UoS is essential for it to be useful in graphics. Ideally, rendering of models at interactive speeds is desirable. Currently, we use the "blobby" display of models for smooth display but it is slow. Polygonization alternatives described in Chapter 2 and in [I.2-Fournier and Ranjan 96] should be considered as several current machines support polygon rendering in hardware. • A.2 Texture mapping on UoS model: There are two issues here: (i) traditional texture mapping, i.e., mapping a rectangular image onto the surface of the UoS model, and (ii) preserving the texture or color during the simplification of the UoS model. Traditional texture mapping is extremely difficult considering that mapping a parametrically square texture onto a single sphere without singularity Chapter 7. Conclusion and Future Work 152 is not possible. Preservation of textures during simplification may be possible to some extent by mapping the texture of the union of all spheres enclosed by the cluster sphere on to the new cluster sphere. • B , l Segmentation of objects i n volumetr ic data: Our algorithms work under the assumption that the object has been adequately segmented with respect to the threshold. The use of the representation will be greatly enhanced if it is possible to meaningfully segment the data using the representation. The number of spheres in the UoS model of an object (after simplification) may provide a useful clue in this regard. • B . 2 Segmentation of a 3D object into components: We believe that using the features, it is possible to specify the critical regions in the objects (such as branching region). By following these branches, it should be possible to segment an object in a stable fashion into regions where it is cylindrical, planar, conical, and the like. This will be useful in decomposing a complex object into parts for better shape interpolation algorithms, providing "skeletal" handles for animation, or for object recognition. We have done some work in this regard and the preliminary results are promising. • B .3 Ver i fy ing "difference" i n shapes: One interesting aspect of the UoS repre-sentation is that it gives a measure of "difference" between shapes. It remains to be seen if it actually matches with the intuitive notion of shape. It will be interesting to conduct experiments with humans to measure the "difference" in the shape of different objects and see if the residual distance that we get from matching their UoS models has any correlation with the human's perception of shape. Chapter 7. Conclusion and Future Work 153 • B .4 Fractal dimension from UoS model: Because the UoS representation is stable (i.e., is related to the shape of the object), it should be possible to get a rough estimate of the fractal dimension of an object from the UoS representation by plotting the number of spheres (after simplification) in the model of an object at different scales. Bibl iography A . Surveys [A.l-Badler and Bajcsy 78] N.I. Badler and R. Bajcsy. "Three-dimensional representa-tions for computer graphics and computer vision," Computer Graphics (SIG-GRAPH '78 Proceedings), Vol. 12, No. 3, August 1978, pp. 153-160. [A.2-Bolle and Vemuri 91] R . M . Bolle and B.C. Vemuri. "On three-dimensional surface reconstruction methods," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 13, No. 1, January 1991, pp. 1-13. [A.3-Brown 92] L .G.Brown. "A survey of image registration techniques," ACM Com-puting Surveys, Vol. 24, No. 4, December 1992, pp. 325-376. [A.4-Elvins 92] T.T. Elvins. "A survey of algorithms for volume visualization," Com-puter Graphics, Vol. 26, No. 3, August 1992, pp. 194-201. [A.5-Hebert et al. 94] M . Hebert, J. Ponce, T .E . Boult, A . Gross, D. Forsyth. "Report of the N F S / A R P A workshop on 3D object representation for computer vision," C U N Y Graduate Center, New York, December 5-7, 1994. [A.6-McCormick et al. 87] B. McCormick, T.S. Defanti, and M.D. Brown. "Visualization in scientific computing," Computer Graphics, Vol. 21, No. 6, November 1987. [A.7-Requicha 80] A . A . G . Requicha. "Representations for rigid solids: theory, methods, and systems," ACM Computing Surveys, Vol. 12, No. 4, 1980, pp. 437-464. [A.8-Srihari 81] S.N. Srihari. "Representation of three-dimensional digital images," ACM Computing Surveys, Vol. 13, 1981, pp. 399-424. [A.9-Stytzi et al. 91] M.R. Stytzi, G. Frieder, and O. Frieder. "Three-dimensional med-ical imaging: Algorithms and computer systems," ACM Computing Surveys, Vol. 23, No. 4, December 1991, pp. 421-499. B . 3D data acquisi t ion i n medicine [B.l-Bates et al. 83] R.H.T. Bates, K . L . Garden, and T . M . Peters. "Overview of com-puterized tomography with emphasis on future developments," Proceedings of the IEEE, Vol. 71, No. 3, March 1983, pp. 356-372. 154 Bibliography 155 [B.2-Jones 90] T. Jones. "Positron emission tomography," Clinical Physics and Physio-logical Measurement, Vol. 11, No. A , 1990, pp. 27-36. [B.3-Knoll 83] G.F. Knoll . "Single photon emission computed tomography," Proceedings of the IEEE, Vol. 71, No. 3, March 1983, pp. 320-329. [B.4-Maguire et al. 91] G.Q. Maguire Jr., M . E . Noz, H. Rusinek, J. Jaeger, E .L . Kramer, J.J. Sanger, and G. Smith. "Graphics applied to medical image registration," IEEE Computer Graphics and Applications, Vol. 11, No. 2, March 1991, pp. 20-28. [B.5-Paddpck 91] S.W. Paddock. "The laser-scanning confocal microscope in biomedical research," Proceedings of the Society For Experimental Biology & Medicine., Vol. 198, No. 3, December 1991, pp. 772-780. [B.6-Riederer 89] S.J. Riederer. "Recent advances in magnetic resonance imaging," Pro-ceedings of the IEEE, Vol. 76, No. 9, September 1989, pp. 1095-1105. [B.7-Wells 88] P.N.T. Wells. "Ultrasound imaging," Journal of Biomedical Engineering, Vol. 10, No. 6, November 1988, pp. 548-554. C. 3D representations C . l Surface methods [C. l . l -Eck et al. 95] M . Eck, T. DeRose, T. Duchamp, H . Hoppe, M . Lounsbery, and W. Stuetzle. "Multiresolution analysis of arbitrary meshes," Computer Graphics (SIGGRAPH '95 Proceedings), August 1995, pp. 173-182. [C.1.2-Fuchs et al. 77] H. Fuchs, Z . M . Kedem, and S.P. Uselton. "Optimal surface recon-struction from planar contours," Communications of the ACM, Vol. 20, No. 10, October 1977, pp. 693-702. [C.1.3-Hoppe et al. 92] H. Hoppe, T. DeRose, T. Duchamp, J . MacDonald, and W. Stuetzle. "Surface reconstruction from unorganized points," Computer Graphics (SIGGRAPH '92 Proceedings), Vol. 26, No. 2, July 1992, pp. 71-78. [C.1.4-Hoppe et al. 93] H. Hoppe, T. DeRose, T. Duchamp, J . MacDonald, and W. Stuetzle. "Mesh optimization," Computer Graphics (SIGGRAPH '93. Pro-ceedings), Vol. 27, August 1993, pp. 19-26. Bibliography 156 [C.1.5-Lin et al..89] W.C. Lin, S.Y. Chen, and C.T. Chen. "A new surface interpolation technique for reconstructing 3D objects from serial cross-sections," Computer Vision, Graphics and Image Processing, Vol. 48, No. 1, October 1989, pp. 124-143. [C.1.6-Lorensen and Cline 87] W . E . Lorensen and H.E. Cline. "Marching cubes: A high resolution 3D surface construction algorithm," Computer Graphics (SIGGRAPH '87 Proceedings), Vol. 21, No. 4, July 1987, pp. 163-169. [C.1.7-Miller et al. 91] J .V. Miller, D.E. Breen, W.E . Lorensen, R . M . O'Bara, and M . J . Wozny. "Geometrically deformed models: A method for extracting closed geometric models from volume data," Computer Graphics (SIGGRAPH '91 Pro-ceedings), Vol. 25, No. 4, July 1991, pp. 217-226. [C.1.8-Reissell 93] L . M . Reissell. "Multiresolution geometric algorithms using wavelets I: Representations for parametric curves and surfaces," Technical Report 93-17, Department of Computer Science, U B C , May 1993. [C.1.9-Schroeder et al. 92] W.J . Schroeder, J A . Zarge, and W.E . Lorensen. "Decimation of triangle meshes," Computer Graphics (SIGGRAPH '92 Proceedings), Vol. 26, No. 2, July 1992, pp. 65-70 [C.l.lO-Meyers et al. 91] D. Meyers, S. Skinner, and K . Sloan. "Surfaces from contours: The correspondence and branching problems," Proceedings of Graphics Interface '91, June 1991, pp. 246-254. [C . l . l l -Tu rk 92] G. Turk. "Re-tiling polygonal surfaces," Computer Graphics (SIG-GRAPH '92 Proceedings), Vol. 26, No. 2, July 1992, pp. 55-64. C.2 Voxel methods [C.2.1-Meagher 82] D. Meagher. "Geometric modeling using octree encoding," Com-puter Graphics and Image Processing, Vol. 19, June 1982, pp. 129-147. [C.2.2-Montani and Scopigno 90] C. Montani and R. Scopigno. "Rendering volumetric data using the STICKS representation scheme," Computer Graphics (San Diego Workshop on Volume Visualization), Vol. 24, No. 5, November 1990, pp. 87-93. [C.2.3-Muraki 93] S. Muraki. "Volume data and wavelet transform," IEEE Computer Graphics and Applications, Vol. 13, No. 4, July 1993, pp. 50-56. [C.2.4-Udupa and Odhner 91] J .K. Udupa and D. Odhner. "Fast visualization, manip-ulation, and analysis of binary volumetric objects," IEEE Computer Graphics and Applications, Vol. 11, No. 6, November 1991, pp. 53-62. Bibliography 157 [C.2.5-Udupa and Odhner 93]. J .K. Udupa and D. Odhner, "Shell rendering," IEEE Computer Graphics and Applications, Vol. 13, No. 6, November 1993, pp. 58-67. [C.2.6-Wilhelms and Gelder 90] J. Wilhelms and A . Van Gelder. "Octrees for faster iso-surface generation extended abstract," Computer Graphics (San Diego Workshop on Volume Visualization), Vol. 24, No. 5, November 1990, pp. 57-62. C.3 Volume methods [C.3.1-Boissonnat 84] J. Boissonnat. "Geometric structures for three-dimensional shape representation," ACM Transactions on Graphics, Vol. 3, No. 4, October 1984, pp. 266-286. [C.3.2-Edelsbrunner et al. 94] H . Edelsbrunner and E.P. Miicke. "Three-dimensional al-pha shapes," ACM Transactions on Graphics, Vol. 13, No. 1, January 1994, pp. 43-72. [C.3.3-Edelsbrunner 92] H. Edelsbrunner. "Weighted alpha shapes," Technical Report UIUCDCS-R-92-1760, Department of Computer Science, UIUC, Urbana, Illinois, 1992. [C.3.4-Naylor 92a] B .F . Naylor. "Interactive solid geometry via partitioning trees," Pro-ceedings of Graphics Interface '92, May 1992, pp. 11-18. [C.3.5-Naylor 92b] B .F . Naylor. "Partitioning tree image representation and generation from 3D geometric models," Proceedings of Graphics Interface '92, May 1992, pp. 201-212. [C.3.6-Solina and Bajcsy 90] F. Solina and R. Bajcsy. "Recovery of parametric mod-els from range images: The case for superquadrics with global deformations," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 12, No. 2, February 1990, pp. 131-147. [C.3.7-Terzopoulos and Metaxas 91] D. Terzopoulos and D. Metaxas. "Dynamic 3D models with local and global deformations: Deformable superquadrics," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 13, No. 7, July 1991, pp. 703-714. Bibliography 158 D . U n i o n of Spheres models D . l Generation [D.l. l-Badler et al. 79] N.I. Badler, J. O'Rourke, and H. Toltzis. "A spherical represen-tation of a human body for visualizing movement," Proceedings of the IEEE, Vol. 67, No. 10, October 1979, pp. 1397-1403. [D.1.2-Hubbard 96] P .M. Hubbard. "Approximating polyhedra with spheres for time-critical collision detection," ACM Transactions on Graphics, Vol. 15, No. 3, July 1996, pp. 179-210. [D.1.3-Hubbard 95] P .M. Hubbard. "Collision detection for interactive graphics appli-cations," Ph.D. Thesis, Department of Computer Science, Brown University, Providence, U.S.A. [D.1.4-Mohr and Bajcsy 83] R. Mohr and R. Bajcsy. "Packing volumes by spheres," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 3, No. 1, January 1983, pp. 111-116. [D.1.5-0'Rourke and Badler 79] J. O'Rourke and N.I. Badler. "Decomposition of three-dimensional objects into spheres," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 1, No. 3, July 1979, pp. 295-305. D.2 Display [D.2.1-Bloomenthal 88] J . Bloomenthal. "Polygonisation of implicit surfaces," Computer Aided Geometric Design, Vol. 5, No. 4, 1988, pp. 341-355. [D.2.2-Connolly 83a] M . Connolly. "Solvent-accessible surfaces of proteins and nucleic acids," Science, August 1983, pp. 709-713. [D.2.3-Connolly 85] M.Connolly. "Molecular surface triangulation," Journal of Applied Crystallography, Vol. 18, 1985, pp. 499-505. [D.2.4-Franklin 81] W.R. Franklin. "An exact hidden sphere algorithm that operates in linear time," Computer Graphics and Image Processing, Vol. 15, Apri l 1981, pp. 364-379. [D.2.5-Fuchs et al. 85] H. Fuch's, J. Goldfeather, J.P. Hulrquist, S. Spach, J.D. Austin, F.P. Brooks, J .G. Eyles, and J . Poulton. "Fast spheres, shadows, textures, trans-parencies, and image enhancements in pixel-planes," Computer Graphics (SIG-GRAPH '85 Proceedings), July 1985, pp. 111-120. Bibliography 159 [D.2.6-Halperin and Overmars 94] D. Halperin, and M . H . Overmars. "Spheres, molecules, and hidden surface removal," Proceedings of the Tenth Annual Sym-posium on Computational Geometry, 1994, pp. 113-122. [D.2.7-Knowlton 81] K . Knowlton. "Computer-aided definition, manipulation and de-piction of objects composed of spheres," Computer Graphics, Vol. 15, No. 4, December 1981, pp. 352-375. [D.2.8-Knowlton and Cherry 77] K . Knowlton and L. Cherry. "ATOMS, a three-d opaque molecule system," Computers and Chemistry, Vol. 1, No. 3, 1977, pp. 161-166. [D.2.9-Max 79] N.L . Max. " A T O M L L L : A T O M S with shading and highlights," Com-puter Graphics (SIGGRAPH '79 Proceedings), Vol. 13, No. 2, August 1979, pp. 165-173. [D.2.10-Overweld and Wyvil l 93] C . W . A . M . van Overweld and B. Wyvil l . "Potentials, polygons and penguins: An efficient adaptive algorithm for triangulating an equi-potential surface," Proceedings of the Fifth Annual Western Computer Graphics Symposium, Vernon, Canada, March 1993. [D.2.11-Sanner et al. 95] M.F . Sanner, A . J . Olson, and J. Spehner. "Fast and robust computation of molecular surfaces," Short communication, Eleventh Annual Symposium on Computational Geometry, Vancouver, Canada, June 1995. D.3 Computing mass properties [D.3.1-Aurenhammer 87] F. Aurenhammer. "Power diagrams: Properties, algorithms and applications," SIAM Journal of Computing, Vol. 16, 1987, pp. 78-96. [D.3.2-Aurenhammer 88] F. Aurenhammer. "Improved algorithms for disks and balls using power diagrams," Journal of Algorithms, Vol. 9, 1988, pp. 151-161. [D.3.3-Avis et al. 88] D. Avis, B . K . Bhattacharya, and H. Imai. "Computing the volume of the union of spheres," The Visual Computer, Vol. 3, No. 6, May 1988, pp. 323-328. [D.3.4-Connolly 83b] M . Connolly. "Analytical molecular surface calculation," Journal of Applied Crystallography, Vol. 16, 1983, pp. 548-558. [D.3.5 •-Edelsbrunner 93] H. Edelsbrunner. "The union of balls and its dual shape," Proceedings of the 9th Annual Symposium on Computational Geometry, 1993, pp. 218-229. Bibliography 160 [D.3.6-Lee and Requicha 82] Y . T . Lee and A . A . G . Requicha. "Algorithms for computing the volume and other integral properties of solids. II. A family of algorithms based on representation conversion and cellular approximation," Communications of the ACM, Vol. 25, September 1982, pp. 642-650. D.4 Blobby models of objects [D.4.1-Blinn 81] J.F. Blinn. "A generalization of algebraic surface drawing," ACM Transactions on Graphics, Vol. 1, No. 3, August 1981, pp. 227-234. ) [D.4.2-Kacic-Alesic 91] Z. Kacic-Alesic and B. Wyvi l l . "Controlled blending of pro-cedural implicit surfaces," Proceedings of Graphics Interface '91, June 1991, pp. 236-245. [D.4.3-Wyvill et al. 86] G. Wyvil l , C. McPheeters, and B. Wyvi l l . "Data structure for soft objects," The Visual Computer, Vol. 2, No. 4, August 1986, pp. 227-234. D.5 Skeletal models [D.5.1-Attali and Montanvert 94] D. Attali and A . Montanvert. "Semicontinuous skele-tons of 2D and 3D shapes," Proceedings of Visual Form, Analysis, and Recog-nition, Capri, 1994. [D.5.2-Brandt 92] J.W. Brandt. "Describing a solid with a three-dimensional skeleton," Proceedings of the SPIE, Vol. 1830, Curves and Surfaces in Computer Vision and Graphics III, 1992, pp. 258-268. [D.5.3-Brandt and Algazi 92] J.W. Brandt and V.R . Algazi. "Continuous skeleton com-putation by Voronoi diagram," CVGIP: Image Understanding, Vol. 55, No. 3, May 1992, pp. 329-338. [D.5.4-Goldak et al. 91] J .A. Goldak, X . Yu, A . Knight, and L. Dong. "Constructing discrete medial axis of 3-D objects," Internat. J. Comput. Geom. Appl.,\o\. 1, No. 3, 1991, pp. 327-339. [D.5.5-Kirkpatrick 79] D.G. Kirkpatrick. "Efficient computation of continuous skele-tons," Proceedings of the 20th Annual IEEE Symposium on Foundations of Computer Science, 1979, pp. 18-27. [D.5.6-Morse et al. 94] B.S. Morse, S.M. Pizer, and D.S. Fritsch. "Robust object repre-sentation through object-relevant use of scale," Proceedings of SPIE Conference on Medical Imaging, SPIE Vol. 2167, 1994, pp. 104-114. Bibliography 161 [D.5.7-Tsao and Fu 84] Y . F . . Tsao and K.S. Fu. "Stochastic skeleton modeling of ob-jects," Computer Vision, Graphics, and Image Processing, Vol. 25,1984, pp. 348-370. E . Reconstruction from range data [E.l-Besl and McKay 92] P.J. Besl and N.D. McKay. "A method of registration of 3-D shapes," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 14, No. 2, February 1992, pp. 239-256. [E.2-Chen and Medioni 92] Y . Chen and G. Medioni. "Object modelling by registration of multiple range images," Image and Vision Computing, Vol. 10, No. 3, Apri l 1992, pp. 145-155. [E.3-Curless and Levoy 96] B . Curless and M . Levoy. "A volumetric method for building complex models from range images," To appear in, Computer Graphics (SIG-GRAPH '96 Proceedings). [E.4-Muraki 91] S. Muraki. "Volumetric shape description of range data using "Blobby Model"," Computer Graphics (SIGGRAPH '91 Proceedings), Vol. 25, No. 4, July 1991, pp. 227-235. [E.5-Schmitt et al. 91] F. Schmitt, X . Chen, and W.-H. Du. "Geometric modelling from range image data," In W. Purgathofer, editor, Eurographics '91, North-Holland, ; September 1991, pp. 317-328. [E.6-Soucy and Laurendeau 92] M . Soucy and D. Laurendeau. "Multi-resolution surface modeling from multiple range views," Proc. IEEE Conf. Computer Vision and Pattern Recognition, Illinois, 1992, pp. 348-353. [E.7-Turk and Levoy 94] G. Turk and M . Levoy. "Zippered polygon meshes from range images," Computer Graphics (SIGGRAPH '94 Proceedings), pp. 311-318. [E:8-Vemuri and Aggarwal 86] B .C . Vemuri and J .K. Aggarwal. "3-D model construc-tion from multiple views using range and intensity data," Proc. IEEE Conf. Com-puter Vision and Pattern Recognition, 1986. F. Shape transformation [F.l-Beier and Neely 92] T. Beier and S. Neely. "Feature-based image metamorphosis," Computer Graphics (SIGGRAPH '92 Proceedings), Vol. 26, No. 2, July 1992, pp. 35-42. Bibliography 162 [F.2-Burtnyk and Wein 76] N . Burtnyk and M . Wein. "Interactive skeleton techniques for enhancing motion dynamics in key frame animation," Communications of the ACM, Vol. 19, 1976, pp. 564-569. [F.3-He et al. 94] T. He, S. Wang, and A . Kaufman. "Wavelet-based volume morphing," Proceedings of IEEE Visualization '94, 1994, pp. 85-92. [F.4-Hong et al. 88] T . M . Hong, N . Magnenat-Thalmann, and D. Thalmann. "A general algorithm for 3-D shape interpolation in a facet-based representation," Proceed-ings of Graphics Interface '88, June 1988, pp. 229-235. [F.5-Hughes 92], J. Hughes. "Scheduled Fourier volume morphing," Computer Graphics (SIGGRAPH '92 Proceedings), Vol. 26, No. 2, July 1992, pp. 43-46. [F.6-Kaul and Rossignac 91] A . Kaul and J. Rossignac. "Solid-interpolating deforma-tions: construction and animation of PIPs," Proceedings of Eurographics '91, pp. 493-505. [F.7-Kent et al.] J.R. Kent, W.E . Carlson, and R.E. Parent. "Shape transformation for polyhedral objects," Computer Graphics (SIGGRAPH '92 Proceedings), Vol. 26, No. 2, July 1992, pp. 47-54. [F.8-Lerios et al. 95] A . Lerios, C D . Garfinkle, and M . Levoy. "Feature-based volume metamorphosis," Computer Graphics (SIGGRAPH '95 Proceedings), Vol. 29, No. 2, August 1995, pp. 449-456. [F.9-Payne and Toga 92] B. Payne and A . Toga. "Distance field manipulation of surface models," IEEE Computer Graphics and Applications, Vol. 12, No. 1, January 1992, pp. 65-71. [F.lO-Reeves 81] W.T. Reeves. "Inbetweening for computer animation utilizing moving point constraints," Computer Graphics (SIGGRAPH '81 Proceedings), August 1981, pp. 263-269. [F.ll-Schmitt 91] F. Schmitt, X . Chen, and W.-H. Du. "Geometric modelling from range image data," In W. Purgathofer, editor, Eurographics '91, North-Holland, September 1991, pp. 317-328. [F.12-Sederberg et al. 93] T .W. Sederberg, P. Gao, G. Wang, and H . Mu. "2D shape blending: An intrinsic solution to the vertex path problem," Computer Graphics (SIGGRAPH '93 Proceedings), August 1993, pp. 15-18. [F.13-Sederberg and Greenwood 92] T .W. Sederberg, and E. Greenwood. "A physically based approach to 2-D shape blending," Computer Graphics (SIGGRAPH '92 Proceedings), July 1992, pp. 25-34. Bibliography 163 [F.14-Shapira and Rappoport 95] M . Shapira and A . Rappoport . "Shape blending using the star-skeleton representation," IEEE Computer Graphics and Applications, V o l . 15, No. 2, March 1995, pp. 44-50. G . Miscellaneous [G. l -Amanat ides 87] J . Amanatides and A . Woo. " A fast voxel traversal algorithm for ray tracing," Eurographics '87, August 1987, pp. 1-10. [G.2-Barr 81] A . H . Barr . "Superquadrics and angle-preserving transformations," IEEE Computer Graphics and Applications, V o l . 1, No. 1, January 1981, pp. 11-23. [G .3 -Blum 77] H . B l u m and R. Nagel. "Shape description using weighted symmetric axes features," Proceedings of the IEEE Conference on Pattern Recognition and Image Processing, June 1977, pp. 203-215. [G.4-Borgefors 84] G . Borgefors. "Distance transformations in arbitrary dimensions," Computer Vision, Graphics, and Image Processing, V o l . 27, 1984, pp. 321-345. [G.5-Catmul l 78] E . Ca tmul l . "The problems of computer-assisted animation," Com-puter Graphics (SIGGRAPH '78 Proceedings), August 1978, pp. 348-353. [G.6-Drebin et al. 88] R . A . Drebin, L . Carpenter, and R Hanrahan. "Volume render-ing," Computer Graphics (SIGGRAPH '88 Proceedings), V o l . 22, No. 4, August 1988, pp. 65-74. [G.7-Edelsbrunner and Miicke 90] H . Edelsbrunner and E .P . Miicke. "Simulation of sim-plicity: A technique to cope wi th degenerate cases in geometric algorithms," IEEE Transactions on Graphics, V o l . 9, No. 1, January 1990, pp. 66-104. [G.8-Edelsbrunner et al. 92] H . Edelsbrunner, P. Fu , N . Akki ra ju , and C . Varela. "Shape2D v l . l , " Department of Computer Science, and Nat ional Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, 1993. [G.9-Forsey 91] D . Forsey. " A surface model for skeleton-based character animation," Eurographics Workshop on Animation and Simulation, September 1991, pp. 55-73. [G. lO-Ghosh and Jain 93] P . K . Ghosh and P K . Jain. " A n algebra of geometric shapes," IEEE Computer Graphics and Applications, V o l . 13, No. 5, September 1993, ' pp. 50-59. Bibliography 164 [G.ll-Heckbert 82] P. Heckbert. "Color image quantization for frame buffer display," Computer Graphics Proceedings of S I G G R A P H '82, Vol. 16, No. 3, July 1982, pp. 297-307. [G.12-Herbison-Evans 82] D. Herbison-Evans. "Real-time animation of human figure drawings with hidden lines omitted," IEEE Computer Graphics and Applica-tions, Vol. 2, No. 9, November 1982, pp. 27-33. [G.l3-0vermars 92] M . H . Overmars. "A graphical user interface toolkit for Silicon Graphics workstations, v2.1," Department of Computer Science, Utrecht Uni-versity, the Netherlands, 1992. [G.14-Sabella 88] P. Sabella. "A rendering algorithm for visualizing 3D scalar fields," Computer Graphics (SIGGRAPH '88 Proceedings), Vol. 22, No. .4, August 1988, pp. 51-58. [G.15-Sederberg and Parry 86] T .W. Sederberg and S.R. Parry. "Free-form deformation of solid geometric models," Computer Graphics (SIGGRAPH '86 Proceedings), Vol. 20, No. 4, August 1986, pp. 151-160. [G.16-Welzl 91] E. Welzl. "Smallest enclosing disks, balls and ellipsoids," Tech. Report No. B 91-09, Fachbereich Mathematik, Freie Universitat Berlin, Berlin, Germany, 1991. H . Books J [H.l-Bookstein 91] F .L . Bookstein. Morphometric Tools for Landmark Data: Geometry and Biology, Cambridge University Press, Cambridge, England, 1991. [H.2-Chui 92] O K . Chui. An Introduction to Wavelets, Academic Press, Boston, 1992. [H.3-Duda and Hart 73] R.O. Duda and P.E. Hart. Pattern Classification and Scene Analysis, Wiley, New York, 1973. [H.4-Foley et al. 90] J.D. Foley, A . van Dam, S.K. Feiner, and J.F. Hughes. Computer Graphics: Principles and Practice, Addison-Wesley, 1990. [H.5-Kaufman 91] A . Kaufman. Volume Visualization, I E E E Computer Society Press, 1991. [H.6-Lord and Wilson 84] E .A. Lord and O B . Wilson. The Mathematical Description of Shape and Form, Ellis Horwood Limited, Chichester, England, 1984. [H.7-Marr 82] D. Marr. Vision, W . H . Freeman and Company, New York, 1982. Bibliography 165 [H.8-Mortenson 85] M . E . Mortenson. Geometric Modeling, John Wiley and Sons, 1985. [H.9-Paul 81] R. Paul. Robot Manipulators, M.I.T. Press, Cambridge, M A , 1981. [H.lO-Pavlidis 82] T. Pavlidis. Algorithms For Graphics And Image Processing, Com-puter Science Press, Rockville, M D , 1982. [H.ll-Preparata and Shamos 85] F.P. Preparata and M.I. Shamos. Computational Ge-ometry: An Introduction, Springer-Verlag, New York, 1985. [H.12-Rosenfeld and Kak 82] A . Rosenfeld and A . C . Kak. Digital Picture Processing, 2nd edition, Academic Press, New York, 1982. [H.13-Serra 82] J. Serra. Image Analysis and Mathematical Morphology, Academic Press, London, 1982. [H.14-Tarjan 83] R.E . Tarjan. Data-structures and network algorithms, S IAM, Philadel-phia, 1983. I. Publications arising from this thesis [I.l-Ranjan and Fournier 96] V . Ranjan and A . Fournier. "Matching and interpolation of shapes using unions of circles," to appear in: Proceedings of Eurographics '96.. [I.2-Fournier and Ranjan 96] A . Fournier and V . Ranjan. "Tessellation of union of spheres with spherical triangles," manuscript for Technical Report, Department of Computer Science, University of British Columbia. [I.3-Fournier and Ranjan 95a] A. Fournier and V . Ranjan. "Drawing curves with arcs of circles," Technical Report, Department of Computer Science, University of British Columbia. . [I.4-Ranjan and Fournier 95b] V . Ranjan and A . Fournier. "Creating union of spheres models from multiple views," to appear in: Journal of Graphics Tools. [I.5-Ranjan and Fournier 95c] V . Ranjan and A. Fournier. "Union of spheres (UoS) model for volumetric data," Short communication, Eleventh Annual Symposium on Computational Geometry, Vancouver, Canada, June 1995. [I.6-Ranjan and Fournier 94] V . Ranjan and A . Fournier. "Volume models for volumetric data," IEEE Computer, special issue on volume visualization, Vol. 27, No. 7, July 1994, pp. 28-36. Notation O The abstract concept of an object B(0) Boundary'of O 1Z The abstract concept of a model for the object B{1Z) Boundary of 1Z V A set of points on the boundary of an object T A transformation Instances of objects and models Model for object A using R representation. For example, Model for A using boundary points representation Union of spheres model for A derived from PA constants scalar quantities vector quantities A matrix p, q points, voxels, spheres pq line c curve S = {si, s2,.. •, sn} A set of spheres U(S) Union of spheres in S A, B RA PA UoS(P A ) e, S, Tyo l,d n. d M 166 Notation 167 D Diameter A Area V Volume FQ, /() ' Functions 0() Order d(,) Distance between two arguments, heavily overloaded; meaning is clear from the arguments and/or subscript. Sometimes, pseudo-distance, i.e., does not maintain the triangle inequality. djs(,) Euclidean distance between two points dc{i) Length of curve between two points on the curve UoT Union of Triangles UoC Union Of Circles UoTet Union of Tetrahedra UoS Union of Spheres Glossary B Boolean operations — Set operations such as union, intersection, and difference. boundary points — Points on the surface of an object or its model. boundary cell — A cell through which the surface of the object passes. boundary point density — A measure of the denseness of sampled boundary points. A boundary point density of 5 means that a sphere with any boundary point as its center and radius greater than S necessarily contains some sampled boundary point. 168 Glossary 169 C cell — A region of volumetric data defined by voxels (see Figure 63). D Delaunay triangulation — Triangulation which is a dual of the Voronoi diagram. distance — A distance measure d(x, y) means: d(x, y) > 0, d(x, y) = 0, if x = y d(x,y) = d(y,x), and d{x, y) < d(x, z) + d(z, y) (triangle inequality) We will call pseudo-distance a quantity that verifies everything except the triangle in-equality. distance transform — A n image analysis tool for analyzing binary voxel data; it gives the approximate distance from every interior voxel to the nearest exterior voxel. E empty circle property (ECP) — The property of Delaunay triangulations which says that the circumscribing circle (sphere) of any triangle (tetrahedron) cannot contain any other vertex. exterior grid-point — A grid-point that is outside the defined object. exterior voxel — A voxel that is outside the defined object. G general position — In 2D, the condition that no 3 points are collinear, and no 4 points lie on a circle. In 3D, the condition that no 4 points lie in a plane, and no 5 points lie on a sphere. grid-point — Point in 3-space for which we have data. Refer to volumetric data. Glossary 170 H hierarchical — In this thesis, it means a model or display method organized in a manner such that the model/display with coarser details is a subset of the model/display with finer details. I interior grid-point — A grid-point that is inside the defined object. interior voxel — A voxel that is inside the defined object. iso-value — A particular value chosen from the range of values in a volumetric data. iso-surface — A surface defined in the volumetric data by points at the chosen iso-value. M manifold — A manifold surface has the property that, around every one of its points, there exists a neighborhood that is homeomorphic to the plane. A manifold is orientable if we can distinguish its two different sides. Closed, orientable manifolds partition the space into three regions, the inside, the surface, and the outside. manipulation — A vague term meaning any operation that can be performed. Opera-tions can be transformations (translation, rotation, dilation), Boolean operations, defor-mations (skew, bend, stretch etc.), interpolation etc. medial axes The locus of centers of all maximal spheres inside an object; the skeleton or symmetric axes (surfaces in 3D). median point — Center of mass. model — A specific instance of a representation. O object — A closed region in.3-space (in 2-space for 2D) that we want to model and visualize. A n object can have several disconnected components. Glossary 171 P power diagram — A generalization for spheres of the Voronoi diagram for points, primitives — The elements from which a representation is.comprised. R raycasting — A technique to render images by tracing imaginary rays from the eye and through the viewing plane and intersecting it with objects in the scene, registration — The process of alignment of two data sets. These data sets may have been acquired by different imaging modalities or by the same modality at different times, representation — A data-structure and computational method to describe objects in a computer. S skeletonization — See thinning. stability — The property of representations to inhibit large changes in the model for a small change in the data, such as noise or small distortions. T thinning — An image analysis method which derives stick-figure like model (called skeleton) of object from binary voxel data. threshold — A value used for binary classification. For example, in volumetric data any voxel with value greater than the threshold can be treated as inside the object, and with value less than the threshold to be outside. topological properties — Properties of an object that do change with its continuous deformation such as the number of connected components, holes, and cavities. V visualization — The use of graphical media to interpret and understand complex data sets. Often used abusively as synonymous with "display." / Glossary 172 volume models — Models constructed using volume primitives. volume primitives — Primitives in 3D space which have definite interior and exterior regions, e.g., spheres, ellipsoids, superquadrics etc.. volume rendering — A way to display volumetric data without using an intermediate geometric model. volumetric (or 3D) data — Three-dimensional discrete data made of voxels, i.e., (lo-cation, value) tuples. This data can be seen as a 3D array of grid-points denoted by the location field in the tuple with a value available at all grid-points. Examples include med-ical, astrophysical, computational fluid dynamics (CFD), geological, and oceanographic data. ; volumetric object — Object or region of interest defined in volumetric data. Voronoi diagram — Given a set of vertices, their Voronoi diagram is a set of cells so that every point within a cell is closer to one vertex than to any other vertex, voxel — A volume element. Generally, it means a (location, value) tuple in the volu-metric data. Each voxel is usually representative of properties of some volume around the grid-point denoted by the location field in the tuple. ( Index alignment of models, 100 alpha-shape, 22 approximation properties general approach, 76 of union of arcs, 78 of union of circles, 81 of union of spheres, 86 of union of tetrahedra, 86 of union of triangles, 82 blobbies, 24, 66 Boolean operations, 169 boundary cell, 169 point density, 7, 169 points, 169 BSP tree representation, 21 cell, 170 clustering, 58 contributions, 12, 149 core representation, 26 data 2D calf, 135 cow, 48 fish, 91 giraffe, 66 M R I brain slice, 48, 65 rectangle, 133 3D C T femur, 48 hand, 53, 73 head, 91 knot, 48, 66 P E T brain, 48 rabbit, 48, 66 sources, 4 Delaunay triangulation, 170 display blended spheres, 31 union of spheres, 30 distance, 170 between models, 104 between objects, 121 between spheres, 104 feature, 104, 105 normalization, 122 position, 104 size, 104 transform, 170 weights in, 104 empty circle property, 35, 170 error Conditions A , B , 9 definition, 8 in polylines, 80 in union of arcs, 81 in union of circles, 83 in union of spheres, 86 in union of tetrahedra, 86 in union of triangles, 82 feature definition, 108 distance, 105, 110 examples, 112 general position, 170 grid-point, 170 exterior, 170 interior, 171 hierarchical, 171 histogram 173 Index use in display, 48 inside-outside test for polygon-mesh data, 47 for range data, 43 for ray-cast data, 47 for volumetric data, 38 interpolation examples, 132 problem definition, 97 using UoS, 132 iso-surface, 171 iso-value, 171 manifold, 171 manipulation, 171 marching cubes, 16 mass properties from UoS, 27 matching applications, 114 examples 2 D , 135 3 D , 141 general strategies, 95 problems, 94 spheres, 112 using UoS, 91, 99 medial axes, 171 median point, 171 model, 171 definition, 5 moments, 100 object, 5, 171 recognition, 94 shape, 6 volumetric, 173 octrees, 18 polylines, 80 power cell, 28 power diagram, 28, 64, 171 primitives, 172 properties desirable, list of, 9 intrinsic, 9 ' non-intrinsic, 9 radius field, 105 radius gradient, 106, 107 raycasting, 172 registration, 94, 172 examples, 117 problem definition, 96 using UoS, 115 representation, 172 BSP tree, 21 core, 26 definition, 5 octree, 18 shell, 19 skeletal, 25 surface-based, 15 tetrahedra, 21 volume primitives, 20 volume-based, 18 voxel-based, 18 wavelet, 19 r-regular, 76 sampling density, 7 segmentation using UoC, 101 using UoS, 103 shape comparison, 94 interpolation, 94 of an object, 6 shell representation, 19 simplification algorithm, 62 examples, 65 motivation, 58 Index of UoS, 27 skeleton representation, 25 skeletonization, 172 sphericity definition, 60 selection, 65 stability definition, 9, 172 example, 59 of UoS, 126 surface models from contours, 17 from range data, 18 from volumetric data, 16 tetrahedra representation, 21 thinning, 172 threshold, 172 topological properties, 172 transformation non-shape-preserving, 6 shape-preserving, 6 Union of arcs, 75 approximation properties, 80 Union of circles, 75 approximation properties, 82 generation, 39 simplification, 65 Union of spheres approximation properties, 86 blended, 23, 31, 66 definition, 10 display methods, 30 . generation, 23 approach, 34 examples, 47 from computer models, 46 from polygon-mesh, 47 from range data, 42 from volumetric data, 37 motivation, 10 pipeline, 12 polygonization, 32 simplification, 27 Union of tetrahedra, 75 approximation properties, 86 Union of triangles, 75 approximation properties, 82 generation, 39 verification of inside spheres, 45 of simplified spheres, 64 visualization, 172 volume models, 172 primitives, 172 rendering, 173 volumetric data, 173 object, 173 Voronoi diagram, 173 voxel, 173 exterior, 170 interior, 17.1 wavelet representation, 19
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A union of spheres representation for 3D objects
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A union of spheres representation for 3D objects Ranjan, Vishwa 1996
pdf
Page Metadata
Item Metadata
Title | A union of spheres representation for 3D objects |
Creator |
Ranjan, Vishwa |
Date Issued | 1996 |
Description | Reconstruction of an object from a set of points sampled from its boundary is an important problem in graphics and vision. Several methods exist to compute and display surface (e.g., polygonal) and volumetric (e.g., polyhedral) models of objects from the boundary points. In order to display, transform, and compare objects, it is often convenient or necessary to use different representations of objects. Basic desired properties of representations of objects are efficiency of computation, storage, and display. Other important properties include stability (small changes in the data, such as noise or small distortions, cause small changes in the model), the ability to determine the similarities between two data sets, and the computation of simplified models. A survey of common representations of objects (e.g., surface, octrees, etc.) shows that some important properties are lacking in these representations. In this thesis we study a union of spheres representation (UoS) for the volume enclosed by an object's boundary. We present algorithms to generate stable union of spheres models of objects from various sources of data, such as volumetric data (e.g., data from CT or MRI scanners), range data, and other existing models. The spheres can be simplified to obtain multi-scale models. We present a method to establish correspondence between two objects using their union of spheres models and use this to calculate distances between objects, to register objects, and to interpolate between objects. This establishes a measure to study and compare such models. Examples with simple and complex objects show how this measure corresponds closely to the intuitive human understanding of shape. |
Extent | 23619288 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051598 |
URI | http://hdl.handle.net/2429/6140 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1996-148181.pdf [ 22.53MB ]
- Metadata
- JSON: 831-1.0051598.json
- JSON-LD: 831-1.0051598-ld.json
- RDF/XML (Pretty): 831-1.0051598-rdf.xml
- RDF/JSON: 831-1.0051598-rdf.json
- Turtle: 831-1.0051598-turtle.txt
- N-Triples: 831-1.0051598-rdf-ntriples.txt
- Original Record: 831-1.0051598-source.json
- Full Text
- 831-1.0051598-fulltext.txt
- Citation
- 831-1.0051598.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0051598/manifest