Multi-resolution Surface Approximation for AnimationByLIFENG WANGB.Sc., University of Science and Technology of ChinaChina, 1990A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE STUDIES(DEPARTMENT OF COMPUTER SCIENCE)We accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAApril, 1993© Lifeng Wang, 1993In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)Department of ^Computer ScienceThe University of British ColumbiaVancouver, CanadaDate^April 28th, 1993DE-6 (2/88)AbstractThis thesis addresses the problem of approximating a set of gridded data points obtainedfrom a three-dimensional digitizing system to create a representation with a hierarchical bicu-bic B-spline surface that is suitable for further manipulation and animation. Chord lengthparameterization is obtained using 2-D deformation technique.A full multigrid (FMG) numerical method is used to solve the surface approximation andthe multi-resolution elements created are used directly to define the overlays in a hierarchicalB-spline surface. The direct use of FMG multi-resolution data offers reasonable surface shapebehaviour, but the number of non-zero offsets is large. Storage cost is reduced either by elimi-nating offsets whose magnitude is below a certain tolerance or by reducing all offsets in a givenlevel by a user specified amount. The resulting spline surface is modifiable, both locally andglobally while retaining surface details of the digitized data.An interactive system based on these methods was created and the results of approximatingtwo large data sets are presented.iiContentsAbstract^ iiContents iiiList of Figures^ viAcknowledgement ix1 Introduction 11.1 Overall Goal ^ 11.2 Surface Approximation ^ 21.2.1^Mathematical Formulation ^ 21.2.2^Discrete Data 21.2.3^Interpolation and Approximation 31.3 Problem Definition ^ 41.3.1^Test Data 41.3.2^Data Parameterization ^ 61.4 Related Work ^ 61.5 Overview 82 Mathematical Background 102.1 Uniform Bicubic B-Spline Surfaces ^ 102.2 Subdivision Scheme ^ 142.3 Hierarchical B-splines 162.3.1^Practical Consideration of Tensor-Product Surfaces ^ 162.3.2^Localized Subdivision-Hierarchical B-splines 172.4 Parameterization ^ 192.4.1^Parameterization Schemes ^ 202.4.2^Summary 222.5 Multigrid Numerical Method 232.5.1^Iterative Methods on Linear Algebraic Equations ^ 232.5.2^Multigrid Philosophy^ 25iii3^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 Parametric Space Mapping2627273031343.1 Parameterization ^ 343.2 Initial Parameterization ^ 363.3 Deforming Parametric Space 383.3.1^Deformation Function 383.3.2^Refinement and Local Control ^ 403.3.3^Deformation ^ 403.4 Sampling Data Points 473.5 Summary ^ 484 Multi-resolution Surface Fitting Using Multigrid 494.1 Deriving Equations and Constraints ^ 494.1.1^Data Point Distribution 494.1.2^Deriving the Positional Interpolation Equations ^ 514.1.3^Deriving the Boundary Conditions ^ 524.2 Solving For Control Vertices ^ 534.2.1^Full Multigrid Method 534.2.2^Full Multigrid Operators 544.3 Performance Evaluation ^ 595 Hierarchy Building 625.1 Inverse Subdivision ^ 625.1.1^Hierarchical Surface File Format ^ 625.1.2^Inverse Subdivision Scheme 645.1.3^Remarks: Solution to Triangular Linear Equation^ 675.1.4^Problems with Inverse Subdivision ^ 685.2 Use of MG Data ^ 685.3 Zeroing Non-zero Offsets ^ 705.4 Smoothing Algorithm 715.5 Manipulation of Resulting Surface ^ 765.6 Implementation and Results ^ 766 Conclusions and Future Work 786.1 Conclusions ^ 786.2 Future Work 79ivBibliography^ 80A Plates 83Appendix^ 83List of Figures1.1 Digitized Data Topology. ^ 52.1 A bicubic B-spline surface with a single patch defined by 16 control vertices.. .^112.2 A B-spline surface is a mosaic of surface patches^ 112.3 The parametric range for a surface is the union of the parametric range for eachpatch. ^ 132.4 Subdivision splits a single patch into four patches. ^ 152.5 Subdivision splits the parametric space at the midpoint 162.6 Subdivision through non-local knot insertion. 172.7 A wave on Slh (N = 16) is projected onto C2 2h (N = 8). On the coarse grid, thewave seems more oscillatory than on the fine grid. ^ 252.8 The schedule for grid visiting when the multigrid V-cycle is executed. Each levelis identified with its spacing. Multigrid V-cycle starts from the finest level andproceeds to the coarsest level, then returns to the finest level^ 282.9 The schedule for grid visiting of W-cycle scheme^ 302.10 Full multigrid grid visiting schedule. ^ 303.1 Labels for the tensor-product B-spline surface. ^ 353.2 Two approaches to approximating the head data: a) wrap a cylindrical surfacearound the head, b) fold the four edges of a rectangular surface to cover up thehead with the opening at the neck^ 353.3 The initial mapping from data space to parametric space. Each row of data ismapped onto a square in parametric space starting from the center to the edge. 363.4 Two coordinate frames imposed on the parametric space. 383.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 asin ST space. ^ 393.6 The control grid in UV space is uniformly refined from a 2 x 2 grid into a 3 x 3grid. The control points inside brackets are the indices of the 4 control points asthey are indexed at level 1. ^ 423.7 The three cases of the sample point falling into a data quadrilateral. ^ 453.8 A control point is displaced from its original position. ^ 46vi^3.9 The UV space is evenly sampled for produce the parametric pair set 484.1 The (X +1) x (Y +1) data points are mapped onto X x Y patches. Every patchcorner corresponds to a data point^ 504.2 A data point Z23 is mapped onto the upper right corner of a single patch. At thesame time, it is the corner point of three adjacent patches. ^ 514.3 A illustration of the four-colour relaxation scheme on a uniform grid. First the0 points are swept, then the 0 points, A points and finally ^ points. ^ 545.1 The inverse subdivision process generate a lower level curve from a upper levelcurve in a way that reduces the number of non-zero offsets while the hierarchyis built using subdivision. 655.2 A portion of hierarchical B-pline data file without zeroing non-zero offsets. . . .^705.3 Offset number generated in hierarchical surface with given tolerance value. Rawdata ranges from -100 to +100. ^ 725.4 Residual per data point caused by offset zeroing with given tolerance value. Rawdata ranges from -100 to +100. 725.5 Maximum Residual caused by offset zeroing with given tolerance value. Rawdata ranges from -100 to +100^ 735.6 Offset zeroing side effect. a) Before the offset 2 is zeroed. b) The zeroing of offset2 results in a bump. ^ 745.7 Offset zeroing with smoothing function. a) Before the offset 2 is zeroed. b) Afteroffset 2 is zeroed with smoothing function. ^ 75A.1 Victor Hugo raw data^ 84A.2 Deformed UV isoparametric lines in ST space with 5 level refinement. ^ 84A.3 Deformed UV isoparametric lines in ST space with 5 level refinement. Datapoints are plotted over the space^ 85A.4 FMG fit at level 1 and level 2. 86A.5 FMG fit at level 3 and level 4. 86A.6 FMG fit at level 5 and level 6. ^ 87A.7 FMG fit at level 7 and level 8. 87A.8 Left: FMG fit at level 9. Right: Hierarchical B-spline surface without zeroingoffsets. Offset number is 90949^ 88A.9 Left: Hierarchical B-spline surface with offset zeroing threshold 0.5 and smooth-ing algorithm. Offset Number is 3913. Right: Hierarchical B-spline surface withoffset zeroing threshold 0.5 without smoothing algorithm. 88A.10 Left: Hierarchical B-spline surface with offset zeroing threshold 0.1 and smooth-ing algorithm. Offset Number is 23270. Right: Hierarchical B-spline surface withoffset zeroing threshold 1.0 and smoothing algorithm. Offset Number is 1918. . . 89viiA.11 Left: Hierarchical B-spline surface with offset zeroing threshold 2.0 and smooth-ing algorithm. Offset Number is 930. Right: Hierarchical B-spline surface withoffset zeroing threshold 5.0 and smoothing algorithm. Offset Number is 318. . . . 89A.12 Left: Surface manipulation - Smile 1. Right: Surface manipulation - Smile 2. . . 90A.13 Left: Surface manipulation - Smile 3. Right: Surface manipulation - Nean-derthal man. ^ 90A.14 Another fit to a digitized face. ^ 91viiiAcknowledgementsI would like to thank Dr. David R. Forsey, my thesis supervisor, for his guidance andencouragement throughout my work on this thesis. The valuable ideas and suggestions fromhim 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 readingthrough 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 grammarand 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 onthis thesis. His comments have helped a lot on improving my work. Many other people in thedepartment offered time and effort my work. It is difficult to list all the individuals who hasbeen involved with me. But those people who have provided help know that they deserved anote of THANKS.If this is the chance, I would like to direct the most important acknowledgement to myparents, who have been providing priceless love and support over the past. Without theirencouragement 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.ixChapter 1Introduction1.1 Overall GoalThis research is a part of a long term project to create a system for the animation and renderingof human animal forms. An integral part of this system is a mechanism whereby digitized surfaceinformation is used to determine the surface characteristics ( shape and colour) of an animatedfigure. The underlying surface formulation for this system is hierarchical B-spline [17] [20], amulti-resolution approach to modelling with tensor product surfaces. Thus the mechanisms forapproximating 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 approxi-mating regular gridded data with a hierarchical spline. The contribution of this work is the useof 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. Surfaceapproximation for animation adds one additional complication. Because the approximatingsurface may only be a part of a much larger surface, the question of what part of the surfaceapproximates what part of the data must also be addressed.1Chapter 1. Introduction^ 21.2 Surface Approximation1.2.1 Mathematical FormulationWe are interested in fitting a tensor-product B-spline surfacem nS(L,V) = E Eto a given set of data points.To generate the equations for surface approximation, each data point (x, y, z) = DA E R3is associated with a domain point (u, v) = SA E R 2 of the spline. This forms the equation:m nE E^DA.i=0 j=0(1.2)Data is gridded ifu E {2/0, • • • , um} (1.3)E {vo, • " , vN}and if the Ss consist of all points in {Ito, • • • ,^x {vo, • • • , vN} the surface approximationequations becomem nE E (1.4)i=0 j=0for s = 0, • ,M and t = 0, • • • , N. Thus the number of equations is equal to the number ofdata points.Hierarchical splines are a multi-resolution approach to splines for use in interactive cre-ation of free-form surfaces [17] [40]. The problem of applying these equations to approximatea digitized surface with a hierarchical spline is complicated by two factors: an unusual param-eterization, and the multi-resolution nature of the hierarchical formulation itself.Readers who are familiar with parameterization and surface fitting can skip the followingsections and go directly to Chapter 3.1.2.2 Discrete DataSystematic measuring devices such as laser scanners, imaging systems and optical scanners areused to measure the positions of points on an object's surface with a resolution high enoughChapter 1. Introduction^ 3to provide an accurate representation of the original object. These devices work in threedimensional space and typically produce a 2-D array of data points. A medical imaging systemlike CT (Computed Tomography) records both exterior and interior information, producing aregularly spaced 3-D lattice.Raw data sets for human figures are usually too large to manipulate in real-time on generalpurpose graphics workstations. Techniques used to approximate the raw data must be veryefficient both in terms of computing time and storage complexity.1.2.3 Interpolation and ApproximationA common technique for recreating the smooth surface from a set of discrete digitized datapoints discrete data set as a continuous shape is a graph. Other techniques include fitting sur-faces (interpolation and approximation), volume rendering and various visualization techniques.There are basically two ways to reconstruct the continuous model from digitized data: inter-polation and approximation. Interpolation requires that the surface passes through every datapoint. Piecewise parametric polynomial interpolation is a common technique [14] and there arevarious other mathematical formulations which meet different requirements [4] [16] [24]. Mostinterpolation methods are simple in concept, have unique solutions and provide an accuraterepresentation of the original data.Approximation is designed to fairly represent the data without necessarily interpolatingthe individual data points and are often useful for data sets comprised of measured valuesincluding errors where their exact interpolation is not desired. Approximate algorithms canbe adaptive or non-adaptive. Adaptive algorithms begin with a simple approximation definedby a small number of coefficients. The number of degrees of freedom in the approximation isadaptively increased in areas of high error until the surface is within a given tolerance of thedata. Adaptive algorithms are efficient because some sub-regions of data have less gradientvariation than other regions. The smoother regions can be approximated with coarse mesheswithout exceeding certain error tolerances, reducing the storage cost. The bumpy regions shouldbe approximated with fine meshes.Chapter 1. Introduction^ 4There are numerous approaches to both approximation and interpolation. Each mathemat-ical representation will produce different approximations for the same set of data. Piecewiseparametric surfaces are one family of tools for approximating shapes [14]. Notable examplesof such formulations include Bezier curves [14], Coons patches [11], and B-splines [3]. Theserepresentations have been used extensively throughout the CAGD field to describe shape incomputer systems. Some mathematical background concerning these representations will begiven in Chapter 2. More details can be found in books by FarM [14] and by Bartels et. al. [3].1.3 Problem Definition1.3.1 Test DataThe data set chosen for use in this thesis is representative of the size and detail of the modelsused in animation, the bust of Victor Hugo (courtesy F. Schmidt, see Plate A.1) obtained froma commercial three-dimensional digitizing system developed by the Laboratoire IMAGE of theEcole Nationale Superieure des Telecommunications [34]. In this integrated device, the micro-processor controls a video-laser data acquisition system capable of measuring approximately200,000 points per minute. The rotation axis (Y axis in Figure 1.1) is usually perpendicularto the beam and along the center line of the object. Figure 1.1 shows the topology of the dataobtained from these kinds of devices.The resulting data are represented as a 2-D array of radii; that is, any non-boundary datapoint 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 approx-imating the surface, the radius values rij are converted to the coordinate triple (xij, yij, zij) inR3 as follows:xij = rij x cos(j),yij = M —^ (1.5)zij = —rij x sin(j),where M is the number of rows of points. In the Victor Hugo data set, there are 264 rowsof 361 data points (95304 points in all). For notational purposes, the data is indexed by rownumber {i, 0 < i < 263} and by column number {j, 0 < j < 360}.Chapter 1. Introduction^ 5Plate 1.1: Digitized Data Topology.Chapter 1. Introduction^ 61.3.2 Data ParameterizationFor 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 splinesurface 1 .For this particular application, the digitized data is used to specify the shape of a portionof 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 singleparametric surface encasing the entire body (This is part of the mechanism used to control thedeformation of the surface around the articulations.). The particular portion of this surfacededicated to the head is not a cylinder, but a rectangular sub-region of the surface where theedges of the region correspond to the opening at the neck and the interior of the region thatencompasses 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 toapproximate cylindrical data.1.4 Related WorkThe literature on surface fitting and approximation is vast. This section will not attempt tosurvey the field but will briefly examine work on surface approximation for animation. Chapter 2will also cover specific issues of data parameterization, numerical methods, and parametricsurface formulations.The work of Nahas et al. [30] addresses the same problem of using digitized data to create ananimated human figure. A raw digitized mesh of data points is used without further processingas the control vertices of a bicubic B-spline surface. For animation, broad-scale changes tothe surface shape are accomplished by embedding the control vertices of the model in a coarserspline 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 istermed parameterization.Chapter 1. Introduction^ 7model of the surface. When the characteristic points are moved, they displace the vertices ofthe 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 ageometrically continuous collection of rectangular bicubic Bezier patches. Bezier patch basisfunctions do not automatically ensure the continuity between adjacent patches, thus duringthe approximation process, extra constraints must be imposed to ensure at least geometriccontinuity (collinear tangents) between patches. The adaptive behaviour is implemented byfitting 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 controlledby using a constrained least squares approximation method where the constraints are imposedby the continuity conditions. A given number of constraints per patch are required to piecethe Bezier patches together to form a continuous composite surface. Because of the propertiesof the basis functions,a Bezier surface will have more control vertices (9 times) than that ofthe equivalent B-spline surface. (Note, however, that with Bezier patch, it has considerablyeasier to control the continuity between patches.) Because each Bezier patch is essentiallyindependent, 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-splinesand their rational counterparts, which provide continuity without the explicit imposition ofconstraints, must refine an entire row or column of patches to increase the number of patchesin the surface.A hierarchical spline [17] [20] allows local refinement of uniform B-splines or Beta-splinesand 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 adaptiveapproach to fitting regular (gridded) data in three dimensions, with a hierarchical bicubic B-spline, and generalized to multivariate B-splines in any dimension [18]. This approach utilizes aleast squares to fit a coarse spline surface (i.e., on the order of 1-16 patches) of the appropriatetopology. 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 ofChapter 1. Introduction^ 8inadequate fit become separated and surrounded by regions within tolerance. These regions areindividually accommodated using local refinement and constrained least squares approximationwithin each separable region. This approach results in a multi-resolution surface formulationin a hierarchical format.The hierarchical approach provides a further contribution by enabling a certain economy ofrepresentation for the final composite surface and produces reasonable results. However, theresulting surface suffers from oscillations when presented with data with high frequency compo-nents. The oscillations were reduced considerably through the use of non-uniform hierarchicalB-splines and an improved parameterization procedure [35].The surfaces resulting from the above approach were not appropriate for further manipula-tion via the hierarchy. Because each level in the hierarchical surface is the result of a separateleast squares approximation of the data, corresponding parametric locations in each overlaymay have widely varying slopes. This plays havoc with the surface when broad-scale changesin shape are made.In contrast to the top-down approach in [33] and [20], Lyche and Morken [26] work bottom-up by first approximating the surface with a fine mesh of spline patches and then removingknots in those regions where they will not cause the surface to move out of tolerance. Thisthesis presents a bottom-up approach to surface approximation using hierarchical splines.1.5 OverviewChapter 2 covers the mathematical background for B-splines, refinement and mid-point sub-division, hierarchical B-spline, parameterization and multigrid numerical method. Chapter 3gives in detail the modified chord-length parameterization method used for the test data, andChapter 4 discusses the construction of the linear equations used for the approximate and theirsolution using the full multigrid method. Details of how the full multigrid method is used toconstruct a hierarchical surface is given in Chapter 5 along with one approach to reducing thestorage requirement of the approximating surface. Conclusions and future work are presentedChapter 1. Introduction^ 9in Chapter 6.Chapter 2Mathematical Background2.1 Uniform Bicubic B-Spline SurfacesThe use of parametric surfaces originated in the field of CAGD. B-spline tensor-product surfacespermit the most supple movements for the purposes of animation [5]. Within this family ofsurfaces, we are interested in the uniform case when knot spacing is constant in both the Uand V parametric directions. A detailed derivation of the uniform B-spline basis functions andefficient algorithms for their evaluation, along with the reasons for their extensive use in CAGDwere presented in the paper by Barsky [5]. The properties presented in this thesis, such aslocal control and control of the degree of continuity at the joints between surface patches, arevery desirable for fitting a spline surface to discrete data [3]. In this section, we will brieflyreview 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 oforder 4 defined by 4 x 4 control vertices. Notice that the surface does not necessarily interpolatethe control vertices.As shown in Figure 2.2, a bicubic B-spline surface is arranged in a piecewise manner, whereeach piece is a segment of the surface called a surface patch. Each surface patch is defined bya subset of the control vertices. The number of control vertices required to define a surfacepatch depends upon the order of the underlying basis functions. The entire surface S(u, v) is a10Chapter 2. Mathematical Background^ 11V12 V22V00 V10 V20 V30Plate 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^ 12mosaic of these patches sewn together and has continuity appropriate to the order of the basisfunctions (i.e., Ck -1 where k is the degree of the B-spline). Tensor-product surfaces have arectangular topology and we will label each individual patch as Pii. A bicubic B-spline surfaceconsists 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 inthree-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 followingway:1^ 1Pii (et,v)= E E vi-F,,J+8Br,k(u)Bs,h(v) (u,u) E RUi_i, Vi_i) x (Ui, (2.1)r=-2 s=-2where k and h are the polynomial orders of the basis functions (for cubic B-spline, they areboth 4.) and the x sign means the span in parametric space and will be explained related toEquation 2.3.The bivariate basis function Bi,k(ui)Bi,h(v .i) is the tensor-product of a set of univariatebasis functions, where (Ui, E R2 is an array of knot positions in parametric space. Foruniform 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 oneimportant property of a B-spline formulation since it means that a control vertex controls thesurface shape locally. A single bicubic B-spline surface patch is fully defined by 4 x 4 controlvertices and is unaffected by any other control vertex in the surface. Conversely, a given controlvertex 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 functionsare:[B_2(u) B_1(u) Bo(u) Bi(u)] = —61 [u3u2u 1 1]—1_ 313—0°4—33311000(2.2)Chapter 2. Mathematical Background^ 13A surface S(u, v) is a smooth composition of all patches. It is defined as:m nS(11 5 V) = E^(u, E RUO, VO) x ( Um , V.11 )]•i=1 j=1(2.3)Note that the parameter ranges for the U and V directions are the union of each patch'sparametric range (Figure 2.3) in Equation 2.1, namelyRUo, Vo ) x (U., V.)] = RUo, vo ) x (U1,14) ] u RU1,^x (U2, V2)] u • • • R^Vn -1 ) x ( Um , VOLA one-to-one mapping provides a unique position in parametric space for each surface pointS(u, v). However, when we employ only the uniform B-spline formulation, separate parameterranges for each patch is possible as long as they are consistent. For example, we can have theparameter range (0, 1) in both U and V directions for all the patches.(U0, Vn) (Um, Vn)(Uo Vo) (Um, Vo)Plate 2.3: The parametric range for a surface is the union of the parametric range for eachpatch.Chapter 2. Mathematical Background^ 142.2 Subdivision SchemeAn important operation that can be performed on spline curve and surfaces is subdivision (alsoknown as refinement) which takes a set of control vertices and basis functions and representsthe surface using more control vertices and basis functions. In this section, we will presenta brief review of this scheme and discuss a specific but important type of subdivision calledmidpoint subdivision.Definition 2.1 Given a B-spline representationm nS(L,V) = E E Bi ,k(u)Bj,h(V)Vi J.i=1 j=1The procedure of describing S(u, v) in terms of a larger set of basis functions Ni ,k and Ni ,h anda larger net of control points Wij is called subdivision, namelyM NS(u, V) = E E Ni,k(u)Nj(v)Wi ,j.i=1 j=1M and N are the number of control vertices in U and V dimension respectively after subdi-vision.0In Figure 2.4, a surface consisting of a single patch is refined into four patches. Note thatafter subdivision, the patch does not change shape, there are simply more control vertices inthe control graph. Subdivision is achieved by inserting new knot positions into the existingknot sequences defining Bi,k(U) and Bj,h(v). The newly generated basis functions are expressedin terms of discrete splines cxj,k(j) [9] in the following way:Bi ,k(U) = i ai,k(i)N j,k(u)^ (2.4)j=1where Su is the number of knots inserted in the U parametric direction. The subdivision in theV direction is similar.1 patch4 patchescontrol graphsubdivisionpatch subdivisionChapter 2. Mathematical Background^ 15From Equation 2.4, the net of control points Wii can be derived in the same way:m nWij = E E ar ,k(i) a 3 ,i(i)Vr 0 •^ (2.5)r=1 s=1We can express Equation 2.5 in matrix form as:[W] = [ai][V][a2]•^ (2.6)Plate 2.4: Subdivision splits a single patch into four patches.As stated before, the most straightforward subdivision splits each parametric range in Uand V at its midpoint. This converts each patch determined by the control vertices [V] intofour patches (Figure 2.4). In parametric space, the region for the original patch is split intofour regions along the midlines as shown in Figure 2.5. For uniform cubic B-splines:[al] =^[M1] for regions I and III{ [M2] for regions II and IV,for regions I and II[a2] = { [Mi l[M2] for regions III and IV,(2.7)/ 1 3 1 0o I ^oo^3 1 ) .0 0(2.8)where100og 1[m2] =[mi] =Chapter 2. Mathematical Background^ 16Plate 2.5: Subdivision splits the parametric space at the midpoint.2.3 Hierarchical B-splines2.3.1 Practical Consideration of Tensor-Product SurfacesSubdivision for tensor-product surfaces is non-local because the nature of knot insertion is non-local. For example, if patch Pii in Figure 2.6a is the place where we want to refine, this canonly be done by inserting two knots. One sequence of knot positions is along the U parametricdirection { Uik 11 < k < n} and the other sequence of knot positions {Vkj I 1 < k < m} is alongthe 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 columnChapter 2. Mathematical Background^ 17j 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 patchesversus 3 patches). When modelling large complex surfaces, this property could add hundredsof unnecessary patches.Plate 2.6: Subdivision through non-local knot insertion.2.3.2 Localized Subdivision—Hierarchical B-splinesThe hierarchical B-spline [17] is designed to localize the effect of subdivision through the useof overlay surfaces. Overlay surfaces are hierarchically controlled subdivisions. It is importantto keep in mind that the W surface (Equation 2.5) is an exact re-representation of its parentV surface. If one of the control vertices of W is moved, then the W surface departs from its Vparent, 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 Vdefinition as our basic description of the surface, and save only those control vertices of W which1 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^ 18define the position of the displaced surface. The child surface defined by only those retainedcontrol vertices of W is called an overlay. This portion is regarded as a separate surface to bemanipulated. 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 onlythose 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 editedarea are expected to remain embedded in that area. They will follow editing changes only ifthey can be dynamically tied to that area. This means that their control vertices must movein accordance with the movement of the section of the surface being edited.Hierarchical surface representation meets this requirement by introducing the reference andoffset concepts. At level 1, each control vertex Vii of any part of the surface is of the form of:Vlj = Ri i + Oij (2.9)where R lij is the reference control vertex and (4 is the offset. Usually, Rij is derived fromthe parent surface definition through subdivision, and 0!_ i is the change to the original Wirsvectors resulting from modification of the overlay surface. Thus the surface representation inEquation 2.3 at level 1 can be rewritten as:= EME7Li= n..7!. 1. BRui)g. (v j)(R I• • +^-) (2.10)=^E7=1 ni (^\ PP12'3J- \-'7271311Ein_i M(Iii).Bi(Vi)0! J.1 J.J + z_ai=Here the superscript 1 is used to indicate the parameter associated with a certain hierarchylevel 1. Now we can see that each overlay surface is composed of two components. One is thereference surface:SIR(U, V) = E E /3! (no B((v)•.R(3^3^2,3i=1 j=1and the other is the offset surface:nzt ni4(7.1, V) = E E ki! (n i )B( (v •)(:)( •3 ^JOi=1 j=1(2.11)(2.12)S1 (u, v)Chapter 2. Mathematical Background^ 19Editing changes to the V-defined surface will automatically cause revisions in the W-definedsurface by directly changing the R's, and through them revisions are made dynamically to theW's. Editing changes to the W-defined surface are recorded as changes to the O's. Thus anychanges in coarser levels will change all of the finer level R's and result in a global change insurface shape. Changes to the offsets at the current level affect the current level of refinementand all of the finer levels.Since the reference information R i i at level 1 is always directly derivable from its parentsurface, the crucial information needed to store the upper levels are the non-zero offsets. Re-peated substitution of Equation 2.11 into Equation 2.10 results in the recursive formulation ofEquation 2.3 with L being the finest level:S(u, v) =E7i.1-21E11BP(Ui)BY(Vi)ViCijE rin=11 E 7.11 1 /31. (Ui)BI(VAOti+^Ilf'(ui)Bt (vi)0 '41i .(2.13)Thus, a full representation of a hierarchical surface requires only the control nodes of the rootlevel, 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 acomplex shape because the storage for the surface is proportional to the number of editedpoints.2.4 ParameterizationDefinition 2.2 Given a data set {Dij^(xii,yii,zij)1Dii E R3 , (1,1) < (i, j) < (M, N)} anda parametric position set^= (ujj,vii)1Uii E [(Uoo, Voo) x (Umn, limn)] C R2 }. Find theparametric surface S(u, v) E R3 such thatS(ujj, vii) = Dij^ (2.14)wherem nS(Uij, vij) = E E Bi(uii)Bi(vij)Vrs•r=0 s=0Chapter 2. Mathematical Background^ 20If Equation 2.14 contains all of the M x N data points, surface S is an interpolation of thedata D.In order to generate the equations for the surface approximation problem, each data pointDid must be associated with a corresponding parametric domain point (uii, vii). Defining thisassociation 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 calleda parametric problem if we are attempting to approximate data from a physical object wheredomain information is unknown and therefore must be estimated.2.4.1 Parameterization SchemesParameterization is a difficult problem, but a number of different algorithms exist which havebeen proven effective for univariate and bivariate problems [16] [32]. Different parameterizationsof 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 usdiscuss the different parameterization schemes in the context of curve fitting using a singleparameter. Given a set of data points D2 = (x2 , yi) and a knot sequence {to, ti, , t n } withincreasing order, we fit the data to a curve C(t) = (x(t), y(t)). Without loss of generality, we canassume to = 0, and hi = ti+i where hi forms the knot spacing vector H = (ho, h1 . , hn_1).Uniform ParameterizationThe easiest way to determine hi is to simply set all hi's to be the same constant value, namelyhi = Const, for i = 0, 1, , n — 1. This is called uniform or equidistant parameterization andwas proposed by de Boor [12] and Ahlberg [1]. It is sometimes used because the computationsare simple since time is not needed to compute each hi. However, the amount of computationaltime saved is minor, and curves with loops and cusps can occur even though the definingequations would appear to make that impossible. Examples are cited by de Boor [12] andChapter 2. Mathematical Background^ 21Epstein [13]. This method is too simplistic to cope with most practical situations, with oneexception given by Foley [15], where the uniform parameterization is superior to chord-lengthparameterization. The overall poor performance of the uniform parameterization is attributedto the fact that it "ignores" the geometry of the data points. If we interpret the parameter tof the curve as a time variable, as time passes from time to to time t ri , the point C(t) tracesout the curve from point C(to) to C(tn). With uniform parameterization, point C(t) spendsthe same amount of time between any two adjacent data points, irrespective of their relativedistances. A good analogy is that of a car driving along the interpolating curve. We have tospend the same amount of time between any two data points. If the distance between twopoints is large, we must move with a high velocity. If the next two data points are close to eachother, 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 ParameterizationGiven the physical analogy above, it is perhaps more reasonable to adjust the velocity accordingto the geometric distribution of the data points. One way of achieving this is to have the knotspacing 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 thanuniform knot spacing does, although not in all cases. Foley [15] cites an example where chordlength parameterization does not give the shape users have in mind and wild oscillations occurbecause 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 chordlength will not produce curves with cusps at the data points, thus giving some theoreticaladvantage to the chord length parameterization over the uniform one.Chapter 2. Mathematical Background^ 22Centripetal ParameterizationAnother similar parameterization based on the distance norm was developed by Lee [24] andis called the centripetal model. This parameterization is based on the same paradigm as chord-length parameterization. If we sethi = VtII Did-1 — Di lf,^ (2.16)the resulting motion of a point on the curve will minimize the centripetal force acting on it. Sofor trajectory based designs where we care about machine wear, this is the right candidate.Affine Invariant ParameterizationThe uniform, chord length and centripetal methods are each invariant under rotations, trans-lations 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 thex and y coordinates of the data are scaled differently, the resulting fit is not a simple scaledversion of the original fit. The geometric properties of scale, rotation and translation invari-ance are important in character font applications [31] Affine invariant chord length and angleparameterizations are designed to meet these kinds of practical use [31].2.4.2 SummaryParameterization is a difficult problem that is very application dependent. Usually the inputdata points decide which parameterization method to choose and there is no general methodwhich is suitable for most practical case. (e.g., the case when scale, data point spacing, andangles between adjacent data points all vary largely.) Norm based parameterizations suchas chord length and centripetal parameterization have certain important geometric propertieswhich make them superior to the uniform case. Chord length and centripetal parameterizationare satisfactory for most geometric applications. The problem with these two schemes is thatthey 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^ 23Affine invariant parameterization is a useful property. Uniform parameterization has thisproperty but it does not produce good results when the data point spacings are large andvariable. Affine invariant chord length and angle parameterizations generally provide pleas-ing fitting results. However, both will fail if the data points fall in the one dimension-lesshyperplane [35].2.5 Multigrid Numerical Method2.5.1 Iterative Methods on Linear Algebraic EquationsThe 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 vectorarray 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 theyare extremely inefficient with computational costs of 0(N3). When the matrix A is sparse, theseelimination methods fill up those entries which are originally zero and a lot of storage space iswasted.For many practical applications, iterative methods are used to provide faster solutions. Thespeed of the algorithm depends on the rate of convergence, and since each iteration costs moreor less the same amount of computational time, a rapid convergence rate is crucial. Given aninitial guess of u°, at the r th iteration the iterative procedure calculates the new approximationusing:Ur+1 = CUT ± C (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 alge-braic equations. Among these, Gauss-Seidel and Jacobi relaxation methods are the two typesChapter 2. Mathematical Background^ 24commonly used. The Jacobi method is very similar to the Gauss-Seidel method. The only dif-ference is that the Jacobi method keeps using the solution from iteration r for all variables whilecomputing new values for iteration r + 1. Gauss-Seidel makes use of a new value immediatelyafter it is calculated.Here, we give a brief overview of the formulation and function of the Gauss-Seidel relaxationscheme which is employed throughout this thesis. We denote the matrices representing thediagonal, lower-triangular and upper-triangular parts of A by D, L and U respectively, withL and U having zero elements on the diagonal, thus A = L + D + U. For the case of Gauss-Seidel relaxation, where each entry xi in ur is replaced immediately when 4+ 1 is available, theiteration equation is:(D L)ur+ 1 = d Uu r .The corresponding matrix C and the vector c can be computed as:C =^+ L) -1 U, c = (D L) -1 b(2.19)(2.20)where v is the exact solution and u is the approximate solution to equation 2.17. The errorvector 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 todamp the rapid fluctuations in the residuals is by Gauss-Seidel relaxation scheme. During thefirst few iterations, Gauss-Seidel has very fast convergence, with residuals rapidly decreasingfrom one iteration to the next. But once the high frequency components in the error aresmoothed out, the convergence rate levels off. This inherent computational sluggishness can bestudied from a spatial frequency perspective [37]. A local Fourier analysis [7] shows that theconvergence is fast as long as the residuals have high frequencies compared to the scale of thegrid spacing. As soon as the residuals are smoothed out, convergence slows down. The overallconvergence is asymptotical, and thus slow overall if high accuracy in the solution is required.0^1^2^3^5^6^7^8^9^10 II I^13 14 15^6Chapter 2. Mathematical Background^ 25One solution to this is through the use of multigrid methods.2.5.2 Multigrid PhilosophyOne 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 ofrelaxation on a coarse grid and then use interpolation to get the resulting approximation as aninitial guess on a finer grid. A coarse grid is chosen because relaxation on a coarse grid is lessexpensive than on a fine grid since there are fewer unknowns to be updated. Another reasona coarse grid relaxation is faster is that there is a faster convergence rate on a coarse grid. Ifwe examine the smooth components of the residual on a coarse grid, they appear less smooththan on a grid that is twice as fine as shown in Figure 2.7. This suggests that when relaxationstalls, signalling the predominance of smooth error modes, it is advisable to move to a coarsergrid, on which those smooth error modes appear more oscillatory [8]. The relaxation will thenbe more effective. This is the basic idea of multigrid method.a) Wave ea a N = 16 gridob) Wave on a N = 8 gridPlate 2.7: A wave on C2 h (N = 16) is projected onto 12 2h (N = 8). On the coarse grid, the waveseems more oscillatory than on the fine grid.Chapter 2. Mathematical Background^ 262.5.3 Basic StrategiesSuppose one has a set of grids 1th°, Ohl, , O h" over the same domain C2 (hi is the grid spacingat 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 aretwo strategies which can be used:• Nested Iteration.This strategy incorporates the idea of using coarse grids to obtain a better initial guessfor 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 nh to obtain the approximation Uh2. Compute the residual r = d — Auh.3. Relax residual equation Ae = r on I22h to obtain the approximation to the error e2h4. Correct the approximation obtained on I22h uh uh + Interpolate(e2h ).These two strategies will be adopted in our multigrid method discussed below. They coop-erate with each other and demonstrate great efficiency. First, we will give the various methodsused to compose multigrid method as a whole, which will be called operators.Chapter 2. Mathematical Background^ 272.5.4 Multigrid Operators and ParametersNested iteration requires that the original problem be re-defined on the coarser grid [8]. Oneproblem is that there might still be smooth components in the error when we reach the finestgrid 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 questionsunanswered. The relaxation method (Gauss-Seidel) used in the coarse grid correction schemeneeds to be determined. Different relaxation methods might greatly affect the algorithm'sperformance. We have to define the restriction operator used to transfer the residuals from Ohto C/2h and the prolongation operator which is used to transfer the error estimate back to C2hfrom C22h . The restriction operator from a fine grid to a coarse grid is denoted as gh and theprolongation operator from a coarse grid to a fine grid 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 optimalperformance.In multigrid method, the intralevel processes are usually basic relaxation methods (such asGauss-Seidel relaxation), the interlevel prolongation processes involve local Lagrange interpola-tions, and the restriction processes involve local averaging operations. The exact form of theseoperations are problem-dependent. Fortunately, there are well established standards that havebeen obtained from numerical experiment. For general analysis of the theory, there are detailedreferences by Hackbusch [21] and Brandt [7]. Further discussion concerning the choices for theconcrete problem will be presented in later chapters.2.5.5 Multigrid V-cycleUsing the intergrid transfer operators and necessary parameters, the two-grid coarse grid cor-rection scheme can be refined into a multigrid method. The idea of residual correction on acoarse grid is applied recursively and the correction process is repeated on successively coarserChapter 2. Mathematical Background^ 28grids until a direct solution of the residual equation at the coarsest level is inexpensive. Thecoarse grid correction scheme has been incorporated into the multigrid V-cycle to provide fastnumerical performance.h2h4h8h16hPlate 2.8: The schedule for grid visiting when the multigrid V-cycle is executed. Each level isidentified with its spacing. Multigrid V-cycle starts from the finest level and proceeds to thecoarsest level, then returns to the finest level.Figure 2.8 shows the schedule for the visiting order of the grids. Because of the pattern ofthis diagram, this algorithm is called the V-cycle. The V-cycle scheme has a compact recursivedefinition which is given by:Algorithm 2.1 Multigrid V-Cycle Operator MGVUh 4- MGVh (Uh ,dh )1. Relax vi times on Ahvh = dh with a given initial guess Uh .2. IF Clh = coarsest grid, THEN go to step 4ELSE d2h 4- igh (dh - Ah uh )ugh ,_ 0ugh 4_ m Gv2h (u2h , d2h) .Chapter 2. Mathematical Background^ 293. Correct uh 4- uh + 41 112h,_ m Rh (uh , dh \)4. REPEAT uh ^vo times.5. END.The MR operator is a exact solution operator which relaxes the equations vo times and returnsthe exact solution at the coarsest level.The V-cycle is just one of a family of multigrid cycling schemes. The entire family is calledthe v-cycle method and is defined recursively by the following algorithm:Algorithm 2.2 Multigrid v-Cycleuh <— MG Vh (Uh , dh )1. Relax vi times on A hvh = dh with a given initial guess uh .2. IF Oh = coarsest grid, THEN go to step 4ELSE d2h ,_ gh(dh — Ahuh)u2h +._ 0REPEAT u2h 4_ mGv2h(u2h , d2h) v times.3. Correct uh 4___ uh + p2u2h4. REPEAT uh 4-- MRh (uh , dh) vo times.5. END.In practice, only v = 1 (which gives a V-cycle) and v = 2 are used. Figure 2.9 shows theschedule of grid visiting for v = 2, which is called the W-cycle scheme.Chapter 2. Mathematical Background^ 30h2h4h8hPlate 2.9: The schedule for grid visiting of W-cycle scheme.2.5.6 Full Multigrid V-cycleIn Section 2.5.3, we proposed two distinct, basic multigrid strategies; nested iteration and coarsegrid correction. In the multigrid V-cycle, the coarse grid correction scheme was employed tosmooth out the low frequency error components. Nested iteration involves solving a small scaleproblem on the coarser grid and transferring the solution back to the finer grid to provide theinitial 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 Grid2h4h8h Coarse GridPlate 2.10: Full multigrid grid visiting schedule.Here we have the compact algorithm:Chapter 2. Mathematical Background^ 31Algorithm 2.3 Full Multigrid V-Cycle Operator FMGVvh 4-- FMG Vh (uh , dh )I. IF 12h = coarsest grid, THEN go to step 3.ELSE d2h 4_ gh oh — Ahuh )u2h 4_ 0u2h <-- FMGV2h (u2h, d2h).2. Correct uh ,_ uh + 4ihu2h3. REPEAT uh 4_ mRh(uh , dh \) vo times.4. END.FMG method start from the coarsest grid and proceed to the finest grid. Each V-cycleinside a full multigrid iteration is preceded by a smaller V-cycle designed to provide the bestpossible initial guess. As shown in Section 2.5.7, the extra work done in these preliminaryV-cycles is not only inexpensive, but it generally pays for itself by providing better initial guessfor the finest grid.Full multigrid method is the complete knot into which the many threads of numericalmethods are tied. It is a remarkable synthesis of ideas and techniques which individuallyhave been well known and used for a long time. Taken alone, many of these methods haveserious 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 toremove the limitations. The result is a simple but powerful algorithm.2.5.7 Performance EvaluationMultigrid is a very efficient method for solving a system of linear equations. Both the storageand computational costs are of optimal order.Chapter 2. Mathematical Background^ 32Storage CostConsider 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 sideof Equation 2.17 and the unknowns from the left hand side. So for the storage requirements,the finest grid 1h needs 2MN units, the next coarser grid C2 2h needs 4 as many, etc. If weaccumulate all of the storage units at each level, we have the upper bound on storage which is:S = 2MN(1 + 2 -2 2 -4 + + 2-2n) < 2MN (2.22)1 — 2 -2 — 8 MN.3^*Since we only store the components necessary for carrying out the relaxation, namely theinput and output of the equations, and have no intermediate storage, the storage cost is optimal.Computational CostNow let us consider the computational cost. It is convenient to measure these costs in termsof work units (WU). A work unit is the cost of performing one relaxation sweep on the finestgrid. Equivalently, it is the cost of expressing the fine grid operator. It is customary to neglectthe cost of intergrid transfer operations which could only amount to 15-20 percent of the costof the entire cycle.First consider a multigrid V-cycle with one pre-relaxation and one post-relaxation, namelyvl = v2 = 1. Since the coarsest level has a constant number of relaxations and is only visitedonce, we can simplify by assigning vo = 2. Adding these costs and again using the geometricseries for an upper bound gives:2wu 8 wu^TMGV = 2 ( 1 + 2-2 + 2-4 + + 2-2n) < ^2-2 WU ^•For a full multigrid V-cycle, we just sum up the costs for each single V cycle. Starting fromthe finest grid, the computational cost is bounded by 3 WU as just calculated above. The nextcoarser grid costs 4as much as the finer grid. Then the computational cost for full multigridV-cycle is:^2 „ ^32T tFMGV =,^1+ Z 2 + 2-4 + + 2-2n) < ^^1 — 2 -2)^(1 — 2-2)2 WU = T WU.^(2.24)(2.23)Chapter 2. Mathematical Background^ 33Since each sweep at the finest level is O(MN), the total computational cost of full multigridV-cycle is also O(MN) for a single cycle, which is of optimal order.Convergence RateIn order to determine the final cost for solving a system of equations, we have to know theconvergence rate. The total cost is the cost for each cycle multiplied by the number of cy-cles required to achieve a particular accuracy. The number of cycles depends totally on theconvergence rate for a fixed error tolerance.Convergence rate analysis is a complex area [21] [7]. Local mode analysis [21] provides anon-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:1r = ,^(vi + v2) 2 .(2.25)Due to the preliminary cycling through coarser grids, only 0(1) V-cycles are needed by thetime the algorithm reaches the finest grid. Therefore, the total computational cost of FMGmethod is still 0(MN) [8].Chapter 3Parametric Space Mapping3.1 ParameterizationTensor-product B-splines (the surface formulation used in this thesis) are rectangular. Fig-ure 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 Westedges are mapped to the first column of data and the N and S edges are mapped to the lastrow of data.As stated previously, we wish to approximate a cylindrically digitized object with a B-splinesurface that the boundary of the surface corresponds to the neck of the figure (i.e., the last rowof data). This goal can not be met by mapping a cylindrical surface onto cylindrical data asshown in Figure 3.2a, the mapping from the surface shown in Figure 3.2b to the head will meetthis requirement.After the initial mapping, we also wish to more closely approximate those regions which arearticulated, such as the mouth and the eyes. These areas contain the most interesting featuresand are important for modelling and editing. Finally, we would like the parameterizationalgorithm to be as fast as possible.This goal is achieved in a three step process. First the data is mapped into an intermediateuniform parametric space ST. Then the mapping function M(ui, vi) = (s2, ti) is defined tomap each parametric pair (ui, vj) in UV space into ST space. This mapping function is defined34a) Cylindrical Surface b) Folded SurfaceChapter 3. Parametric Space Mapping^ 35Plate 3.1: Labels for the tensor-product B-spline surface.Plate 3.2: Two approaches to approximating the head data: a) wrap a cylindrical surface aroundthe head, b) fold the four edges of a rectangular surface to cover up the head with the openingat the neck.Chapter 3. Parametric Space Mapping^ 36through 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 mea-sure. At last, the UV space is uniformly sampled and mapped back to data space throughfunction M(ui, vi) in order to form a set of gridded data points D23 in Equation reffittingdef.3.2 Initial ParameterizationImagine that we put a square rubber sheet over Victor Hugo's bust and let the edge of thesheet match the opening at the neck. The first row of the data will be mapped to the center ofthe sheet with the rest of the rows spreading out to the edge (Figure 3.3).last row of data'^r+1r(q+2, r),Neste.(264, 264)q^,q+1q+21originSPlate 3.3: The initial mapping from data space to parametric space. Each row of data ismapped 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 initialparameterization of the data. The ST space sits under the UV space which is used for the finaldata parameterization. The ST space acts as the rubber sheet here. The range of the region isChapter 3. Parametric Space Mapping^ 37from 0 to 2M where M is the number of rows of data. In the case of Victor Hugo, there are 264rows. 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 centerof S and ended at the point where it intersects with the boundary of S (Figure 3.3). Each rayspans 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 segmentdivide it into M sub-segments with equal spacing. Here is the algorithm for calculating theinitial parametric value pair (s, t) in parametric space ST(s, t) for each data point Dad E R2 :Algorithm 3.1 Procedure InitialMapanglestep = 360.0 / N;for column = 1 to N doangle = column * anglestep;for row = 1 to M doif(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^ 38end.end.The (s, t) pair is not the final parametric value used to form the fitting equations. The STspace is an intermediate space where the mapping deformation occurs.3.3 Deforming Parametric Space3.3.1 Deformation FunctionThe deformation takes place in UV space which sits on top of the ST space introduced inSection 3.2. Note that we currently have three different spaces to manipulate, the data spacewhich belongs to R3 , the intermediate space ST and parametric space UV.last row of dataii+1i+21(0, 0)^.1UPlate 3.4: Two coordinate frames imposed on the parametric space.The deformation is defined in terms of a tensor-product bivariate Bernstein polynomial. Weimpose a grid of control points P23 on the square region in ST space (a control point is definedto be the point lying on the control grid and should be distinguished from control vertex whichT(V)^•1^1 1i+2Chapter 3. Parametric Space Mapping^ 39is a spline control vertex). They form m 1 lines in the S direction, and n 1 lines in the Tdirection (Figure 3.5). Their locations are defined to be:xix 528. j 528 Paj Xo "4-^S+^T.m^nwhere 528 is the length and width of the parametric region.S(U)Plate 3.5: A control grid imposed on the parametric space with grid size of m = n = 7. Eachcontrol 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 truefor any data point in R2 . The deformation is achieved by moving the control points from theirundisplaced grid positions. The deformed position X' in ST of an arbitrary point X = (u, v) inUV space is found by first computing its (s, t) coordinate when there is no deformation. Underthis 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 —^ / (528 —i =0 3 =0(3.2)where P23 is a vector containing the Cartesian coordinates of the control point in ST space.(3.1)Chapter 3. Parametric Space Mapping^ 403.3.2 Refinement and Local ControlThe control points here act as the coefficients of a Bernstein polynomial. Since the Bernsteinpolynomial is the basis function for Bezier curves and surfaces, there are meaningful relation-ships between the deformation and the placement of the control points. The parametric regionin ST space can be viewed as a Bezier surface patch defined by these control points. Theorder of the basis function is defined by m and n. Whenever we move a control point, we arereshaping the Bezier surface by moving a control vertex.A Bezier surface has the property of interpolating the boundary points. This is importantbecause we can constrain the boundary of our control grid to the boundary of our parametricregion. Consequently, the movement of interior control points or boundary control points willonly effect the parametric mapping inside the parametric region. Otherwise, the boundaryof our parametric region will not be stable and would not provide a constant size parametricspace for mapping the data. For a B-spline formulation, it is known that the surface doesnot 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 theboundary of a surface patch and constrain the surface points by interpolating the boundarycontrol vertices automatically.Bilinear interpolation is adequate for our particular application which means both m and nEquation 3.2 are chosen as 1. In order to form the chord length parameterization, the controlgrid is refined and the control points are displaced. The following section provides the criteriafor grid refinement and control point displacement.3.3.3 DeformationDiscrete GradientTo construct a chord-length parameterization, the first step is to calculate how much variationthere is in the digitized data.Chapter 3. Parametric Space Mapping^ 41Definition 3.1 The discrete gradient with respect to i in parametric direction gi at data pointDid 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 magnitudeof point Dii is the radial distance from the central axis to the surface. Given the distance zbetween two adjacent sample points, then the distance between two adjacent data points Di_Liand Dig is d = I Z2 + g2 1 . Because of the fine sampling grid from the digitizer, the magnitudeof z is constant and is usually very small compared to the magnitude of the non-zero gi's. Sowe can approximate the above formula as d = gi.This implies that if we direct more parametric lines to those high gradient areas, the datapoints that are farther apart will also be farther apart in parametric space. This is the generalidea of chord length parameterization. If the refinement of the control grid over parametricspace is fine enough, more control points will be generated and displaced. Thus we can achievethe chord length parameterization within a certain tolerance.Refinement and Free Control PointsThe deformation is built starting from a single bilinear patch defined by four control points410, Pd1 and Pill making up the control quadrilateral. These four points are not displacedsince they bound the region being approximated. Thus all points P ik have the same (s, t) pairsin ST space as the corresponding (u, v) pairs in UV space.Each refinement locally subdivides the 1-level control quadrilateral and generates foursmaller control quadrilaterals residing at level 2. At level 1, the criterion to decide whether1 Since the sample points are really fine, we neglect that the surface is curved and take it as fiat.itit itII.Chapter 3. Parametric Space Mapping^ 42or not to refine a control quadrilateral is the discrete gradient sum inside it. First, the overalldiscrete gradient sum SG for all of the data points which cover the whole parametric area Eare calculated. Second, the tolerance threshold T for each control quadrilateral with area A isset asT = SG —AE . (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 thenclicking on the refine button.The refinement is performed in UV space. The mid-point subdivision technique is adoptedfor simplicity. As shown in Figure 3.6 in UV space, the original control quadrilateral is evenlysplit into four equal sub-control quadrilaterals.e2 tab^it ^Ft (11 b10Ur?.(8.)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 atlevel 1.u^Poe( ) Di02^'01Po i( V )Poo^'1)0p2^U )P i01 00- 1- 2^2 ^ PiP10 00^P11—^) plU11^ 00Chapter 3. Parametric Space Mapping^ 43So we have the default coordinates for newly generated control points asp2^U) Poi12 01)p2 = (v ) . 1 '^P 20^'10) p2 _ ( ( U ^±—P11. \ ' 21^ ' 0=( V^1II2^)) p p22 '' 11vuThe (s, t) pairs for these newly generated control points can be derived from Equation 3.2by substituting the (u, v) coordinates calculated above.As shown in Figure 3.6, after refinement there are 9 control points altogether at the secondlevel. The four corner 134's stay at the same positions as the four Pi's and are not be reposi-tioned any more. This leaves 5 new control points that can be moved to deform the parametricspace. We will refer these control points as free control points.Free Control Point DisplacementWe will now discuss how to deform the parametric space by moving the free control point. Asshown in Figure 3.6, we have five control points which can be displaced, namely pro , Pot , ph, PIand P.Chapter 3. Parametric Space Mapping^ 44To displace a free control point, for example Pro , the first thing is to calculate the discretegradient sum along the chord Pj0P1-0• This is done by sampling along the chord PcloPilo andobtaining the discrete gradient value at each point and then summing them up. The samplingrate is constant. This means that the number of points sampled (Ms) is proportional to thelength of the chord PiloPI-0 in ST space. Each sampled point (s, t) is mapped back into discretedata space. This can be assured to be a one-to-one mapping by "inversely" applying Algorithm3.1 since this algorithm is one-to-one. We sum up the discrete gradient values at each sampleddata point Si by:Mss 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 fourdata points Di , Di+i Di+i,j+i, Dij±i. There are three cases for this point to be located inthe 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 thediscrete gradient at the sampled position. For example, we get Di_ f_si j +63 when we map Si todata 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 thediscrete gradient So = 1S is found. Then control point pro is moved to the position 0Chapter 3. Parametric Space Mapping^ 45Dija) Case 1: The sample point is one of the corner points.Dij+ 1DijDijjb) Case 2: The sample point is on an edge ofthe data quadrilateral.c) Case 3: The sample point is in the dataquadrilateral.Plate 3.7: The three cases of the sample point falling into a data quadrilateral.Chapter 3. Parametric Space Mapping^ 46to balance the discrete gradient sum on either side of Ph. Note the this displacement ofPP happens in ST space. Figure 3.8 shows that after Ph is moved toward Plo , the sub-quadrilateral P4PloPliPh will cover less discrete data points in data space than the adjacentsub-quadrilateral noPhilini . In UV space, these two quadrilaterals cover the same amountof parametric space. This means we have more parametric space to cover the high gradientareas 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 balancethe gradient sum along that edge on which it sits.Control point Ph is displaced based on the same idea of discrete gradient sum balance. Aniterative process is designed to compute the approximate point position B = (s, t) such thatthe discrete gradient sums inside each of the four sub-quadrilaterals surrounding B are veryclose to each other. The area discrete grid sum is calculated by expanding Equation 3.6 fromone dimension to two dimensions:Ms NsSarea = EE 95, 3i=1 j=1(3.9)Chapter 3. Parametric Space Mapping^ 47where Ms and Ng are the number of sample points in the U and V parametric directionsrespectively, gs,, is the discrete gradient value at sample point Sij and Sarea is the total discretegradient sum of the area.If we refine the grid of control points enough times, a near chord-length parameterizationis obtained for those areas with a fine control grid. Color plate 2 shows the control grid afterfive levels of refinement.3.4 Sampling Data PointsThe FMG method requires gridded data, but because of the odd configuration of the splinesurface, 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 parametricdomain of the spline (UV space) is sampled on a grid and mapped back onto data spacethrough the intermediate ST space.There are X and Y sampling points in the U and V directions respectively (Figure 3.9). Inpractice the values of X and Y depend on the density of the original data. The step size waschosen as 257 in both U and V to provide a sufficiently dense sampling to catch most of hedetails of the surface and to provide the ratio of 2 : 1 for full multigrid methods. Thus we have_528X = Y = 257 and ustep = vstep = 264'We now have the parameter set {(ujj, vij) = (i x ustep, j x vstep)10 < i < X, 0 < j < Y} tofind 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 usingEquation 3.2. This is a one-to-one mapping. In the same way that we calculated the discretegradient, 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.7and Equation 3.8) was used to calculate the exact value at Di+oi ,j+63 = Zqr from the fourexisting 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 inV^X samplesustepN\Y samplesChapter 3. Parametric Space Mapping^ 48UPlate 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 newdata 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 controlvertices which define the parametric B-spline surface. This will be discussed in the followingchapter.3.5 SummaryA parametric space deformation technique based on image warping was introduced. The tech-nique deforms the isoparametric lines to follow rough regions in the raw data. The deformationwas controlled by control points that form a control grid in parametric space that can be re-fined and displaced. The gradient distribution criterion used to move a control point generatesa near chord-length parameterization. The parametric space after deformation was sampledevenly and then mapped onto the data space. Bilinear interpolation is used to obtain the exactposition of the data point associated with parameter pair.Chapter 4Multi-resolution Surface FittingUsing Multigrid4.1 Deriving Equations and ConstraintsThe remaining problem is to solve the set of linear equations to obtain a surface that interpolatesthe gridded data points {Z12 E R310 < i _< X, 0 < j < Y} The control vertices V are theunknowns and the data points Z are knowns. For each data point Zii, the interpolationcondition can be expressed as an equation in terms of the unknown control vertices. Thenumber of equations should be equal to the number of data points.4.1.1 Data Point DistributionThe data points to be interpolated must now be distributed among the patches to make acorrespondence from a data point to a surface point. As shown in Figure 4.1, if we sampleat the four corners of each patch, X x Y patches provide enough freedom for interpolatingthe 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 boundaryconditions to determine the unique solution of Equation 2.14.The additional equations will be covered in Section 4.1.3.49zon ZinnPmnzoo^Zlo^ zino202 ...112zoiD22 222rZ11/ Ili / P21 z21/Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 50Plate 4.1: The (X + 1) x (Y + 1) data points are mapped onto X x Y patches. Every patchcorner corresponds to a data point.Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 514.1.2 Deriving the Positional Interpolation EquationsFor 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 thesame 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) = [ oUse 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,.; +^ (4.1)i = 1, , m — 1 and j = 1, ,n — 1.According to the mapping from data points to surface points, Pij(1,1) = Zii. We canexpress Equation 4.1 in matrix form as:Av = d^ (4.2)Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 52where v is the vector containing the control vertices from V11 to Vx_iy_i and d is the vectorcontaining the data points transformed by certain control vertices (certain unknown controlvertices are moved to right hand side to make the matrix A banded). The correspondingmatrix A is/ 16 4 4 14 16 4 1 44 16 4 1 44 164 1 16 41 4 1 4 16 41 4 1 4 16 41 4 4 16• /Note this matrix is symmetric, positive definite and thus it is well-conditioned. It is alsodiagonally strong which makes it possible for relaxation scheme to converge. The matrix issparse 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. Theseequations are expressed in (X + 1) x (Y + 1) unknown control vertices.4.1.3 Deriving the Boundary ConditionsIn order to completely determine the unknown control vertices, an additional 2X +2Y equationsare required. Each of these equations corresponds to an extra condition imposed on each of the2X + 2Y boundary data points. The boundary data points are:Zoj j = O , . , n — 1;Zi n i = 0, , M — 1;Zind j = 1, . . . , n;Z, 7 0^i = 1, ... ,7n.The remaining degrees of freedom are determined by specifying the tangents and curvatureat the edge of the surface which interpolates the boundary data points. At each boundary datapoint, the curvature is set to be zero and the surface tangent to the control graph.Chapter 4. Multi-resolution Surface Fitting Using Multigrid^53The 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 = (4Vo.; + 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)The corresponding matrix for these equations is:forl<j<Y-1forl<i<X-1forl<i<X-1forl<j<Y-1(4.4)/ 4 11 4 11 4 11 4 11 4 /Up until now, for the 265 x 265 data points we derived by sampling the parametric space, wehad 265 x 265 linear equations, one for each data point. These 265 x 265 linear equations involve265 x 265 unknown B-spline control vertices which define the final approximating surface. Inthe following section, we will discuss the solution of these equations.4.2 Solving For Control Vertices4.2.1 Full Multigrid MethodThe system of linear equations defining the approximation is solved using the full multigridmethod. The benefits are two-fold:• The excellent numerical behaviour of FMG provides a fast and stable way of solving theseequations.• FMG method provides a multi-resolution surface which can be used to directly create thelevels in a hierarchical spline surface.The FMG method was described in Section 2.5.7 showing that it is an excellent candidatefor a well conditioned sparse matrix on a uniform grid.Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 544.2.2 Full Multigrid OperatorsIn the following several sections, we will discuss the relaxation schemes, coarsest level solver,intergrid operators and the interpolation operator used in our implementation.Relaxation OperatorFrom the matrix corresponding to Equation 4.1, we know that it is a nine point formula 1on 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 colourGauss-Seidel relaxation scheme for the nine point formula because one quarter of the residualscalculated after post-relaxation are zero [21].co o o c o 0v,, ^ a a ^ a0 0^o 0 0 0D ^ • ^ '^^ •o 0 0 C 0 0^ a ^ a a^0 0o^0^0^0^0UPlate 4.3: A illustration of the four-colour relaxation scheme on a uniform grid. First the 0points 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, shownin Figure 4.3. We divide the whole set of data points on the uniform grid into four sub-setsand relax them in turn. The four sub-sets of data points are represented with O, 0, ^ and A1 If we take the control vertices as unknowns, we call an equation composing a matrix a nine point formula ifeach linear equation contains nine unknowns, which are adjacent to each other on a grid.Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 55respectively. First all the 0 points (which contain every other original data point on the gridin 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 datapoints. 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 whichgreatly speeds up the whole multigrid cycle since most of the computation cost is spent onrelaxation.Another relaxation scheme implemented is line Gauss-Seidel relaxation based on the fol-lowing analogy. In parametric space, these gridded control vertices in both u and v directionhave very clear geometric meaning. In the u direction, each "line" of control vertices defines aB-spline curve residing inside the surface which is a blending of curves along u isoparametricdirection with the curves along the v isoparametric direction. The surface is defined by thewhole set of control vertices which can also be considered to be a blending. The surface is calleda tensor product of those single curves. The line Gauss-Seidel relaxation method relaxes the udirection unknowns and then the v direction unknowns "line" by "line". Because the unknownsare coupled closely in both u and v directions, this method is very effective in reducing theresiduals.Coarsest Level OperatorThe coarsest level operator MR (Algorithm 2.1) is Gauss-Seidel relaxation. The explicit equa-tion for relaxation is:vi 1= —liidi — Ek<i+1 - E likvik>i(4.5)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.Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 56Prolongation OperatorThe most common prolongation method is interpolation. A common prolongation operator isbilinear 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 122h into^usereverse square brackets. A stencil defined on the fine grid is centered at a point present incoarse 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 theformula forRestriction OperatorThere are two types of restriction operators that are commonly used, injection and full weight-ing. (In 2-D, injection is defined by vg = v1 2i .) In 2-D the full-weighting operator. In 2-D, itis expressed as2 (v2i 2j+1 ,h^,h^„h_^(„,h .^„,h. „,h^„,h^16 k"2z-1-1,23-1-1 "2i-12j+1^'2i-1,2j-1+^"2i,2j-1^u2i+1,2j -F "2i-12j+1)^4V2hi,2j)^ (4.6)The stencil forms of these operators are as follows:0 0 00 1 00 0 0andrespectively.Thus, for injection, the value at the coarse grid points is just the value at the correspondingfine grid point. Full-weighting uses a weighted average of the values at the corresponding fineChapter 4. Multi-resolution Surface Fitting Using Multigrid^ 57grid point and its nearest neighbours. In our implementation, the full-weighting operator waschosen.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, sothe stencil in this case is:cio 1 t20For a corner point, there are only three neighbours, thus the stencil is:00 s -§- [.0 0 0Interpolation OperatorThe interpolation operator takes a solution at the coarser level and interpolates it to a finerlevel between adjacent V-cycles. As suggested by Brandt [7], bicubic interpolation is a bettercandidate for interpolation operator than bilinear interpolation. Because it introduces less highfrequency errors in the initial guess at a finer level than bilinear interpolation does. Thus, abetter 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:vi(i,j). _ 371 vi-i(i - 3,j) j79 vi--i(i - + 179 vi-i(i - 161 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]:/ 109- 162- 100000009_ 1^9- 167^16 - 172-0^0^081 9^81196 2^162^196271^19^186171 17 TV0^0^0_ 9 _ 1 _ 97.7^16^167000000117709- 162—10 0WChapter 4. Multi-resolution Surface Fitting Using Multigrid^ 58For 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 whichyields171 (ih, jh) = 471-1 (ih — 3h, jh) + 713 V I-1 (ih — h, jh) +171-1 (ih + h, jh)^(4.8)which leads to matrix:/ 0 0^3 _ 1^3 0 164^8^32^64^0 0 0^0^0^0 00 0 9^3^9 30 -0 0320 02^14^16t 06^I- -Q3^719 ' '30 -64^8^32^T Z0 0 0^0^0^0 00 0 0^0^0^0 0Suppose we have X x Y solutions after we complete the V-cycle for the coarser level. Theinterpolation operator will generate (2X — 1) x (2Y — 1) "control vertices" as an initial guess forupper level from the X x Y solutions. Using interpolation, the geometric shape of the surfacedefined by the X x Y control vertices at the coarser level is different than the surface definedby the (2X — 1) x (2Y — 1) control vertices. This suggests that a more suitable alternativefor interpolation is B-spline surface subdivision. (see Chapter 2) Midpoint subdivision willgenerates (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 surfacedefined 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 approximationshape at the finest level, we prefer using the exact shape generated at coarser level as the initialguess at finer level to using the shape interpolated to finer level by bicubic interpolation (thelatter one could add high frequency components). The user will not experience unusual changesin shape which can be introduced by the bicubic interpolation.Now we give a general analysis of the computation cost of both bicubic interpolation andmid-point subdivision. For each new control vertex created by bicubic interpolation, nineneighbouring points are involved. For the subdivision scheme, an even indexed control vertex iscalculated from two neighbouring data points and an odd indexed control vertex is calculatedChapter 4. Multi-resolution Surface Fitting Using Multigrid^ 59from three neighbouring points. Based on the computation cost of each generated controlvertex, subdivision is less expensive. However, bicubic interpolation makes use of all the existingcoarser level points and transfers them directly to finer level in constant time without furthercalculation. Thus bilinear interpolation only needs to compute 4of the total control pointsfor the finer level, decreasing computation cost. The disadvantages and advantages of thismethodology even out in practice. For example, when implementing with a 265 x 265 grid, theinterpolation process costs almost the same amount of time for both methods.4.3 Performance EvaluationFor the Victor Hugo data set, the FMG method was used to solve a system of 257 x 257 linearequations using a nine level grid scheme scheme. The first level has 2 x 2 data points. Atthe finest level, all of the sampled data points Zip are interpolated by a surface consisting of256 x 256 patches.The final result must be a multi-resolution surface definition, thus the approximation resultat 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] tostore the fit result after each smaller v-cycle. This backup process is done each time that thefitting at a certain level is done. This means of a cost of V ( 2 5 7 x 257) storage units. Still, itis 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 inapproximately 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 showsthe time cost of each iteration at every level. Here time t is in hundredths of seconds, and Land I indicate the level index (the smaller, the coarser) and iteration number respectively. Theresidual tolerance is set to be 0.0001, which is very small compared to the general magnitudeof a control vertex which was 200.Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 60L I r t1 11 0.0000013258 152 1 0.4906180920 222 2 0.0276220148 232 3 0.0012658379 232 4 0.0000636812 232 5 0.0000032423 233 1 0.5924855008 983 2 0.0340934001 973 3 0.0018144016 993 4 0.0001093170 973 5 0.0000067578 974 1 0.6936782689 4154 2 0.0331101110 4174 3 0.0018249213 4164 4 0.0001103976 4164 5 0.0000068922 4165 1 0.5606160961 17055 2 0.0383129169 16875 3 0.0019832604 17045 4 0.0001116628 16955 5 0.0000069643 1704Table 4.1: The computation cost at each level in FMG.Chapter 4. Multi-resolution Surface Fitting Using Multigrid^ 61Since FMG pursues a direct solution to the linear equations instead of an approximatesolution like least squares method, noise components in the original data set do not effect theefficiency. However, the algorithm does not have the ability to remove or reduce the effect ofnoise on the approximation result.As for the different relaxation and interpolation techniques used inside the FMG cycle, abrief analysis will be provided for comparison. At the finest level, the time used for line Gauss-Seidel is almost the same as for four-colour Gauss-Seidel relaxation methods, which is about 2seconds for each sweep. Then line Gauss-Seidel method does not smooth the error componentsas fast as the four-colour Gauss-Seidel. When we used line Gauss-Seidel, more sweeps wererequired to meet the same residual tolerance. The reason for this phenomena requires moreinvestigation. As expected, the subdivision interpolation method performs better than thebicubic interpolation method because the former one produces a smaller residual at the finerlevel.Chapter 5Hierarchy BuildingThe FMG approach produces approximation to the data having multiple resolutions. The nextstep is to generate hierarchical B-spline from the results.Lyche and Morken [26] work bottom up by first approximating the surface with a fine meshof spline patches then removing knots in those regions where knot removal will not cause thesurface to move out of tolerance.In the following sections, we are going to present a bottom-up approach for generating thehierarchical representation of the surface that we obtained from the FMG solver.5.1 Inverse Subdivision5.1.1 Hierarchical Surface File FormatThe hierarchical B-spline surface is generated by generating a file defining a hierarchical B-splinesurface from the multi-resolution representation. Recall from Chapter 2 that a hierarchical B-spline surface is very economic in terms of space, we only need to store those non-zero offsetvectors and the first level control vertices as reference vectors. The following script file is atypical 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.62Chapter 5. Hierarchy Building^ 63Surface: 'Plane 0Root Frame: 0 'LinkONumber of Frames: 1Number of Offsets: 35Topology: UO VOSymmetry: 1.500000 VEuler Angle Orientation: 0.000000 0.000000 0.000000Position: 0.000000 0.000000 0.000000Offset DataLevel 0 - spacing 1 4x4UO 0:F 0:-60 0 -601:F 0:-20 0 -602:F 0:20 0 -603:F 0:60 0 -60U1 0:F 0:-60 0 -201:F 0:-20 -25.5 -202:F 0:20 -25.5 -203:F 0:60 0 -20U2 0:F 0:-60 0 201:F 0:-20 -25.5 202:F 0:20 -25.5 203:F 0:60 0 20U3 0:F 0:-60 0 601:F 0:-20 0 602:F 0:20 0 603:F 0:60 0 60Level 4 - spacing 0.0625U24 20:T :29.0296 23.0948 47.4123Chapter 5. Hierarchy Building^ 6428:T :29.0296 —23.0948 47.4123This file contains the surface definition as well as a definition for the coordinate frames foran articulated figure. The header file contains the statistical data about the surface such as thenumber of frames and number of non-zero offsets.5.1.2 Inverse Subdivision SchemeThis section describes an approach that greatly reduces the number of offsets, but does notproduce a hierarchical B-spline surface that is animatable.This scheme is designed to minimize the number of non-zero offsets. We begin from thefinest level / in FMG. At this level, we have the entire definition of the surface and the finestsurface 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-zerooffsets. This idea can be described using a spline curve. Suppose we have a cubic B-splinecurve defined by 11 control vertices (Figure 5.1a) containing 8 curve segments, and we let thiscurve be at level / in the hierarchy. Imagine a spline curve obtained by subdividing a 4-segmentB-spline curve (Figure 5.1b) with the same order, and adjusting the control vertices by offsetreferencing.This 4-segment curve is at level / — 1. Calculating the 7 control vertices defining the 4-2 This minimization is not proven in theory, it is a minimum number easily obtained.V6 •IV7• 1VsI V^I ^V10^11^V9 •••1V4•1V5•1- 11- 11-1V3^V 4•-t1-1V71-1^•V6•I--------------V21-1^••Vi •• ---------------1-Chapter 5. Hierarchy Building^ 65a) A curve defined by 11 control vertices and having 8 segments.• 1-4V5b) 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 levelcurve in a way that reduces the number of non-zero offsets while the hierarchy is built usingsubdivision.0 0 0 0 0 2If a unique solution is desired, weo \ / v l0 "/210 14-1V3,10/- 1V3 V50 1-1V4 0, 1 (5.1)0 v51-1 VC0 /-1V6 "/0 \ V 17-1"/u91V10to the following 7 x 7 matrix:/can reduce the matrixChapter 5. Hierarchy Building^ 66segment curve is an overdetermined problem. The equations look like this:/^1^17^3g 10 80^00^00^0\0^00 0 0 0 0\ / 1-1 \i vi1^ ,/1^0 0 0 0^u23 1 ,1-1R 0 0 0^u3f 3 1^/-10 0^v4R 1 § 1 0^i—i0 g V5 1g 1 /-10 0 1 ^V60 0 0 29: 7 1 \ v 17-1 Iv lV10V11 /(5.2)This matrix is well-conditioned and thus we know that the solution is unique. The resultingseven 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 level1, namely vi, 74, v4 , vs , 43 , vim , vir We term this technique inverse subdivision. So if we usethis idea to build the hierarchy in our bottom-up approach and generate each lower level usinginverse 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 Subdivisionfor level = highestLevel to lowestLevel+1 doInverseSubdivide(level);0 \0/ vil-1 \v21-1/v2/ -10 V3 V40 V4I- 1 V I6 (5.3)0 / -1V5 7/806 V61-1 V101^i k^v71-1 jv Vii(5.4)= 2 ,•••, 7 ./1)3100000 0 101)320000001,34000001)35000001)33 0 0000000where 13i is defined by the recurrence relation:= 1A = 1(6 -Chapter 5. Hierarchy Building^ 67end.SaveBaseLevel(lowestLevel);for level = lowestLevel to finestLevel-1 doSubdivideLevel(level);foreach CV at level+1if(Offset(CV) != 0)SaveNonzeroOffset(Offset(CV));end.end.5.1.3 Remarks: Solution to Triangular Linear EquationA simple symbolic manipulation proves that the inverse of the subdivision matrix is filledeverywhere and all the entries change when the dimension changes.Exploiting the banded structure of the matrix, we can reduce it to a unit upper triangularform having an upper triangle consisting of all zeros excepts for the superdiagonal. This canbe done by elementary row operations. That is, the system is transformed to:Chapter 5. Hierarchy Building^ 68The ilis can be calculated as follows:= 2v 1 14v;-v:_ 1= 2, ... , 7.Now that the transformation of the system of equations is complete, the equations can besolved by back substitution from the bottom up by performing the following computations:1-1V7 = V774-1 = V i -^* Vi+i i = 6, ... , 1.(5.6)5.1.4 Problems with Inverse SubdivisionBecause of the nature of inverse subdivision, the offsets generated by Algorithm 5.1 are notwell behaved. The magnitude of offsets become large and the base level control vertices donot 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 controlvertices 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 whichbehaves well at each level.5.2 Use of MG DataRemember that the multi-resolution data is available from the FMG solver. These multi-levelcontrol vertices define a reasonable shape at each resolution level. Instead of just making useof the finest level result to build the hierarchy, we take advantage of the availability of all levelsof data to build the overlays.We begin from the lowest level, a flat single-patch surface (Plate A.4). The 16 controlvertices VBA4 (Algorithm 2.3) , which are solutions for the coarsest level in FMG, definethe base level of our hierarchical B-spline having a parametric spacing of one. Then the baselevel is subdivided and level 1 control vertices VBAKii are generated. Since we use uniformknot spacing, the parametric spacing for subsequent levels is halved recursively. The FMG data(5.5)Chapter 5. Hierarchy Building^ 69V. at level 1 with compatible dimension is used as the final surface definition at level 1. Thusoffsets at level / are calculated byC4i = VBA4 — Vili^(5.7)where we use Vh to indicate the control vertices generated by refining VBAK ii l at level / — 1.The algorithm is given below:Algorithm 5.2^ Hierarchy BuildingAssignBaseLevel (0) ;for level = 0 to finestLevel doSubdivideLevel (level) ;for = 0 to mfor j = 0 tonoffset = 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 directionsrespectively 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 arewell behaved. Algorithm 5.2 is not intended to reduce the non-zero offsets, as Algorithm 5.1was. Adjacent levels in FMG are just interpolations of points at different resolutions. Nothingin the formulation determines the magnitude of the offset. It turns out that almost all of theoffsets generated with algorithm 5.2 are non-zero. If we store this as a hierarchical B-splineChapter 5. Hierarchy Building^ 70surface, the cost is approximately the same as the cost of storing all of the levels of the FMGsolution. The following section shows how to reduce this cost.5.3 Zeroing Non-zero OffsetsThe FMG approach works well because the solution at one level is a good initial guess for thesolution at the next finer level. When converted into a hierarchical B-spline surface, most ofthe offsets are non-zero but their magnitude is small (Figure 5.2). The challenge is to reducethe 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 tolerancebeginning at the coarsest resolution. Since this operation changes the reference control verticesRid of the next level, the offsets at that level are recalculated to restore the approximation atthat level.100:A : 0.0946581752 0.0206951908 -0.7312056357101:A : 0.1082406374 0.0363397667 -0.1349227312102:A : -0.0510964493 -0.0067193945 0.2410037345103:A : -0.4924907012 -0.0094809132 1.4855188643104:A : -0.0812187087 0.0016742268 0.1157502732105:A : 0.2453602217 0.0028018600 -0.5366043224106:A : 1.0591058826 -0.0002960616 -2.5841572814107:A : 0.2264022796 -0.0016349198 -0.5210445247108:A : -0.1049558748 -0.0005351651 0.3200479296109:A : -1.0981894014 0.0036999245 2.1388922244110:A : 0.2164414669 0.0024895610 -0.4276355013Plate 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 ZeroingAssignBaseLevel(0);Chapter 5. Hierarchy Building^ 71for level = 0 to finestLevel doSubdivideLevel (level) ;for i = 0 to mfor j = 0 to noffset = VBAK Clevel+1] [i] [j] - V [level] [i] [j] ;if (abs (offset) > tolerance * level)SaveNonZeroOff set (offset) ;elseVBAK 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 thezeroing process and the tolerance used. If a relatively large tolerance is used, a large number ofzero 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 causedby offset zeroing at the finest level divided by the number of data points interpolated.) is plottedin Figure 5.4. And the maximum residual is plotted in Figure 5.5.5.4 Smoothing AlgorithmPlate A.9 shows the surface with a tolerance of 0.5 out of 200. Because offset are blindly trun-cated, 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 thesame amount the surface recedes back to its reference shape. A reasonable amount is that of theoffset tolerance, however it can also be specified explicitly by user to meet concrete applicationsChapter 5. Hierarchy Building^ 72Plate 5.3: Offset number generated in hierarchical surface with given tolerance value. Raw dataranges from -100 to +100.Plate 5.4: Residual per data point caused by offset zeroing with given tolerance value. Rawdata ranges from -100 to +100.Chapter 5. Hierarchy Building^ 73Plate 5.5: Maximum Residual caused by offset zeroing with given tolerance value. Raw dataranges from -100 to +100.(Figure 5.7).Here is the complete algorithm incorporating smoothing techniques with the offset zeroingalgorithm.Algorithm 5.4 Offset Zeroing with Smoothing FunctionAssignBaseLevel(0);for level = 0 to finestLevel doSubdivideLevel (level) ;for i = 0 to mfor j = 0 to noffset = VBAK [level+1] [i] [j] - V [level] [i] [j] ;if (abs (offset) > tolerance * level){offset -= off setTolerance ;SaveNonZeroOff set (offset) ;Chapter 5. Hierarchy Building^ 74b) After zeroing offset 2.Plate 5.6: Offset zeroing side effect. a) Before the offset 2 is zeroed. b) The zeroing of offset 2results in a bump.Chapter 5. Hierarchy Building^ 75b) After zeroing offset 2 with smoothing function.Plate 5.7: Offset zeroing with smoothing function. a) Before the offset 2 is zeroed. b) Afteroffset 2 is zeroed with smoothing function.Chapter 5. Hierarchy Building^ 76}elseVBAK [level+1] [i] [j] = V [level] [i] [j] ;end.5.5 Manipulation of Resulting SurfaceBy interactive adjustment of the tolerance and smoothing factor, a compact hierarchical rep-resentation of the original surface is produced. The resulting surface can be used a input tothe Dragon editor for further manipulation and animation. Within the editor, the absoluteoffsets generated by the approximation code are converted into the tangent offsets that allowfine 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 ResultsThis 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 normalmagnitude around 200. The value of each data point originally represented the radius from acentral 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 datapoints) and approximated using the FMG method described above. A 9-level hierarchicalsurface interpolating all the sampled points contains 90494 offsets. With a tolerance of 0.1, thisnumber is reduced to 23270 offsets, about 25% of the sample data set. Smoothing would furtherreduce storage, but the amount of reduction is highly dependent upon how much smoothing isrequired, which differs from situation to situation.end.end.Chapter 5. Hierarchy Building^ 77Another example is a data set of a digitized face. With 66049 data points, we approximateit with 256 x 256 patches at the finest level. With a tolerance of 0.05, the offset number in itshierarchical representation is 19235 (Plate A.14).Chapter 6Conclusions and Future Work6.1 ConclusionsIn this thesis, the digitized Victor Hugo's bust has been used as a testbed of the surfaceapproximation techniques designed to solve the problem stated above. For this specific case,the data parameterization challenge is met by a parametric space deformation scheme. Thisscheme 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 equa-tions have to be solved and a fast algorithm used to solve these equations is desired. Full multi-grid methods (FMG) was chosen to carry out this task because of its optimal performance. TheFMG 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 an-imation, a hierarchical B-spline surface was used to represent the fitted surface. The surfacehierarchy was built using two different techniques. The inverse subdivision offered a good com-pression of the surface storage in terms of non-zero offset number, but the resulting surface wasnot suitable for multi-resolution editing. The direct use of multi-resolution data from the FMGfitter produced a well-behaved shape using slightly more storage. To further compact the re-sulting surface, offset zeroing techniques were employed to remove those non-zero offsets withina 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 offset78Chapter 6. Conclusions and Future Work^ 79zeroing threshold and smoothing factor. The hierarchical B-spline approximation is suitablefor 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 WorkLike 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 adjacentdata points in parametric space. Subsequently, these angles resulted in notches in the surfaceunder chord-length parameterization. Better initial parameterization is required and a para-metric space deforming scheme such as the one presented by Waters and Terzopoulos [39] inimproving the final surface fit. The chosen parameterization scheme, although adequate forthis particular application, is not general enough. Other methods, such as used in [10], will beinvestigated with the goal of using scanned data as a mold that can be applied to an existingsurface to give it shape.Future work will also look into better restriction methods for FMG and at the possibilityof 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 andSurface 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 ComputerGraphics and Geometric Modelling, Morgan Kaufmann Publishers Inc., Palo Alto, Cali-fornia, 1987.[4] Barsky, B. and Greenberg, D. "Determining a Set of B-spline Control Vertices to Generatean 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", Mathematicsof 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 Tech-niques in Computer Aided Geometric Design and Computer Graphics", Computer Graphicsand 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. Lycheand L. Schumaker, Academic Press, Boston, 1989.[11] Coons, S., "Surface Patches and B-spline Approximation", Computer Aided GeometricDesign: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.80Bibliography^ 81[13] Epstein, M., "On the Influence of Parameterization in Parametric Interpolation", SIAMJournal of Numerical Analysis, 13(2):261-268, 1976.[14] Farin, G., Curves and Surfaces for Computer Aided Geometric Design, Academic PressInc., 1988.[15] Foley, T., "Interpolation with Interval and Point Tension Controls Using Cubic Weightedv-Splines", ACM Transactions on Mathematical Software, V.13(1):68-96, 1987.[16] Foley, T. and Nielson G., "Knot Selection for Parametric Spline Interpolation", Mathe-matical 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 An-imation, Univ. of Waterloo Technical Report CS-90-28, (Ph.D. Thesis) 1990.[18] Forsey, D. and Bartels, R. "Tensor Products and Hierarchical Fitting", Proceedings IEEECurves 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 Trans-actions on Graphics.[20] Forsey, D., and Bartels, R. H., "Hierarchical B-spline Refinement", Proceedings of SIG-GRAPH '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 GeometricDesign, V.5:27-31, 1988.[23] Kass, M., Witkin, A. and Terzopolous, D., "Snakes: Active Contour Models", InternationalJournal of Computer Vision, V.4:321-331, 1988.[24] Lee, E., "On Choosing Nodes in Parametric curve Interpolation", SIAM Conference onApplied Geometry, Albany, New York, 1985.[25] Lin, C., Chang, P. and Luh, J., "Formulation and Optimization of Cubic PolynomialJoint Trajectories for Mechanical Manipulation", Proceedings of 21st IEEE Conference onDecision 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 Interpola-tion 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 ofApplied Mathematics, V.29(2):203-220, 1982.[30] Nahas, M., Huitric, H. and Saintourens, M., "Animation of a B-Spline Figure", The VisualComputer, V.3(4), 1987.[31] Nielson, G., "Coordinate Free Scattered Data Interpolation", Topics in MultivariateApproximation: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", Math-ematical Methods in Computer Aided Geometric Design:445-467, edited by T. Lyche andL. Schumaker, Academic Press, Boston, 1989.[33] Schmitt, F., Barsky, B. and Du, W. "An Adaptive Subdivision Method for Surface-Fittingfrom Sampled Data", Proceedings of SIGGRAPH '86:179-188, Dallas, Texas, 1986.[34] Schmitt, F., Maitre, H., Clainchard, A. and Lopez-Krahm, J., "Acquisition and Represen-tation of Real Object Surface Data", SPIE Proceedings Vol. 602 of Biostereometrics '85Conference, 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", Proceed-ings of SIGGRAPH '86:151-160, Dallas, Texas, 1986.[37] Terzopoulos, D., "Image Analysis Using Multigrid Relaxation Methods", IEEE Transac-tions on Pattern Analysis and Machine Intelligence, V.8(2):129-139, 1986.[38] Terzopoulos, D.,"Multi-level Surface Reconstruction", Computer Vision, Graphics and Im-age 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 APlatesThis appendix includes the pictures generated by the implementation as results.Appendix A. Plates^ 84Plate A.1: Victor Hugo raw data.Plate A.2: Deformed UV isoparametric lines in ST space with 5 level refinement.Appendix A. Plates^ 85Plate A.3: Deformed UV isoparametric lines in ST space with 5 level refinement. Data pointsare plotted over the space.Appendix A. Plates^ 86Plate A.4: FMG fit at level 1 and level 2.Plate A.5: FMG fit at level 3 and level 4.Appendix A. Plates^ 87Plate A.6: FMG fit at level 5 and level 6.Plate A.7: FMG fit at level 7 and level 8.Appendix A. Plates^ 88Plate 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 smoothingalgorithm. Offset Number is 3913. Right: Hierarchical B-spline surface with offset zeroingthreshold 0.5 without smoothing algorithm.Appendix A. Plates^ 89Plate A.10: Left: Hierarchical B-spline surface with offset zeroing threshold 0.1 and smoothingalgorithm Offset Number is 23270. Right: Hierarchical B-spline surface with offset zeroingthreshold 1.0 and smoothing algorithm. Offset Number is 1918.Plate A.11: Left: Hierarchical B-spline surface with offset zeroing threshold 2.0 and smoothingalgorithm. Offset Number is 930. Right: Hierarchical B-spline surface with offset zeroingthreshold 5.0 and smoothing algorithm. Offset Number is 318.Appendix A. Plates^ 90Plate A.12: Left: Surface manipulation - Smile 1. Right: Surface manipulation - Smile 2.Plate A.13: Left: Surface manipulation - Smile 3. Right: Surface manipulation - Neanderthalman.Appendix A. Plates^ 91Plate A.14: Another fit to a digitized face.
- 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 |
FileFormat | 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 |
GraduationDate | 1993-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | 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