Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Multiresolution surface construction for hierarchical B-splines Wong, David 1995

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1995-0246.pdf [ 8.09MB ]
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

Full Text

M U L T I R E S Q L U T I O N S U R F A C E C O N S T R U C T I O N FOR H I E R A R C H I C A L B-SPLINES by D A V I D W O N G B.Sc , Simon Fraser University, 1994 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E IN T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F C O M P U T E R S C I E N C E We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA Apr i l , 1995 Copyright & © David Wong, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Ca*npt>t&-r S c . The University of British Columbia Vancouver, Canada D a t e 1 1 ^ A r . i , mi* DE-6 (2/88) A B S T R A C T 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[ f c + 1^ of control points, a mesh with half the resolution M ^ , 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 ^ 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. ii T A B L E O F C O N T E N T S 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.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 2.3 Surface approximation 12 2.3.1 Regularization 12 2.3.2 Outliers . . . . . . . . . . . . . 13 iii Chapter 3 Related Work 14 3.1 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 3.2 Multiresolution approximations 17 3.2.1 Multilevel visual surfaces with discontinuities 17 3.2.2 Wavelets 18 3.2.3 Hierarchical B-splines 19 3.3 Outliers handling 20 3.3.1 The outlier process 20 3.3.2 Robust statistics 21 3.4 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.2.1 The S matrix . - 42 5.2.2 The P matrix 43 5.2.3 The W t f c + 1 l matrix 44 5.3 Solving the system 44 v Chapter 6 Results 45 6.1 Hierarchical B-spline surfaces 45 6.1.1 Simple geometries 45 6.1.2 Complex geometries 49 6.1.3 Use in animation ; 51 6.2 B-spline surface 53 6.2.1 Fixed tolerance versus Variable tolerance 53 6.3 Digitized data 56 Chapter 7 Conclusions 58 Bibliography 59 Appendix A The multi-grid method 63 A . l A brief introduction to multi-grid method 63 A. 1.1 The multi-grid method 65 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 76 Appendix C Solving the overdetermined system 78 C. l A and A* A , 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 L I S T O F F I G U R E S 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 fit. 39 5.1 Three cases of the use of subdivision template 43 viii 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. . . . 54 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 A1 A with control mesh size 11x11 79 C.2 The matrices of L + L1 of A1 A with control mesh size 67 X 67 81 ix A C K N O W L E D G M E N T S I'd like to thank everyone who helped me with this research. My 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. D E D I C A T I O N I dedicate this thesis to my Mom and Dad for their life of hard working xi C H A P T E R O N E I N T R O D U C T I O N 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 broad-scale 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 B-spline 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 1 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 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 M O D E L L I N G P R O B L E M S A D D R E S S E D B Y T H I S T H E S I S 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 sur-face 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 hierar-chical 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 essen-tially 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 repre-sentation 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 (a) Perspective view (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 5 (e) Smooth - Level 0 (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 7 1.2 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 back-ground 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 implemen-tation. 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. C H A P T E R T W O B A C K G R O U N D 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 ex-port 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 re-quired. 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 sur-faces [22]. Chapter 2: Background 10 2.2.1 B A S I C FORMULATION Let, H[k](u,v) = J2Y. vi3Bffi(«)Bi*i(w) (2-i) t=0 j=o where surface H ^ ( t t , v) is defined by an n X m mesh of control vertices vj*- , and piecewise polynomial basis functions, B ^ ( i ) and Bjfc].(<) , of degree K and r respectively. The parameters u and v vary independently from some minimum umin, u T O t n to some maximum u m a x , vmax. 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,..., vn+T}. The basis 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 R E F I N E M E N T 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: R[k+i] = sv[k] (2 .2) where R i f c + 1 ] and denote column vectors of their corresponding mesh points. 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. 2.2.3 O F F S E T 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*^"1' in mesh M^" 1 " 1 ! is represented in offset form as: V[fc+i] = R[fc+i]0o[fc+1] (2.3) where Ri f c + 1 ] 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,-**1^ directly after a new level is created is equal to Rfj 1"^ and, by definition, o[-*^ = 0. Any change to the position of a control node VJ'j*"1^5 1S represented in the offset vector O J - * * 1 ^ as a relative change in position from the reference point R [ - * ^ . Changes to V,-^"1^ (and thus to pj^"1') do hot change the position of the reference point. On the other hand, changes in the shape of a parent surface alter RJ** 1 ' , result in a procedurally change of vj** 1 ' as defined by the reference-offset operator. 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 m i n ^ \ \ S ( u i , V i ) - 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 ' f c + 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. C H A P T E R T H R E E R E L A T E D 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. Al 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 Catmull-Clark 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(n2) where n is the 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 geomet-ric 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 MVC (Minimum Variation Curves) func-tionals [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 MVC functional produces very good result even in the parametric case. Since the MVC functional is compu-tationally 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 MVC functional respectively. These data depen-dent 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 dif-ferential 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 17 3 . 1 . 3 U S I N G M E S H 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 ex-tended 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 multi-grid 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 s u r f a c e discontinuities as one of distributed parameter iden-tification [52] [14] [43] within a variational f o r m u l a t i o n of v i s i b l e T S u r f a c e 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 con-version. 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(n2) 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 non-orthogonal 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 esti-mating 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 regu-larization 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 mea-sured 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. General-ized 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 trans-forms 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 con-taining 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 min-imum of the generalized cross-validation function for the smooth contour curve. The method consists of a linear search in log 1 0(A) scale to,find two points on the opposite sides of the 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 1 0(A) scale by evaluating the function at a number of equally spaced grid points to find the 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 M U L T I - R E S O L U T I O N B - S P L I N E A P P R O X I M A T I O N 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 ' f c + 1 l at level k + 1 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 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[/ , n e s'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 M.Mtnest\ 4.2 I N I T I A L A P P R O X I M A T I O N Given a control mesh, M' f c + 1l, with control points^ V^-k+1\ denning a B-spline surface, an initial 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 -no rm 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 deter-ministic 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 ] _ v^ 1! Since R ^ 4 " 1 ! = SV^, we can rewrite the equation as which leads to solving the linear system (4.1) 2 (4.2) (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 approxi-mation of the spike example is shown in Figure 4.2. 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 not the parametric surface H ^ l , thus, we treat the mesh points 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]. A Figure 4.2: Coarser level approximation by least squares. 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(||vifc]||2 + ||vifc]||2yWi; (4.4) 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 + \(u2x + u2y)dxdy (4.5) where (u — d)2 measures the deviation of the approximation surface from the data. By the calculus of variation, the unique solution to this minimization problem is by solving the functional d^}^ — 0. An extremely efficient multi-grid method [30] can be applied to find 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 hier-archy, 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 29 (e) A = 0.1 (f) A = 1.0 Figure 4.4: The effect of varying A on the membrane model. Chapter 4- Multi-resolution B-spline Approximation 30 4 . 5 T H E T H I N P L A T E 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 K?y[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 + )dudv (4.7) where A is the regularization parameter that controls the stiffness of the plate, and the sub-scripts 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 MVC 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 - 1 . 1 1 -1 (4-9) 14.10) Chapter 4'- Multi-resolution B-spline Approximation 31 (VJj)v„ = second order uv-derivative computed by central differencing 1 -2 1 (4.11) 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 + \(u2xx + 2u2xy + u2yy)dxdy (4.12) with the necessary condition of minimum by calculus of variation is \A2u + u = d (4.13) where the A2u — uxxxx + 2uxxyy + uyyyy is the biharmonic operator in 2D. There are many 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 Fig-ure 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 lin-ear 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 33 (e) A = 0.01 (f)A = 1.0 Figure 4.5: The effect of varying A on the thin plate model. Chapter 4- Multi-resolution B-spline Approximation 34 4.5.3 S P A R S E D A T A S E T F O R S M O O T H 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 35 (e) 30% data distribution (f) Approximate surface for 30% data Figure 4.6: The thin plate model on sparse data set. Chapter 4- Multi-resolution B-spline Approximation 36 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 S M O O T H N E S S C O N S T R A I N T Assuming that V^"1"1] will be given, we want to find V"W that minimizes the number of offsets 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 f c + 1J and replace the approximation surface points u(x, y) with the refined mesh points R i f c + 1 l = S V ^ . 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 mm svw - v[f c + 1i + A p'vM (4.14) where S is the subdivision matr ixof 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^k+1^ but not in R[* + 11 (the refined control mesh from V W ) . To get better result, 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. 4.7 W E I G H T E D 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 V t f c + 1 l and not R j * + 1 ] , and which belong to R [ f c + 1 l , we can estimate these geometric outliers by computing the second order partial derivatives in the control mesh V ^ 1 ) with respect 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 wt f c + 1 l . We then solve this'over-determined system using least-squares. / W [ * + i ] S \ / W [ f c + l ] y [ f c + l ] 0 V v /AP j V Which corresponds to the following minimization, problem (4.16) 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 39 Figure 4.9: Coarser level approximation, using weighted fit. 4.8 C A L C U L A T I N G 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. Iter-atively 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] _ sy[o] ( 4 - 1 8 ) To reduce storage space, we can threshold T(.) the offsets as we compute each level, only keeping those with large magnitudes. Note that, by so doing, we will change the for k > 0. The modified is denoted by v^. They are calculated as v[k] = Sv[k-i] + . r ( Q W ) (4.19) where is V ^ . 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 I M P L E M E N T A T I O N 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 A L G O R I T H M The basic algorithm for constructing of hierarchical B-spline surfaces is 1. Input the finest level control mesh M [ / , n e s l l of a B-spline surface with control points ylfmest] a n ( j v a m e f o r t n e 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 ^ * 1 ' to the approximation (see Section 4.7 and 5.2.3). (d) Combine the matrices, control mesh points and the regularization parameter A to form the Equation 4.16. 1 wi*+1is ^ , v A P , / wtfc+1]vffc+1J \ \ 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 ' F C ' for each level k during the construction of hierarchical B-spline. We will use a lexicographical order for the indexing in Y J 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 4 6 4 1 1 4 6 4 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, V r ^ 1 " 1 ' , 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 [k] Chapter 5: Implementation 43 vjj they apply to, and in the row corresponding to v{.^ +1^ of the S matrix. Figure 5.1 shows 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' f c ] . • • • v t - i , j - i yM V « , J + 1 • • • v i + l j - l v t + l , j vt fc ] • 1 v « + l , j + l • • •. 1 1 1 3 1 1 J _ 3 J _ q 16 2 1 6 3 1 3 1 1 1 3 1 1 16 4 16 2 16 4 16 2 16 4 16 2 16 4 1 16 1 16 &6 1 16 1 16 16^ 16^ 16^ 16^ (5.2) 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 E P 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: Implementation 44 for a particular control point V*|*- is v » - i , j v i - i , j + i ••• v » , j - i v i,j v t',j 1 V2 1 - 2 1 v/2 -v^ (5.3) 5.2.3 T H E wtfc+1l MATRIX The way to assign weights to the approximation is described in Section 4.7. The weight matrix Wf f c + 1 ] is a diagonal matrix, with weights Wij on the diagonal corresponding to each of the control point V?**1^. u>o,o W0,2 w. (5-4) 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(n2), where n is the total number of control 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{n2). If memory is a concern, the 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. C H A P T E R SIX R E S U L T S 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 S U R F A C E S Surface Original No. No. of offsets No. of B-spline of offsets obtained patches equivalent Spike (Figure 6.1) 17 17 34 ' 49 Hump (Figure 6.2) 19 19 40 143 Plop (Figure 6.3) 19 19 28 121 Eye (Figure 6.4) 422 253 565 2625 Dog (Figure 6.5) 257 302 609 1176 Mask (Figure 6.6) 478 467 1383 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 geome-tries. 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 4$ 6 .1 .2 C O M P L E X G E O M E T R I E S 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) Al l levels Figure 6.5: Dog head (609 patches, 257 offsets) Chapter 6: Results 51 6.1.3 U S E I N 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: Results 53 6 . 2 B - S P L I N E 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 Reboot (Figure 6.7) 648 658 600 576 Reboot w / var tol (Figure 6.8) 648 324 369 576 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 55 Figure 6.8: Reboot head: Generated from B-spline with variable tolerance. Chapter 6: Results 56 6.3 D I G I T I Z E D 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 No. of offsets obtained No. of patches B-spline equivalent Spock (Figure 6.9) 43520 6492 19567 40376 Table 6.3: Result for a digitized scan 3D data Level No. of patches No. of offsets 0 40 64 1 160 195 2 640 505 3 2551 1250 4 9178 2650 5 19567 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 57 fd) Level 3 (e) Level 4 (f) Level5 (g) All levels (h) Smiling Spock (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 interme-diate 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 representa-tion 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. Introduction to optimization theory. Wiley, New York, 1974. G. Alexander. Kronecker Products and Matrix Calculus with Applications. Halsted Press, 1981. D. M . Bates, M . J . Lindstrom, G. Wahba, and B. S. Yandel. Gcvpack - routines for gen-eralized 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. Technical report, Department of Mathematics, Linkoping University, 1967. 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. Visual reconstruction. The MIT Press, Cambridge, Mas-sachusetts, 1987. T. Boult. What is regular in regularization. In Proceedings of the First International Conference on Computer Vison, London, pages 457-462, 1987. 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 identifica-tion 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 IAM. 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 Analysis of Regularization and Constrained Least Squares Methods. Linkoping, 1977. [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, Apri l 1993. [25] D. R. Forsey and L. Wang. Multi-resolution surface approximation for animation. In Proceedings of Graphics Interface '93, pages 192-200, May 1993. [26] D. Geiger and R. A . M. Pereira. The outlier process. In Neural networks for signal processing, pages 60-69, 1991. [27] S. Geisser. The predictive sample reuse method with applications. In Journal American Statistic Association, pages 320-328, 1975. [28] B. R. Gossick. Hamilton's Principle and Physical Systems. Academic Press, New York and London, 1967. [29] G. Greiner. Variational design and fairing of spline surfaces. In Proceedings of Eurograph-ics '94, pages C143-C154, 1994. [30] W. Hackbusch. Multi-grid methods and applications. Computational mathematics. Springer-Verlag, 1985. 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. J in, 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. Regularization Methods for Large-Scale Problems, pages 253-315. Springer-Verlag, 1993. W. J . Kennedy and J . E. Gentle. 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 rep-resentation 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 Proceed-ings 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 Math-ematik, 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 syn-thesizing overview. In IEEE Transactions on pattern analysis and machine intelligence, pages 45-61, 1976. H. Qin and D. Terzopoulos. D-NURBS: Dynamic Non-Uniform Rational B-splines. In ACM Transactions on Graphics, pages 103-136, Apr i l 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, Apr i l 1993. Bibliography 62 [46] R. T. Rockafellar. The theory of subgradients and its application to problems of optimiza-tion: 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 210-218, June 1986. [48] B. Shahraray and D. J . Anderson. Optimal estimation of contour properties by cross-validated regularization. In IEEE Transactions on pattern analysis and machine intelli-gence, pages 600-610, June 1989. [49] M . Stone. Cross-validatory choice and assessment of statistical prediction. In Journal Royal Statistic Society Series B, pages 111-147, 1974. [50] R. Szeliski and D. Tonnesen. Surface modeling with oriented particle systems. In Pro-ceedings 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 Transac-tions 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 regularization method. In Soviet Mathematics Doklady, pages 1035-1038, 1963. [57] G. Wahba. A survey of some smoothing problems and the method of generalized cross-validation 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 cross-validation for determining the degree of smoothing. In Communications in statistics, pages 125-141, 1975. A P P E N D I X A T H E M U L T I - G R I D M E T H O D A. l A B R I E F I N T R O D U C T I O N T O M U L T I - G R I D 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^ finest — /finest on the finest grid for U f i n e s t , given Ljinest and /finest are defined in the following. A sequence of grids HQ, ill, • • •, Hh • • ••> Hfinest corresponding grid spacing ho, hi, . . . ,h\, . . . , h f i n e s t We assume hj = 2hj+i, ilj C Hj+i for all j C {0,1, . . .,/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 - c y c l e Figure A . l : One multi-grid iteration for V-cycle and W-cycle on a 5 levels grid Appendix A 65 A . 1 . 1 T H E 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~1(fi — Liii[) (restrict the computed defect) vi-i := 0 (initialize) for j:=l to 7 mgm(l — 1, vi-i, d/_i) (coarse grid correction) ui := ui + l\_xvi-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 dF}^ = 0. In the multi-grid language, the fine and coarser grid operator, X/, for all / in valid range, is a discretization of the functional d F } ^ at different level of grid. And f f i 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. Appendix A 66 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 THIN P L A T E M O D E L First, we discretize Equation 4.12, and then use the necessary condition of minimum XA2u+u = 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) + (ui+1j - di+1j) + (uj+ij+i - dj+ij+i) + (ui:j+i - dij+i) } and Sij is defined depending on where the cell is in the domain. ( i+ l j+1) Cell ( i + l j ) Appendix A 67 A . 3 . 1 D E F I N I N G Si. Before defining Sij, lets define some intermediate labels: ul :\U2 +U2 + U2 +U2 1 A r r / 2 + r / 2 + U 2 + U 2 i B o t h the Uxx and the Uyy are average over the four gr id points on the cell. The big circle in the following pictures denote the cell of interest, and the small circles on the gr id crossings denotes the gr id points involved. A . 3 . 2 I N T E R I O R C E L L ( ~\ ( r \ ) \ ) -\ \ ( J \ \ ( b \ ! ) \ \ ( ) - i \ ) \ ( J V \ ( ) \ \ ) \ J \ ) We get the following: Si3 = h2{U2x + 2U2y + U2y} Appendix A 68 A.3.3 B O U N D A R Y C E L L c 1 ( ( K •\ ( J \ "\ ( t =4 J \ N M b J t N M N = By the natural boundary conditions, Uyy = 0 on the grid point located on the horizontal boundary. Su = h2{Ulx + Wly + \[U^+l + U^i+l.+l +0 + 0]} For the other three boundary cases, the appropriate rotation is made. A.3.4 C O R N E R C E L L o ( D O Q Similarly, using the natural boundary condition again. Both Uxx = 0 and Uyy = 0 on the corner grid point. Su = h\wly + \[Ult+lj + U^i+1.+1 + 0 + 0] + -4[U^]+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 69 A.3.5 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 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. dF dUij T E M P L A T E 0 • • • p • • h • • • (Uij dij) + ^ 2 - 8 2 1 - 8 20 - 8 1 2 - 8 2 1 T E M P L A T E 1 • • • 4 p • • • (Uij dij) + ^4 1 -1 2 8 19 2 - 6 8 2 Appendix A 70 T E M P L A T E 2 • • • • (Ui j dij) + ^ 1 -I 2 -12 4 16 - 8 1 T E M P L A T E 3 • • p • • (Uij — dij) + ^ 2 - 8 -6 18 2 - 6 2 - 8 1 2 = 0 T E M P L A T E 4 A N D 5 • • • 0 • (Uij - dij) + 4 - 1 2 4 -6 15 - 8 1 Template 4 is the above and template 5 is just the reflection of template 4 along the diagonal. T E M P L A T E 6 • p • (Uij dij) + ^4 -12 8 12 -12 2 = 0 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 O F S O L V I N G M E T H O D S The approximation problem can be written as a single over-determined system, we can solve it by least squares techniques. Recall Equation 4.16. Let ' wi*+1's N ( w[ f c + 1M f c + 1] ^ A = W [ f c + i ] s b = I Wlk+l]v[k+l] \ V 0 where A G » M x J V , x G $lNx3, and b G 3 £ M x 3 . The dimensions M and N are the number of control points in level k + 1 and k: M - m[fc+1W fc+1l and N = m ^ n W . 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 Appendix B 72 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, A1 Ax = Alb (B.I) This is the classical method of solving least squares problem which dates back to Gauss. Since our system has full rank, A1 A is symmetric and positive definite. Therefore, we can use Cholesky factorization to solve it. Perform Cholesky factorization, A* A = LL1 (B.2) followed by solving two triangular systems, Ly = Alb Llx = y Gains in performance may be obtained by using various ordering of the columns of A . For ex-ample, minimum degree and reverse Cuthil l-McKee ordering methods. However, the condition number, K ( A ' A ) = K\A) and the formation of matrix A1 A may loose information. Thus, we need to compensate the unstable property of normal equations approach by iterative refinement. A correction term 6x is computed from AtA6x = A ' r for the residual r = b — Ax. A computed solution x is refined by r = b - Ax R*y = R6x = Alr 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 (CNE) (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 A1 0 V If the M first pivots of Gaussian elimination are chosen along the leading identity block, then the resulting system / I A 0 -A* A W r X ' b ^ v x I -A% reduces to the system of normal equations. 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 an is the smallest eigenvalue of A. Since we don't know an, we can not find aopt-Therefore, we adapt the heuristic used in M A T L A B with value of a equals: a = max _ " i,j 1000 Appendix B B . 3 T H E Q R M E T H O D The QR factorization decomposes a matrix A G sflM*N into a product of an orthogonal matrix Q G $lMxM and an upper triangular matrix R G $lNxN. A = QR Since the Z/2-norm is invariant under orthogonal transformations, we can multiply Ql through 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 \\Qt(QR)x-Qtb\\2 ( V R \\RX - Cl\\2 + \\c22\ (B.3) Thus, the objective function \\Ax — b\\2 can be minimized by solving just Rx = C\. The second 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 75 B . 4 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*0 be an initial approximation. Take r° = b- Ax° p° = s° = A * r ° 7o = l W and for A; = 0,1,2,... compute qk = Apk, ak = fk/\\qk\\2, xk+1 = xk + ctkpk, r k + 1 = rk - akqk, Pk = fk+ihk, Pk+1 = sk+1 + /3kpk 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;0 be an initial approximation. Take r° = b - Ax° 70 = INT and for k — 0,1,2,... compute Appendix B 76 x qk = AM~1pk, •*+* = xk + akM.-1pk, sk+i _ M - ^ V ^ 1 , Pk = lk+i/lk, « f c = 7 . / l k " l l 2 , rk+i _ rk _ akqk^ Thus, the extra cost per iteration is used to solve the system y = M~xp and s = M~t(Atr). Solving for y = M_1p is not bad, we can use similar technique as in ordinary preconditioned C G to solve it without forming M-1. However, the s = M~t(Atr) involves M~l which is not 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 Mart in Hanke et. al. [35] developed ways to convert general regularization problem into the standard form, which can then be solved efficiently. Given our problem where var relate to original var by the transformation. The standard form can then be solved by bi-diagonalizations of W[ f c + 1^S. Convert by QR factorization and others in 0(M2) to the standard form w [ f c + 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 = V£. Thus, the time to repeatedly solve this system is very small, only 0(M). -APPENDIX C S O L V I N G T H E O V E R D E T E R M I N E D S Y S T E M 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 A1 A Let's examine the matrix A and A1 A for inputs of different size control meshes. The A1 A matrix is of interest, because the CNE (Corrected Normal Equations) method needs it. Normally, we would expect that A1 A will have more fills (non-zeros) than the original sparse A. However, in our case, this turns out to be the opposite. We actually reduce the number of non-zeros by forming AlA\ This is true for different input control meshes and different A values. Figure C.l shows an example of the A and A1 A matrices. And Table C.l shows the number of non-zeros of A and At A for difference size control meshes. Size of Control mesh 5 x 5 7 x 7 11 x 11 19 x 19 35 x 35 67 x 67 A 228 443 1083 3203 10803 39443 A1 A 196 361 841 2401 7921 28561 Table C . l : Number of non-zeros in A and A1 A 78 Appendix C 79 ( a ) A (b) A'A Figure C . l : The matrices A and AtA with control mesh size 1 1 x 1 1 Appendix C 80 C.2 K R O N E C K E R P R O D U C T The interesting thing to observe is that the A1 A matrices of all sizes have a self similar structure, 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[k] Cl = V t f c + 1 l " non—vec w where V , [fc] non—vec one coordinate of the 2D mW x control points matrix of kth level without vectorization (i.e. a X matrix). y[fc+i] = one coordinate of the k+1th level control mesh points after vectorization (i.e. m^k+1^n^k+1^ X 1 matrix), the B-spline refinement matrix for ID (i.e. mJ-k+1^ X mW matrix), the B-spline refinement matrix for the other dimension (i.e. n^k+1^ X matrix). B C The vectorization view of the same thing is: ( C\\B C\iB ... ^ \ : ••• Cn[k]n[k+i]B J v[fc] _ v[fc+i] 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 ) V W = (C ® Byv^+V Appendix C 81 rewritten as another Kronecker product, [ ( C ' C ) ® ( f l*f l ) ]VM = (C* ® ^ j v f ^ 1 ] Thus, this is why we have the self similarity structure of A1 A. This way of viewing the fit, 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 A1 A still 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 A1 A in Equation B.2. We will compare the amount of fill-in with no ordering, reverse Cuthil l-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 A1 A with control mesh size 67 x 67 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 1 A in Section C . l . Size of Control mesh 5x5 7x7 11 x 11 19 x 19 35 x 35 67 x 67 No ordering 216 445 1281 5137 26961 170065 R C M ordering 212 457 1353 5777 32321 213537 M D ordering 204 445 1295 5991 25183 144473 Table C.2: Number of non-zeros in L + L* C . 4 C H O L E S K Y FACTORIZATION The number of fill-ins for minimum degree is usually the least, while reverse Cuthil l-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 5x5 7x7 11 X 11 19 x 19 35 x 35 67 x 67 No ordering 3387 7316 22020 97628 664812 6602636 M D ordering 3267 7338 22615 126859 646241 6139044 R C M ordering 3347 7492 23536 118864 971056 11170672 Table C.3: Number of floating point operations in Cholesky factorization with ordering Appendix C 83 Size of Control mesh 5x5 7x7 11 x 11 19 x 19 35 x 35 67 x 67 No ordering 0.0333 0.0500 0.0667 0.1333 0.5333 3.7667 M D ordering 0.0500 0.0667 0.1000 0.2500 1.2333 9.5333 R C M ordering 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 DIFFERENT 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 QR and Regularization solving method have not been tested, since they both require the sparse QR 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 Atb — A*Ax. The C G applied directly to (scaled) normal 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 5 x 5 7 x 7 k l l x 11 19 x 19 35 x 35 67 x 67 C N E w/o ordering 0.0167 0.0500 0.0833 0.4000 1.8167 7.8167 C N E w/ M D 0.0500 0.0500 0.1000 0.4500 2.6000 12.7500 C N E w/ R C M 0.0500- 0.0833 0.1167 0.4500 1.8667 10.2167 Augmented System 0.1000 0.4333 2.0167 8.6833 110.9000 749.3500 C G L S 0.3000 0.7333 2.8500 13.4500 62.7333 220.6167 P C G L S w/ diag. scaling 0.5667 1.0333 3.5333 12.0000 6.6.6500 263.8667 C G directly on N E 0.0667 0.0833 0.2.500 0.5333 4.3167 20.3333 C G directly on scaled N E 0.0667 0.1167 0.3333 0.9667 12.2833 61.8333 Table C .5: C P U time (in seconds) to solve the overdetermined system Size of Control mesh 5x5 7x7 11 x 11 19 x 19 35 x 35 67 x 67 C N E w/o ordering 13.804 29743 80736 2949,68 1502380 10946124 C N E w/ M D 13468 29765 81583 339571 1451805 10021876 C N E w/ R C M 13771 29757 83548 327724 1905104 16296656 Augmented System 37357 288861 2113828 9848333 139560047 1161712834 C G L S 111430 324397 1625416 7284865 38541559 126965372 P C G L S w/ diag. scaling 147612 315968 1433746 6092994 33786798 126083194 C G directly on N E 19276 44923 146134 545230 3664677 16602967 C G directly on scaled N E 16438 39658 163530 758668 10049005 40620965 Table C.6: Number of floating point operations to solve the overdetermined system Appendix C 85 C M 5 x 5 7 x 7 11 x 11 19 x 19 35 x 35 67 x 67 C N E w/o ordering X 1.8849e-14 2.1181e-14 3.5899e-14 8.8837e-14 1.3941e-13 2.6942e-13 y 2.8517e-15 5.4582e-15 1.9939e-14 4.5710e-14 8.9880e-14 2.2623e-13 z 5.7370e-15 1.9325e-14 2.7661e-14 6.1171e-14 1.2439e-13 2.6463e-13 C N E w/ M D X 1.0391e-14 1.3758e-14 4.3982e-14 9.2037e-14 1.4612e-13 2.8227e-13 y 1.8420e-15 8.0413e-15 2.0580e-14 4.3827e-14 8.9761e-14 2.2543e-13 z 4.4991e-15 2.0116e-14 3.1522e-14 6.0382e-14 1.2793e-13 2.6038e-13 C N E w/ R C M X 1.1096e-14 1.9160e-14 3.4788e-14 9.1271e-14 1.3635e-13 2.7497e-13 y 4.9110e-15 4.4878e-15 1.9771e-14 4.5480e-14 8.8926e-14 2.3096e-13 z 4.6975e-15 2.2004e-14 3.2738e-14 5.9113e-14 1.2534e-13 2.6003e-13 Augmented System X 2.2248e-14 3.5994e-14 5.0306e-14 9.9019e-14 2.0749e-13 4.6143e-13 y 7.6193e-15 1.2949e-14 6.1030e-14 8.3130e-14 3.8978e-13 4.2786e-13 z 1.0249e-14 2.8360e-14 6.2684e-14 9.5297e-14 2.7171e-13 4.6292e-13 C G L S X 1.0026e-07 2.5655e-ll 2.1449e-07 5.3940e-07 6.7814e-07 8.8677e-07 y 3.6455e-09 2.0846e-07 2.0132e-07 6.0962e-07 7.5690e-07 8.5071e-07 z 1.2218e-09 3.8156e-09 6.6629e-08 3.8487e-07 9.0999e-07 9.0097e-07 P C G L S w/ diagonal scaling X 9.4027e-10 2.9077e-09 4.0308e-07 7.5398e-07 8.1299e-07 8.5652e-07 y 6.6658e-08 8.6759e-09 3.9354e-08 6.8864e-07 1.1341e-06 1.0070e-06 z 1.8293e-07 2.6246e-08 1.9201e-07 1.0424e-06 1.2101e-06 9.2589e-07 C G c irectly on N E X 1.8043e+01 1.5345e+01 7.2108e+00 1.7080e+00 2.2390e-01 3.0970e-01 y 2.7159e+00 4.0115e+00 7.2950e-01 6.5885e+00 4.3588e+00 1.1940e-01 z 8.7654e-01 1.5279e+00 3.3378e+00 2.4584e+00 6.3486e-01 3.1590e-01 C G c irectly on sea led N E X 3.3559e+01 1.4981e+01 2.0632e+00 5.8871e-02 1.3529e-05 8.6802e-07 y 8.1622e+00 3.4391e+00 1.1915e+00 5.3942e-01 9.7705e-07 9.2334e-07 z 4.5260e+00 3.0764e+00 5.2185e-01 2.5363e-01 2.4031e-03 8.5101e-07 Table C.7: The error residual Atb - A% Ax Glossary Animation: a motion picture made by photographing successive positions of inanimate ob-jects 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 filter-ing 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. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items