MULTIRESQLUTION SURFACE CONSTRUCTION FOR HIERARCHICAL B-SPLINES by DAVID W O N G B . S c , Simon Fraser University, 1994 A THESIS SUBMITTED IN P A R T I A L F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE IN T H E F A C U L T Y O F G R A D U A T E DEPARTMENT OF COMPUTER STUDIES SCIENCE We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH A p r i l , 1995 Copyright & © David Wong, 1995 COLUMBIA In presenting degree freely thesis in partial fulfilment at the University of British Columbia, available copying this for reference of this department publication or thesis by for scholarly his of this and study. thesis or her for financial of Ca*npt>t&-r Sc. The University of British Columbia Vancouver, Canada D DE-6 a t e (2/88) 11^ A .i , r mi* the I agree I further agree requirements that that gain shall It is for an advanced the Library shall make it permission for extensive purposes may b e granted representatives. permission. Department of by the head understood not be allowed that without of my copying or my written ABSTRACT This thesis presents a method for automatically generating a hierarchical B-spline surface suitable for animation and interactive modification from an initial set of control points. Given an existing mesh M[ ^ of control points, a mesh with half the resolution M ^ , is fc+1 constructed by simultaneously approximating the finer mesh and minimizing a smoothness or fairness constraint using weighted least squares. The membrane and thin plate models have been tried as smoothness constraints. Various numerical solving techniques have also been tested. Curvature measures of M ^ ! are used to identify features that need only be 1 represented in the finer mesh. Some hierarchical B-spline surfaces, a B-spline surface and a digitized surface have been tested by this method. The resulting hierarchical surface generated accurately and economically reproduces the original mesh, is free from excessive undulations in the intermediate levels, and is well behave for editing and animation. ii T A B L E OF CONTENTS Abstract ii Table of contents iii List of figures viii Acknowledgments x Dedication xi Chapter 1 Introduction 1 1.1 Modelling problems addressed by this thesis 2 1.2 Goal 7 1.3 Thesis overview 7 Chapter 2 Background 8 2.1 Problems with tensor product B-spline 8 2.2 Hierarchical B-splines 9 2.3 2.2.1 Basic formulation 10 2.2.2 Refinement 10 2.2.3 Offset from parent surface 11 2.2.4 Types of offset operators 11 Surface approximation 12 2.3.1 Regularization 12 2.3.2 Outliers . . . . . . . . . . . . iii . 13 Chapter 3 Related Work 3.1 3.2 3.3 3.4 14 Smooth approximation surfaces 14 3.1.1 Using fairness energy norms 14 3.1.2 Using particle systems 16 3.1.3 Using mesh optimization 17 Multiresolution approximations 17 3.2.1 Multilevel visual surfaces with discontinuities 17 3.2.2 Wavelets 18 3.2.3 Hierarchical B-splines 19 Outliers handling 20 3.3.1 The outlier process 20 3.3.2 Robust statistics 21 Optimal choice of the regularization parameters 21 3.4.1 Generalized cross-validation 21 3.4.2 Search methods 22 iv Chapter 4 Multi-resolution B-spline Approximation 24 4.1 Divide and conquer 24 4.2 Initial approximation 25 4.3 Smoothness constraints 26 4.4 The membrane model 27 4.4.1 The membrane model for noise removal 27 4.4.2 The effect of varying A on the membrane model 27 4.5 The thin plate model 30 4.5.1 The thin plate model for noise removal 31 4.5.2 The effect of varying A on the thin plate model 32 4.5.3 Sparse data set for smooth approximation 34 4.6 Least squares approximation with smoothness constraint 36 4.7 Weighted approximation 37 4.8 Calculating and thresholding the offsets 39 Chapter 5 Implementation 41 5.1 Algorithm 41 5.2 Setting up the matrices 42 5.3 5.2.1 The S matrix . 5.2.2 The P matrix 43 5.2.3 The W t 44 f c + 1 - l matrix Solving the system 42 44 v Chapter 6 6.1 6.2 6.3 Results 45 Hierarchical B-spline surfaces 45 6.1.1 Simple geometries 45 6.1.2 Complex geometries 49 6.1.3 Use in animation ; B-spline surface 53 6.2.1 53 Fixed tolerance versus Variable tolerance Digitized data Chapter 7 56 Conclusions 58 Bibliography Appendix A A.l 51 59 The multi-grid method 63 A brief introduction to multi-grid method 63 A. 1.1 65 The multi-grid method A.2 A multi-grid solver for noise removal using the membrane model 65 A.3 A multi-grid solver for noise removal using the thin plate model 66 A.3.1 Defining Stj 67 A.3.2 Interior cell 67 A.3.3 Boundary cell 68 A.3.4 Corner cell 68 A.3.5 Minimization 69 vi Appendix B Survey of solving methods 71 B.I The normal equations method 72 B.2 The augmented system method 73 B.3 The QR method 74 B.4 The conjugate gradient method 75 B.4.1 With preconditioning 75 B. 5 The regularization method Appendix C C. l 76 Solving the overdetermined system A and A* A 78 , 78 C.2 Kronecker product 80 C.3 FiU-in 81 C.4 Cholesky factorization 82 C.5 Comparison of the different solving methods 83 Glossary 86 vii LIST OF FIGURES 1.1 The original test surface at level 3 3 1.2 Comparison between smooth and least squares hierarchy 5 1.3 Comparison between smooth and least squares hierarchy under deformation by adjusting the same single control point by the same amount in level 1 6 2.1 Spiny B-spline surface, 16524 control points 8 4.1 The spike (19 X 19 control points) 24 4.2 Coarser level approximation by least squares 26 4.3 The original smooth surface data 28 4.4 The effect of varying A on the membrane model 29 4.5 The effect of varying A on the thin plate model 33 4.6 The thin plate model on sparse data set 35 4.7 Coarser level approximation, with smoothing A = 0.1 37 4.8 Weighting of control mesh 38 4.9 Coarser level approximation, using weighted 5.1 Three cases of the use of subdivision template viii fit. 39 43 6.1 Spike (28 patches, 17 offsets) 46 6.2 Hump (40 patches, 19 offsets) 47 6.3 Plop (28 patches, 19 offsets) 48 6.4 Eye (565 patches, 253 offsets) 49 6.5 Dog head (609 patches, 257 offsets) 50 6.6 Mask: Comparison of original and generated hierarchies 52 6.7 Reboot head: Generated from a B-spline surface with uniform tolerance. 6.8 Reboot head: Generated from B-spline with variable tolerance 55 6.9 Spock head: Generated from Cyberware data, 43520 data points 57 A.l One multi-grid iteration for V-cycle and W-cycle on a 5 levels grid 64 C.l The matrices A and A A with control mesh size 11x11 79 C.2 The matrices of L + L of A A with control mesh size 67 X 67 81 1 1 1 ix ... 54 ACKNOWLEDGMENTS I'd like to thank everyone who helped me with this research. M y thesis supervisor Dave Forsey who introduced me to this interesting research topic, continuously giving me good advice. He also donated a huge amount of his very valuable time for making sure that this thesis is well written and well organized. Moreover, he is the person that made me work overnight for the first time to submit our paper. Thanks to Alain Fournier, as my second reader who gave me very valuable advice on my thesis draft, for his valuable family time over the weekend to proof read my draft. He also taught me a lot of important concepts in the field of computer graphics. Thanks to my other graphics professor Kelly Booth, who widened my view in the field of computer graphics, for convincing me stay here at U B C to do my Masters. Uri Asher and Jim Varah taught me very useful numerical analysis in their courses. My friendly office mate, Jason Harrison, as my student reader for this thesis, gave me very valuable suggestions on organizing my thesis, and painstakingly corrected my grammar in English. My friend, David Moiilton, gave me good advice and helped me out with my first multi-grid method on the thin plate for noise removal. All my undergraduate professors at S F U taught me a solid background in the field of computing science. Last but not least, personal thanks to all my friends and relatives for their support. DEDICATION I dedicate this thesis to my M o m and Dad for their life of hard working xi CHAPTER O N E INTRODUCTION The geometry of the objects used in computer generated imagery and animation are generated either procedurally (by some algorithm), interactively (using a surface editing application), or by digitizing real objects. These 3D geometric models take many forms, usually represented in polygonal, Bezier or B-spline form. These representations are "single-resolution" in that there is only one scale relating all the data points forming the surface. A hierarchical B-spline (see Section 2.2) developed by Forsey [23] represents a surface using a series of surfaces definitions each at a different resolution. This provides a compact surface representation allows multiresolution surface editing that retains surface detail during broadscale surface manipulation, a property that is very useful during the creation and animation of free-form surfaces. However, hierarchical B-spline surfaces are currently constructed using Dragon, an interactive surface modelling system developed by Forsey [23] and the ability to construct hierarchical surfaces from other sources of 3D data would considerably widen the applicability of the hierarchical approach. This thesis presents a scheme that, given an initial control mesh for a tensor-product Bspline surface, will construct a multiresolution hierarchical B-spline surface with smooth 1 intermediate level surfaces that are appropriate'for use in modelling and animation. We assume that the initial mesh already sufficiently represents the desired surface and the parameterization of the underlying surface is'given. We are concerned solely with constructing the multi-resolution hierarchical B-spline structure. To create this initial mesh from polyhedral surface data, it must be first converted into a Smooth or fair in the context of this thesis is .with respect to the underlying control mesh not to the surface which the mesh defined, and it is measured with energy norms like the thin plate or the membrane model. 1 1 Chapter 1: Introduction 2 B-spline surface by a B-spline approximation routines [24]. Similarly, a preprocessing step is needed to first convert other spline-based representations (like Bezier) to B-spline form. The problems associated with these conversions are beyond the scope of this thesis. 1.1 MODELLING PROBLEMS ADDRESSED B Y THIS THESIS Surface approximation is the process of approximating a surface in a particular representation given a set of data points. The most common approximation technique is finding the surface that minimizes the sum of squares of all the errors that measure the distances from the approximation surface to the data points. This is called a least squares approximation. Surface approximation using a hierarchical B-spline is the process of not only determining the coefficients that produces a surface that adequately represents the data, but also for the intermediate surfaces levels as well. This process is complicated by the fact that there is an infinite number of intermediate level surfaces that would produce identical approximation. This is similar to the basic.approximation problem where there are fewer number of degree of freedom in representation than in the data. In such cases, least squares is often used. We can also use least squares approximation to generate a hierarchical B-spline representation. However, a multi-resolution representation created by least squares approximation for hierarchical B-splines [24] or wavelet decomposition [19] produces filtered surfaces that may undulate to accommodate high frequency components. For example, the surface in Figure 1.1 is essentially flat with two spikes protruding from it. Figure 1.2 shows the result of approximated this surface using least squares at multiple resolutions. This approach is quite satisfactory if the goal is to simply represent a static surface, given a surface of fine enough resolution (i.e. sufficient degrees of freedom) and the approximating surface with fit the data exactly. However one of the main strengths of the hierarchical representation is to have a surface where small scale details follow the broad changes in shape of the surface (see Figure 1.3a). If a simple least squares approach is used to create the intermediate Chapter 1: Introduction 8 (a) Perspective view A (b) Orthogonal view Figure 1.1: The original test surface at level 3. levels of the hierarchy, the additional undulation in the coarser level representations cause non-predictable behaviour of the surface when it is being manipulated (see Figure 1.3b) For hierarchical B-spline surfaces constructed by hand, the coarser level surfaces are typically not "blurred" versions of the final form, but are built up in much the same way as a clay sculpture is; a basic shape is first generated to which details are added. The hierarchy for Figure 1.1 would be modelled as a flat plane at all intermediate levels with the two spikes added at the finest resolution surface. This not only reduces the number of offsets needed to store the representation but a hierarchy with smooth intermediate level surfaces is desirable for multi-scale modelling and animation purposes because broad-scale editing of the surface produces reasonable and predictive results, where details follow the global change in shape. (How this is accomplished is covered in Section 2.2.) However, if we apply the same technique to the undulating hierarchy generated by least squares, unexpected distortions of the surface arise. For example, in Figure 1.3 are two representations, one created by hand and the other using a least squares approximation for the multiresolution representation. If the same left edge middle control point in each representation is moved by the same amount the results are quite different. A reasonably predictable result is obtained from the smooth hierarchy where the left spike follows the change in global shape, Chapter 1: Introduction 4 whereas an unexpected distortion occurs in the other representation: the left spike became much thinner, the unrelated right spike got much thicker and another small spike in the middle pointing downward was introduced. This type of behaviour is unacceptable for a representation that must be edited or animated. Thus, the notion of a surface that minimizes the positional error is too simplistic. Another metric must be found that can more fully recreate the type of hierarchical B-spline surfaces that are created by hand. It is unlikely to quantitatively determine how an expert might build these surface, however this thesis demonstrates that a hierarchy with intermediate levels that contain no extra undulations will produce a better behaved surface. Chapter 1: Introduction (e) Smooth - Level 0 5 (f) Least squares - Level 0 t Figure 1.2: Comparison between smooth and least squares hierarchy Chapter 1: Introduction 6 Figure 1.3: Comparison between smooth and least squares hierarchy under deformation by adjusting the same single control point by the same amount in level 1. Chapter 1: Introduction 1.2 7 G O A L The goal of this thesis is, given an initial control mesh M ^ m e s ' l of a B-spline surface, which may be generated from data, to automatically generate a hierarchical B-spline representation that, if possible, satisfies the following: e The representation should accurately approximate the original mesh. • Intermediate levels are "smooth" according to some measure. • The storage requirement is reduced. 1.3 T H E S I S O V E R V I E W The remainder of this thesis is structured as follows: -Chapter 2 provides the necessary background on hierarchical B-spline surfaces and surface approximation. Chapter 3 presents a literature search on related work. Chapter 4 describes the basic approach for automatically constructing hierarchical B-spline surfaces. Chapter 5 details the algorithm and implementation. Chapter 6 shows the results of applying the algorithm to a variety of inputs. The conclusions and future work follows in Chapter 7. Various numerical techniques tested or surveyed for the solution method and a glossary are at the end of this thesis. CHAPTER T W O BACKGROUND This chapter describes the problems with tensor product B-splines, introduces the hierarchical B-spline surface representation, and presents the surface approximation problem. 2.1 P R O B L E M S W I T H T E N S O R P R O D U C T B - S P L I N E Figure 2.1: Spiny B-spline surface, 16524 control points. For example, the surface in Figure 2.1 as a B-spline surface, requires 16524 control points to define its shape because the non-local property of tensor product B-spline knot insertion forces the addition of knots along two bands in the two parametric directions on the surface in order to represent small scale local details. This non-local refinement property tremendously increases the number of degree of freedom (DOF) of a B-spline surface, resulting in increased time and space overhead when using splines for complicated multi-scale detailed surface approximation, rendering, or polygonizing for export to other non-spline based applications. 8 Chapter 2: Background 9 A tensor product B-spline surface is a single scaled representation. It requires manipulation of a large number of control points in order to make a broad-scale change to the surface shape. It is also cognitively very hard for a user to maintain the small-scale surface details without distortion while editing at a broader scale. The unavailability of multi-scale editing greatly increases the time and cost of manipulating free-form shapes for tasks such as facial animation. Hierarchical B-spline [23] is a better way of representing multi-scaled detailed- surfaces which allows for multiresolution surface editing and displaying, and at the same time, economical storaging. 2.2 H I E R A R C H I C A L B - S P L I N E S In 1988 Forsey and Bartels [23] introduced the concept of locally refining a B-spline surface. In this scheme a surface is represented as a series of levels called overlays, each of which consists of a collection of B-spline surfaces with twice the resolution of the parent surface at level k — 1. It is a multi-resolution surface representation that provides: 1. Local refinement of tensor-product spline surfaces that localizes detail where it is required. 2. Multi-resolution surface editing that retains surface detail during broad-scale surface manipulation. 3. Multi-resolution surface displaying. 4. Economical surface representation. (The surface in Figure 2.1 requires only 178 data points as a hierarchical spline.) 5. Multi-resolution animation capabilities that make it easier to animate complicated surfaces [22]. Chapter 2: Background 2.2.1 BASIC 10 FORMULATION Let, H (u,v) [k] J2Y. i3 ffi(«) i*i( ) v = t=0 where surface H^(tt, v) is defined by an n B B w (2-i) j=o X m mesh of control vertices vj*- , and piecewise polynomial basis functions, B ^ ( i ) and Bj ].(<) , of degree K and r respectively. The parameters fc u and v vary independently from some minimum u i , m n u T O t n to some maximum u max , v . max Certain values of u and v, called knots, correspond to the joints between polynomial patches and form a knot sequence in each parameter {uo,..., u m + K } , and {vo,..., v + }. The basis n T functions are assumed to have local support that is, each basis function affects only a restricted area of the spline. To be used in a hierarchical surface representation the bases are required to be refinable in the sense that each one can b.e re-expressed as a linear combination of one or more "smaller" basis functions with additional knots without affecting the shape of the surface. 2.2.2 REFINEMENT Let be the initial mesh at level 0, denned by control vertices V ^ . Points R.M in the refined mesh are computed by linear combinations of the' V l ° l . Expressed in matrix form for general level k: [k+i] R where R i f c + 1 ] and = [k] sv denote column vectors of their corresponding mesh points. (2.2) These notations are overloaded, but readers should be able to deduce whether they refer to a mesh of points or column vectors of those points within context. The components of the subdivision matrix S depend upon the order and type of basis functions, the number and location of the inserted knots and the topology of the surface. In tensor- Chapter 2: Background 11 product surfaces, refinement is non-local in that knots are added in the parametric domain resulting in the addition of an entire new row or column of patches across the entire surface. O F F S E T 2.2.3 F R O M P A R E N T S U R F A C E In a hierarchical B-spline surface representation, each control vertex vj*^" ' in mesh M^" " ! is 1 1 1 represented in offset form as: [fc+i] V where R i fc+1 [fc+i]0o[ ] fc+1 = (2.3) R ] is derived by mid-point refinement of V ' * l as in Equation 2.2, O ^ 1 ! is the offset vector and "0" is the reference-offset operator. The position of V,-** ^ directly after a new level is created is equal to Rfj "^ and, by definition, 1 1 o[-*^ = 0. Any change to the position of a control node VJ'j*" ^5 1 vector O J - * * 1 ^ 1S represented in the offset as a relative change in position from the reference point R [ - * ^ . Changes to V,-^" ^ (and thus to pj^" ') do hot change the position of the reference point. On the other 1 1 hand, changes in the shape of a parent surface alter RJ** ', result in a procedurally change of 1 vj** ' as defined by the reference-offset operator. 1 The combined reference plus offset representation with the locally supported B-spline basis functions, allows a hierarchical B-spline surface to be locally refined so that patches can be added to a restricted region of the surface [23]. Additionally, this decreases the storage required in comparison to a B-spline surface because the non-zero offsets completely define the shape of the surface. 2.2.4 T Y P E S O F O F F S E T O P E R A T O R S Offsets correspond to the difference between the shape of the surface at two neighbouring levels of representation. The offset operator .determines how the difference is interpreted. In the simplest case the operator is simply vector addition and is, in operation, similar to wavelet Chapter 2: Background 12 decomposition. However the offset operator can be arbitrarily complex. For example, it could use a local frame of reference derived from the tangent and normal of the parent surface. The use of the offset operator enhances the behaviour of the surface so that detailed features follow broad scale deformations of the surface [23]. Other offset operators include deformations induced by an underlying skeleton [22]. The hierarchical approach is not restricted to B-splines but is applicable to any representation where the basis functions have local support and are refinable. For example, Finkelstein and Salesin [21] uses a wavelet basis while incorporating the offset form for the wavelet coefficients. 2 . 3 S U R F A C E A P P R O X I M A T I O N Surface approximation is the task of approximating a surface with a particular representation with a particular norm ||.||, given a set V of data points d{. That is finding a surface S(u, v) in the given representation, and a parameterization where U{ and V{ is associated with each data point d{ such that it satisfies min^\\S(ui,Vi) - d \\ { ten The data may be obtained by physical measurements of an existing surface. The set of data points may be scattered around in space with no explicit connectivity information. The surface representation may be just a mesh of points in space or a surface defined by some bases functions with control points. In our case, we consider the finer mesh M [ f c + 1 l as the set of data points given and the approximation surface is the coarser level mesh M ^ . Therefore, we focus on grid data approximation, and ignore the problem of parameterizing the resulting surface. 2 . 3 . 1 R E G U L A R I Z A T I O N The surface approximation problem is generally ill-posed in the sense that the solution to the problem may not be unique (ie., given a different parameterization), or the solution is very sensitive to a change in the data (ie., the solution surface can wildly fluctuate with a small Chapter 2: Background 13 change in data). The non-uniqueness problem is avoided by specifying a parameterization. The choice of a norm can affect the degree of sensitivity. However, in general it still remains. The standard way [55] [56] to handle ill-posed problems is to regularize them by introducing a stabilizer function with a regularization parameter usually denoted by A . The regularization parameter is used to control the degree of regularization. The stabilizer functions for surface approximation are usually some kind of energy functional of surface stretching or blending, like the membrane and the thin plate models [52]. These energy functipnals are also used extensively in the modelling of smooth or fair deformable surfaces [44] [16] [31] in computer graphics. A deformable surface is one which alters its shape in response to stress or force. 2.3.2 O U T L I E R S A n outlier in a statistical data set is a data point which deviate significantly from the main body of the data set. Outliers are usually arose from noise in measurements. The presence of outliers in data is a serious problem in surface approximation. The outliers can bias most norms, like the Lz-noxm. This resultsin fitting the outliers most at the expense of other data. In our case of automatic construction for hierarchical B-spline surface, we can think of the set of control points V ' fc+1 ^ representing the appropriate frequency components in the level k + 1 that can not be adequately represented in the coarse level k where only lower frequencies can be represented, as outliers in the point of view from the coarse level k. Having methods of eliminating the outliers from data set used for surface approximation is important, and the handling of outliers in statistics may give us insight on how to handle them in our geometric case. CHAPTER RELATED THREE W O R K The chapter presents related work in the following areas: smooth approximation surfaces, multiresolution approximations, outliers and how to find the optimal choice of a regularization parameter. 3 . 1 S M O O T H A P P R O X I M A T I O N S U R F A C E S The desire for smooth approximation surfaces in Computer Graphics is related to free-form deformable surfaces. The models for deformable surfaces have generally used fairness energy norm, like a derivative of the membrane model, the thin plate model, or the minimum variation of curvature model. Others like oriented particle systems or mesh optimization exist. A l l of these methods are used to model fair surfaces. 3 . 1 . 1 U S I N G F A I R N E S S E N E R G Y N O R M S Qin and Terzopoulous [44] developed a dynamics version of NURBS model, called D-NURBS. D-NURBS is a physics-based model that incorporates mass distribution, internal deformation energies, and other physical quantities into the popular NURBS geometric representation. It adopts the thin-plate under tension energy model [52] for deformation energy. The model produces a set of dynamic differential equations derived from Lagrangian dynamics [28], which are integrated numerically, and automatically evolve the control points and weights such that the D-NURBS move and deform in a "physically intuitive" manner. Sculpting tools adjust control points and weights, or apply forces and specify local or global shape constraints. For a complex surface with hundreds of control points, it can not achieve interactive rates on SGI Indigo machines. 14 Chapter 3: Related Work 15 Clough and Tocher [16] developed finite element stiffness matrices for analysis of plate blending. The matrices can be used to model a deformed plate. The advantage of the Clough-Tocher interpolant is that the construction is local. However, not all the remaining degrees of freedom are directly constrained by the data. The unconstrained degrees of freedom are often set using local heuristics and typically result in surfaces that are not fair. Celniker and Gossard [13] extended the Clough-Tocher interpolation method to set the remaining degrees of freedom to minimize a fairness norm. The fairness norm they use is based on a linear combination of thin-plate and membrane energies, and is quadratic, so it can be minimized by solving a sparse linear system. Halstead, Kass and DeRose [31] developed an efficient and fair interpolation using CatmullClark surfaces. It allows the construction of a smooth surface that interpolates the vertices of a mesh of arbitrary topology type. Normal vectors can also be interpolated at an arbitrary subset of the vertices. The basic idea is to compute a control mesh whose Catmull-Clark subdivision surface interpolates the given data and minimizes a smoothness or fairness measure of the surface. It follows Celniker and Gossard's work [13] in the use of a fairness norm based on a linear combination of the thin-plate and membrane energies. The approach taken is first to solve the system of interpolation constraints to yield a new control mesh, subdivide this new mesh once to add new degrees of freedom, and use these new DOFs to minimize the fairness norm subject to the interpolation constraints. The algorithm is 0(n ) where n is the 2 total number of control points. One disadvantage of this approach is it introduces too many extra degrees of freedom, roughly four times the number of original DOFs, making further manipulation of the surface not easys. Moreton and Sequin [40] developed a simple-to-use mechanism for the creation of complex smoothly shaped surface of any genus or topological type, by specifying interpolated geometric constraints consisting of positions, and optionally, surface normals and surface curvature. Non-linear optimization techniques are used to minimize fairness functionals based on the variation of curvature, to produce minimum variation curves (MVC) and minimum variation surfaces (MVS). Their methods produces very high quality surfaces with predictable, intuitive Chapter 3: Related Work 16 behaviour, while generating, where possible, simple shapes, such as cylinders, sphere and tori. The only drawback is the computationally expensive non-linear optimization that is needed. Giinther Greiner [29] discussed the choice of a fairness functional to tradeoff between high quality and computational effort. In particular, the comparison between the simplified thin plate energy functionals for parametric surfaces and M V C (Minimum Variation Curves) functionals [40] is presented. He pointed out that the simplified thin plate functional is only good in the case of bivariated function, but not in the parametric case, while the M V C functional produces very good result even in the parametric case. Since the M V C functional is computationally expensive, this suggests the use of data dependent fairness functionals instead of a universal one, and introduced two fairing functionals which are a good approximation to the minimal energy thin plate functional and the M V C functional respectively. These data dependent functionals are quadratic and positive definite, thus efficient to solve. However, the data dependence requires some knowledge of the geometry of the final surface like surface topology. 3.1.2 U S I N G P A R T I C L E S Y S T E M S Szeliski and Tonnesen [50] proposed a very different technique to model elastic surfaces based on interacting particle systems, which unlike previous techniques, can be used to split, join, or extend surfaces without the need for manual' intervention. The particles have long-range attraction forces and short-range repulsion forces and follow Newtonian dynamics. To model surface elements instead of point masses or volume elements, an orientation is added to each particle's state. They devise new interaction potentials for oriented particles which favour locally planar or spherical arrangements. Surface stretching and growing is modelled with automatically generation of new particles. It can be used for interpolation of 3D point set of unknown topology. The resulting surfaces is obtained through integrating system of differential equations through time, and the execution time is O(ralogn) with hierarchical data organization. One disadvantage is the lack of precise control over the mathematical form of the surface, needed in engineering applications. Chapter 3: Related Work 3 . 1 . 3 U S I N G M E S H 17 O P T I M I Z A T I O N Hoppe et al. [34] tackled the problem of surface reconstruction of a set of scattered 3D data ; points and mesh simplification by a technique called mesh optimization. It involves minimizing an energy function that captures the competing desire of tight geometric fit and compact representation over the mesh. The energy functional consists of three terms. First, a distance energy that measures the closeness of fit. Second, a representation energy that penalizes meshes with a large number of vertices. Finally, a regularizatipn term that conceptually places springs of rest-length zero on the edges of the mesh. This results in a smooth mesh whose edges align themselves along directions of low curvature, and whose vertices concentrates in areas of high Gaussian curvature. It also recovers sharp edges and corners. The minimization algorithm partitions the problem into two nested subproblems: an inner continuous minimization and an outer discrete minimization. The search space consists of all meshes homeomorphic to the given initial mesh of arbitrary topology. A method for the topological type determination of the initial mesh is presented in an earlier paper by Hoppe et al. [33]. 3 . 2 M U L T I R E S O L U T I O N A P P R O X I M A T I O N S Multiresolution representations of curves or surfaces can reduce the required storage space compared to traditional unit scale representations. Multiresolution representations also allow the hiding of unnecessary details when viewing or manipulating at a particular scale. They allow operations to be performed on different scales like multi-scale editing, which can make interactive manipulation of the resulting surfaces or curves easy. 3 . 2 . 1 M U L T I L E V E L V I S U A L S U R F A C E S W I T H D I S C O N T I N U I T I E S Terzopoulos [51] presented a computational theory of visual surface reconstruction that extended the multilevel structure of the earliest processing stages in vision to later stages that reconstruct full surface representations from local information about surface shape at scat- Chapter 3: Related Work 18 tered locations. The approach is to formulate the problem as a variational principle describing the equilibrium of a thin flexible plate. It uses an efficient multilevel algorithm, the multigrid method, to solve the hierarchy of full visual surface representation. In a related paper Terzopoulos [53] presented techniques for handling discontinuities on the visual surfaces. He views the problem of detecting surface discontinuities as one of distributed parameter iden- tification [52] [14] [43] within a variational formulation of visibleTSurface reconstruction. He proposed two estimation procedures which dynamically adjust the controlled-continuity model during surface reconstruction such that its continuity becomes consistent with discontinuities implied by data. The first approach is the use of local validation, which detects discontinuities by locally monitoring sharp deflections in the evolving surface. The second approach is the use of non-convex optimization [15] [46] for variational continuity control. 3 . 2 . 2 W A V E L E T S Finkelstein and Salesin [21] proposed a multiresolution curve representation, based on wavelets [18] that allows a curve to be edited at any scale while preserving its details, to be smoothed at continuous levels, and to be approximated to within, any given error tolerance for scan conversion. The storage requirement• is 0 ( n ) , where n is the original number of data points. If required, compression can be done by thresholding the resulted wavelet coefficients. The wavelets basis used are the endpoint-interpolating cubic B-splines. The reconstruction filters have a banded structure, allowing construction'in O(n) time. However, the analysis filters are dense, which would seem to imply an 0(n ) 2 time decomposition algorithm. Fortunately, Quak and Weyrich [45] have developed a clever trick for performing the decomposition in linear time. The result is an efficient algorithm both in terms'of storage and speed to represent curves by wavelets. The resulted wavelet representations allow changes in the overall "sweep" of a curve while maintaining its fine details, or "character", and vice versa. DeRose, Lounsbery, and Warren [19] extended multiresolution analysis to surfaces of arbitrary genus or topological types using techniques from subdivision surfaces. Currently, no one knows Chapter 3: Related Work 19 of a construction for subdivision surfaces leading to locally support pre-wavelets [18] (that span the orthogonal complement but are not mutually orthogonal), nor if such a construction exists. The approach they took is to obtain locally supported quasi-wavelets [18] which span a nonorthogonal complement. It is possible to make the quasi-wavelets span as close as needed to the orthogonal complement, at the expense of increasing their support. These locally supported quasi-wavelets imply a linear speed decomposition algorithm. However, truncation of a wavelet expansion leads to parametric least-squares best approximation; truncation of a quasi-wavelet expansion does not. The example application they gave is multiresolution rendering, where details are rendered only when visible. 3 . 2 . 3 H I E R A R C H I C A L B - S P L I N E S Forsey and Bartels [24] presented a way to decompose a surface fitting problem into a sequence of curve fitting processes, making the computation particularly efficient, when tensor product splines are .used with grid data. They employed the Kronecker product algebra [2] to derive this curve fitting processes. The use of the hierarchical bicubic B-spline representation [23] adds further efficiency by adaptively decomposing the fitting process into subproblems involving only a portion of the data, and focusing on out-of-tolerance residual points. The curve fitting processes are done by least squares. Linear constraints can also be imposed without affecting the efficiency, as long as the linear constraints retain the Kronecker structure. Comparison with multiresolution analysis is also made. Since only a least squares fit applied, the resulting surfaces could have unwanted wrinkles. We will address this particular issue in Section 4.3. Forsey and Wang [25] considered the problem of approximating a 3D digitized surface with a hierarchical bicubic B-spline [23] to produce a surface that can be further modelled or animated. A modified full multi-grid (FMG) technique is employed to obtain a high resolution B-spline approximation. The intermediate results of the F M G calculations generate the component overlays of a hierarchical B-spline surface representation. The fitting and the variation of chord length parameterization used are both adaptive, and concentrate the surface into those Chapter 3: Related Work 20 regions with a high gradient. The resulting fitting process may be more effective. However, the adaptive parameterization yields a non-intuitive parameterization that is hard to manipulate by hand. Since the usefulness of a parameterization depends on how it will be used, automatic generation would unlikely to fit all the needs of different user. 3 . 3 O U T L I E R S H A N D L I N G Some effective ways of how to handle outliers in statistics are presented below. However, they are generally inefficient as they are trying to solve a non-convex minimization problem. 3 . 3 . 1 T H E O U T L I E R P R O C E S S Geiger and Pereira [26] extended the weak membrane model by introducing an additional binary process, the outlier process, interacting with the data and capable of suppressing the noise data points for fitting. This outlier process is similar to the line process [53] [8] which detects surface discontinuities. Both the line process and the outlier process are used in their work. These two processes detect changes from the smoothness assumption and compete against each other. The outlier process detects outliers while the line process detects actual discontinuities. The resulting non-convex system is solved by applying deterministic continuous methods, such as deterministic annealing [9] and graduated non-convexity [8]. These methods first start by exactly solving for a specific, set of parameters, then continuously change the specific set of parameters and at the same time solve the set of equations for the updated parameters. Of course, as pointed out earlier, these methods are generally slow to converge. Black and Rangarajan [7] generalized the concept of a line process to that of an analog outlier process, and showed that a problem formulated in terms of outlier process can be viewed in terms of robust statistics [32]. They also characterized a class of robust statistical problems that can be converted into an outlier process formulation. Thus, they unified the line process, or outlier process approaches for regularization with discontinuities, and robust statistical Chapter 3: Related Work 21 techniques. The power of this result lies in the ability to recover an outlier process from a robust estimation problem, thus permitting the straightforward extension of results in robust statistics to problems in surface reconstruction. 3 . 3 . 2 R O B U S T S T A T I S T I C S Hampel etal. [32] identified the main goals of robust statistics. The first is the description of the structure best fitting the bulk of the data. Second, to identify deviating data points (outliers) or deviating substructures for further treatment, if desired. Least squares estimate is notoriously sensitive to outliers; the problem being that the outliers unduly influence the quadratic estimator. They proposed the use of influence functions to analyze the behaviour of an estimator. The influence function characterizes the bias that a particular measurement has on the solution and is proportional to the derivative of the estimator. To increase robustness, an estimator must be more tolerant of outlying measurements. A n estimator where the influence of outliers that tends to zero is called a redescending estimator. The Lorentzian estimator [32] is one example. 3 . 4 O P T I M A L C H O I C E O F T H E R E G U L A R I Z A T I O N P A R A M E T E R S The optimal choice of the regularization parameter A for the regularized surface approximation problem can be found using a technique • called generalized cross-validation [49] (see below). However, this method usually results in a search in a space with multiple local extrema, and requires repeated evaluations of the computationally expensive generalized cross-validation function. 3 . 4 . 1 G E N E R A L I Z E D C R O S S - V A L I D A T I O N Shahraray and Anderson [48] discussed the use of generalized cross-validation [49] for estimating the regularization parameter for smooth contour curves. The basic idea behind cross- Chapter 3: Related Work 22 validation is independent of whether we are handling a curve, a surface or even an other regularization problem. Predictive sample reuse [27] is another name for ordinary cross-validation. In the predictive sample reuse, the "goodness" of the parameter (or set of parameters) is measured in terms of the ability of the model to predict some of the observations based on the other observations. The optimum value of the parameter (or parameters) for m samples is found by minimizing an average discrepancy between the observed values of n data samples and the predicted values for these samples obtained by applying the predictive function to the remaining m — n samples. This average discrepancy is measured by the cross-validation mean square error, thus minimization is done over this error. Wahba and Wold [58] have applied this method to the determination of the smoothing parameter for the smoothing splines. Generalized cross-validation is a modified version of ordinary cross-validation that has been shown both theoretically and experimentally to have very desirable properties [17]. Performing generalized cross-validation is equivalent to first finding a rotation in the Euclidean n-space that transforms the problem into the symmetric case, and then applying ordinary cross-validation [57]. A minimum seeking algorithm is used, which iteratively reduces the size of the interval containing the minimum of the function by repeatedly evaluates the generalized cross-validation function. The generalized cross-validation function is related to the ordinary cross-validation mean square error. 3 . 4 . 2 S E A R C H M E T H O D S Shahrary and Anderson [47] presented an derivative based search method for finding the minimum of the generalized cross-validation function for the smooth contour curve. The method consists of a linear search in log (A) scale to,find two points on the opposite sides of the 10 minimum, followed by a binary search to converge to the minimum. Notice that the search space may have many local minima, in addition to the global minimum. In that case, this kind of gradient descent method may become trapped' in local minimum. Bates et al. [3] presented a set of routines for generalized cross-validation. A nonderivative Chapter 3: Related Work 23 based method is presented consisting of two stages. The first stage, performs a grid search-in log (A) scale by evaluating the function at a number of equally spaced grid points to find the 10 approximation of the minimum. The second stage, performs a golden section search [1] [36] in the neighbourhood of the minimizing grid point. The golden section search is preferable to the binary search, since it only requires one function evaluation for each iteration, compared to two for a binary search, to approximate the derivative of the function needed in the second stage. This method is less susceptible to finding a local minimum than a simple derivative based method. Other searching methods can be found in most standard texts on optimization theory [1] [42] [15] [46]. CHAPTER FOUR MULTI-RESOLUTION B-SPLINE APPROXIMATION The various stages of the construction method for hierarchical B-spline surfaces is described in this chapter, and tested with the spike example shown in Figure 4.1. The spike is a hard problem for hierarchical surface approximation. The single high frequency outlier causes most surface approximation algorithms to produce surfaces with undesirable oscillations. Figure 4.1: The spike (19 x 19 control points). 4 . 1 D I V I D E A N D C O N Q U E R The problem of generating a hierarchical B-spline representation from a given B-spline control mesh is decomposed into a smaller set of sub-problems, each of which dealing with just a single level of the hierarchy. The sub-problem is, given the current control mesh M ' of the hierarchy, find the next coarser level control mesh such that it: • is a close approximation in some norms to the surface of current level. • is smooth according to some smoothness measure. 24 f c + 1 l at level k + 1 Chapter 4'- Multi-resolution B-spline Approximation 25 • minimizes the number of non-zero offsets needed to store, in order to recover the current level from the coarser level representation. where the finest level will be just the control mesh M[/ ,nes 'l of the original B-spline surface. We note that we are operating on the control" mesh, not on any data points that may have been used to create 4.2 INITIAL M.M \ tnest APPROXIMATION Given a control mesh, M' fc+1 l, with control points^ V^- \ denning a B-spline surface, an initial k+1 approximation to the next coarser mesh, M^, with control points, V^, is to minimize the magnitude of the offsets at level k + 1. Minimizing the offsets implies choosing a norm for minimization. The Xi-norm may be good for approximation, since it is insensitive to outliers. The search space for using i j - n o r m is non-differentiable, thus no derivative methods can be used to solve it. We are left with very inefficient direct search methods [42]. The Lorentzian estimator [32] which is differentiable can also be used, as shown in related work Section 3.3. Unfortunately the problem is still non-convex and can only be solved by inefficient deterministic continuous methods, like deterministic annealing [9] and graduated non-convexity [8]. Therefore, we choose to use the simple but efficient Z2-norm for minimizing the magnitude of the offsets R[*+I] _ Since R ^ 4 " 1 ! = SV^, (4.1) v^ ! 1 we can rewrite the equation as 2 (4.2) which leads to solving the linear system (4.3) Chapter 4' Multi-resolution B-spline Approximation 26 in least squares. No smoothness constraint imposed yet. The resulting multi-resolution representation of this least squares approximation 1 or wavelet decomposition produces filtered surfaces can undulate. For example, a coarser level approximation of the spike example is shown in Figure 4.2. A Figure 4.2: Coarser level approximation by least squares. 4.3 S M O O T H N E S S C O N S T R A I N T S One way to avoid unwanted undulation in the solution of Equation 4.3, is to impose a penalty on the smoothness or variation of curvature on the resulting surface, similar to the use of the fairness energy norm for free form deformable surface shown in Section 3.1.1. The underlying assumption is that most objects are smooth over most of the surface except at some areas. Therefore we can impose a smoothness constraint upon the coarser approximation to reduce the undulation. Various fairness energy norm or smoothness energy models are available. In our case we want to smooth out the control mesh points thus, we treat the mesh points not the parametric surface H ^ l , as mesh of sampled points of an underlying "control sur- face". Notice that the smoothness constraint is simultaneously imposed with the interpolation constraint, not after the approximation mesh is generated. This is because the smoothness S o l v i n g this least squares approximation can be sped up tremendously by employing the tensor product property of B-splines to decompose the problem into a sequence of curve fitting processes [24]. Chapter 4- Multi-resolution B-spline Approximation 27 constraint is used as a regularization. 4 . 4 T H E M E M B R A N E M O D E L One possible smoothness energy model is the membrane model [13], that models the energy of a membrane undergoing stretching. The energy involved is computed from the first order partial derivatives of the surface, JJ A(||vi || + ||vi || yWi; fc] 2 fc] (4.4) 2 where subscripts u and v are the parametric partial derivatives, and A is a regularization parameter that controls how stiff the membrane is., 4 . 4 . 1 T H E M E M B R A N E M O D E L F O R N O I S E R E M O V A L The membrane model is quite commonly used for noise removal over a height field surface. The process is to find a smooth surface u(x,y) which approximates the data d(x,y) by minimizing F(u) = JJ(u - df + \(u 2 + u )dxdy (4.5) 2 x y where (u — d) measures the deviation of the approximation surface from the data. 2 By the calculus of variation, the unique solution to this minimization problem is by solving the functional ^}^ — 0. A n extremely efficient multi-grid method [30] can be applied to find d the solution (see Appendix A ) . 4 . 4 . 2 T H E E F F E C T O F V A R Y I N G A O N T H E M E M B R A N E M O D E L Given the smooth surface in Figure 4.3, we added random noise to it with mean 0 to generate the noisy test data in Figure 4.4a. Different values of A yield different amounts of smoothing. Chapter 4'- Multi-resolution B-spline Approximation 28 Figure 4.3: The original smooth surface data Figure 4.4 presents the resulting surfaces with respect to different values of A . In Figure 4.4d, although all the noise has been/removed, the computed surface deviates quite a bit from the original data in Figure 4.3. Moreover, the membrane model has a tendency of making the computed surface flat (ie. constant height), because of penalizes large first order partial derivatives. This tendency is not suitable for smoothing out intermediate levels within the surface hierarchy, because most 3D surfaces can have large first order partial derivatives with respect to the parameter space for each coordinate, that is they can be slanted with respect to the uv parameterization. Chapter 4- Multi-resolution B-spline Approximation (e) A = 0.1 29 (f) A = 1.0 Figure 4.4: The effect of varying A on the membrane model. Chapter 4- Multi-resolution T H E T H I N 4.5 B-spline P L A T E 30 Approximation M O D E L The thin plate [29] model models the bending energy of a very thin plate of material. The energy is calculated by (4-6) where a n d ?y[k] K denote the principle curvatures of the surface V ^ . du> denotes the surface element, and Q denotes the entire surface domain. This functional is highly nonlinear and the numerical minimization of the function is very computationally expensive, thus often one sees the simplified thin plate energy [12][13] [54], that minimizes IIq vt*i +2 + (4.7) )dudv where A is the regularization parameter that controls the stiffness of the plate, and the subscripts uu, uv, and vv are the second order partial derivatives. This simplification may not be a very accurate approximation in the parametric discussed by Giinthef Greiner [29] (see Section 3.1.1). However, a better approximation called M V C functional proposed by Giinther Greiner [29] is still computationally expensive. Thus, in terms of tradeoff between high quality and computational effort, this simplified thin plate energy (Equation 4.7) is still worth using. Discretizing the energy measure for the control mesh yields t'=0 j=0 + 2 (vg) (4.8) + where (y\j)uu = second order derivative computed by central differencing (v\j)uv — second order uv-derivative computed by forward differencing 1 -2 1 (4-9) -1. 1 1 -1 14.10) Chapter 4'- Multi-resolution B-spline Approximation 31 1 (VJj) „ = second order uv-derivative computed by central differencing v -2 (4.11) 1 The above three matrices for the second order derivatives are called templates, and they are collected into the smoothness energy matrix P , which can be de-coupled into a linear system for each of the x, y, and z coordinates. Although the actual second order partial derivatives of the underlying surface that the mesh defined can be used to measure the smoothness, we are more interested in fairing out the control mesh. 4.5.1 T H E T H I N P L A T E M O D E L F O R N O I S E R E M O V A L Like the membrane model, it is quite common to use the thin plate model for noise removal over a height field surface. This use of the thin plate model demonstrates some of the properties of the functional. The computational requirement for the thin plate is more complex than the membrane model. However, it produces smoother approximations of noisy surfaces with smaller mean squares error. The formulation finds a smooth surface u(x,y) which approximate the data d(x,y) by minimizing F(u) = Jj{u - df + \(u 2 xx + 2u 2 xy + u )dxdy 2 yy (4.12) with the necessary condition of minimum by calculus of variation is \A u + u = d 2 where the A u — u 2 xxxx + 2u xxyy (4.13) + u yy is the biharmonic operator in 2D. There are many yy references on how to solve biharmonic equations (see [11] [41] [10] and references therein). One discretization of Equation 4.12 and solving method is given by Terzopoulos, for computing visible-surface representation in the field of robotics vision [51] [53]. We have implemented a. Chapter 4- Multi-resolution B-spline Approximation 32 different efficient multi-grid method for solving this noise removal problem (see Appendix A). 4.5.2 T H E E F F E C T O F V A R Y I N G A O N T H E T H I N P L A T E M O D E L Like the membrane model, different A values yield different amounts of smoothing (see Figure 4.5). The original smooth surface without noise is in Figure 4.3. From Figure 4.5 we notice that the thin plate model is very effective at removing noise. In the membrane model, a large enough A to remove all noise produces a flat surface that does not represent the data anymore. Moreover, in the extreme case when A = 1.0, the thin plate model computes a slanted linear plane, whereas the membrane model produces a horizontal plane. This is because we are minimizing a higher order derivative, at the cost of dealing with a slightly more complicated system of equations and a slightly higher computational cost. We will use the thin plate model for the smoothness constraint. Chapter 4- Multi-resolution B-spline Approximation (e) A = 0.01 33 (f)A = Figure 4.5: The effect of varying A on the thin plate model. 1.0 Chapter 4- Multi-resolution 4.5.3 S P A R S E D A T A B-spline S E T F O R Approximation S M O O T H 34 A P P R O X I M A T I O N Both the membrane and the thin plate model work reasonably well at computing a smooth approximation of a sparse data set. Figure 4.6 c, shows using only 30% of the randomly distributed noise data available, the thin plate model for smooth approximation can still recover the original surface quite well. This is true for the membrane model, except the membrane model produces a less smooth surface. Chapter 4- Multi-resolution B-spline Approximation (e) 30% data distribution (f) Approximate surface for 30% data Figure 4.6: The thin plate model on sparse data set. 35 Chapter 4- Multi-resolution B-spline Approximation 4.6 L E A S T S Q U A R E S A P P R O X I M A T I O N W I T H 36 S M O O T H N E S S C O N S T R A I N T Assuming that V^" " ] will be given, we want to find V"W that minimizes the number of offsets 1 1 to represent the hierarchical B-spline surface and also be smooth according to the thin plate functional. That is, we want to replace the data d(x,y) in the noise removal problem with the given mesh points vf mesh points R i f c + 1 fc+1 J and replace the approximation surface points u(x, y) with the refined l =SV^. A n efficient multi-grid method may be desirable for this problem. However, due to the way B-spline refinement works, the coarser boundary control points may be throwing away at finer levels. The smoothness constraints and the interpolation constraints are on two adjacent levels. These two factors make the derivation of a multi-grid method non-trivial needing many special cases. We may want other means to solve this problem. One way is to consider it as a least square problem sv w - v[ mm fc+1 i + A p'vM (4.14) where S is the subdivision m a t r i x o f the B-spline, and P is the thin plate smoothness energy matrix. This can be written as an over-determined least square problem S \ (4.15) However, there are problems with just-doing a planar least squares. It Will attempt to fit the outlier most, since they will have much higher residuals due to the squaring. Thus, we will have a fit that is a closer approximation to the outliers at the expense of a non-exact fit and some small amplitude wrinkles on the remaining control points. The magnitude of wrinkles is small due to smoothing. For the spike example, we will have the coarser level surface as shown in Figure 4.7. This problem is due to the fact that we are trying to fit those control points that are supposed to be in V^ ^ but not in R[* 1 (the refined control mesh from V W ) . To get better result, k+1 +1 Chapter 4' Multi-resolution B-spline Approximation Figure 4.7: Coarser level approximation, with smoothing A = 0.1. we need a way to estimate which control points should be fitted and which should not. This requires using some sort of weighting on the approximation. W E I G H T E D 4.7 A P P R O X I M A T I O N Although there does not seem to exist a way of finding out which control points belong to Vt f c + 1 l and not R j * ] , and which belong to R [ + 1 f c + 1 l , we can estimate these geometric outliers by computing the second order partial derivatives in the control mesh V ^ ) with respect 1 to the uv parameter space. We then sum the resulting magnitudes together, and made the following observation. If a control point has a very high magnitude of second order partial derivations (which give crude approximation of curvature), it is very likely that this point will not be in R ^ ' + 1 l We can use actual curvature for this outlier estimation. However, the computational cost for second order partial derivatives are lower than that for curvature. The basic approach is to first eliminate the areas of high curvature by assigning weights to the control vertices, according to the above second order partial estimates. By subtracting all the average partials from the maximum partial, and then normalizing them by the maximum, we get weights for each control point in the range [0,1]. The computed weights will control to what degree the thin plate will approximate those control vertices. A 0 weight means don't fit through the control point at all, while weight 1 means try to fit to the control point as much Chapter 4'- Multi-resolution B-spline Approximation. 38 as we can without violating the smoothness constraint (see Figure 4.8 for an example). As we noted earlier in Section 4.5.3, the thin.plate model is quite robust in terms of incomplete data set, therefore this weighted system will also be robust. The approximation constraints are weighted by the assigned weights. We denote the weight matrix by w t fc+1 l . We then solve this'over-determined system using least-squares. / W [ * + i ] S \ / W [fc+l]y[fc+l] (4.16) V v/AP j 0 V Which corresponds to the following minimization, problem mm (4.17) Solving these equations for the spike example, produces a flat coarser level surface as shown in Figure 4.9. Chapter 4- Multi-resolution B-spline Approximation Figure 4.9: 4.8 C A L C U L A T I N G 39 Coarser level approximation, using weighted fit. A N D T H R E S H O L D I N G T H E O F F S E T S By solving, V ^ , in Equation 4.16, we get the control mesh M ^ , for the coarser level. Iteratively solving for until we reach k = 0, where restriction can no longer applied. This gives us the multi-resolution representation of the original surface. In the hierarchical form, the surface shape is stored as offsets (see Equation 2.3). Thus, we calculate the offsets for each level k + 1, starting from level 1, where level 0 is the coarsest level. Q[l] = V [i] _ y[o] s ( 4 - 1 8 ) To reduce storage space, we can threshold T(.) the offsets as we compute each level, only keeping those with large magnitudes. k > 0. The modified is denoted by v^. [k] v where is V ^ . Note that, by so doing, we will change the = They are calculated as [k-i] Sv for + .r(QW) (4.19) Chapter 4- Multi-resolution B-spline Approximation 40 Thus, the offsets for the remaining k > 1 levels are computed by Q[*+I] = y[k+i] _ .§„[*] (4.20) C H A P T E R FIVE IMPLEMENTATION This chapter presents the basic algorithm and the way to set up the appropriate matrices within the system of equations for the construction of hierarchical B-spline surfaces. Various solving methods for the system of equations are presented in the Appendix B and C. 5.1 ALGORITHM The basic algorithm for constructing of hierarchical B-spline surfaces is 1. Input the finest level control mesh M [ / ylfmest] a n ( j v a m e f or t n e , n e s l l of a B-spline surface with control points smoothing parameter A. 2. Repeatedly compute the coarser level k control mesh until the resolution of the mesh can no longer be reduced. (a) Set up the subdivision matrix S for the coarse level, according to the control mesh size and topology (see Section 2.2.2 and 5.2.1). (b) Set up the thin plate energy measure matrix P (see Section 4.5 and 5.2.2). (c) Assign weights W ^ * ' to the approximation (see Section 4.7 and 5.2.3). 1 (d) Combine the matrices, control mesh points and the regularization parameter A to form the Equation 4.16. 1 , wi* is ^ +1 vAP , / wt ]vf J \ fc+1 \ fc+1 0 (e) Solve this overdetermined system using least squares. 41 / Chapter 5: Implementation 42 3. Given the above hierarchy, compute the offsets for each level from coarsest to finest, with thresholding applied to the magnitude of the offsets fitting (see Section 4.8). 4. (Optional) To force exact fit, keep all the offsets for the finest level without thresholding. 5 . 2 S E T T I N G U P T H E M A T R I C E S We have to set up the matrices, S , P , and W ' ' for each level k during the construction of F C hierarchical B-spline. We will use a lexicographical order for the indexing in YJ F C J in their vector form. 5 . 2 . 1 T H E S M A T R I X For a B-spline subdivision surface, the subdivision template is 16 1 1 4 4 6 6 4 4 1 1 (5.1) To use this template, the spacing of the uw-space grid corresponding to the control mesh, M ^ , is divided by half. The subdivision template is then used to calculate the new control points, Vr^ " ', 1 1 [k] as a weighted sum of the original control points V)"- that are underneath a weight of the template. The template is centered at the control point, v[*s^, the value is to be computed. Note that this template handles the boundaries and corners as well as the interior nodes. The S matrix is constructed by pulling out the necessary weights from the template for each v{*s^ and by inserting those weights into their corresponding column according to the nodes Chapter 5: Implementation 43 v j j they apply to, and in the row corresponding to v{.^ ^ of the S matrix. Figure 5.1 shows +1 some cases. • \ c^ 1 J \ ! J 1 1 1 Case 1 Case 2 Case 3 Figure 5.1: Three cases of the use of subdivision template. The corresponding rows in the S matrix, where row numbers (from top to bottom) correspond to case number. r ••• v't - i , j. - i y M fc] v 11 16 4 V 1 3 16 2 1 1 16 4 1 16 1 16 16^ 16^ J_3 16 2 J_q 16 3 «,J+1 • • • 1 3 16 2 i+lj-l 11 16 4 &6 16^ v v t+l,j vt fc] • 13 16 2 «+l,j+l 1 1 16 4 1 16 1 16 v 1 • • •. (5.2) 16^ This yields Equation 2.2: the refinement matrix for B spline surfaces. For other subdivision : surfaces, we just need to change the subdivision template, and for different topology, like cylindrical, the template is wrapped across boundaries. 5 . 2 . 2 T H EP M A T R I X According to the Equation 4.8, the central and forward differencing templates for second order partials, are used to compute the uu, uv, and vv partial derivatives of V ^ L Each of these partials for a particular control point V , - * - is set up in a row of the matrix P in a manner similar to the way we did for the S. matrix. The square of the norm in Equation 4.8 is included in the least squares that will be applied to the system. Thus, the three rows of the P matrix Chapter 5: 44 Implementation for a particular control point V*|*- is v » - i , j v i - i , j + i ••• » , j - i v v 1 i,j v t',j -2 V2 v/2 1 -v^ (5.3) 1 5.2.3 THE wt l fc+1 MATRIX The way to assign weights to the approximation is described in Section 4.7. The weight matrix Wf fc+1 ] is a diagonal matrix, with weights Wij on the diagonal corresponding to each of the control point V?** ^. 1 u>o,o (5-4) W ,2 0 w. 5.3 SOLVING T H E SYSTEM Comparisons of various least squares solving methods for the system (Equation 4.16) are presented in the Appendices. The method chosen is the Corrected Normal Equation (see [38]) without any row ordering of the matrix. This is the most accurate and time efficient method among all those tested. The running time is 0(n ), where n is the total number of control 2 points in the coarser level mesh. Typical execution time for the generation of the whole hierarchy is under 10 seconds for a 67 X 67 mesh and less than 3 minutes for 128 X 195 on a MIPS RS4000 processor. The space requirement is also 0{n ). If memory is a concern, the 2 Conjugate Gradient iterative method (see Appendix B and C) which is linear in storage, but it is slower than the Corrected Normal Equation, can be used. CHAPTER SIX RESULTS The method was tested on a variety of input meshes. Most were created using Dragon, an interactive editor for hierarchical B-splines. Other test surfaces included a B-spline surface created using Softimage (a commercial- modelling and animation package), and a digitized surface data created using a Cyberware laser scanning system. We would like to thank Reboot for the B-spline surface, and to thank Tony D. DeRose, Micheal Lounsbery and Joe Warren for the Cyberware data. The results are summarized in Table 6.1,Table 6.2 and Table 6.3. A l l of the resulting surfaces approximate the initial input mesh within 0.5% tolerance of its overall range. The final column in the table lists the number of control points required to represent the surface as a traditional full tensor-product B-spline. 6.1 H I E R A R C H I C A L B - S P L I N E Surface S U R F A C E S Original No. of offsets No. of offsets obtained No. of patches 17 19 19 422 257 478 17 19 19 253 302 467 34 40 28 565 609 1383 Spike (Figure 6.1) Hump (Figure 6.2) Plop (Figure 6.3) Eye (Figure 6.4) Dog (Figure 6.5) Mask (Figure 6.6) B-spline equivalent ' 49 143 121 2625 1176 24960 Table 6.1: Results for surfaces created using an interactive hierarchical B-spline editor 6 . 1 . 1 S I M P L E G E O M E T R I E S Figures 6.1 through 6.3 demonstrate the hierarchical surface produced for some simple geometries. The intermediate levels very closely-resemble those produced with the interactive editing 45 Chapter 6: Results >,>> system, have no extra undulations, and add no extra offsets. In fact they are so close to the original hierarchies, that one can not tell from their rendered images. Therefore, the original hierarchies are not displayed. Note that the colour images in this section show the distribution of patches of different parametric size. Patches defined at the finest level of refinement are yellow, while patches at coarser levels are reddish. (c) All levels Figure 6.1: Spike (28 patches, 17 offsets) Chapter 6: Results (e) All levels Figure 6.2 Hump (40 patches, 19 offsets) Chapter 6: Results 48 (e) All levels Figure 6.3: Plop (28 patches, 19 offsets) Chapter 6: Results 6.1.2 COMPLEX 4$ GEOMETRIES The human eye in Figure 6.4 was created by a novice user of the interactive editor. The resulting approximation actually reduces the number of offsets used to represent the surface while retaining the final shape. For a later effort by the same individual modelling a dog (see Figure 6.5), the algorithm was unable to improve the representation. Figure 6.4: Eye (565 patches, 253 offsets) Chapter 6: Results r>0 (e) All levels Figure 6.5: Dog head (609 patches, 257 offsets) Chapter 6: Results 6.1.3 USE I N 51 A N I M A T I O N The mask in Figure 6.6 was also created by a novice user of the interactive editor, who wanted to build a mask model for animation. The model was not easy to animate, because part of the control mesh in the intermediate levels were self intersecting making the effect of broad scale deformations difficult to control. Using the finest mesh of this surface as the basic shape a hierarchical B-spline representation was reconstructed. The result was promising, all intermediate levels are smoothed out and the resulting model was used for further editing and in a 30 seconds animation. Chapter 6: Results 52 (a) Original (b) Generated Figure 6.6: Mask: Comparison of original and generated hierarchies Chapter 6 . 2 6: Results B - S P L I N E 53 S U R F A C E One goal of the algorithm is to widen the applicability of hierarchical B-splines by creating them from existing data. The head in Figure 6.7 was created using Softimage for Reboot, a weekly television series executed entirely using computer generated images. Because of the time constraints of the production schedule, this surface is about at the limit of the complexity practical to use. The surface is tedious to manipulate and there are too many patches in places where they are not needed due to the non-local knot insertion property of B-splines (ie., extra knots inserted along the top and back of the head). Although the hierarchical B-spline surface we obtained with a single tolerance did not reduce the storage space, the smooth hierarchy it generated is good for building complicated facial expressions (see Figure 6.7d). Surface Initial No. of control points No. of offsets obtained No. of patches B-spline equivalent 648 648 658 324 600 369 576 576 Reboot (Figure 6.7) Reboot w / var tol (Figure 6.8) Table 6.2: Results for B-spline surfaces 6 . 2 . 1 F I X E D T O L E R A N C E V E R S U S V A R I A B L E T O L E R A N C E The initial mesh of 27 X 24 control points of the B-spline surface produced the hierarchical surface in Figure 6.7. Using a single offset tolerance for the entire surface generate 658 offsets. By selectively decreasing the tolerance in those regions where detail is not necessary (such as the top and back of the head), the number of offsets can be reduced to 324, and the total number of patches required dropped from 658 to 369, because only the surface in the region of interest is refined (see Figure 6.8). Chapter 6: Results (b) Level 1 (c) Level 2 (d) Hilarious Figure 6.7: Reboot head: Generated from a B-spline surface with uniform tolerance. Chapter 6: Results Figure 6.8: Reboot head: Generated from B-spline with variable tolerance. 55 Chapter 6: Results 6.3 D I G I T I Z E D 56 D A T A Finally, a hierarchical surface was generated from digitized surface information. The original data was 170 X 256 points representing distance (see Figure 6.9i), evenly distributed around a vertical axis. These were converted to cartesian coordinates in 3 dimensions and used directly as the B-spline control points for the initial mesh. To produce a more accurate representation, we could have ran a B-spline approximation routine to the initial data to get an initial B-spline control mesh. However, the data samples are so dense, that the control mesh of a B-spline surface approximation with the same number of DOFs is close to the original data. The mesh was approximated using an offset tolerance of 0.5% producing a hierarchical surface (see Figure 6.9a-g) with 5 levels using a total of 6492 offsets. The surface decomposition is show in Table 6.4. Surface Initial No. of data points Spock (Figure 6.9) 43520 No. of offsets obtained 6492 No. of patches 19567 B-spline equivalent 40376 Table 6.3: Result for a digitized scan 3D data Level 0 1 2 3 4 5 No. of patches No. of offsets 40 160 640 2551 9178 19567 64 195 505 1250 2650 1828 Table 6.4: Distribution of offsets for the Spock surface With the surface in hierarchical form, control points at lower levels of detail in the hierarchy can be used to make broad scale changes in surface shape. Figure 6.9h shows the result of moving just 5 control points at level 2 in the hierarchy: three for the eyebrow and two for the corner of the mouth. Chapter 6: Results fd) Level 3 (g) All levels 57 (e) Level 4 (h) Smiling Spock (f) Level5 (i) Original data Figure 6.9: Spock head: Generated from Cyberware data, 43520 data points. CHAPTER SEVEN CONCLUSIONS We have described an efficient method for constructing a hierarchical B-spline surface from an initial mesh of control points. This method uses the thin plate model as a fairness constraint on the solution, making the intermediate level approximation surfaces fair. Various numerical techniques have been examined or employed to speed up the solution process. On input meshes derived from hierarchical B-spline surfaces, the algorithm generates intermediate levels that, in simple cases, are very close to the original and for more complex cases may improve upon the representation built interactively. The improvement is both in terms of reduction in the number of offsets required, and smoother intermediate control meshes that ease the task of animation. On input meshes derived from external sources, the results are also good. For the digitized 3D scan data, the storage compression ratio is high. Moreover, it produces a compact representation that is useful for subsequent modelling and animation. In the related work section, we already gave some pointers to relevant materials for future research. This will include to examine additional efficient methods to improve the identification of features on a control mesh, to introduce derivative discontinuities into the method to facility modelling of even more complicated objects, to determine the optimal regularization parameter value, and to develop better numerical solving methods for meshes of other topologies. 5 8 Bibliography P. R. Adby and M . A . H. Dempster. York, 1974. G . Alexander. 1981. Introduction to optimization theory. Wiley, New Kronecker Products and Matrix Calculus with Applications. Halsted Press, D. M . Bates, M . J . Lindstrom, G . Wahba, and B. S. Yandel. Gcvpack - routines for generalized cross-validation. In Communications in statistics: Simulation and computation, pages 263-297, 1987. A . Bjorck. Iterative refinement of linear least squares solutions I. In BIT, pages 257-278, 1967. A . Bjorck. Least squares methods. Linkoping University, 1967. Technical report, Department of Mathematics, A . Bjorck. SSOR preconditioning methods for sparse least squares problems. In Computer Science and Statistics 12th Annual Symposium on the Interface, pages 21-25, 1979. M . J . Black and A . Rangarajan. The outlier process: Unifying line processes and robust statistics. In Proceedings of IEEE Conference on Computer Vision and Pattern Recogni- tion, pages 15-22, 1994. A . Blake and A . Zisserman. sachusetts, 1987. T. Boult. Visual reconstruction. The M I T Press, Cambridge, Mas- What is regular in regularization. Conference on Computer Vison, London, In Proceedings of the pages 457-462, 1987. First International D. Braess and P. Peisker. On the numerical solution of the biharmonic equation and the role of squaring matrices for preconditioning. In IMA Journal on Numerical Analysis, pages 393-404, 1986. S. C. Brenner. A n optimal-order nonconforming multigrid method for the biharmonic equation. In SI AM Journal on Numerical Analysis, pages 1124-1138, 1989. ' M . Carignan, Y . Yang, N. Magnenat-Thalmann, and D. Thalmann. Dressing animated synthetic actors with complex deformable clothes. In Proceedings of SIGGRAPH '92, Computer Graphics, pages 99-104, July 1992. G . Celniker and D. Gossard. Deformable curve and surface finite-elements for free-form shape design. In Proceedings of SIGGRAPH '91, Computer Graphics, pages 257-266, July 1991. 59 Bibliography 60 [14] G . Chavent. Identification of distributed parameters. In IFAC Symposium on identification and system parameter estimation. New York: American Elsevier, 1973. [15] F. H. Clarke. Optimization and nonsmooth analysis. Wiley-Interscience, New York, 1983. [16] R. Clough and J . Tocher. Finite element stiffness matrices for analysis of plate bending. In Matrix Methods in Structural Mechanics, pages 515-545, 1966. [17] P. Craven and G . Wahba. Smoothing noisy data with spline functions - estimating the correct degree of smoothing by the method of generalized cross-validation. In Numerische Mathematik, pages 377-403, 1979. [18] I. Daubechies. Ten lectures on wavelets. S I A M . Philadelphia, 1992. [19] T. DeRose, M . Lounsbery, and J . Warren. Multiresolution analysis for surfaces of arbitrary topological type. Technical Report 93-10-05, University of Washington, 1993. [20] L. Elden. Numerical Linkoping, 1977. Analysis of Regularization and Constrained Least Squares Methods. [21] A . Finkelstein and D. H. Salesin. Multiresolution curves. In Proceedings of SIGGRAPH '94, Computer Graphics, July 1994. [22] D. R. Forsey. A surface model for skeleton-based character animation. In Proceedings of the Second Eurographics Workshop on Animation and Simulation, September 1991. [23] D. R. Forsey and R. H. Bartels. Hierarchical B-spline refinement. In Proceedings of SIGGRAPH '88, Computer Graphics, pages 205-212, August 1988. [24] D. R. Forsey and R. H. Bartels. Surface fitting with hierarchical splines. Technical report, Computer Graphics Laboratory, Computer Science Department, University of Waterloo, April 1993. [25] D. R. Forsey and L. Wang. Multi-resolution surface approximation for animation. Proceedings of Graphics Interface '93, pages 192-200, May 1993. [26] D. Geiger and R. A . M . Pereira. The outlier process. In In Neural networks for signal processing, pages 60-69, 1991. [27] S. Geisser. The predictive sample reuse method with applications. In Statistic Association, pages 320-328, 1975. [28] B. R. Gossick. Hamilton's and London, 1967. Journal American Principle and Physical Systems. Academic Press, New York [29] G . Greiner. Variational design and fairing of spline surfaces. In ics '94, pages C143-C154, 1994. [30] W . Hackbusch. Multi-grid Springer-Verlag, 1985. Proceedings of Eurograph- methods and applications. Computational mathematics. Bibliography 61 M . Halstead, M . Kass, and T. DeRose. Efficient, fair interpolation using catmull-clark surfaces. In Proceedings of SIGGRAPH '93, Computer Graphics, pages 35-44, August 1993. F. R. Hampel, E. M . Ronchetti, P. J . Rousseeuw, and W . A . Stahel. Robust statistics: The approach based on influence functions. John Wiley and Sons, New York, N Y , 1986. H. Hoppe, T. DeRose, T . Duchamp, M . Halstead, H. Jin, J . McDonald, J . Schweitzer, and W . Stuetzle. Piecewise smooth surface reconstruction. In Proceedings of SIGGRAPH '94, Computer Graphics, pages 295-302, 1994. H. Hoppe, T. DeRose, T. Duchamp, J . McDonald, and W . Stuetzle. Mesh optimization. In Proceedings of SIGGRAPH '93, Computer Graphics, pages 19-26, 1993. M . H. Karlsruhe and P. C. H. Lyngby. pages 253-315. Springer-Verlag, 1993. W. J . Kennedy and J . E. Gentle. Regularization Methods for Large-Scale Problems, Statistical computing. Dekker, New York, 1980. T. Lyche and M . Knut. Spline-wavelets of minimal support. Technical report, Institutt for Informatikk, Universitetei i Oslo, 1992. P. Matstoms. Sparse QR Factorization with Applications to Linear Least Squares Prob- lems. Linkoping, 1994. F. Mokhatarian and A . Mackworth. A theory of multiscale, curvature-based shape representation for planar curves. In IEEE Transactions of Pattern Analysis and Machine Intelligence, pages 789-805, 1992. H. P. Moreton and C. Sequin. Functional optimization for fair surface design. In Proceedings of SIGGRAPH '92, Computer Graphics, pages 167-176, 1992. P. Peisker and D. Braess. A conjugate gradient method and a multigrid algorithm for morley's finite element approximation of the biharmonic equation. In Numerische Mathematik, pages 567-586, 1987. D. A . Pierre. Optimization theory with applications. Dover, New York, 1986. M . P. Polis and R. E. Goodson. Parameter identification in distributed systems: A synthesizing overview. In IEEE Transactions on pattern analysis and machine intelligence, pages 45-61, 1976. H. Qin and D. Terzopoulos. D - N U R B S : Dynamic Non-Uniform Rational B-splines. In ACM Transactions on Graphics, pages 103-136, April 1994. E. Quak and N. Weyrich. Decomposition and reconstruction algorithms for spline wavelets on a bounded interval. C A T Report 294, Center for approximation theory, Texas A & M University, April 1993. Bibliography 62 [46] R. T. Rockafellar. The theory of subgradients and its application to problems of optimization: Convex and nonconvex functions. Helderman Verlag, Berlin, 1981. [47] B. Shahraray and D. J . Anderson. Optimal smoothing of digitized contours. In IEEE Computer Society Conference on Computer Vison and Pattern Recognition, pages 210218, June 1986. [48] B. Shahraray and D. J . Anderson. Optimal estimation of contour properties by crossvalidated regularization. In IEEE Transactions on pattern analysis and machine intelligence, pages 600-610, June 1989. [49] M . Stone. Cross-validatory choice and assessment of statistical prediction. Royal Statistic Society Series B, pages 111-147, 1974. In Journal [50] R. Szeliski and D. Tonnesen. Surface modeling with oriented particle systems. In Proceedings of SIGGRAPH '92, Computer Graphics, pages 185-194, July 1992. [51] D. Terzopoulos. Multilevel computational processes for visual surface reconstruction. In Computer Vison, Graphics, Image Processing, pages 52-96, 1983. [52] D. Terzopoulos. Regularization of inverse visual problems involving discontinuities. In IEEE Transactions on Pattern Analysis and Machine Intelligence, pages 413-424, 1986. [53] D. Terzopoulos. The computation of visible-surface representations. In IEEE Transactions on pattern analysis and machine intelligence, pages 416-438, July 1988. [54] D. Terzopoulos, J . Piatt, A . Barr, and K. Fleischer. Elastically deformable models. In Proceedings of SIGGRAPH '87, Computer Graphics, pages 205-214, July 1987. [55] A . N. Tikhonov. Regularization of incorrectly posed problems. In Soviet Mathematics Doklady, pages 1624-1627, 1963. [56] A . N. Tikhonov. Solution of incorrectly formulated problems and the method. In Soviet Mathematics Doklady, pages 1035-1038, 1963. regularization [57] G . Wahba. A survey of some smoothing problems and the method of generalized crossvalidation for solving them. In Applications of Statistics, pages 507-523, 1977. [58] G . Wahba and S. Wold. Periodic splines for spectral density estimation: The use of crossvalidation for determining the degree of smoothing. In Communications in statistics, pages 125-141, 1975. APPENDIX THE A MULTI-GRID A.l METHOD A BRIEF I N T R O D U C T I O N TO MULTI-GRID M E T H O D A multi-grid method [30] is a numerical solving technique that utilizes different scales of grids to achieve very efficient 0{n) solving time for a given error decrease ratio. The notations for multi-grid method that solve the equation L finest^ on the finest grid for Ufi n e s t , given Lji nest finest and — /finest /finest are defined in the following. A sequence of grids HQ, ill, • • •, Hh • • ••> Hfinest . . . ,h\, . . ., h fi corresponding grid spacing ho, hi, We assume hj = 2hj i, ilj C Hj+i for all j C { 0 , 1 , . . + n e s t .,/inest}. The operators involved are • Relax(l, ui, Ll, fl) is the iterative solver (also called a smoother) that relax on level /, the equation L\u\ = /;. • is the restriction operator that transfers data from a fine grid at level / to a coarser grid at level / — 1. 63 Appendix A • 64 is the prolongation operator that transfers data from a coarse grid at level / — 1 to a finer grid at level /. • L\ is the coarse grid operator where Lfinest is the original system to be solved. The parameters are • v\ is the number of relaxations or smoothing before coarse grid correction. • V2 is the number of relaxations after coarse grid correction. • 7 is the number of coarse grid corrections. If 7 = 1, it is called a V-cycle. If 7 = 2 , it is called a W-cycle (see Figure A . l ) . W-cycle Figure A . l : One multi-grid iteration for V-cycle and W-cycle on a 5 levels grid Appendix A T H E A.1.1 65 M U L T I - G R I D M E T H O D Procedure mgm(l,ui, //) if / = 0 then solve exactly LQUO = fo else Relax(l,ui, Li, fi) for v\ times di-i := lj~ (fi — Liii[) (restrict the computed defect) 1 vi-i := 0 (initialize) for j:=l to 7 mgm(l — 1, vi-i, d/_i) (coarse grid correction) ui := ui + l\_ vi-x x Relax(l,ui, Li, fi) for v>2 times A.2 A M U L T I - G R I D S O L V E R F O R N O I S E R E M O V A L U S I N G T H E M E M B R A N E M O D E L For noise removal we are solving for u that satisfies the functional }^ = 0. In the multi-grid dF language, the fine and coarser grid operator, X/, for all / in valid range, is a discretization of the functional } ^ at different level of grid. And f f i d F n e s t = 0 in this case. Red-black Gauss-Seidel iteration can be used for relaxation. For the grid transfer operations, we can use full weighting with reflection boundary for restriction, iff, and bilinear interpolation for prolongation, Ijj. A v-cycle will do the job. A multi-grid solver has been implemented to test out the behaviour of the membrane model. 66 Appendix A A.3 A M U L T I - G R I D SOLVER F O R NOISE R E M O V A L USING T H E T H I N P L A T E MODEL First, we discretize Equation 4.12, and then use the necessary condition of minimum XA u+u = 2 d (Equation 4.13) to define the fine and coarser grid operators, L\ for all / of different grid levels. These operators are the most complex and tedious to derive, thus we will do them last. The remaining operators, like the smoother and grid transfer operators can be chosen from standard techniques. For the smoother, we use lexicographical Gauss-Seidel iteration. For the grid transfer operators, we use full weighting or bicubic weighting for the restriction operator iff, and bilinear or bicubic interpolation for the prolongation operator Ijj. We are free to choose any pair of the grid transfer operators as long as the total order is greater than or equal to the order of the necessary condition, which is of order 4. Minimizing Equation 4.12 yields mm with discretization N-l i,j=Q where dij) + (u j i+1 - di j) +1 + (uj+ij+i - dj+ij+i) + (u j i i: + and Sij is defined depending on where the cell is in the domain. (i+lj+1) Cell (i+lj) - dij i) + } Appendix A A . 3 . 1 67 D E F I N I N G Si. Before defining Sij, lets define some intermediate labels: :\U 2 ul Arr/2 B o t h the U xx a n d the U yy +U +U 2 + r / +U 2 2 + U 2 2 + U 2 1 i are average over the four grid points on the cell. T h e b i g circle i n the following pictures denote the cell of interest, a n d the small circles o n the g r i d crossings denotes the g r i d points involved. A . 3 . 2 I N T E R I O R C E L L ( ~\ ( \) \ ) r \J (\ \) \ ) (b \ J\ V! \) (\ \J (\ \) -\ \) (\ -i) W e get the following: S = h {U 2 i3 2 + 2U + U } 2 x 2 y y Appendix A 68 B O U N D A R Y A.3.3 c1 KJ ( J•\ ( "\ t \ =4 N M b C E L L ( \ ( J t NM N= By the natural boundary conditions, U yy boundary. Su = h {Ulx + Wly + 2 = 0 on the grid point located on the horizontal \[U^ +l + U^ . +0 + 0]} i+l +l For the other three boundary cases, the appropriate rotation is made. A.3.4 C O R N E R o (DO C E L L Q Similarly, using the natural boundary condition again. Both U xx = 0 and U = 0 on the yy corner grid point. Su = h\wl y + \[Ul t+lj + U^ . i+1 +1 + 0 + 0] + - [U^ 4 ]+1 + + 0 + 0]} For the other three corners, similarly results are obtained by rotating the equation with respect to the center of the grid. Appendix A A.3.5 69 M I N I M I Z A T I O N N-l min F(Uh) = min ^ Dij + XSij i,j=o The necessary condition is dF dUij for all 0 < i,j < N. The following set of templates are obtained which define Lh- Squares indicate cells that involved in the minimization, and circle is the point (i,j). Note: Each of the template represents four actual templates which are the rotation variations of it. T E M P L A T E 0 • • • • p h • • • • 2 - 8 (Uij dij) + ^ 1 -8 2 20 -8 2 - 8 2 1 T E M P L A T E 1 1 • • • • 4p • • (Uij dij) + ^ 2 8 1 -8 19 2 -6 4 2 1 Appendix A TEMPLATE 70 2 2 • • • TEMPLATE (Ui j dij) + ^ -12 • 1 -I 16 2 -8 2 -6 18 2 -6 4 -8 1 3 • • p TEMPLATE 4 AND • • • • (Uij — dij) + ^ -8 1 =0 2 5 • 0 (Uij - dij) + 4 • -12 -6 4 15 -8 1 Template 4 is the above and template 5 is just the reflection of template 4 along the diagonal. TEMPLATE 6 • p • (Uij dij) + ^4 -12 12 =0 8 -12 2 After defining all these templates, we have indirectly defined the Li for all levels of grid. Thus, the multi-grid algorithm is completely defined and ready for execution. APPENDIX B S U R V E Y OF SOLVING METHODS The approximation problem can be written as a single over-determined system, we can solve it by least squares techniques. Recall Equation 4.16. ' wi* 's +1 ( w[ N fc+1 M fc+1 ]^ Let A= W [ f c + i ] s I Wlk+l]v[k+l] b= 0 V where A G » M x J V , x G $l , and b G 3 £ Nx3 control points in level k + 1 and k: M - M x 3 \ . The dimensions M and N are the number of m[ W l and N = m ^ n W . fc+1 fc+1 Now we can consider various way of solving the least square problem Ax = b. Since the system can be decoupled into three systems for each of the x, y, and z coordinates, we can solve them independently. Moreover, the A matrix remains unchanged, therefore, we need to perform factorization only once. 71 72 Appendix B B.I T H E N O R M A L E Q U A T I O N S M E T H O D The solution can be obtained by forming and solving the normal equations, A Ax = A b 1 l (B.I) This is the classical method of solving least squares problem which dates back to Gauss. Since our system has full rank, A A is symmetric and positive definite. 1 Therefore, we can use Cholesky factorization to solve it. Perform Cholesky factorization, A* A = LL 1 (B.2) followed by solving two triangular systems, Ly = A b l Lx = y l Gains in performance may be obtained by using various ordering of the columns of A . For example, minimum degree and reverse Cuthill-McKee ordering methods. However, the condition number, K(A'A) = K\A) and the formation of matrix A A may loose information. Thus, we need to compensate the 1 unstable property of normal equations approach by iterative refinement. A correction term 6x is computed from A A6x = A ' r for the residual r = b — Ax. A computed solution x is refined t by r = b - Ax R*y = A r l R6x = y Appendix B 73 x <— x + 6x Repeat the iteration until we are satisfied with the accuracy of x. This method is called the Corrected Normal Equations ( C N E ) (see [38]). B.2 T H E A U G M E N T E D S Y S T E M M E T H O D A least squares solution x with residual r — b - Ax must satisfy the augmented system I A A 0 1 V If the M first pivots of Gaussian elimination are chosen along the leading identity block, then the resulting system / r W I A 0 -A* A X ' v I x reduces to the system of normal equations. b ^ -A% The same instability properties are expected. However, if we introduce a scaling parameter a, we get a r x V / 0 This allows other more favorable pivot sequences. Bojorck [4] found an optimal choice of a with respect to numerical stability. 1 where a n is the smallest eigenvalue of A. Since we don't know a , n we can not find a t- Therefore, we adapt the heuristic used in M A T L A B with value of a equals: a = max _ " i,j 1000 op Appendix B B.3 T H E Q R M E T H O D The QR factorization decomposes a matrix A G sfl * M Q G $l MxM and an upper triangular matrix R G N into a product of an orthogonal matrix $l . NxN A = QR Since the Z/2-norm is invariant under orthogonal transformations, we can multiply Q through l Ax — 6, and the X2-norm of the resulting equation remains the same. Let Q*b then we have Ull2 {Ax-bW^WQ'Ax-Q'bW \\Q (QR)x-Q b\\ t t 2 ( R V \\R - \\ + \\c \ X 2 Cl 2 2 (B.3) Thus, the objective function \\Ax — b\\ can be minimized by solving just Rx = C\. The second 2 term is the squared residual norm, ||r|| = ||c2||. Pontus Matstoms [38] developed an algorithm for QR factorization of general large and sparse matrices based on multi-frontal techniques. Appendix B B.4 75 T H E C O N J U G A T E G R A D I E N T M E T H O D The conjugate gradient method of Hestenes and Stiefel with slight algebraic rearrangement can be used to solve least squares problem (see [5]). Let a* be an initial approximation. Take 0 r° = b- Ax° p° = s° = A*r° 7o =l W and for A; = 0,1,2,... compute q = Ap , k x k a = = x + ctkp , r Pk = fk+ihk, P k+1 k k f /\\q \\ , k k k+1 k+1 2 k =r k =s k+1 aq, k k + /3 p k k This is called the CGLS method. It requires only one extra matrix multiplication of A ' per iteration, compare with the original CG method. B.4.1 W I T H P R E C O N D I T I O N I N G The CGLS can also be preconditioned with preconditioner M as traditional CG method (see [5]). However, the cost imposed is not trivial. Let a; be an initial approximation. Take 0 r° = b - Ax° 70 and for k — 0,1,2,... compute = INT Appendix B 76 q = 1 x•*+* = x + k k+i s « = 7./lk"ll , a M.- p , k+i _ k _ qk^ AM~ p , k 2 k 1 k fc k r r ak _ M-^V^ , 1 Pk = lk+i/lk, Thus, the extra cost per iteration is used to solve the system y = M~ p and s = x Solving for y = M p _1 M~ (A r). t t is not bad, we can use similar technique as in ordinary preconditioned C G to solve it without forming M . However, the s = M~ (A r) -1 t t involves M~ which is not l easy to solve. For the case M corresponds to SSOR, Bjorck [6] developed efficient ways to solve it. B.5 T H E R E G U L A R I Z A T I O N M E T H O D Another view of the thin plate smoothness constraint is as a regulation technique applied to the ill-posed problem S V ^ = V t f c + 1 l . Lars Elden [20] and Martin Hanke et. al. [35] developed ways to convert general regularization problem into the standard form, which can then be solved efficiently. Given our problem Convert by Q R factorization and others in 0(M ) 2 to the standard form where var relate to original var by the transformation. The standard form can then be solved by bi-diagonalizations of W[ fc+1 ^S. w [fc+i] s _ u (A Appendix B 77 where U and V are orthogonal matrices, and J is the bi-diagonal matrix. The real advantage of this method is that, in the end it will produce an over-determined system with bi-diagonals and diagonal entries. This system can easily be solved again and again, when we change the value of A to try out the other solutions. min where J and c\ are pre-computed, and system is very small, only 0(M). = V£. Thus, the time to repeatedly solve this APPENDIX C SOLVING T H E OVERDETERMINED SYSTEM The overdetermined system (Equation 4.16), will be solved by various least squares methods described in Appendix B for comparison. Moreover, the sparsity of the matrix will be exploited. C.l A AND A A 1 Let's examine the matrix A and A A for inputs of different size control meshes. The A A matrix 1 1 is of interest, because the CNE (Corrected Normal Equations) method needs it. Normally, we would expect that A A will have more fills (non-zeros) than the original sparse A. However, 1 in our case, this turns out to be the opposite. We actually reduce the number of non-zeros by forming A A\ This is true for different input control meshes and different A values. Figure C.l l shows an example of the A and A A matrices. And Table C.l shows the number of non-zeros 1 of A and A A for difference size control meshes. t Size of Control mesh A A A 1 5x5 7x7 11 x 11 19 x 19 35 x 35 67 x 67 228 443 1083 3203 10803 39443 196 361 841 2401 7921 28561 Table C . l : Number of non-zeros in A and A A 1 78 Appendix 79 C (a) A (b) Figure C . l : The matrices A and A A t A'A with control mesh size 1 1 x 1 1 80 Appendix C C.2 KRONECKER PRODUCT The interesting thing to observe is that the A A matrices of all sizes have a self similar structure, 1 where if we treat each block as a single element, the whole matrix resembles the block itself (see Figure C.l). This is due to the fact, that we are dealing with a tensor product surface (B-spline) with limited support. Ignoring the weighting and smoothing part for a second, the tensor product view for each coordinate x, y, and z is BV C = Vt [k] l " non—vec f c + 1 l w where V, [fc] non—vec one coordinate of the 2D mW x control points matrix X of kth level without vectorization (i.e. a y[fc+i] = matrix). one coordinate of the k+1th level control mesh points after vectorization (i.e. m^ ^n^ ^ X k+1 k+1 1 matrix), B the B-spline refinement matrix for ID (i.e. C the B-spline refinement matrix for the other dimension (i.e. mJ- ^ X k+1 mW matrix), matrix). n^ ^ X k+1 The vectorization view of the same thing is: ( C\\B C\iB ... ^ v \ : ••• C [k] [k+i]B n n [fc] _ [fc+i] v J The intermixing B and C matrix above is the Kronecker product (see [2]) of C and B, denoted by C ® B. The Kronecker product algebra will have the normal equation (Equation B.l), (C ® B)\C ® 5)VW = (C ® Byv^+V Appendix C 81 rewritten as another Kronecker product, [ ( C ' C ) ® (fl*fl)]VM = (C* ® ^ j v f ^ ] 1 Thus, this is why we have the self similarity structure of A A. This way of viewing the fit, 1 leads to efficient fitting algorithm by decomposing the fitting problem into a sequence of curve fitting processes. This has been explored by David R. Forsey and Richard H. Bartels [24]. However, for our case here with smoothing introduced, although the structure of the A A still 1 has the Kronecker product structure, the smoothing part can not be represented as a tensor product. Thus, we can not take advantage of the dual views of tensor product. C.3 F I L L - I N Let's examine the fill-in of L + L* of Cholesky factorization of A A in Equation B.2. We will 1 compare the amount of fill-in with no ordering, reverse Cuthill-McKee, and minimum degree ordering (see Figure C.2). (a) No ordering (b) R C M (c) MD ordering Figure C.2: The matrices of L + L of A A with control mesh size 67 x 67 % 1 The minimum degree ordering is expected to produce the least fill-in in general. This is the case here for most of the time. The reverse Cuthill-Mckee is expected to perform much better then without ordering. However, it is the reverse for our case. This kind of pattern holds true Appendix C 82 even for different input control meshes and different A values, which could change the entries in A , or the sparsity of A due to zero weights. Table C.2 shows the fill-in of L + IS on the A A 1 in Section C . l . Size of Control mesh No ordering R C M ordering M D ordering 5 x 5 7 x 7 11 x 11 19 x 19 35 x 35 216 445 1281 5137 26961 212 457 1353 5777 32321 204 445 1295 5991 25183 67 x 67 170065 213537 144473 Table C.2: Number of non-zeros in L + L* C.4 CHOLESKY FACTORIZATION The number of fill-ins for minimum degree is usually the least, while reverse Cuthill-McKee is the worst. However, if we measure the number of floating point operations and C P U times to do the Cholesky factorization on the resulting ordered system plus the overhead of ordering, Cholesky factorization applied to no ordering is the fastest. In order words, when we applied the solving techniques, we are more interested in minimizing the computation time rather than minimizing the fill-ins. Of course, the total solving time and memory space depends on the fill-in of L + IS too. The next section will compare the actual total operations and computation time for solving the entire system for all 3 coordinates of the coarser control mesh points The following is Cholesky factorization on the set of different sizes matrix A from Section C . l . The machine it ran on is a SGI 4DVGX. Size of Control mesh No ordering M D ordering R C M ordering 5 x 5 7 x 7 11 X 11 19 x 19 3387 7316 22020 97628 3267 7338 22615 126859 3347 7492 23536 118864 35 x 35 67 x 67 664812 6602636 646241 6139044 971056 11170672 Table C.3: Number of floating point operations in Cholesky factorization with ordering Appendix C Size of Control mesh No ordering M D ordering R C M ordering 83 5x5 7 x 7 11 x 11 19 x 19 35 x 35 67 x 67 0.0333 0.0500 0.0667 0.1333 0.5333 3.7667 0.0500 0.0667 0.1000 0.2500 1.2333 9.5333 0.0733 0.0833 0.1000 0.2167 0.8500 5.6167 Table C.4: C P U time (in seconds) for Cholesky factorization with ordering C.5 C O M P A R I S O N OF T H E D I F F E R E N T SOLVING METHODS The followings test some of the solving methods in Section B on different size control meshes for the human head test data for a particular A. Tests were run for all the test data and with different A in reasonable useful range. However, only the results with respect to the human head data are shown here for conciseness. The Q R and Regularization solving method have not been tested, since they both require the sparse Q R factorization, which M A T L A B does not support, yet, and it is non-trivial to implement. For the other solving methods, we are more interested in terms of the user time to solve the system on all 3 coordinates of the control mesh. This is estimated by C P U time as shown in Table C.5. Table C.6 shows the number of floating point operations. This operation count ignores the work to do column ordering and other overheads. Errors are shown for each coordinate for each solving method. The C N E used are with two steps correction. Two steps are enough for all the test cases tried. The time for C G L S is large if require errors to be very small, so the stopping criteria is to solve up to 10e-6 change in residual of A b — A*Ax. The C G applied directly to (scaled) normal t equation is not very stable, so the stopping criteria-is the above for C G L S plus if the residual start to increase (in this case, it may diverge due to accumulation error). The residual error is shown in Table C.7. Appendix C 84 Size of Control mesh 5x5 7x7 C N E w / o ordering CNE w/ M D CNE w/ R C M Augmented System CGLS P C G L S w / diag. scaling C G directly on N E C G directly on scaled N E 0.0167 0.0500 0.05000.1000 0.3000 0.5667 0.0667 0.0667 0.0500 0.0500 0.0833 0.4333 0.7333 1.0333 0.0833 0.1167 k l l x 11 19 x 19 35 x 35 67 x 67 0.0833 0.1000 0.1167 2.0167 2.8500 3.5333 0.2.500 0.3333 0.4000 0.4500 0.4500 8.6833 13.4500 12.0000 0.5333 0.9667 1.8167 2.6000 1.8667 110.9000 62.7333 6.6.6500 4.3167 12.2833 7.8167 12.7500 10.2167 749.3500 220.6167 263.8667 20.3333 61.8333 Table C.5: C P U time (in seconds) to solve the overdetermined system Size of Control mesh C N E w / o ordering CNE w/ M D CNE w/ R C M Augmented System CGLS P C G L S w / diag. scaling C G directly on N E C G directly on scaled N E 5x5 7 x 7 11 x 11 19 x 19 35 x 35 67 x 67 13.804 29743 10946124 80736 2949,68 1502380 13468 29765 81583 339571 1451805 10021876 13771 29757 1905104 83548 327724 16296656 37357 288861 2113828 9848333 139560047 1161712834 111430 324397 1625416 7284865 38541559 126965372 147612 315968 1433746 6092994 33786798 126083194 19276 44923 146134 545230 3664677 16602967 16438 39658 163530 758668 10049005 40620965 Table C.6: Number of floating point operations to solve the overdetermined system Appendix CM 85 C 5x5 7x7 C N E w/o ordering X 1.8849e-14 2.1181e-14 2.8517e-15 5.4582e-15 y z 5.7370e-15 1.9325e-14 C N E w/ M D X 1.0391e-14 1.3758e-14 1.8420e-15 8.0413e-15 y z 4.4991e-15 2.0116e-14 C N E w/ R C M X 1.1096e-14 1.9160e-14 4.9110e-15 4.4878e-15 y z 4.6975e-15 2.2004e-14 Augmented System X 2.2248e-14 3.5994e-14 7.6193e-15 1.2949e-14 y z 1.0249e-14 2.8360e-14 CGLS X 1.0026e-07 2.5655e-ll 3.6455e-09 2.0846e-07 y z 1.2218e-09 3.8156e-09 P C G L S w / diagonal scaling X 9.4027e-10 2.9077e-09 6.6658e-08 8.6759e-09 y z 1.8293e-07 2.6246e-08 C G c irectly on N E X 1.8043e+01 1.5345e+01 2.7159e+00 4.0115e+00 y z 8.7654e-01 1.5279e+00 C G c irectly on sealed N E X 3.3559e+01 1.4981e+01 8.1622e+00 3.4391e+00 y z 4.5260e+00 3.0764e+00 11 x 11 19 x 19 35 x 35 67 x 67 3.5899e-14 1.9939e-14 2.7661e-14 8.8837e-14 4.5710e-14 6.1171e-14 1.3941e-13 8.9880e-14 1.2439e-13 2.6942e-13 2.2623e-13 2.6463e-13 4.3982e-14 2.0580e-14 3.1522e-14 9.2037e-14 4.3827e-14 6.0382e-14 1.4612e-13 8.9761e-14 1.2793e-13 2.8227e-13 2.2543e-13 2.6038e-13 3.4788e-14 1.9771e-14 3.2738e-14 9.1271e-14 4.5480e-14 5.9113e-14 1.3635e-13 8.8926e-14 1.2534e-13 2.7497e-13 2.3096e-13 2.6003e-13 5.0306e-14 6.1030e-14 6.2684e-14 9.9019e-14 8.3130e-14 9.5297e-14 2.0749e-13 3.8978e-13 2.7171e-13 4.6143e-13 4.2786e-13 4.6292e-13 2.1449e-07 2.0132e-07 6.6629e-08 5.3940e-07 6.0962e-07 3.8487e-07 6.7814e-07 7.5690e-07 9.0999e-07 8.8677e-07 8.5071e-07 9.0097e-07 4.0308e-07 3.9354e-08 1.9201e-07 7.5398e-07 6.8864e-07 1.0424e-06 8.1299e-07 1.1341e-06 1.2101e-06 8.5652e-07 1.0070e-06 9.2589e-07 7.2108e+00 7.2950e-01 3.3378e+00 1.7080e+00 6.5885e+00 2.4584e+00 2.2390e-01 4.3588e+00 6.3486e-01 3.0970e-01 1.1940e-01 3.1590e-01 2.0632e+00 1.1915e+00 5.2185e-01 5.8871e-02 5.3942e-01 2.5363e-01 1.3529e-05 9.7705e-07 2.4031e-03 8.6802e-07 9.2334e-07 8.5101e-07 Table C.7: The error residual A b - A Ax t % Glossary Animation: a motion picture made by photographing successive positions of inanimate objects or rendering successive positions of 3D objects in computer graphics. Deformation: alteration of form or shape by stress or external influence. Dynamics: the aspects of motion requiring the consideration of forces or accelerations. Discontinuity: lack of continuity or cohesion; for example, an edge. Fair: to smooth a surface, more uniformly distributing the changes in curvature of the surface. Force: an agency or influence that if applied to a free body results chiefly in an acceleration of the body and sometimes in deformation and other effects. Functional: a mapping from a function to a single real-number. Interpolate: to estimate values of a function between two known values. Level: a specific resolution of a multi-resolution surface. Mesh: a grid of points. Multi-resolution: analysis and synthesis using multiple spatial resolutions or levels. Non-convex: not convex in space possibly with many local extrema. Norm: a real-valued nonnegative function defined on a vector space to facilitate the measure of distance. Offset: a stored vector between the reference position of a vertex and it's actual position. Optimization: a process of making something as fully perfect, functional, or effective as possible. Outlier: something that statistically situated away from a main body of data or the control points that belong to a finer level of the B-spline hierarchical surface and not the current coarser level. Prolongation: to compute a higher-resolution representation of a multi-resolution system by interpolating the known values, the operator 86 Glossary 87 Refinement: to add additional degrees of freedom to a system by interpolating the known degrees of freedom. Regularization: a standard technique for solving ill-posed problem. Relaxation: a iterative procedure applied to solve a system of equations. Restrict: to compute a lower-resolution representation of a multi-resolution system by filtering a higher-resolution, the operator • Smooth: smooth with respect to the underlying control mesh of the surface and it is measured with energy norms like the thin plate or the membrane model. Sparse: of few and scattered elements (non-zeros for matrices). Threshold: a level, point, or value above which something is true or will take place and below which it is not or will not. Wavelet: a base which allows for transformation that cuts up data or functions or operators into different frequency components, and then studies each component with a resolution matched to its scale.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multiresolution surface construction for hierarchical...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multiresolution surface construction for hierarchical B-splines Wong, David 1995
pdf
Page Metadata
Item Metadata
Title | Multiresolution surface construction for hierarchical B-splines |
Creator |
Wong, David |
Date Issued | 1995 |
Description | This thesis presents a method for automatically generating a hierarchical B-spline surface suitable for animation and interactive modification from an initial set of control points. Given an existing mesh M[sup k+1] of control points, a mesh with half the resolution M[sup k+1], is constructed by simultaneously approximating the finer mesh and minimizing a smoothness or fairness constraint using weighted least squares. The membrane and thin plate models have been tried as smoothness constraints. Various numerical solving techniques have also been tested. Curvature measures of M[sup k+1] are used to identify features that need only be represented in the finer mesh. Some hierarchical B-spline surfaces, a B-spline surface and a digitized surface have been tested by this method. The resulting hierarchical surface generated accurately and economically reproduces the original mesh, is free from excessive undulations in the intermediate levels, and is well behave for editing and animation. |
Extent | 8477874 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-01-15 |
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.0051305 |
URI | http://hdl.handle.net/2429/3684 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1994-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1995-0246.pdf [ 8.09MB ]
- Metadata
- JSON: 831-1.0051305.json
- JSON-LD: 831-1.0051305-ld.json
- RDF/XML (Pretty): 831-1.0051305-rdf.xml
- RDF/JSON: 831-1.0051305-rdf.json
- Turtle: 831-1.0051305-turtle.txt
- N-Triples: 831-1.0051305-rdf-ntriples.txt
- Original Record: 831-1.0051305-source.json
- Full Text
- 831-1.0051305-fulltext.txt
- Citation
- 831-1.0051305.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-0051305/manifest