Multi-resolution Surface Approximation for Animation By LIFENG WANG B.Sc., University of Science and Technology of China China, 1990 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (DEPARTMENT OF COMPUTER SCIENCE) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April, 1993 © Lifeng Wang, 1993 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. (Signature) Department of ^Computer Science The University of British Columbia Vancouver, Canada ^April 28th, 1993 Date DE-6 (2/88) Abstract This thesis addresses the problem of approximating a set of gridded data points obtained from a three-dimensional digitizing system to create a representation with a hierarchical bicubic B-spline surface that is suitable for further manipulation and animation. Chord length parameterization is obtained using 2-D deformation technique. A full multigrid (FMG) numerical method is used to solve the surface approximation and the multi-resolution elements created are used directly to define the overlays in a hierarchical B-spline surface. The direct use of FMG multi-resolution data offers reasonable surface shape behaviour, but the number of non-zero offsets is large. Storage cost is reduced either by eliminating offsets whose magnitude is below a certain tolerance or by reducing all offsets in a given level by a user specified amount. The resulting spline surface is modifiable, both locally and globally while retaining surface details of the digitized data. An interactive system based on these methods was created and the results of approximating two large data sets are presented. ii Contents Abstract^ ii Contents^ iii List of Figures^ vi Acknowledgement^ ix 1 Introduction 1.1 1.2 1.3 1.4 1.5 2 1 Overall Goal ^ Surface Approximation ^ 1.2.1^Mathematical Formulation ^ 1.2.2^Discrete Data ^ 1.2.3^Interpolation and Approximation ^ Problem Definition ^ 1.3.1^Test Data ^ 1.3.2^Data Parameterization ^ Related Work ^ Overview ^ 1 2 2 2 3 4 4 6 6 8 Mathematical Background 10 2.1 2.2 2.3 10 14 16 16 17 19 20 22 23 23 25 2.4 2.5 Uniform Bicubic B-Spline Surfaces ^ Subdivision Scheme ^ Hierarchical B-splines ^ 2.3.1^Practical Consideration of Tensor-Product Surfaces ^ 2.3.2^Localized Subdivision-Hierarchical B-splines ^ Parameterization ^ 2.4.1^Parameterization Schemes ^ 2.4.2^Summary ^ Multigrid Numerical Method ^ 2.5.1^Iterative Methods on Linear Algebraic Equations ^ 2.5.2^Multigrid Philosophy ^ iii ^ ^2.5.3^Basic Strategies ^ 2.5.4^Multigrid Operators and Parameters ^ 2.5.5^Multigrid V-cycle ^ 2.5.6^Full Multigrid V-cycle ^ 2.5.7^Performance Evaluation ^ 3 Parametric Space Mapping 3.1 3.2 3.3 3.4 3.5 4 Multi-resolution Surface Fitting Using Multigrid 4.1 4.2 4.3 5 Deriving Equations and Constraints ^ 4.1.1^Data Point Distribution ^ 4.1.2^Deriving the Positional Interpolation Equations ^ 4.1.3^Deriving the Boundary Conditions ^ Solving For Control Vertices ^ 4.2.1^Full Multigrid Method ^ 4.2.2^Full Multigrid Operators ^ Performance Evaluation ^ 34 34 36 38 38 40 40 47 48 49 49 49 51 52 53 53 54 59 Hierarchy Building 62 5.1 62 62 64 67 68 68 70 71 76 76 5.2 5.3 5.4 5.5 5.6 6 Parameterization ^ Initial Parameterization ^ Deforming Parametric Space ^ 3.3.1^Deformation Function ^ 3.3.2^Refinement and Local Control ^ 3.3.3^Deformation ^ Sampling Data Points ^ Summary ^ 26 27 27 30 31 Inverse Subdivision ^ 5.1.1^Hierarchical Surface File Format ^ 5.1.2^Inverse Subdivision Scheme ^ 5.1.3^Remarks: Solution to Triangular Linear Equation ^ 5.1.4^Problems with Inverse Subdivision ^ Use of MG Data ^ Zeroing Non-zero Offsets ^ Smoothing Algorithm ^ Manipulation of Resulting Surface ^ Implementation and Results ^ Conclusions and Future Work 6.1 Conclusions ^ 6.2 Future Work ^ 78 78 79 iv Bibliography^ 80 A Plates^ 83 Appendix^ 83 List of Figures 1.1 Digitized Data Topology. ^ 5 2.1 A bicubic B-spline surface with a single patch defined by 16 control vertices.. .^11 2.2 A B-spline surface is a mosaic of surface patches ^ 11 2.3 The parametric range for a surface is the union of the parametric range for each patch. ^ 13 2.4 Subdivision splits a single patch into four patches. ^ 15 2.5 Subdivision splits the parametric space at the midpoint ^ 16 2.6 Subdivision through non-local knot insertion. ^ 17 2.7 A wave on Sl h (N = 16) is projected onto C2 2h (N = 8). On the coarse grid, the wave seems more oscillatory than on the fine grid. ^ 25 2.8 The schedule for grid visiting when the multigrid V-cycle is executed. Each level is identified with its spacing. Multigrid V-cycle starts from the finest level and proceeds to the coarsest level, then returns to the finest level ^ 28 2.9 The schedule for grid visiting of W-cycle scheme ^ 30 2.10 Full multigrid grid visiting schedule. ^ 30 3.1 Labels for the tensor-product B-spline surface. ^ 35 3.2 Two approaches to approximating the head data: a) wrap a cylindrical surface around the head, b) fold the four edges of a rectangular surface to cover up the head with the opening at the neck ^ 35 3.3 The initial mapping from data space to parametric space. Each row of data is mapped onto a square in parametric space starting from the center to the edge. 36 3.4 Two coordinate frames imposed on the parametric space. 38 3.5 A control grid imposed on the parametric space with grid size of m = n = 7. Each control point is undeformed and has the same coordinates in UV space as 39 in ST space. ^ 3.6 The control grid in UV space is uniformly refined from a 2 x 2 grid into a 3 x 3 grid. The control points inside brackets are the indices of the 4 control points as they are indexed at level 1. ^ 42 3.7 The three cases of the sample point falling into a data quadrilateral. ^ 45 3.8 A control point is displaced from its original position. ^ 46 vi ^ 3.9 The UV space is evenly sampled for produce the parametric pair set 48 4.1 The (X +1) x (Y +1) data points are mapped onto X x Y patches. Every patch corner corresponds to a data point ^ 50 4.2 A data point Z23 is mapped onto the upper right corner of a single patch. At the same time, it is the corner point of three adjacent patches. ^ 51 4.3 A illustration of the four-colour relaxation scheme on a uniform grid. First the 0 points are swept, then the 0 points, A points and finally ^ points. ^ 54 5.1 The inverse subdivision process generate a lower level curve from a upper level curve in a way that reduces the number of non-zero offsets while the hierarchy is built using subdivision. 65 5.2 A portion of hierarchical B-pline data file without zeroing non-zero offsets. . . . ^70 5.3 Offset number generated in hierarchical surface with given tolerance value. Raw data ranges from -100 to +100. ^ 72 5.4 Residual per data point caused by offset zeroing with given tolerance value. Raw data ranges from -100 to +100. ^ 72 5.5 Maximum Residual caused by offset zeroing with given tolerance value. Raw data ranges from -100 to +100 ^ 73 5.6 Offset zeroing side effect. a) Before the offset 2 is zeroed. b) The zeroing of offset 2 results in a bump. ^ 74 5.7 Offset zeroing with smoothing function. a) Before the offset 2 is zeroed. b) After offset 2 is zeroed with smoothing function. ^ 75 A.1 Victor Hugo raw data ^ 84 A.2 Deformed UV isoparametric lines in ST space with 5 level refinement. ^ 84 A.3 Deformed UV isoparametric lines in ST space with 5 level refinement. Data points are plotted over the space ^ 85 A.4 FMG fit at level 1 and level 2. ^ 86 A.5 FMG fit at level 3 and level 4. ^ 86 A.6 FMG fit at level 5 and level 6. ^ 87 87 A.7 FMG fit at level 7 and level 8. ^ A.8 Left: FMG fit at level 9. Right: Hierarchical B-spline surface without zeroing offsets. Offset number is 90949 ^ 88 A.9 Left: Hierarchical B-spline surface with offset zeroing threshold 0.5 and smoothing algorithm. Offset Number is 3913. Right: Hierarchical B-spline surface with offset zeroing threshold 0.5 without smoothing algorithm. 88 A.10 Left: Hierarchical B-spline surface with offset zeroing threshold 0.1 and smoothing algorithm. Offset Number is 23270. Right: Hierarchical B-spline surface with offset zeroing threshold 1.0 and smoothing algorithm. Offset Number is 1918. . . 89 vii A.11 Left: Hierarchical B-spline surface with offset zeroing threshold 2.0 and smoothing algorithm. Offset Number is 930. Right: Hierarchical B-spline surface with offset zeroing threshold 5.0 and smoothing algorithm. Offset Number is 318. . . . 89 A.12 Left: Surface manipulation - Smile 1. Right: Surface manipulation - Smile 2. . . 90 A.13 Left: Surface manipulation - Smile 3. Right: Surface manipulation - Nean90 derthal man. ^ 91 A.14 Another fit to a digitized face. ^ viii Acknowledgements I would like to thank Dr. David R. Forsey, my thesis supervisor, for his guidance and encouragement throughout my work on this thesis. The valuable ideas and suggestions from him improved the content and presentation of this thesis significantly. I would also like to thank the second reader of my thesis, Dr.Peter Cahoon for reading through the draft of this thesis, his comments and help made the thesis complete. Karen Kuder read the first draft of my thesis and spent a lot of time changing my grammar and spelling error. I feel much appreciated for her very great care on reading my thesis. Dr. Pierre Poulin volunteered to read the draft paper presented on GI93 which is based on this thesis. His comments have helped a lot on improving my work. Many other people in the department offered time and effort my work. It is difficult to list all the individuals who has been involved with me. But those people who have provided help know that they deserved a note of THANKS. If this is the chance, I would like to direct the most important acknowledgement to my parents, who have been providing priceless love and support over the past. Without their encouragement and comforts, most of the opportunities I enjoyed would not have been possible, including the opportunity to come to Canada to complete the work I present here. ix Chapter 1 Introduction 1.1 Overall Goal This research is a part of a long term project to create a system for the animation and rendering of human animal forms. An integral part of this system is a mechanism whereby digitized surface information is used to determine the surface characteristics ( shape and colour) of an animated figure. The underlying surface formulation for this system is hierarchical B-spline [17] [20], a multi-resolution approach to modelling with tensor product surfaces. Thus the mechanisms for approximating an arbitrary surface definition with a hierarchical B-spline surface are required. This thesis takes the first step toward this goal by addressing the sub-problem of approximating regular gridded data with a hierarchical spline. The contribution of this work is the use of a multi-resolution approach to solving system of linear equations (the full multigrid method, i.e. FMG) to directly generate the multi-resolution surface of hierarchical B-spline. Surface approximation for animation adds one additional complication. Because the approximating surface may only be a part of a much larger surface, the question of what part of the surface approximates what part of the data must also be addressed. 1 Chapter 1. Introduction^ 2 1.2 Surface Approximation 1.2.1 Mathematical Formulation We are interested in fitting a tensor-product B-spline surface =E E m n S(L,V) to a given set of data points. To generate the equations for surface approximation, each data point (x, y, z) = DA E R 3 is associated with a domain point (u, v) = SA E R 2 of the spline. This forms the equation: E E^ m n DA. (1.2) i=0 j=0 Data is gridded if u E {2/0, • • • , um} E {vo, • " , vN} (1.3) and if the Ss consist of all points in {Ito, • • • , ^x {vo, • • • , vN} the surface approximation equations become E E^ m n (1.4) i=0 j=0 for s = 0, • ,M and t = 0, • • • , N. Thus the number of equations is equal to the number of data points. Hierarchical splines are a multi-resolution approach to splines for use in interactive creation of free-form surfaces [17] [40]. The problem of applying these equations to approximate a digitized surface with a hierarchical spline is complicated by two factors: an unusual parameterization, and the multi-resolution nature of the hierarchical formulation itself. Readers who are familiar with parameterization and surface fitting can skip the following sections and go directly to Chapter 3. 1.2.2 Discrete Data Systematic measuring devices such as laser scanners, imaging systems and optical scanners are used to measure the positions of points on an object's surface with a resolution high enough Chapter 1. Introduction^ 3 to provide an accurate representation of the original object. These devices work in three dimensional space and typically produce a 2-D array of data points. A medical imaging system like CT (Computed Tomography) records both exterior and interior information, producing a regularly spaced 3-D lattice. Raw data sets for human figures are usually too large to manipulate in real-time on general purpose graphics workstations. Techniques used to approximate the raw data must be very efficient both in terms of computing time and storage complexity. 1.2.3 Interpolation and Approximation A common technique for recreating the smooth surface from a set of discrete digitized data points discrete data set as a continuous shape is a graph. Other techniques include fitting surfaces (interpolation and approximation), volume rendering and various visualization techniques. There are basically two ways to reconstruct the continuous model from digitized data: interpolation and approximation. Interpolation requires that the surface passes through every data point. Piecewise parametric polynomial interpolation is a common technique [14] and there are various other mathematical formulations which meet different requirements [4] [16] [24]. Most interpolation methods are simple in concept, have unique solutions and provide an accurate representation of the original data. Approximation is designed to fairly represent the data without necessarily interpolating the individual data points and are often useful for data sets comprised of measured values including errors where their exact interpolation is not desired. Approximate algorithms can be adaptive or non-adaptive. Adaptive algorithms begin with a simple approximation defined by a small number of coefficients. The number of degrees of freedom in the approximation is adaptively increased in areas of high error until the surface is within a given tolerance of the data. Adaptive algorithms are efficient because some sub-regions of data have less gradient variation than other regions. The smoother regions can be approximated with coarse meshes without exceeding certain error tolerances, reducing the storage cost. The bumpy regions should be approximated with fine meshes. Chapter 1. Introduction ^ 4 There are numerous approaches to both approximation and interpolation. Each mathematical representation will produce different approximations for the same set of data. Piecewise parametric surfaces are one family of tools for approximating shapes [14]. Notable examples of such formulations include Bezier curves [14], Coons patches [11], and B-splines [3]. These representations have been used extensively throughout the CAGD field to describe shape in computer systems. Some mathematical background concerning these representations will be given in Chapter 2. More details can be found in books by FarM [14] and by Bartels et. al. [3]. 1.3 Problem Definition 1.3.1 Test Data The data set chosen for use in this thesis is representative of the size and detail of the models used in animation, the bust of Victor Hugo (courtesy F. Schmidt, see Plate A.1) obtained from a commercial three-dimensional digitizing system developed by the Laboratoire IMAGE of the Ecole Nationale Superieure des Telecommunications [34]. In this integrated device, the microprocessor controls a video-laser data acquisition system capable of measuring approximately 200,000 points per minute. The rotation axis (Y axis in Figure 1.1) is usually perpendicular to the beam and along the center line of the object. Figure 1.1 shows the topology of the data obtained from these kinds of devices. The resulting data are represented as a 2-D array of radii; that is, any non-boundary data point is connected to exactly four neighbouring points, forming a cylindrical contour(Figure 1.1). The device measures the distance from the surface point to the central axis, but before approximating the surface, the radius values rij are converted to the coordinate triple (xij, yij, zij) in R 3 as follows: xij = rij x cos(j), yij = M —^ zij = —rij x sin(j), (1.5) where M is the number of rows of points. In the Victor Hugo data set, there are 264 rows of 361 data points (95304 points in all). For notational purposes, the data is indexed by row number {i, 0 < i < 263} and by column number {j, 0 < j < 360}. Chapter 1. Introduction ^ Plate 1.1: Digitized Data Topology. 5 Chapter 1. Introduction ^ 6 1.3.2 Data Parameterization For a cylindrical data set, the natural choice for an approximating surface is also cylindrical. This simplifies the problem of finding the correspondence between the data and the spline surface 1 . For this particular application, the digitized data is used to specify the shape of a portion of an existing spline figure with its own notion of which part of the surface is what (head, foot, etc.) For example, the animated figures created with hierarchical B-splines are typically a single parametric surface encasing the entire body (This is part of the mechanism used to control the deformation of the surface around the articulations.). The particular portion of this surface dedicated to the head is not a cylinder, but a rectangular sub-region of the surface where the edges of the region correspond to the opening at the neck and the interior of the region that encompasses the rest of the head. In this case, parameterization that maps a data point in R 3 to a point in parametric space (in R2 ), is not straightforward. Chapter 3 describes the parameterization method used to approximate cylindrical data. 1.4 Related Work The literature on surface fitting and approximation is vast. This section will not attempt to survey the field but will briefly examine work on surface approximation for animation. Chapter 2 will also cover specific issues of data parameterization, numerical methods, and parametric surface formulations. The work of Nahas et al. [30] addresses the same problem of using digitized data to create an animated human figure. A raw digitized mesh of data points is used without further processing as the control vertices of a bicubic B-spline surface. For animation, broad-scale changes to the surface shape are accomplished by embedding the control vertices of the model in a coarser spline mesh, created by selecting a relatively few characteristic points from the original B-spline 'The operation that identifies which part of the data set matches which part of the approximate surface is termed parameterization. Chapter 1. Introduction^ 7 model of the surface. When the characteristic points are moved, they displace the vertices of the model itself. The data set used was sparse and no attempt was made to reduce its size. The paper by Schmitt et al. [33] presented an adaptive process to fit surface data with a geometrically continuous collection of rectangular bicubic Bezier patches. Bezier patch basis functions do not automatically ensure the continuity between adjacent patches, thus during the approximation process, extra constraints must be imposed to ensure at least geometric continuity (collinear tangents) between patches. The adaptive behaviour is implemented by fitting a portion of the data with a single patch, testing the fit for errors within a given tolerance, and then subdividing the patch if the tolerance is exceeded. Geometric continuity is controlled by using a constrained least squares approximation method where the constraints are imposed by the continuity conditions. A given number of constraints per patch are required to piece the Bezier patches together to form a continuous composite surface. Because of the properties of the basis functions,a Bezier surface will have more control vertices (9 times) than that of the equivalent B-spline surface. (Note, however, that with Bezier patch, it has considerably easier to control the continuity between patches.) Because each Bezier patch is essentially independent, each patch in the surface can be split into two or more subpatches (refinement) independently. The broad class of tensor-product spline surfaces, e.g., B-spline, Beta-splines and their rational counterparts, which provide continuity without the explicit imposition of constraints, must refine an entire row or column of patches to increase the number of patches in the surface. A hierarchical spline [17] [20] allows local refinement of uniform B-splines or Beta-splines and their rational counterparts (A brief summary of this approach is provided in Section 2.3). In the paper by Forsey and Bartels [19], this property was exploited to provide an adaptive approach to fitting regular (gridded) data in three dimensions, with a hierarchical bicubic Bspline, and generalized to multivariate B-splines in any dimension [18]. This approach utilizes a least squares to fit a coarse spline surface (i.e., on the order of 1-16 patches) of the appropriate topology. Regions that are out of tolerance, i.e., more than a given distance from the data, are refined and the fitting process is repeated for the refined surface. Eventually, regions of Chapter 1. Introduction^ 8 inadequate fit become separated and surrounded by regions within tolerance. These regions are individually accommodated using local refinement and constrained least squares approximation within each separable region. This approach results in a multi-resolution surface formulation in a hierarchical format. The hierarchical approach provides a further contribution by enabling a certain economy of representation for the final composite surface and produces reasonable results. However, the resulting surface suffers from oscillations when presented with data with high frequency components. The oscillations were reduced considerably through the use of non-uniform hierarchical B-splines and an improved parameterization procedure [35]. The surfaces resulting from the above approach were not appropriate for further manipulation via the hierarchy. Because each level in the hierarchical surface is the result of a separate least squares approximation of the data, corresponding parametric locations in each overlay may have widely varying slopes. This plays havoc with the surface when broad-scale changes in shape are made. In contrast to the top-down approach in [33] and [20], Lyche and Morken [26] work bottomup by first approximating the surface with a fine mesh of spline patches and then removing knots in those regions where they will not cause the surface to move out of tolerance. This thesis presents a bottom-up approach to surface approximation using hierarchical splines. 1.5 Overview Chapter 2 covers the mathematical background for B-splines, refinement and mid-point subdivision, hierarchical B-spline, parameterization and multigrid numerical method. Chapter 3 gives in detail the modified chord-length parameterization method used for the test data, and Chapter 4 discusses the construction of the linear equations used for the approximate and their solution using the full multigrid method. Details of how the full multigrid method is used to construct a hierarchical surface is given in Chapter 5 along with one approach to reducing the storage requirement of the approximating surface. Conclusions and future work are presented Chapter 1. Introduction in Chapter 6. ^ 9 Chapter 2 Mathematical Background 2.1 Uniform Bicubic B-Spline Surfaces The use of parametric surfaces originated in the field of CAGD. B-spline tensor-product surfaces permit the most supple movements for the purposes of animation [5]. Within this family of surfaces, we are interested in the uniform case when knot spacing is constant in both the U and V parametric directions. A detailed derivation of the uniform B-spline basis functions and efficient algorithms for their evaluation, along with the reasons for their extensive use in CAGD were presented in the paper by Barsky [5]. The properties presented in this thesis, such as local control and control of the degree of continuity at the joints between surface patches, are very desirable for fitting a spline surface to discrete data [3]. In this section, we will briefly review the B-spline surface formulation and present the notation used in subsequent chapters. Detailed explanations can be found in [3] and [14]. A B-spline surface is defined by a set of control vertices. Figure 2.1 displays a single patch of order 4 defined by 4 x 4 control vertices. Notice that the surface does not necessarily interpolate the control vertices. As shown in Figure 2.2, a bicubic B-spline surface is arranged in a piecewise manner, where each piece is a segment of the surface called a surface patch. Each surface patch is defined by a subset of the control vertices. The number of control vertices required to define a surface patch depends upon the order of the underlying basis functions. The entire surface S(u, v) is a 10 11 Chapter 2. Mathematical Background^ V12 V00 V10 V22 V20 V30 Plate 2.1: A bicubic B-spline surface with a single patch defined by 16 control vertices. Plate 2.2: A B-spline surface is a mosaic of surface patches. Chapter 2. Mathematical Background ^ 12 mosaic of these patches sewn together and has continuity appropriate to the order of the basis functions (i.e., Ck -1 where k is the degree of the B-spline). Tensor-product surfaces have a rectangular topology and we will label each individual patch as Pii. A bicubic B-spline surface consists of patches that are cubic in both the U and V parametric directions. The control vertices form a two-dimensional array although each vertex Vii is a vector in three-dimensional space (i.e., Vii = (x, y, z) E R3 ). The mathematical formulation for a patch Pii is defined as a tensor-product in the following way: 1^ 1 Pii (et,v)= E E vi-F,, +8Br,k(u)Bs,h(v) (u,u) E RUi_i, Vi_i) x (Ui, J (2.1) r=-2 s=-2 where k and h are the polynomial orders of the basis functions (for cubic B-spline, they are both 4.) and the x sign means the span in parametric space and will be explained related to Equation 2.3. The bivariate basis function Bi,k(ui)Bi,h(v i) is the tensor-product of a set of univariate . basis functions, where (Ui, E R 2 is an array of knot positions in parametric space. For uniform spline patches, the knot spacings in both U and V are constant. Each point Pii(u, v) is an average of the 16 control vertices weighted by its underlying basis functions. This is one important property of a B-spline formulation since it means that a control vertex controls the surface shape locally. A single bicubic B-spline surface patch is fully defined by 4 x 4 control vertices and is unaffected by any other control vertex in the surface. Conversely, a given control vertex exerts influence over only 16 surface patches and has no effect on the remaining patches. B-spline also has other great properties such as variation diminishing, convex hull and etc. [14]. The uniform bicubic B-spline formulation is used throughout this thesis. The basis functions are: [B_2(u) B_1(u) Bo(u) Bi(u)] = — 61 [u 3 u 2 u 1 1] —1 3 —3 1 _3 — 0° 33 00 1 4 1 0 (2.2) Chapter 2. Mathematical Background ^ 13 A surface S(u, v) is a smooth composition of all patches. It is defined as: m n S(11 5 V) = E^(u, E RUO, VO) x ( U m i=1 j=1 , V.11 )]• (2.3) Note that the parameter ranges for the U and V directions are the union of each patch's parametric range (Figure 2.3) in Equation 2.1, namely RUo, Vo ) x (U., V.)] = RUo, vo ) x (U1,14) ] u RU1,^x (U2, V2)] u • • • R^Vn -1 ) x ( Um , VO L A one-to-one mapping provides a unique position in parametric space for each surface point S(u, v). However, when we employ only the uniform B-spline formulation, separate parameter ranges for each patch is possible as long as they are consistent. For example, we can have the parameter range (0, 1) in both U and V directions for all the patches. (U0, Vn) ( Uo Vo ) (Um, Vn) (Um, Vo) Plate 2.3: The parametric range for a surface is the union of the parametric range for each patch. Chapter 2. Mathematical Background^ 14 2.2 Subdivision Scheme An important operation that can be performed on spline curve and surfaces is subdivision (also known as refinement) which takes a set of control vertices and basis functions and represents the surface using more control vertices and basis functions. In this section, we will present a brief review of this scheme and discuss a specific but important type of subdivision called midpoint subdivision. Definition 2.1 Given a B-spline representation m n S(L,V) =E E Bi k(u)Bj,h(V)Vi J. , i=1 j=1 The procedure of describing S(u, v) in terms of a larger set of basis functions Ni k and Ni h and , , a larger net of control points Wij is called subdivision, namely M N S(u, V) =E E i=1 j=1 Ni,k(u)Nj(v)Wi ,j. M and N are the number of control vertices in U and V dimension respectively after subdivision. 0 In Figure 2.4, a surface consisting of a single patch is refined into four patches. Note that after subdivision, the patch does not change shape, there are simply more control vertices in the control graph. Subdivision is achieved by inserting new knot positions into the existing knot sequences defining Bi,k(U) and Bj,h(v). The newly generated basis functions are expressed in terms of discrete splines cxj,k(j) [9] in the following way: Bi k(U) = , i j=1 ai,k(i)N j,k(u) ^ (2.4) where Su is the number of knots inserted in the U parametric direction. The subdivision in the V direction is similar. Chapter 2. Mathematical Background^ 15 From Equation 2.4, the net of control points Wii can be derived in the same way: mn Wij =E E ar ,k(i) a r=1 s=1 3 ,i(i)Vr 0 • ^ (2.5) We can express Equation 2.5 in matrix form as: [W] = [ai][V][a2]•^ (2.6) patch subdivision 1 patch 4 patches control graph subdivision Plate 2.4: Subdivision splits a single patch into four patches. As stated before, the most straightforward subdivision splits each parametric range in U and V at its midpoint. This converts each patch determined by the control vertices [V] into four patches (Figure 2.4). In parametric space, the region for the original patch is split into four regions along the midlines as shown in Figure 2.5. For uniform cubic B-splines: {[M1] for regions I and III [al] =^ [M2] for regions II and IV, for regions I and II [a2] = { [Mi l [M2] for regions III and IV, (2.7) Chapter 2. Mathematical Background ^ where 1 [mi] = o g1 0 0 16 / 1 3 1 0 [m2] = o I ^o o^3 1 ) . (2.8) 00 Plate 2.5: Subdivision splits the parametric space at the midpoint. 2.3 Hierarchical B-splines 2.3.1 Practical Consideration of Tensor-Product Surfaces Subdivision for tensor-product surfaces is non-local because the nature of knot insertion is nonlocal. For example, if patch Pii in Figure 2.6a is the place where we want to refine, this can only be done by inserting two knots. One sequence of knot positions is along the U parametric direction { Uik 11 < k < n} and the other sequence of knot positions {Vkj I 1 < k < m} is along the V direction. However, because of the way the basis functions are defined and evaluated, subdivision does not simply split patch P13 into four patches. Each patch in row i and column Chapter 2. Mathematical Background^ 17 j are also split resulting in many more patches and control vertices than desired or necessary 1. In Figure 2.6b, we can see more surface patches are generated than we expected (9 patches versus 3 patches). When modelling large complex surfaces, this property could add hundreds of unnecessary patches. Plate 2.6: Subdivision through non-local knot insertion. 2.3.2 Localized Subdivision—Hierarchical B-splines The hierarchical B-spline [17] is designed to localize the effect of subdivision through the use of overlay surfaces. Overlay surfaces are hierarchically controlled subdivisions. It is important to keep in mind that the W surface (Equation 2.5) is an exact re-representation of its parent V surface. If one of the control vertices of W is moved, then the W surface departs from its V parent, but only in the area influenced by the displaced control vertex. Outside of that region, the two representations are identical. Because of this correspondence, we could just use the V definition as our basic description of the surface, and save only those control vertices of W which 1 For a bi-cubic B-spline one-patch surface, to split the patch into four, three new knots have to be inserted. Chapter 2. Mathematical Background^ 18 define the position of the displaced surface. The child surface defined by only those retained control vertices of W is called an overlay. This portion is regarded as a separate surface to be manipulated. To maintain the continuity between the overlay surface and its parent surface (the V representation), more than one control vertex of W must be kept and these are only those control vertices that do not disturb the continuity of the surface. When editing takes place at one level of the surface, any overlays resting within the edited area are expected to remain embedded in that area. They will follow editing changes only if they can be dynamically tied to that area. This means that their control vertices must move in accordance with the movement of the section of the surface being edited. Hierarchical surface representation meets this requirement by introducing the reference and offset concepts. At level 1, each control vertex Vii of any part of the surface is of the form of: Vlj = Riii + Oij (2.9) where R lij is the reference control vertex and (4 is the offset. Usually, Ri j is derived from the parent surface definition through subdivision, and 0!_ i is the change to the original Wirs vectors resulting from modification of the overlay surface. Thus the surface representation in Equation 2.3 at level 1 can be rewritten as: S 1 (u, v) = EME7Li 7!. 1. BRui)g. (v j)(R I• • +^-)^ 1 J.Ji^ ni (^\ PP12'3J-+\-'727131 1 i M(Iii).Bi(Vi)0! J. z_ai= 1Ein_ =^E7= =n .. (2.10) Here the superscript 1 is used to indicate the parameter associated with a certain hierarchy level 1. Now we can see that each overlay surface is composed of two components. One is the reference surface: SR I (U, V) = EE /3! (no B((v) •.R(2,3 3^3^ (2.11) E E ki! (n i )B(3(v •)(:)( ^JO• (2.12) i=1 j=1 and the other is the offset surface: nzt ni 4(7.1, V) = i=1 j=1 Chapter 2. Mathematical Background ^ 19 Editing changes to the V-defined surface will automatically cause revisions in the W-defined surface by directly changing the R's, and through them revisions are made dynamically to the W's. Editing changes to the W-defined surface are recorded as changes to the O's. Thus any changes in coarser levels will change all of the finer level R's and result in a global change in surface shape. Changes to the offsets at the current level affect the current level of refinement and all of the finer levels. Since the reference information R iii at level 1 is always directly derivable from its parent surface, the crucial information needed to store the upper levels are the non-zero offsets. Repeated substitution of Equation 2.11 into Equation 2.10 results in the recursive formulation of Equation 2.3 with L being the finest level: S(u, v) =E7i.1-21E11BP(Ui)BY(Vi)ViCij E irn=11 E 711 1 /31. (Ui)BI(VAOti +^Ilf'(ui)Bt (vi)0 4' 1i . . (2.13) Thus, a full representation of a hierarchical surface requires only the control nodes of the root level, the knot spacing at each level and the non-zero vector offset values C4i within each level. This offers a considerable amount of storage saving. This is especially suitable for editing a complex shape because the storage for the surface is proportional to the number of edited points. 2.4 Parameterization Definition 2.2 Given a data set {Dij^(xii,yii,zij)1Dii E R 3 , (1,1) < (i, j) < (M, N)} and a parametric position set^= (ujj,vii)1Uii E [(Uoo, Voo) x (Umn, limn)] C R2 }. Find the parametric surface S(u, v) E R3 such that S(ujj, vii) = Dij^ where mn S(Uij, vij) = E E Bi(uii)Bi(vij)Vrs• r=0 s=0 (2.14) Chapter 2. Mathematical Background^ 20 If Equation 2.14 contains all of the M x N data points, surface S is an interpolation of the data D. In order to generate the equations for the surface approximation problem, each data point Did must be associated with a corresponding parametric domain point (uii, vii). Defining this association for a particular data set and a particular surface is referred to as parameterization. If the association is already known, this is referred to as a non-parametric problem. It is called a parametric problem if we are attempting to approximate data from a physical object where domain information is unknown and therefore must be estimated. 2.4.1 Parameterization Schemes Parameterization is a difficult problem, but a number of different algorithms exist which have been proven effective for univariate and bivariate problems [16] [32]. Different parameterizations of the same data set will result in very different surfaces in terms of their shape and continuity. Since the parameterization scheme used in the U and V direction is independent, let us discuss the different parameterization schemes in the context of curve fitting using a single parameter. Given a set of data points D2 = (x2 yi) and a knot sequence {to, ti, , t n } with , increasing order, we fit the data to a curve C(t) = (x(t), y(t)). Without loss of generality, we can assume to = 0, and hi = t i+i where hi forms the knot spacing vector H = (ho, h1 . , h n _1). Uniform Parameterization The easiest way to determine hi is to simply set all hi's to be the same constant value, namely hi = Const, for i = 0, 1, , n — 1. This is called uniform or equidistant parameterization and was proposed by de Boor [12] and Ahlberg [1]. It is sometimes used because the computations are simple since time is not needed to compute each hi. However, the amount of computational time saved is minor, and curves with loops and cusps can occur even though the defining equations would appear to make that impossible. Examples are cited by de Boor [12] and Chapter 2. Mathematical Background^ 21 Epstein [13]. This method is too simplistic to cope with most practical situations, with one exception given by Foley [15], where the uniform parameterization is superior to chord-length parameterization. The overall poor performance of the uniform parameterization is attributed to the fact that it "ignores" the geometry of the data points. If we interpret the parameter t of the curve as a time variable, as time passes from time to to time t ri , the point out the curve from point C(to) to C(t) traces C(t n ). With uniform parameterization, point C(t) spends the same amount of time between any two adjacent data points, irrespective of their relative distances. A good analogy is that of a car driving along the interpolating curve. We have to spend the same amount of time between any two data points. If the distance between two points is large, we must move with a high velocity. If the next two data points are close to each other, we will overshoot since we can not abruptly change our velocity and acceleration (i.e., the first derivative and second derivative respectively in the physical counterparts) [16]. Chord-length Parameterization Given the physical analogy above, it is perhaps more reasonable to adjust the velocity according to the geometric distribution of the data points. One way of achieving this is to have the knot spacing proportional to the distances between data points as calculated in Euclidean space. The equation is as follows: hi =II Di+i — Di (2.15) A parameterization of the form Equation 2.15 is called a chord length parameterization. This is a commonly used parameterization scheme and usually produces better results than uniform knot spacing does, although not in all cases. Foley [15] cites an example where chord length parameterization does not give the shape users have in mind and wild oscillations occur because the data is poorly scaled or the direction of data changes abruptly. Epstein [13] proves that the parameterization based on a vector-norm distance such as chord length will not produce curves with cusps at the data points, thus giving some theoretical advantage to the chord length parameterization over the uniform one. Chapter 2. Mathematical Background^ 22 Centripetal Parameterization Another similar parameterization based on the distance norm was developed by Lee [24] and is called the centripetal model. This parameterization is based on the same paradigm as chordlength parameterization. If we set hi = VtII Did-1 — Di lf,^ (2.16) the resulting motion of a point on the curve will minimize the centripetal force acting on it. So for trajectory based designs where we care about machine wear, this is the right candidate. Affine Invariant Parameterization The uniform, chord length and centripetal methods are each invariant under rotations, translations and equal scaling in the x and y coordinates (consider data fitting in 2-D). However, the chord length and centripetal methods are not scale invariant. This means that when the x and y coordinates of the data are scaled differently, the resulting fit is not a simple scaled version of the original fit. The geometric properties of scale, rotation and translation invariance are important in character font applications [31] Affine invariant chord length and angle parameterizations are designed to meet these kinds of practical use [31]. 2.4.2 Summary Parameterization is a difficult problem that is very application dependent. Usually the input data points decide which parameterization method to choose and there is no general method which is suitable for most practical case. (e.g., the case when scale, data point spacing, and angles between adjacent data points all vary largely.) Norm based parameterizations such as chord length and centripetal parameterization have certain important geometric properties which make them superior to the uniform case. Chord length and centripetal parameterization are satisfactory for most geometric applications. The problem with these two schemes is that they are not scale invariant, that is, if the data is scaled differently in the x and y directions, the interpolation result will be different. Chapter 2. Mathematical Background^ 23 Affine invariant parameterization is a useful property. Uniform parameterization has this property but it does not produce good results when the data point spacings are large and variable. Affine invariant chord length and angle parameterizations generally provide pleasing fitting results. However, both will fail if the data points fall in the one dimension-less hyperplane [35]. 2.5 Multigrid Numerical Method 2.5.1 Iterative Methods on Linear Algebraic Equations The surface construction problem discussed here involves a large system of linear equations. Let the matrix form of our system of equations be: Av = d^ (2.17) where A represents the coefficient matrix, the vector v contains the unknowns and the vector array d contains the known data point triples. The existing methods can be classified into two types, direct methods and iterative methods. Direct methods, such as Gaussian elimination, provide unique and exact solutions, however they are extremely inefficient with computational costs of 0(N 3 ). When the matrix A is sparse, these elimination methods fill up those entries which are originally zero and a lot of storage space is wasted. For many practical applications, iterative methods are used to provide faster solutions. The speed of the algorithm depends on the rate of convergence, and since each iteration costs more or less the same amount of computational time, a rapid convergence rate is crucial. Given an initial guess of u°, at the r th iteration the iterative procedure calculates the new approximation using: Ur+1 = UT ± C CU (2.18) where C is the iteration matrix and c is a constant vector. Many iterative methods have to been developed in the past for numerically solving algebraic equations. Among these, Gauss-Seidel and Jacobi relaxation methods are the two types Chapter 2. Mathematical Background^ 24 commonly used. The Jacobi method is very similar to the Gauss-Seidel method. The only difference is that the Jacobi method keeps using the solution from iteration r for all variables while computing new values for iteration r + 1. Gauss-Seidel makes use of a new value immediately after it is calculated. Here, we give a brief overview of the formulation and function of the Gauss-Seidel relaxation scheme which is employed throughout this thesis. We denote the matrices representing the diagonal, lower-triangular and upper-triangular parts of A by D, L and U respectively, with L and U having zero elements on the diagonal, thus A = L + D + U. For the case of GaussSeidel relaxation, where each entry xi in ur is replaced immediately when 4+ 1 is available, the iteration equation is: (D L)u r + 1 = d Uu r . (2.19) The corresponding matrix C and the vector c can be computed as: C =^+ L) -1 U, c = (D L) -1 b (2.20) where v is the exact solution and u is the approximate solution to equation 2.17. The error vector e = v — u, satisfies the residual function: Ae = r^ (2.21) where r = d — Au and is called the residual. The residual vector contains error components of different frequencies. An efficient way to damp the rapid fluctuations in the residuals is by Gauss-Seidel relaxation scheme. During the first few iterations, Gauss-Seidel has very fast convergence, with residuals rapidly decreasing from one iteration to the next. But once the high frequency components in the error are smoothed out, the convergence rate levels off. This inherent computational sluggishness can be studied from a spatial frequency perspective [37]. A local Fourier analysis [7] shows that the convergence is fast as long as the residuals have high frequencies compared to the scale of the grid spacing. As soon as the residuals are smoothed out, convergence slows down. The overall convergence is asymptotical, and thus slow overall if high accuracy in the solution is required. Chapter 2. Mathematical Background^ 25 One solution to this is through the use of multigrid methods. 2.5.2 Multigrid Philosophy One way to improve a relaxation scheme, at least in its early stages, is to use a good initial guess. A technique for obtaining an improved initial guess is to perform some preliminary iterations of relaxation on a coarse grid and then use interpolation to get the resulting approximation as an initial guess on a finer grid. A coarse grid is chosen because relaxation on a coarse grid is less expensive than on a fine grid since there are fewer unknowns to be updated. Another reason a coarse grid relaxation is faster is that there is a faster convergence rate on a coarse grid. If we examine the smooth components of the residual on a coarse grid, they appear less smooth than on a grid that is twice as fine as shown in Figure 2.7. This suggests that when relaxation stalls, signalling the predominance of smooth error modes, it is advisable to move to a coarser grid, on which those smooth error modes appear more oscillatory [8]. The relaxation will then be more effective. This is the basic idea of multigrid method. 0^1^2^3^5^6^7^8^9^10 II I^13 14 15^6 a) Wave ea a N = 16 grid o b) Wave on a N = 8 grid Plate 2.7: A wave on C2 h (N = 16) is projected onto 12 2h (N = 8). On the coarse grid, the wave seems more oscillatory than on the fine grid. Chapter 2. Mathematical Background ^ 26 2.5.3 Basic Strategies Suppose one has a set of grids 1th°, Ohl, , O h " over the same domain C2 (hi is the grid spacing at the ith grid). The grid spacing ratio is defined as hk +i : hk, which is usually chosen as 2 : 1. The reason for choosing this value as the grid spacing ratio will be explained later. There are two strategies which can be used: • Nested Iteration. This strategy incorporates the idea of using coarse grids to obtain a better initial guess for fine grids. 1. Relax Av = d on level 0 and interpolate to level 1 providing initial guess. 2. Relax Av = d on level 1 and interpolate to level 2 providing initial guess. n. Relax Av = d on level n-1 and interpolate to level n providing initial guess. n+1. Relax Av = d on level n to obtain the final solution. • Coarse Grid Correction. This strategy incorporates the idea of using the residual equation to relax on the error. 1. Relax Av = d on n h to obtain the approximation U h 2. Compute the residual r = d — Auh. 3. Relax residual equation Ae = r on I2 2h to obtain the approximation to the error e 2h 4. Correct the approximation obtained on I22h uh uh + Interpolate(e 2h ). These two strategies will be adopted in our multigrid method discussed below. They cooperate with each other and demonstrate great efficiency. First, we will give the various methods used to compose multigrid method as a whole, which will be called operators. Chapter 2. Mathematical Background^ 27 2.5.4 Multigrid Operators and Parameters Nested iteration requires that the original problem be re-defined on the coarser grid [8]. One problem is that there might still be smooth components in the error when we reach the finest grid due to the process of transferring from a coarser grid to a finer grid. As we will see below, the scheme designed to get rid of this defect is another pathway to multigrid method. There is a rationale for using the coarse grid correction. However, it leaves a lot of questions unanswered. The relaxation method (Gauss-Seidel) used in the coarse grid correction scheme needs to be determined. Different relaxation methods might greatly affect the algorithm's performance. We have to define the restriction operator used to transfer the residuals from Oh to C/ 2 h and the prolongation operator which is used to transfer the error estimate back to C2h . from C2 2 h The restriction operator from a fine grid to a coarse grid is denoted as prolongation operator from a coarse grid to a fine grid gh and the 4h . In the coarse grid correction scheme, the number of pre-relaxations, before restricting the residuals onto the next coarser grid, is vi, and the number of post-relaxations after prolongating the error solution onto a finer grid is v2. Also, the grid spacing ratio p should be chosen according to provide multigrid method optimal performance. In multigrid method, the intralevel processes are usually basic relaxation methods (such as Gauss-Seidel relaxation), the interlevel prolongation processes involve local Lagrange interpolations, and the restriction processes involve local averaging operations. The exact form of these operations are problem-dependent. Fortunately, there are well established standards that have been obtained from numerical experiment. For general analysis of the theory, there are detailed references by Hackbusch [21] and Brandt [7]. Further discussion concerning the choices for the concrete problem will be presented in later chapters. 2.5.5 Multigrid V cycle - Using the intergrid transfer operators and necessary parameters, the two-grid coarse grid correction scheme can be refined into a multigrid method. The idea of residual correction on a coarse grid is applied recursively and the correction process is repeated on successively coarser Chapter 2. Mathematical Background^ 28 grids until a direct solution of the residual equation at the coarsest level is inexpensive. The coarse grid correction scheme has been incorporated into the multigrid V-cycle to provide fast numerical performance. h 2h 4h 8h 16h Plate 2.8: The schedule for grid visiting when the multigrid V-cycle is executed. Each level is identified with its spacing. Multigrid V-cycle starts from the finest level and proceeds to the coarsest level, then returns to the finest level. Figure 2.8 shows the schedule for the visiting order of the grids. Because of the pattern of this diagram, this algorithm is called the V-cycle. The V-cycle scheme has a compact recursive definition which is given by: Algorithm 2.1 Multigrid V-Cycle Operator MGV U h 4- MGV h (U h ,d h ) 1. Relax vi times on A h vh = d h with a given initial guess U h . 2. IF Cl h = coarsest grid, THEN go to step ELSE d 2h 4- ig h (d h - A h u h ) ugh ,_ 0 ugh 4_ m Gv2h ( u 2h , d2h) . 4 Chapter 2. Mathematical Background ^ 3. Correct uh 4- 29 u h + 41 11 2h ,_ m Rh ( u h , dh \^vo times. 4. REPEAT u h ) 5. END. The MR operator is a exact solution operator which relaxes the equations vo times and returns the exact solution at the coarsest level. The V-cycle is just one of a family of multigrid cycling schemes. The entire family is called the v-cycle method and is defined recursively by the following algorithm: Algorithm 2.2 Multigrid v-Cycle u h <— MG V h (U h , d h ) 1. Relax vi times on A h vh = dh with a given initial guess u h . 2. IF Oh = coarsest grid, THEN go to step 4 ELSE d2h ,_ gh(d h — Ahuh) u2h +._ 0 REPEAT u 2h 4_ mG v 2h( u 2h , d2h) v times. 3. Correct uh 4___ uh + p2u2h 4. REPEAT u h 4-- MR h (uh , d h ) vo times. 5. END. In practice, only v = 1 (which gives a V-cycle) and v = 2 are used. Figure 2.9 shows the schedule of grid visiting for v = 2, which is called the W-cycle scheme. Chapter 2. Mathematical Background ^ 30 h 2h 4h 8h Plate 2.9: The schedule for grid visiting of W-cycle scheme. 2.5.6 Full Multigrid V-cycle In Section 2.5.3, we proposed two distinct, basic multigrid strategies; nested iteration and coarse grid correction. In the multigrid V-cycle, the coarse grid correction scheme was employed to smooth out the low frequency error components. Nested iteration involves solving a small scale problem on the coarser grid and transferring the solution back to the finer grid to provide the initial approximation. Clearly, another recursive path leads to the coarsest grid. The algorithm that joins nested iteration with V-cycle is called the full multigrid V cycle. - The visiting schedule is shown in Figure 2.10. h Fine Grid 2h 4h 8h Coarse Grid Plate 2.10: Full multigrid grid visiting schedule. Here we have the compact algorithm: Chapter 2. Mathematical Background^ 31 Algorithm 2.3 Full Multigrid V-Cycle Operator FMGV v h 4-- FMG V h (u h , d h ) I. IF 12h = coarsest grid, THEN go to step 3. ELSE d2h 4_ gh oh — Ah u h ) u 2h 4_ 0 u2h <-- FMGV2h (u2h, d2h). 2. Correct uh ,_ uh + 4ihu2h 3. REPEAT uh 4_ mRh( u h , dh \) vo times. 4. END. FMG method start from the coarsest grid and proceed to the finest grid. Each V-cycle inside a full multigrid iteration is preceded by a smaller V-cycle designed to provide the best possible initial guess. As shown in Section 2.5.7, the extra work done in these preliminary V-cycles is not only inexpensive, but it generally pays for itself by providing better initial guess for the finest grid. Full multigrid method is the complete knot into which the many threads of numerical methods are tied. It is a remarkable synthesis of ideas and techniques which individually have been well known and used for a long time. Taken alone, many of these methods have serious defects (Gauss-Seidel relaxation, coarser grid correction, nested iteration, and etc.). Full multigrid is a technique for integrating them so that they can work together in a way to remove the limitations. The result is a simple but powerful algorithm. 2.5.7 Performance Evaluation Multigrid is a very efficient method for solving a system of linear equations. Both the storage and computational costs are of optimal order. ^ Chapter 2. Mathematical Background^ Storage Cost Consider a two-dimensional grid with M x N (which we will represent by MN for later reference) grid points at the finest level. At each level, two arrays must be stored, the right hand side of Equation 2.17 and the unknowns from the left hand side. So for the storage requirements, the finest grid 1 h needs 2MN units, the next coarser grid C2 2h needs 4 as many, etc. If we accumulate all of the storage units at each level, we have the upper bound on storage which is: 2MN S = 2MN(1 + 2 -2 2 -4 + + 2-2n) < — 8 MN. 1 — 2 -2 3^* (2.22) Since we only store the components necessary for carrying out the relaxation, namely the input and output of the equations, and have no intermediate storage, the storage cost is optimal. Computational Cost Now let us consider the computational cost. It is convenient to measure these costs in terms of work units (WU). A work unit is the cost of performing one relaxation sweep on the finest grid. Equivalently, it is the cost of expressing the fine grid operator. It is customary to neglect the cost of intergrid transfer operations which could only amount to 15-20 percent of the cost of the entire cycle. First consider a multigrid V-cycle with one pre-relaxation and one post-relaxation, namely vl = v2 = 1. Since the coarsest level has a constant number of relaxations and is only visited once, we can simplify by assigning vo = 2. Adding these costs and again using the geometric series for an upper bound gives: ^TMGV ^1 32 = 2 ( 1 + 2- 2 + 2-4 + + 2-2n) < wu 2wu 8 WU ^• ^2-2 (2.23) For a full multigrid V-cycle, we just sum up the costs for each single V cycle. Starting from the finest grid, the computational cost is bounded by 3 WU as just calculated above. The next coarser grid costs 4as much as the finer grid. Then the computational cost for full multigrid V-cycle is: 2 „ ,^ ^32 ) 1+ Z 2 + 2-4 + + 2-2n) < ^ WU = t — 2 -2 ^(1 T WU.^(2.24) 2-2)2 — T FMGV = Chapter 2. Mathematical Background^ 33 Since each sweep at the finest level is O(MN), the total computational cost of full multigrid V-cycle is also O(MN) for a single cycle, which is of optimal order. Convergence Rate In order to determine the final cost for solving a system of equations, we have to know the convergence rate. The total cost is the cost for each cycle multiplied by the number of cycles required to achieve a particular accuracy. The number of cycles depends totally on the convergence rate for a fixed error tolerance. Convergence rate analysis is a complex area [21] [7]. Local mode analysis [21] provides a non-rigorous way to predict the convergence rate which conforms well to experimental results. A key observation is that when vi = v2 = 1, there is an error reduction rate of 1 for each sweep. The relationship between vi, v2 and the error deduction rate r can be expressed as: 1 r = ,^ (vi + v2) 2 . (2.25) Due to the preliminary cycling through coarser grids, only 0(1) V-cycles are needed by the time the algorithm reaches the finest grid. Therefore, the total computational cost of FMG method is still 0(MN) [8]. Chapter 3 Parametric Space Mapping 3.1 Parameterization Tensor-product B-splines (the surface formulation used in this thesis) are rectangular. Figure 3.1 shows such as surface with the boundary edges marked North, South, East and West. A cylindrical mapping on the data of Figure 3.1 is shown in Figure 3.2a. The East and West edges are mapped to the first column of data and the N and S edges are mapped to the last row of data. As stated previously, we wish to approximate a cylindrically digitized object with a B-spline surface that the boundary of the surface corresponds to the neck of the figure (i.e., the last row of data). This goal can not be met by mapping a cylindrical surface onto cylindrical data as shown in Figure 3.2a, the mapping from the surface shown in Figure 3.2b to the head will meet this requirement. After the initial mapping, we also wish to more closely approximate those regions which are articulated, such as the mouth and the eyes. These areas contain the most interesting features and are important for modelling and editing. Finally, we would like the parameterization algorithm to be as fast as possible. This goal is achieved in a three step process. First the data is mapped into an intermediate uniform parametric space ST. Then the mapping function M(ui, vi) = (s2, ti) is defined to map each parametric pair (ui, vj) in UV space into ST space. This mapping function is defined 34 Chapter 3. Parametric Space Mapping^ 35 Plate 3.1: Labels for the tensor-product B-spline surface. a) Cylindrical Surface b) Folded Surface Plate 3.2: Two approaches to approximating the head data: a) wrap a cylindrical surface around the head, b) fold the four edges of a rectangular surface to cover up the head with the opening at the neck. Chapter 3. Parametric Space Mapping^ 36 through Bernstein polynomials having a given knot spacing. Initially, a single patch is created, then the patch is refined and control vertices are displaced according to the chord-length measure. At last, the UV space is uniformly sampled and mapped back to data space through function M(ui, vi) in order to form a set of gridded data points D23 in Equation reffittingdef. 3.2 Initial Parameterization Imagine that we put a square rubber sheet over Victor Hugo's bust and let the edge of the sheet match the opening at the neck. The first row of the data will be mapped to the center of the sheet with the rest of the rows spreading out to the edge (Figure 3.3). last row of data '^r+1 r (q+2, r) Neste. (264, 264) q^, q+1 , q+2 1 origin S Plate 3.3: The initial mapping from data space to parametric space. Each row of data is mapped onto a square in parametric space starting from the center to the edge. A parametric space ST (same as the UV space used before) is introduced here for initial parameterization of the data. The ST space sits under the UV space which is used for the final data parameterization. The ST space acts as the rubber sheet here. The range of the region is Chapter 3. Parametric Space Mapping ^ 37 from 0 to 2M where M is the number of rows of data. In the case of Victor Hugo, there are 264 rows. We assign the lower left corner of S to be the origin, so the center of S is at (264, 264). If we view the mapping by column, each column of data forms a ray emitted from the center of S and ended at the point where it intersects with the boundary of S (Figure 3.3). Each ray spans a constant of angle anglestep which is determined by the number of columns N. First, we calculate the length of the ray segment. The M data points sitting on the ray segment divide it into M sub-segments with equal spacing. Here is the algorithm for calculating the initial parametric value pair (s, t) in parametric space ST(s, t) for each data point Dad E R 2 : Algorithm 3.1 Procedure InitialMap anglestep = 360.0 / N; for column = 1 to N do angle = column * anglestep; for row = 1 to M do if(angle <= 45.0 II angle > 315.0) s = 264.0 + row; t = 264.0 + row * tan(angle); if (angle > 45.0 && angle <=135.0) s = 264.0 + row / tan(angle); t = 264.0 + row; if(angle > 135.0 && angle <= 225.0) s = 264.0 - row; t = 264.0 - row * tan(angle); if(angle > 225.0 && angle <= 315.0) s = 264.0 - row / tan(angle); t = 264.0 - row; Chapter 3. Parametric Space Mapping ^ 38 end. end. The (s, t) pair is not the final parametric value used to form the fitting equations. The ST space is an intermediate space where the mapping deformation occurs. 3.3 Deforming Parametric Space 3.3.1 Deformation Function The deformation takes place in UV space which sits on top of the ST space introduced in Section 3.2. Note that we currently have three different spaces to manipulate, the data space which belongs to R 3 , the intermediate space ST and parametric space UV. last row of data i i+1 i+2 1 (0, 0)^1. U Plate 3.4: Two coordinate frames imposed on the parametric space. The deformation is defined in terms of a tensor-product bivariate Bernstein polynomial. We impose a grid of control points P23 on the square region in ST space (a control point is defined to be the point lying on the control grid and should be distinguished from control vertex which Chapter 3. Parametric Space Mapping ^ 39 is a spline control vertex). They form m 1 lines in the S direction, and n 1 lines in the T direction (Figure 3.5). Their locations are defined to be: Paj ix 528. jx 528 Xo "4-^S+^T. m^n (3.1) where 528 is the length and width of the parametric region. T(V) 11 1^ ^• i+2 S( U) Plate 3.5: A control grid imposed on the parametric space with grid size of m = n = 7. Each control point is undeformed and has the same coordinates in UV space as in ST space. When first defined, each (u, v) of P23 pair is equal to its (s, t) pair. Actually this is also true for any data point in R 2 . The deformation is achieved by moving the control points from their undisplaced grid positions. The deformed position X' in ST of an arbitrary point X = (u, v) in UV space is found by first computing its (s, t) coordinate when there is no deformation. Under this condition our UV coordinate system is identical to our intermediate ST coordinate system, (s, t) = (u, v). Then the vector valued bivariate Bernstein polynomial is evaluated: X =^ 171) (528 —^ i =0^ / (528 — (3.2) 3 =0 where P23 is a vector containing the Cartesian coordinates of the control point in ST space. Chapter 3. Parametric Space Mapping^ 40 3.3.2 Refinement and Local Control The control points here act as the coefficients of a Bernstein polynomial. Since the Bernstein polynomial is the basis function for Bezier curves and surfaces, there are meaningful relationships between the deformation and the placement of the control points. The parametric region in ST space can be viewed as a Bezier surface patch defined by these control points. The order of the basis function is defined by m and n. Whenever we move a control point, we are reshaping the Bezier surface by moving a control vertex. A Bezier surface has the property of interpolating the boundary points. This is important because we can constrain the boundary of our control grid to the boundary of our parametric region. Consequently, the movement of interior control points or boundary control points will only effect the parametric mapping inside the parametric region. Otherwise, the boundary of our parametric region will not be stable and would not provide a constant size parametric space for mapping the data. For a B-spline formulation, it is known that the surface does not interpolate the boundary control vertices. Compared with the B-spline basis functions, the Bezier surface's Bernstein basis functions save the computational cost of calculating the boundary of a surface patch and constrain the surface points by interpolating the boundary control vertices automatically. Bilinear interpolation is adequate for our particular application which means both m and n Equation 3.2 are chosen as 1. In order to form the chord length parameterization, the control grid is refined and the control points are displaced. The following section provides the criteria for grid refinement and control point displacement. 3.3.3 Deformation Discrete Gradient To construct a chord-length parameterization, the first step is to calculate how much variation there is in the digitized data. ^ Chapter 3. Parametric Space Mapping^ 41 Definition 3.1 The discrete gradient with respect to i in parametric direction gi at data point Did is defined as gi = Dii —^2 < i < M; 1 < j < N ^ (3.3) In the same way, gi is defined as =^—^1 < < M, 2 < j <N ^ (3.4) We define the sum of gi + gi as the discrete gradient at a data point Dii. ^ Without loss of generality, let us consider the i index direction. Recall that the magnitude of point Dii is the radial distance from the central axis to the surface. Given the distance z between two adjacent sample points, then the distance between two adjacent data points Di_Li and Dig is d = I Z 2 + g2 1 . Because of the fine sampling grid from the digitizer, the magnitude of z is constant and is usually very small compared to the magnitude of the non-zero gi's. So we can approximate the above formula as d = gi. This implies that if we direct more parametric lines to those high gradient areas, the data points that are farther apart will also be farther apart in parametric space. This is the general idea of chord length parameterization. If the refinement of the control grid over parametric space is fine enough, more control points will be generated and displaced. Thus we can achieve the chord length parameterization within a certain tolerance. Refinement and Free Control Points The deformation is built starting from a single bilinear patch defined by four control points 410, Pd1 and Pill making up the control quadrilateral. These four points are not displaced since they bound the region being approximated. Thus all points P ik have the same (s, t) pairs in ST space as the corresponding (u, v) pairs in UV space. Each refinement locally subdivides the 1-level control quadrilateral and generates four smaller control quadrilaterals residing at level 2. At level 1, the criterion to decide whether 1 Since the sample points are really fine, we neglect that the surface is curved and take it as fiat. Chapter 3. Parametric Space Mapping ^ 42 or not to refine a control quadrilateral is the discrete gradient sum inside it. First, the overall discrete gradient sum SG for all of the data points which cover the whole parametric area E are calculated. Second, the tolerance threshold T for each control quadrilateral with area A is set as A T = SG — E (3.5) . Whenever the gradient sum inside a control quadrilateral exceeds its tolerance, it is refined; otherwise, it is only refined if the user interactively refines it by selecting the region and then clicking on the refine button. The refinement is performed in UV space. The mid-point subdivision technique is adopted for simplicity. As shown in Figure 3.6 in UV space, the original control quadrilateral is evenly split into four equal sub-control quadrilaterals. e2 tab^it it it ^Ft (11 b it II. 10 r?.(8.) U Plate 3.6: The control grid in UV space is uniformly refined from a 2 x 2 grid into a 3 x 3 grid. The control points inside brackets are the indices of the 4 control points as they are indexed at level 1. ^ Chapter 3. Parametric Space Mapping ^ 43 So we have the default coordinates for newly generated control points as ( V ) Poo^'1)0 U ) p2^ i 01^ P00 u ( ) Di ^Poe^'01 ^02 ^2 - 1 2 ^ P10^ Pi 00 U —^ ) pl P11 11^ 00 Po i U p2^ 12^) Poi 01 ^P )p2 = (v ) . 1 ' 20^'10 ) p2 _ ( ( U — ^± 1 . \ ' 21^ 'P1 10 vu = ( V^ II ) ) p^ 2^ p 1 22^'' 11 The (s, t) pairs for these newly generated control points can be derived from Equation 3.2 by substituting the (u, v) coordinates calculated above. As shown in Figure 3.6, after refinement there are 9 control points altogether at the second level. The four corner 134's stay at the same positions as the four Pi's and are not be repositioned any more. This leaves 5 new control points that can be moved to deform the parametric space. We will refer these control points as free control points. Free Control Point Displacement We will now discuss how to deform the parametric space by moving the free control point. As shown in Figure 3.6, we have five control points which can be displaced, namely and P. pro , Pot , ph, PI Chapter 3. Parametric Space Mapping ^ 44 To displace a free control point, for example Pro , the first thing is to calculate the discrete gradient sum along the chord Pj 0 P1-0 • This is done by sampling along the chord Pclo Pilo and obtaining the discrete gradient value at each point and then summing them up. The sampling rate is constant. This means that the number of points sampled (Ms) is proportional to the length of the chord Pilo PI-0 in ST space. Each sampled point (s, t) is mapped back into discrete data space. This can be assured to be a one-to-one mapping by "inversely" applying Algorithm 3.1 since this algorithm is one-to-one. We sum up the discrete gradient values at each sampled data point Si by: s Ms E gs (3.6) , to obtain the sum S. As shown in Figure 3.7, each sampled point Si will fall into a quadrilateral formed by four data points Di , Di+i Di+i,j+i, Dij±i. There are three cases for this point to be located in the quadrilateral: • At a data point (a corner). • On an edge. • Inside the quadrilateral. In the first case, we take the discrete gradient at the data point. For the other two cases, we assume that the data space is continuous and use bilinear interpolation to calculate the discrete gradient at the sampled position. For example, we get Di_ f_si j + 63 when we map Si to data space. Consider the second case and suppose 6 j = O. Then the gradient will be: . (1 — Si)gii +^0 < Si < 1.^ (3.7) In the third case, when neither Si nor S j is zero, the gradient value should be calculated as: (1 — Si)(1 — j)gii + (1 — 5i)Sjg2,j + 1 + Si(1 — j)gi + Li + bib j^(3.8) After the gradient sum is found along a chord, the gradient midpoint 0 at which the discrete gradient So = 1S is found. Then control point pro is moved to the position 0 Chapter 3. Parametric Space Mapping^ Dij a) Case 1: The sample point is one of the corner points. Di j + 1 Dijj b) Case 2: The sample point is on an edge of the data quadrilateral. Dij c) Case 3: The sample point is in the data quadrilateral. Plate 3.7: The three cases of the sample point falling into a data quadrilateral. 45 Chapter 3. Parametric Space Mapping ^ 46 to balance the discrete gradient sum on either side of Ph. Note the this displacement of PP happens in ST space. Figure 3.8 shows that after Ph is moved toward Plo , the subquadrilateral P4Plo Pli Ph will cover less discrete data points in data space than the adjacent sub-quadrilateral no Phili ni . In UV space, these two quadrilaterals cover the same amount of parametric space. This means we have more parametric space to cover the high gradient areas in data space. Plate 3.8: A control point is displaced from its original position. The other three on edge control points ill , Ph, PA are displaced in the same way to balance - the gradient sum along that edge on which it sits. Control point Ph is displaced based on the same idea of discrete gradient sum balance. An iterative process is designed to compute the approximate point position B = (s, t) such that the discrete gradient sums inside each of the four sub-quadrilaterals surrounding B are very close to each other. The area discrete grid sum is calculated by expanding Equation 3.6 from one dimension to two dimensions: Ms Ns Sarea = EE i=1 j=1 95, 3 (3.9) Chapter 3. Parametric Space Mapping^ 47 where Ms and Ng are the number of sample points in the U and V parametric directions respectively, gs,, is the discrete gradient value at sample point Sij and Sarea is the total discrete gradient sum of the area. If we refine the grid of control points enough times, a near chord-length parameterization is obtained for those areas with a fine control grid. Color plate 2 shows the control grid after five levels of refinement. 3.4 Sampling Data Points The FMG method requires gridded data, but because of the odd configuration of the spline surface, the image of the gridded data in the parametric space is anything but gridded. To obtain gridded data points for the right hand side of Equation 2.14, the parametric domain of the spline (UV space) is sampled on a grid and mapped back onto data space through the intermediate ST space. There are X and Y sampling points in the U and V directions respectively (Figure 3.9). In practice the values of X and Y depend on the density of the original data. The step size was chosen as 257 in both U and V to provide a sufficiently dense sampling to catch most of he details of the surface and to provide the ratio of 2 : 1 for full multigrid methods. Thus we have X = Y = 257 and ustep = vstep = _528 264' We now have the parameter set {(ujj, vij) = (i x ustep, j x vstep)10 < i < X, 0 < j < Y} to find the 3-D data position associated with each (ujj, vij) pair. Each (uij, vij) is mapped onto the ST space to obtain the corresponding (sij, tij) pair using Equation 3.2. This is a one-to-one mapping. In the same way that we calculated the discrete gradient, the (sij, tij) pair is mapped into data space and the real pair (i + Si, j +6j) is returned. Here i and j are the indices for an existing data point Dij. Bilinear interpolation (Equation 3.7 and Equation 3.8) was used to calculate the exact value at Di+oi ,j+6 3 = Zqr from the four existing data points Dij, Di+Lj and Di + i,j + 1. This 3-D coordinate triple (xij,yij, zij) associated with parameter pair (ujj, vij) was used to generate the system of linear equations in Chapter 3. Parametric Space Mapping^ 48 V^X samples N\ Y samples ustep U Plate 3.9: The UV space is evenly sampled for produce the parametric pair set T. Equation 2.14. With respect to the original data set, the parameterization is near chord-length. The new data set also approximates the original surface with approximately the same number of points. Feeding these data points into the FMG solver, we can solve the equations for the control vertices which define the parametric B-spline surface. This will be discussed in the following chapter. 3.5 Summary A parametric space deformation technique based on image warping was introduced. The technique deforms the isoparametric lines to follow rough regions in the raw data. The deformation was controlled by control points that form a control grid in parametric space that can be refined and displaced. The gradient distribution criterion used to move a control point generates a near chord-length parameterization. The parametric space after deformation was sampled evenly and then mapped onto the data space. Bilinear interpolation is used to obtain the exact position of the data point associated with parameter pair. Chapter 4 Multi-resolution Surface Fitting Using Multigrid 4.1 Deriving Equations and Constraints The remaining problem is to solve the set of linear equations to obtain a surface that interpolates the gridded data points {Z12 E R310 < i _< X, 0 < j < Y} The control vertices V are the unknowns and the data points Z are knowns. For each data point Zii, the interpolation condition can be expressed as an equation in terms of the unknown control vertices. The number of equations should be equal to the number of data points. 4.1.1 Data Point Distribution The data points to be interpolated must now be distributed among the patches to make a correspondence from a data point to a surface point. As shown in Figure 4.1, if we sample at the four corners of each patch, X x Y patches provide enough freedom for interpolating the samplings. For bicubic B-spline surface, X x Y patches are defined by (X + 3) x (Y + 3) control vertices. Therefore, additional 2X +2Y equations have to be derived based on boundary conditions to determine the unique solution of Equation 2.14. The additional equations will be covered in Section 4.1.3. 49 Chapter 4. Multi-resolution Surface Fitting Using Multigrid ^ 50 Zinn zon Pmn 202 . . D 112 Z11 r22 222 zoi / Ili / P21 z21/ zoo ^ Zl o ^ zino Plate 4.1: The (X + 1) x (Y + 1) data points are mapped onto X x Y patches. Every patch corner corresponds to a data point. Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 51 4.1.2 Deriving the Positional Interpolation Equations For each data point, the corresponding parametric location in the patch is (1,1) (Figure 4.2). Plate 4.2: A data point Zij is mapped onto the upper right corner of a single patch. At the same time, it is the corner point of three adjacent patches. Substitute u = 1 into Equation 2.2, we have: [ B_2(1) B_1(1) Bo(1) Bi(1) = [ o Use this precalculated basis function value in Equation 2.1, we have: Pij (1, 1) =^ + vi_i, j +1+ 414,1 + +4vi, j +1+ vi+1, _ i +4vi+1,.; +^ ; i = 1, , m — 1 and j = 1, ,n — (4.1) 1. According to the mapping from data points to surface points, Pij(1,1) = Zii. We can express Equation 4.1 in matrix form as: Av = d^ (4.2) Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 52 where v is the vector containing the control vertices from V11 to Vx_iy_i and d is the vector containing the data points transformed by certain control vertices (certain unknown control vertices are moved to right hand side to make the matrix A banded). The corresponding matrix A is / 16 4 4 16 4 4 4 1 1 4 4 1 16 4 4 16 16 4 1 1 4 1 1 4 4 16 1 4 1 4 16 4 4 16 4 4 •/ Note this matrix is symmetric, positive definite and thus it is well-conditioned. It is also diagonally strong which makes it possible for relaxation scheme to converge. The matrix is sparse which is the key condition for multigrid solver to perform well. The above equation is a set of (X — 1) x (Y — 1) equations, one for each data point. These equations are expressed in (X + 1) x (Y + 1) unknown control vertices. 4.1.3 Deriving the Boundary Conditions In order to completely determine the unknown control vertices, an additional 2X +2Y equations are required. Each of these equations corresponds to an extra condition imposed on each of the 2X + 2Y boundary data points. The boundary data points are: Zoj j = O . , n — 1; Zi n i = 0, , M — 1; Zin d j = 1, . . . , n; , Z, 0^i = 1, ... ,7n. 7 The remaining degrees of freedom are determined by specifying the tangents and curvature at the edge of the surface which interpolates the boundary data points. At each boundary data point, the curvature is set to be zero and the surface tangent to the control graph. Chapter 4. Multi-resolution Surface Fitting Using Multigrid ^ 53 The following equations are derived for each corner data point: Voo = ZOO VOn = ZOY VXO = ZmO VXY = ZXY ^ (4.3) Other boundary data points fall into the following equations: Zoj = ( 4 Vo.; + Vo,j—i + Vo,j+i) Zio = t-(4Vio + Vi-1,0 + Vi+i,o) Ziy = 1 (4ViY + Vi-1,Y Vi+1,Y) Zyj = (4Vxj VX,j-1 VX,j+1) forl<j<Y-1 forl<i<X-1 forl<i<X-1 forl<j<Y-1 (4.4) The corresponding matrix for these equations is: /4 1 1 4 1 1 4 1 1 4 1 1 4/ Up until now, for the 265 x 265 data points we derived by sampling the parametric space, we had 265 x 265 linear equations, one for each data point. These 265 x 265 linear equations involve 265 x 265 unknown B-spline control vertices which define the final approximating surface. In the following section, we will discuss the solution of these equations. 4.2 Solving For Control Vertices 4.2.1 Full Multigrid Method The system of linear equations defining the approximation is solved using the full multigrid method. The benefits are two-fold: • The excellent numerical behaviour of FMG provides a fast and stable way of solving these equations. • FMG method provides a multi-resolution surface which can be used to directly create the levels in a hierarchical spline surface. The FMG method was described in Section 2.5.7 showing that it is an excellent candidate for a well conditioned sparse matrix on a uniform grid. Chapter 4. Multi-resolution Surface Fitting Using Multigrid ^ 54 4.2.2 Full Multigrid Operators In the following several sections, we will discuss the relaxation schemes, coarsest level solver, intergrid operators and the interpolation operator used in our implementation. Relaxation Operator From the matrix corresponding to Equation 4.1, we know that it is a nine point formula 1 on the uniform data point grid. The general practice has been to use the Red-Black Gauss- Seidel relaxation scheme for the five point formula. This leads us to choose the four colour Gauss-Seidel relaxation scheme for the nine point formula because one quarter of the residuals calculated after post-relaxation are zero [21]. v o o c o ,, ^ a a 0 0^o D ^ • o 0 0 ^ a 0 c o 0 ^ a 0 0 '^^ • C 0 0 a ^ a 0 ^ ^ o^0^0^ 0 0^0 U Plate 4.3: A illustration of the four-colour relaxation scheme on a uniform grid. First the 0 points are swept, then the Q points, A points and finally ^ points. The four colour relaxation scheme is similar to the Red-Black Gauss-Seidel relaxation, shown in Figure 4.3. We divide the whole set of data points on the uniform grid into four sub-sets and relax them in turn. The four sub-sets of data points are represented with O, 0, ^ and A 1 If we take the control vertices as unknowns, we call an equation composing a matrix a nine point formula if each linear equation contains nine unknowns, which are adjacent to each other on a grid. Chapter 4. Multi-resolution Surface Fitting Using Multigrid ^ 55 respectively. First all the 0 points (which contain every other original data point on the grid in both the u and v parametric directions.) are relaxed (Gauss-Seidel relaxation) one by one. Then the Q points, ^ points and A points are relaxed in the same way in turn. From Figure 4.3, we can see each 0 point being relaxed is only calculated from non-0 data points. Thus before all of the 0 points are totally relaxed, they do not refer to each other. If one has a multi-processor machine, it is possible to parallelize the relaxation process which greatly speeds up the whole multigrid cycle since most of the computation cost is spent on relaxation. Another relaxation scheme implemented is line Gauss-Seidel relaxation based on the following analogy. In parametric space, these gridded control vertices in both u and v direction have very clear geometric meaning. In the u direction, each "line" of control vertices defines a B-spline curve residing inside the surface which is a blending of curves along u isoparametric direction with the curves along the v isoparametric direction. The surface is defined by the whole set of control vertices which can also be considered to be a blending. The surface is called a tensor product of those single curves. The line Gauss-Seidel relaxation method relaxes the u direction unknowns and then the v direction unknowns "line" by "line". Because the unknowns are coupled closely in both u and v directions, this method is very effective in reducing the residuals. Coarsest Level Operator The coarsest level operator MR (Algorithm 2.1) is Gauss-Seidel relaxation. The explicit equation for relaxation is: vi 1 =— l ii di — E k<i +1 - E likvi (4.5) k>i Experiments indicate that to reduce the residual down to be 0.001% of the exact solution, only 11 iterations are required at the coarsest level. So the convergence rate is very fast. ^ ^16 Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 56 Prolongation Operator The most common prolongation method is interpolation. A common prolongation operator is bilinear interpolation. (Others may be found in [21]) The operator is expressed as: which is called the stencil form of the operator. Stencils representing maps of 12 2h into^use reverse square brackets. A stencil defined on the fine grid is centered at a point present in coarse grid. Each entry indicates the contribution of the coarse grid point to the three (or nine) corresponding fine grid points. Specifically the central entry is the coefficient of v1 ,2j in the formula for Restriction Operator There are two types of restriction operators that are commonly used, injection and full weighting. (In 2-D, injection is defined by vg = v1 2i .) In 2-D the full-weighting operator. In 2-D, it is expressed as _^(„,h ^„,h.^„,h^„,h ^ k"2z-1-1,23-1-1^"2i-12j+1^'2i-1,2j-1+ 2 (v 2i 2j+1 ,h^,h^„h "2i,2j-1^u2i+1,2j -F "2i-12j+1)^4V 2hi,2j) . (4.6) The stencil forms of these operators are as follows: 0 0 0 0 1 0 0 0 0 and respectively. Thus, for injection, the value at the coarse grid points is just the value at the corresponding fine grid point. Full-weighting uses a weighted average of the values at the corresponding fine Chapter 4. Multi-resolution Surface Fitting Using Multigrid ^ 57 grid point and its nearest neighbours. In our implementation, the full-weighting operator was chosen. The stencils above are for an interior point which has all non-empty neighbouring points. At a boundary point which is not a corner point, there are only 5 neighbours instead of 8, so the stencil in this case is: o 1 ci t2 0 For a corner point, there are only three neighbours, thus the stencil is: s 0 0 -§- [. 0 0 0 Interpolation Operator The interpolation operator takes a solution at the coarser level and interpolates it to a finer level between adjacent V-cycles. As suggested by Brandt [7], bicubic interpolation is a better candidate for interpolation operator than bilinear interpolation. Because it introduces less high frequency errors in the initial guess at a finer level than bilinear interpolation does. Thus, a better initial guess can be obtained at a finer level from the interpolated result. For an interior grid point (i, j) with grid spacing h where i is odd and j is even, the points (i + 1, j) and (i + 3, j) are coarser grid points. The interpolated value for (i, j) is: v i( i,j ). 1 vi-i(i _ 37 3, j) 9 vi--i(i j7 - + 17 9 vi-i(i - 16 1 vi-i(i + 3, j) (4.7) We can repeat the interpolation process with respect to the v direction at (i, j) with j odd. Applying this interpolation to a vector that is 1 at some interior point and vanishes elsewhere, we obtain the matrix [21]: / 1 0 9 - 162 - 1 0 W 1^9 0 - 169_ 7^16 - 1720 0^0^0 0 181^9^81 96 2^162^1962 1 0 71^ 9^1861 0 7 1 17 TV 0 0^0^0 9 _ 1 _ 9 0 _ 7.7^ 16^167 0 0 0 0 0 0 1 177 0 9 - 16 2 —10 ^ Chapter 4. Multi-resolution Surface Fitting Using Multigrid ^ 58 For a boundary grid point, Equation 4.7 is not be applied since one or two terms are missing. It is more convenient to use the one-sided quadratic interpolation in the interior domain which yields 171 (ih, jh) = 47 1-1 (ih — 3h, jh) + 713 V I-1 (ih — h, jh) +17 1-1 (ih + h, jh) which leads to matrix: ^ (4.8) 3 _ 1^3 0 1 / 0 0^64^8^32^64 0 0 0^0^0^0 0 ^3^9^3 0 0 92^ 0 - 32 4^16 1 6^I 0 0^ t 0 - -Q 3^719^' 0 0 64^8^32^T 0 - '3Z 0 0 0^0^0^0 0 0 0 0^0^0^0 0 Suppose we have X x Y solutions after we complete the V-cycle for the coarser level. The interpolation operator will generate (2X — 1) x (2Y — 1) "control vertices" as an initial guess for upper level from the X x Y solutions. Using interpolation, the geometric shape of the surface defined by the X x Y control vertices at the coarser level is different than the surface defined by the (2X — 1) x (2Y — 1) control vertices. This suggests that a more suitable alternative for interpolation is B-spline surface subdivision. (see Chapter 2) Midpoint subdivision will generates (2X — 1) x (2Y — 1) control vertices from X x Y ones. The central (2X — 3) x (2Y — 3) control vertices define a surface which consists of as four times many patches as the surface defined by the X x Y control vertices. Moreover, both surfaces have exactly the same shape. Since we expect that different levels of the surface are converging to the final approximation shape at the finest level, we prefer using the exact shape generated at coarser level as the initial guess at finer level to using the shape interpolated to finer level by bicubic interpolation (the latter one could add high frequency components). The user will not experience unusual changes in shape which can be introduced by the bicubic interpolation. Now we give a general analysis of the computation cost of both bicubic interpolation and mid-point subdivision. For each new control vertex created by bicubic interpolation, nine neighbouring points are involved. For the subdivision scheme, an even indexed control vertex is calculated from two neighbouring data points and an odd indexed control vertex is calculated Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 59 from three neighbouring points. Based on the computation cost of each generated control vertex, subdivision is less expensive. However, bicubic interpolation makes use of all the existing coarser level points and transfers them directly to finer level in constant time without further calculation. Thus bilinear interpolation only needs to compute 4of the total control points for the finer level, decreasing computation cost. The disadvantages and advantages of this methodology even out in practice. For example, when implementing with a 265 x 265 grid, the interpolation process costs almost the same amount of time for both methods. 4.3 Performance Evaluation For the Victor Hugo data set, the FMG method was used to solve a system of 257 x 257 linear equations using a nine level grid scheme scheme. The first level has 2 x 2 data points. At the finest level, all of the sampled data points Zip are interpolated by a surface consisting of 256 x 256 patches. The final result must be a multi-resolution surface definition, thus the approximation result at each level must be saved before it is interpolated to finer level. This means that for each level, we need not only a matrix [Vij] to store the unknowns, but also a backup matrix [VBAKii] to store the fit result after each smaller v-cycle. This backup process is done each time that the fitting at a certain level is done. This means of a cost of V ( 2 5 7 x 257) storage units. Still, it is 0(N) cost given the number of data points to be interpolated. Running on a Silicon Graphics Crimson, solutions to all levels of equations are found in approximately 30 seconds. According to Hackbusch [21], if we set vi = v2 = 2 and -y = 2 (i.e. use W-cycle instead of V-cycle), the residual r of each iteration is reduced by a factor of (vi + v2) 2 = 16. As shown in the following table, this criteria is met. Table 4.1 also shows the time cost of each iteration at every level. Here time t is in hundredths of seconds, and L and I indicate the level index (the smaller, the coarser) and iteration number respectively. The residual tolerance is set to be 0.0001, which is very small compared to the general magnitude of a control vertex which was 200. Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 60 L 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 r I 11 0.0000013258 1 0.4906180920 2 0.0276220148 3 0.0012658379 4 0.0000636812 5 0.0000032423 1 0.5924855008 2 0.0340934001 3 0.0018144016 4 0.0001093170 5 0.0000067578 1 0.6936782689 2 0.0331101110 3 0.0018249213 4 0.0001103976 5 0.0000068922 1 0.5606160961 2 0.0383129169 3 0.0019832604 4 0.0001116628 5 0.0000069643 t 15 22 23 23 23 23 98 97 99 97 97 415 417 416 416 416 1705 1687 1704 1695 1704 Table 4.1: The computation cost at each level in FMG. Chapter 4. Multi-resolution Surface Fitting Using Multigrid ^ 61 Since FMG pursues a direct solution to the linear equations instead of an approximate solution like least squares method, noise components in the original data set do not effect the efficiency. However, the algorithm does not have the ability to remove or reduce the effect of noise on the approximation result. As for the different relaxation and interpolation techniques used inside the FMG cycle, a brief analysis will be provided for comparison. At the finest level, the time used for line GaussSeidel is almost the same as for four-colour Gauss-Seidel relaxation methods, which is about 2 seconds for each sweep. Then line Gauss-Seidel method does not smooth the error components as fast as the four-colour Gauss-Seidel. When we used line Gauss-Seidel, more sweeps were required to meet the same residual tolerance. The reason for this phenomena requires more investigation. As expected, the subdivision interpolation method performs better than the bicubic interpolation method because the former one produces a smaller residual at the finer level. Chapter 5 Hierarchy Building The FMG approach produces approximation to the data having multiple resolutions. The next step is to generate hierarchical B-spline from the results. Lyche and Morken [26] work bottom up by first approximating the surface with a fine mesh of spline patches then removing knots in those regions where knot removal will not cause the surface to move out of tolerance. In the following sections, we are going to present a bottom-up approach for generating the hierarchical representation of the surface that we obtained from the FMG solver. 5.1 Inverse Subdivision 5.1.1 Hierarchical Surface File Format The hierarchical B-spline surface is generated by generating a file defining a hierarchical B-spline surface from the multi-resolution representation. Recall from Chapter 2 that a hierarchical Bspline surface is very economic in terms of space, we only need to store those non-zero offset vectors and the first level control vertices as reference vectors. The following script file is a typical hierarchical B-spline surface file generated by the Dragon 1 [17]. Hierarchical B—Spline Save file: Version 9 ' Dragon is a Hierarchical B-spline surface editor developed by Dr. David Forsey at University of Waterloo. 62 Chapter 5. Hierarchy Building^ Surface: 'Plane 0 Root Frame: 0 'LinkO Number of Frames: 1 Number of Offsets: 35 Topology: UO VO Symmetry: 1.500000 V Euler Angle Orientation: 0.000000 0.000000 0.000000 Position: 0.000000 0.000000 0.000000 Offset Data Level 0 - spacing 1 4x4 UO 0:F 0:-60 0 -60 1:F 0:-20 0 -60 2:F 0:20 0 -60 3:F 0:60 0 -60 U1 0:F 0:-60 0 -20 1:F 0:-20 -25.5 -20 2:F 0:20 -25.5 -20 3:F 0:60 0 -20 U2 0:F 0:-60 0 20 1:F 0:-20 -25.5 20 2:F 0:20 -25.5 20 3:F 0:60 0 20 U3 0:F 0:-60 0 60 1:F 0:-20 0 60 2:F 0:20 0 60 3:F 0:60 0 60 Level 4 - spacing 0.0625 U24 20:T :29.0296 23.0948 47.4123 63 Chapter 5. Hierarchy Building^ 64 28:T :29.0296 —23.0948 47.4123 This file contains the surface definition as well as a definition for the coordinate frames for an articulated figure. The header file contains the statistical data about the surface such as the number of frames and number of non-zero offsets. 5.1.2 Inverse Subdivision Scheme This section describes an approach that greatly reduces the number of offsets, but does not produce a hierarchical B-spline surface that is animatable. This scheme is designed to minimize the number of non-zero offsets. We begin from the finest level / in FMG. At this level, we have the entire definition of the surface and the finest surface mesh. The next lower level / — 1 is generated in such a way that if we refine level / — 1 and generate the finest level by offset reference, we have s minimum number 2 of non-zero offsets. This idea can be described using a spline curve. Suppose we have a cubic B-spline curve defined by 11 control vertices (Figure 5.1a) containing 8 curve segments, and we let this curve be at level / in the hierarchy. Imagine a spline curve obtained by subdividing a 4-segment B-spline curve (Figure 5.1b) with the same order, and adjusting the control vertices by offset referencing. This 4-segment curve is at level / — 1. Calculating the 7 control vertices defining the 42 This minimization is not proven in theory, it is a minimum number easily obtained. ^ Chapter 5. Hierarchy Building^ 1 V4 ^I • 65 I V 11 ^V10^ • V9^ • 1 V5 • • V6 •I V7 •1 Vs a) A curve defined by 11 control vertices and having 8 segments. 1-1 V7 1-1^• V6 1- 1 1 1-1V31-^ V4 V2 1-1^•• Vi^• • ---------------1- • • -t I -------------- • 1-4 V5 b) A curve obtained by inversely subdividing the curve in a). Plate 5.1: The inverse subdivision process generate a lower level curve from a upper level curve in a way that reduces the number of non-zero offsets while the hierarchy is built using subdivision. ^ ^ ^ Chapter 5. Hierarchy Building^ 66 segment curve is an overdetermined problem. The equations look like this: o \ / 0 2 V31 14-1 0 0 ,1 /- 1 V3 1-1 V4 0 0 0 V5 0, 1 v51-1 (5.1) VC "/ /-1 V6 0 1 0 0 0 0 0 2 vl "/ "/ u9 \ V 71 -1 V 10 / If a unique solution is desired, we can reduce the matrix to the following 7 x 7 matrix: /^1^1 7^3 g 0 81 0^0 0^0 0^0 \0^0 0 0 0 0 0\ i/vi1-1 \ ,/1 1^ 0 0 0 0^u2 3 1^,1-1 R 0 0 0^u3 f 3 1^ 0 0^v4/-1 (5.2) R 1 § 1 0 0 g^ ^V5i—i 1g 1^/-1 0 0 g V6 9: 0 0 0 21 7^ 1 \ v 17-1 vl V10 I V 11 / This matrix is well-conditioned and thus we know that the solution is unique. The resulting seven control vertices define a curve at level 1 with the property that when we refine this curve, seven out of the resulting eleven control vertices are identical to the corresponding ones at level 1, namely vi, 74, v4 , vs , 43 , vim , vi r We term this technique inverse subdivision. So if we use this idea to build the hierarchy in our bottom-up approach and generate each lower level using inverse subdivision, more than half of the offsets at each level will be guaranteed to be zero. This is a big bonus towards an economical storage cost for our final hierarchical surface. The bottom-up hierarchy building process is presented in the following algorithm: Algorithm 5.1 Inverse Subdivision for level = highestLevel to lowestLevel+1 do InverseSubdivide(level); Chapter 5. Hierarchy Building^ 67 end. SaveBaseLevel(lowestLevel); for level = lowestLevel to finestLevel-1 do SubdivideLevel(level); foreach CV at level+1 if(Offset(CV) != 0) SaveNonzeroOffset(Offset(CV)); end. end. 5.1.3 Remarks: Solution to Triangular Linear Equation A simple symbolic manipulation proves that the inverse of the subdivision matrix is filled everywhere and all the entries change when the dimension changes. Exploiting the banded structure of the matrix, we can reduce it to a unit upper triangular form having an upper triangle consisting of all zeros excepts for the superdiagonal. This can be done by elementary row operations. That is, the system is transformed to: /1)310000 01)32000 0 0 1 )33 0 0 0001,340 00001)35 000001 000000 where 0 \ / v il-1 \ v21-1 0 0 0 0 06 1^i / -1 V3 v2 V4 V4I- 1 / -1 V5 V61-1 k^v7 v 1-1 / V 6I (5.3) 7/8 V 10 j V ii 13i is defined by the recurrence relation: =1 A = 1(6 - = 2 ,•••, 7 . (5.4) Chapter 5. Hierarchy Building^ 68 The ili s can be calculated as follows: = 2v 11 1 4v;-v:_ 1 = 2, ... , 7. (5.5) Now that the transformation of the system of equations is complete, the equations can be solved by back substitution from the bottom up by performing the following computations: 1-1= V7 V7 74 -1 = V i -^* Vi+i i = 6, ... , 1. (5.6) 5.1.4 Problems with Inverse Subdivision Because of the nature of inverse subdivision, the offsets generated by Algorithm 5.1 are not well behaved. The magnitude of offsets become large and the base level control vertices do not approximate the original surface at all. It turns out that with each inverse subdivision, the lower level we get, the "wilder" the surface shape becomes. At the base level, the control vertices go off everywhere and the surface is not a reasonable approximation of the finer levels. This instability motivates us to find another way to build a hierarchical surface which behaves well at each level. 5.2 Use of MG Data Remember that the multi-resolution data is available from the FMG solver. These multi-level control vertices define a reasonable shape at each resolution level. Instead of just making use of the finest level result to build the hierarchy, we take advantage of the availability of all levels of data to build the overlays. We begin from the lowest level, a flat single-patch surface (Plate A.4). The 16 control vertices VBA4 (Algorithm 2.3) , which are solutions for the coarsest level in FMG, define the base level of our hierarchical B-spline having a parametric spacing of one. Then the base level is subdivided and level 1 control vertices VBAKii are generated. Since we use uniform knot spacing, the parametric spacing for subsequent levels is halved recursively. The FMG data Chapter 5. Hierarchy Building ^ 69 V. at level 1 with compatible dimension is used as the final surface definition at level 1. Thus offsets at level / are calculated by C4i = VBA4 — Vili^(5.7) where we use Vh to indicate the control vertices generated by refining VBAK iii l at level / — 1. The algorithm is given below: Algorithm 5.2^ Hierarchy Building AssignBaseLevel (0) ; for level = 0 to finestLevel do SubdivideLevel (level) ; for = 0 to m for j = 0 ton offset = VBAK [level-14]^[j] - V [level] [i] [j] ; if (offset !=0) SaveNonZeroOff set (off set) ; end. end. end. where m and n are used to indicate the dimensions of control vertices in the u and v directions respectively at each level. The values vary with the level index number. Because each level is a filtered version of the next higher resolution level, the vertices are well behaved. Algorithm 5.2 is not intended to reduce the non-zero offsets, as Algorithm 5.1 was. Adjacent levels in FMG are just interpolations of points at different resolutions. Nothing in the formulation determines the magnitude of the offset. It turns out that almost all of the offsets generated with algorithm 5.2 are non-zero. If we store this as a hierarchical B-spline Chapter 5. Hierarchy Building^ 70 surface, the cost is approximately the same as the cost of storing all of the levels of the FMG solution. The following section shows how to reduce this cost. 5.3 Zeroing Non-zero Offsets The FMG approach works well because the solution at one level is a good initial guess for the solution at the next finer level. When converted into a hierarchical B-spline surface, most of the offsets are non-zero but their magnitude is small (Figure 5.2). The challenge is to reduce the amount of information required without significantly altering the shape of the surface. The first approach is to zero any offsets whose magnitude is less than a given tolerance beginning at the coarsest resolution. Since this operation changes the reference control vertices Rid of the next level, the offsets at that level are recalculated to restore the approximation at that level. 100:A 101:A 102:A 103:A 104:A 105:A 106:A 107:A 108:A 109:A 110:A : : : : : : : : : : : 0.0946581752 0.1082406374 -0.0510964493 -0.4924907012 -0.0812187087 0.2453602217 1.0591058826 0.2264022796 -0.1049558748 -1.0981894014 0.2164414669 0.0206951908 0.0363397667 -0.0067193945 -0.0094809132 0.0016742268 0.0028018600 -0.0002960616 -0.0016349198 -0.0005351651 0.0036999245 0.0024895610 -0.7312056357 -0.1349227312 0.2410037345 1.4855188643 0.1157502732 -0.5366043224 -2.5841572814 -0.5210445247 0.3200479296 2.1388922244 -0.4276355013 Plate 5.2: A portion of hierarchical B-pline data file without zeroing non-zero offsets. Here is the compact algorithm for zeroing out the offsets within tolerance: Algorithm 5.3 Offset Zeroing AssignBaseLevel(0); Chapter 5. Hierarchy Building^ 71 for level = 0 to finestLevel do SubdivideLevel (level) ; for i = 0 to m for j = 0 to n offset = VBAK Clevel+1] [i] [j] - V [level] [i] [j] ; if (abs (offset) > tolerance * level) SaveNonZeroOff set (offset) ; else VBAK Clevel+1] [i] [j] = V [level] [i] [j] ; end. end. end. Figure 5.3 shows the relationship between the number of non-zero offsets remaining after the zeroing process and the tolerance used. If a relatively large tolerance is used, a large number of zero offsets are created, but the resulting surface is a much poorer approximation of the data. With the same offset zeroing technique, the residual per data point (the total residual caused by offset zeroing at the finest level divided by the number of data points interpolated.) is plotted in Figure 5.4. And the maximum residual is plotted in Figure 5.5. 5.4 Smoothing Algorithm Plate A.9 shows the surface with a tolerance of 0.5 out of 200. Because offset are blindly truncated, the difference between offsets which are just above and below the tolerance is exaggerated (Figure 5.6). Rather than truncating certain offsets, the magnitude of all the offsets are reduced by the same amount the surface recedes back to its reference shape. A reasonable amount is that of the offset tolerance, however it can also be specified explicitly by user to meet concrete applications Chapter 5. Hierarchy Building^ 72 Plate 5.3: Offset number generated in hierarchical surface with given tolerance value. Raw data ranges from -100 to +100. Plate 5.4: Residual per data point caused by offset zeroing with given tolerance value. Raw data ranges from -100 to +100. Chapter 5. Hierarchy Building^ 73 Plate 5.5: Maximum Residual caused by offset zeroing with given tolerance value. Raw data ranges from -100 to +100. (Figure 5.7). Here is the complete algorithm incorporating smoothing techniques with the offset zeroing algorithm. Algorithm 5.4 Offset Zeroing with Smoothing Function AssignBaseLevel(0); for level = 0 to finestLevel do SubdivideLevel (level) ; for i = 0 to m for j = 0 to n offset = VBAK [level+1] [i] [j] - V [level] [i] [j] ; if (abs (offset) > tolerance * level){ offset -= off setTolerance ; SaveNonZeroOff set (offset) ; Chapter 5. Hierarchy Building^ 74 b) After zeroing offset 2. Plate 5.6: Offset zeroing side effect. a) Before the offset 2 is zeroed. b) The zeroing of offset 2 results in a bump. Chapter 5. Hierarchy Building ^ 75 b) After zeroing offset 2 with smoothing function. Plate 5.7: Offset zeroing with smoothing function. a) Before the offset 2 is zeroed. b) After offset 2 is zeroed with smoothing function. Chapter 5. Hierarchy Building^ 76 } else VBAK [level+1] [i] [j] = V [level] [i] [j] ; end. end. end. 5.5 Manipulation of Resulting Surface By interactive adjustment of the tolerance and smoothing factor, a compact hierarchical representation of the original surface is produced. The resulting surface can be used a input to the Dragon editor for further manipulation and animation. Within the editor, the absolute offsets generated by the approximation code are converted into the tangent offsets that allow fine surface details to be retained as broad scale manipulation are made [17]. Plates A.20 — A.23 show some of the modified shapes. 5.6 Implementation and Results This method was tested on the data taken from a bust of Victor Hugo (courtesy F. Schmidt) with 264 rows of 361 data points (95304 points. See Figure 1.1 and Plate A.1) with normal magnitude around 200. The value of each data point originally represented the radius from a central axis. As an approximation this information was converted to points in R 3 . The scanned data was sampled in parametric space with a 257 x 257 grid (66049 data points) and approximated using the FMG method described above. A 9-level hierarchical surface interpolating all the sampled points contains 90494 offsets. With a tolerance of 0.1, this number is reduced to 23270 offsets, about 25% of the sample data set. Smoothing would further reduce storage, but the amount of reduction is highly dependent upon how much smoothing is required, which differs from situation to situation. Chapter 5. Hierarchy Building^ 77 Another example is a data set of a digitized face. With 66049 data points, we approximate it with 256 x 256 patches at the finest level. With a tolerance of 0.05, the offset number in its hierarchical representation is 19235 (Plate A.14). Chapter 6 Conclusions and Future Work 6.1 Conclusions In this thesis, the digitized Victor Hugo's bust has been used as a testbed of the surface approximation techniques designed to solve the problem stated above. For this specific case, the data parameterization challenge is met by a parametric space deformation scheme. This scheme achieves a near chord-length parameterization of the raw data points. For a detailed representation of the original digitized object, a large number of linear equations have to be solved and a fast algorithm used to solve these equations is desired. Full multigrid methods (FMG) was chosen to carry out this task because of its optimal performance. The FMG solver offers us interactive control of the fitting process and a multi-resolution surface. To reduce the storage cost (90K points) and to produce a representation suitable for animation, a hierarchical B-spline surface was used to represent the fitted surface. The surface hierarchy was built using two different techniques. The inverse subdivision offered a good compression of the surface storage in terms of non-zero offset number, but the resulting surface was not suitable for multi-resolution editing. The direct use of multi-resolution data from the FMG fitter produced a well-behaved shape using slightly more storage. To further compact the resulting surface, offset zeroing techniques were employed to remove those non-zero offsets within a certain length, and smoothing was used where offset-zeroing produces unsuitable bumps. The approximation to Victor Hugo's bust was controlled through the selection of the offset 78 Chapter 6. Conclusions and Future Work ^ 79 zeroing threshold and smoothing factor. The hierarchical B-spline approximation is suitable for manipulation with the Dragon editor. Complicated shapes like a human face can be fitted to its very details at a reasonable cost. This is a direct way to achieve realistic models without consuming much interactive manpower. The results are encouraging in that suitable results are possible at a reasonable cost. 6.2 Future Work Like most of the research topics, this thesis generated more problems than it actually solved. The initial parameterization scheme results in some sharp turning angles between adjacent data points in parametric space. Subsequently, these angles resulted in notches in the surface under chord-length parameterization. Better initial parameterization is required and a parametric space deforming scheme such as the one presented by Waters and Terzopoulos [39] in improving the final surface fit. The chosen parameterization scheme, although adequate for this particular application, is not general enough. Other methods, such as used in [10], will be investigated with the goal of using scanned data as a mold that can be applied to an existing surface to give it shape. Future work will also look into better restriction methods for FMG and at the possibility of optimizing the number of zero-offsets in the surface definition. Bibliography [1] Ahlberg, J., Nilson, E. and Walsh, J., The Theory of Splines and Their Applications, Academic Press, New York, 1967. [2] Barsky, B. "End Conditions and Boundary Conditions for Uniform B-spline Curve and Surface Representations", Computers in Industry, 3(1&2):17-29, 1982. [3] Bartels, R., Beatty, J. and Barsky, B. An Introduction to Splines for Use in Computer Graphics and Geometric Modelling, Morgan Kaufmann Publishers Inc., Palo Alto, California, 1987. [4] Barsky, B. and Greenberg, D. "Determining a Set of B-spline Control Vertices to Generate an Interpolating Surface", Computer Graphics and Image Processing, V.14:203-209, 1979. [5] Barsky, B., A study of the Parametric Uniform B-spline Curve and Surface Representation, Univ. of California at Berkeley Technical Report UCB-CSD-83-118, May 1983. [6] Birkhoff, C. and MacLane, S., A Survey of Modern Algebra, The Macmillan Company, New York, 1965. [7] Brandt, A., "Multi-level Adaptive Solutions to Boundary Value Problems", Mathematics of Computation, V.31(138):333-390, 1977. [8] Briggs, W. A Multigrid Tutorial, SIAM, 1987. [9] Cohen, E. and Lyche, T., and Reidenfeld, R., "Discrete B-splines and Subdivision Techniques in Computer Aided Geometric Design and Computer Graphics", Computer Graphics and Image Processing, V.14(2):87-111. [10] Cohen, E., O'Dell, C., "A Data Dependent Parameterization for Spline Approximation", Mathematical Methods in Computer Aided Geometric Design:155-166, edited by T. Lyche and L. Schumaker, Academic Press, Boston, 1989. [11] Coons, S., "Surface Patches and B-spline Approximation", Computer Aided Geometric Design:1 16, edited by Barnhill, R. and Riesenfeld, R., Academic Press, New York, 1974. - [12] de Boor, C., A Practical Guide to Splines, Spring-Verlag, New York, 1978. 80 Bibliography^ 81 [13] Epstein, M., "On the Influence of Parameterization in Parametric Interpolation", SIAM Journal of Numerical Analysis, 13(2):261-268, 1976. [14] Farin, G., Curves and Surfaces for Computer Aided Geometric Design, Academic Press Inc., 1988. [15] Foley, T., "Interpolation with Interval and Point Tension Controls Using Cubic Weighted v-Splines", ACM Transactions on Mathematical Software, V.13(1):68-96, 1987. [16] Foley, T. and Nielson G., "Knot Selection for Parametric Spline Interpolation", Mathematical Methods in Computer Aided Geometric Design:261-271, edited by T. Lyche and L. Schumaker, Academic Press, Boston, 1989. [17] Forsey, D. Motion Control and Surface Modelling of Articulated Figures in Computer Animation, Univ. of Waterloo Technical Report CS-90-28, (Ph.D. Thesis) 1990. [18] Forsey, D. and Bartels, R. "Tensor Products and Hierarchical Fitting", Proceedings IEEE Curves and Surfaces in Computer Vision and Graphics II, Boston, Nov. 1991. [19] Forsey, D. and Bartels, R. "Surface Fitting with Hierarchical Splines" to appear in Transactions on Graphics. [20] Forsey, D., and Bartels, R. H., "Hierarchical B-spline Refinement", Proceedings of SIGGRAPH '88:205-212, Atlanta, Georgia, 1988. [21] Hackbusch, W., Multigrid Methods and Applications, Springer-Verlag, 1985. [22] Hoschek, J., "Intrinsic Parameterization for Approximation", Computer Aided Geometric Design, V.5:27-31, 1988. [23] Kass, M., Witkin, A. and Terzopolous, D., "Snakes: Active Contour Models", International Journal of Computer Vision, V.4:321-331, 1988. [24] Lee, E., "On Choosing Nodes in Parametric curve Interpolation", SIAM Conference on Applied Geometry, Albany, New York, 1985. [25] Lin, C., Chang, P. and Luh, J., "Formulation and Optimization of Cubic Polynomial Joint Trajectories for Mechanical Manipulation", Proceedings of 21st IEEE Conference on Decision and Control, V.1:330-335, Orlando, Florida, December, 1982. [26] Lyche, T. and Morken, V., "Knot removal for parametric B-spline curves and surfaces", Computer Aided Geometric Design, V.4(3):217-230, 1987. [27] Marin, S., "An Approach to Data Parameterization in Parametric Cubic Spline Interpolation Problems", Journal of Approximation Theory, V.41:64-86, 1984. Bibliography^ 82 [28] McCormick, S., "Multilevel Adaptive Methods for Partial Differential Equations", SIAM, 1989. [29] Mullineux, G., "Approximating Shapes Using Parameterized Curves", IMA Journal of Applied Mathematics, V.29(2):203-220, 1982. [30] Nahas, M., Huitric, H. and Saintourens, M., "Animation of a B-Spline Figure", The Visual Computer, V.3(4), 1987. [31] Nielson, G., "Coordinate Free Scattered Data Interpolation", Topics in Multivariate Approximation:175-184, edited by C. Chui, 1. Schumaker and F. Utreras, Academic Press, New York, 1987. [32] Nielson, G. and Foley, T., "A Survey of Applications of an Affine Invariant Norm", Mathematical Methods in Computer Aided Geometric Design:445-467, edited by T. Lyche and L. Schumaker, Academic Press, Boston, 1989. [33] Schmitt, F., Barsky, B. and Du, W. "An Adaptive Subdivision Method for Surface-Fitting from Sampled Data", Proceedings of SIGGRAPH '86:179-188, Dallas, Texas, 1986. [34] Schmitt, F., Maitre, H., Clainchard, A. and Lopez-Krahm, J., "Acquisition and Representation of Real Object Surface Data", SPIE Proceedings Vol. 602 of Biostereometrics '85 Conference, Cannes, France, 2-6 December, 1985. [35] Screckovic, M., Adaptive Hierarchical Fitting of Curves and Surfaces, Master's Thesis, University of Waterloo, 1991. [36] Sederberg, T. and Parry, S. "Free-Form Deformation of Solid Geometric Models", Proceedings of SIGGRAPH '86:151-160, Dallas, Texas, 1986. [37] Terzopoulos, D., "Image Analysis Using Multigrid Relaxation Methods", IEEE Transactions on Pattern Analysis and Machine Intelligence, V.8(2):129-139, 1986. [38] Terzopoulos, D.,"Multi-level Surface Reconstruction", Computer Vision, Graphics and Image Processing, V.24:52-96, 1983. [39] Waters, K. and Terzopoulos, D., "Modelling and Animating Faces Using Scanned Data", The Journal of Visualization and Computer Animation, V2:123-128, 1991. [40] Welch, W. and Witkin, A. "Variational Surface Modelling", Proceedings of SIGGRAPH '92:157-166, Chicago, Ill, 1992. [41] Wolberg, G., "Digital Image Warping", IEEE Computer Society Press, 1990. Appendix A Plates This appendix includes the pictures generated by the implementation as results. Appendix A. Plates^ Plate A.1: Victor Hugo raw data. Plate A.2: Deformed UV isoparametric lines in ST space with 5 level refinement. 84 Appendix A. Plates^ 85 Plate A.3: Deformed UV isoparametric lines in ST space with 5 level refinement. Data points are plotted over the space. Appendix A. Plates^ Plate A.4: FMG fit at level 1 and level 2. Plate A.5: FMG fit at level 3 and level 4. 86 Appendix A. Plates^ Plate A.6: FMG fit at level 5 and level 6. Plate A.7: FMG fit at level 7 and level 8. 87 Appendix A. Plates ^ 88 Plate A.8: Left: FMG fit at level 9. Right: Hierarchical B-spline surface without zeroing offsets. Offset number is 90949. Plate A.9: Left: Hierarchical B-spline surface with offset zeroing threshold 0.5 and smoothing algorithm. Offset Number is 3913. Right: Hierarchical B-spline surface with offset zeroing threshold 0.5 without smoothing algorithm. Appendix A. Plates^ 89 Plate A.10: Left: Hierarchical B-spline surface with offset zeroing threshold 0.1 and smoothing algorithm Offset Number is 23270. Right: Hierarchical B-spline surface with offset zeroing threshold 1.0 and smoothing algorithm. Offset Number is 1918. Plate A.11: Left: Hierarchical B-spline surface with offset zeroing threshold 2.0 and smoothing algorithm. Offset Number is 930. Right: Hierarchical B-spline surface with offset zeroing threshold 5.0 and smoothing algorithm. Offset Number is 318. Appendix A. Plates ^ 90 Plate A.12: Left: Surface manipulation - Smile 1. Right: Surface manipulation - Smile 2. Plate A.13: Left: Surface manipulation - Smile 3. Right: Surface manipulation - Neanderthal man. Appendix A. Plates ^ Plate A.14: Another fit to a digitized face. 91
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multi-resolution surface approximation for animation
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multi-resolution surface approximation for animation Wang, Lifeng 1993
pdf
Page Metadata
Item Metadata
Title | Multi-resolution surface approximation for animation |
Creator |
Wang, Lifeng |
Date Issued | 1993 |
Description | This thesis addresses the problem of approximating a set of gridded data points obtained from a three-dimensional digitizing system to create a representation with a hierarchical bicubic B-spline surface that is suitable for further manipulation and animation. Chord length parameterization is obtained using 2-D deformation technique. A full multigrid (FMG) numerical method is used to solve the surface approximation and the multi-resolution elements created are used directly to define the overlays in a hierarchical B-spline surface. The direct use of FMG multi-resolution data offers reasonable surface shape behaviour, but the number of non-zero offsets is large. Storage cost is reduced either by eliminating offsets whose magnitude is below a certain tolerance or by reducing all offsets in a given level by a user specified amount. The resulting spline surface is modifiable, both locally and globally while retaining surface details of the digitized data. An interactive system based on these methods was created and the results of approximating two large data sets are presented. |
Extent | 5315802 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-09-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051338 |
URI | http://hdl.handle.net/2429/1877 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1993-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1993_spring_wang_lifeng.pdf [ 5.07MB ]
- Metadata
- JSON: 831-1.0051338.json
- JSON-LD: 831-1.0051338-ld.json
- RDF/XML (Pretty): 831-1.0051338-rdf.xml
- RDF/JSON: 831-1.0051338-rdf.json
- Turtle: 831-1.0051338-turtle.txt
- N-Triples: 831-1.0051338-rdf-ntriples.txt
- Original Record: 831-1.0051338-source.json
- Full Text
- 831-1.0051338-fulltext.txt
- Citation
- 831-1.0051338.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0051338/manifest