A dynamic framework for subdivision surfaces

MISSING IMAGE

Material Information

Title:
A dynamic framework for subdivision surfaces
Physical Description:
xv, 168 leaves : ill. ; 29 cm.
Language:
English
Creator:
Mandal, Chhandomay, 1972-
Publication Date:

Subjects

Subjects / Keywords:
Computer vision -- Mathematical models   ( lcsh )
Computer graphics -- Mathematical models   ( lcsh )
Image processing -- Digital techniques   ( lcsh )
Computer and Information Science and Engineering thesis, Ph.D   ( lcsh )
Dissertations, Academic -- Computer and Information Science and Engineering -- UF   ( lcsh )
Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Summary:
Subdivision surfaces are extensively used to model smooth shapes of arbitrary topology. Recursive subdivision on a user-defined initial control mesh generates a visually pleasing smooth surface in the limit. However, users have to carefully select the intiail mesh and.or manipulate the control vertex positions at different levels of subdivision hierarchy to satisfy the functional and aesthetic requirements in the smooth limit surface. This modeling drawback results fromt eh lack of direct manipulation tools for the limit surface. In this dissertation, techniques from physics-based modeling are integrated with geometric subdivision methodology, and a dynamic framework is presented for direct manipulation of the smooth limit surface generated by the subdivision schemes using physics based "force" tools. In a typical subdivision scheme, the user starts with an initial control mesh which is refined recursively using a fixed set of subdivision rules, and a smooth surface is produced in the limit. Most often this limit surface does not have an analytic expression, and hence poses challenging problems in incorporating mass and damping distribution functions, internal deformation energy, forces, and other physical quantities required to develop a physics-based subdivision surface model. In this dissertation, local parameterization techniques suitable for embedding the geometric subdivision surface model in a physics-based modeling framework have been developed. Specific local parameterization techniques have been fully developed for the Catmull-Clark, modified butterfly and the Loop subdivision schemes. Techniques for assigning material properties to geometric subdivision surfaces are derived, and a motion equation for the dynamic model has been formulated using Lagrangian dynamics. Furthermore, advantages of the physics-based deformable models are incorporated into the conventional subdivision schemes, and a dynamic hierarchical control of this model is introduced. Finally, a multiresolution representation of the control mesh is developed and a unified approach for deriving subdivision surface-based finite elements is presented. The proposed dynamic framework enhances the applicability of the subdivision surfaces in modeling applications. It is also useful for hierarchical shape recovery from large range and volume data sets, as well as for non-rigid motion tracking from a temporal sequence of data sets. Multiresolution representation of the initial mesh controlling the smooth limit surface enables global and local editing of the shape as desired by the modeler. This dynamic framework has also been used for synthesizing natural terrain models from sparse elevation data.
Thesis:
Thesis (Ph. D.)--University of Florida, 1998.
Bibliography:
Includes bibliographical references (leaves 159-167).
Statement of Responsibility:
by Chhandomay Mandal.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 030026745
oclc - 41869700
System ID:
AA00012984:00001


This item is only available as the following downloads:


Full Text













A DYNAMIC FRAMEWORK FOR SUBDIVISION SURFACES


By

CHHANDOMAY MANDAL



















A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY


UNIVERSITY OF FLORIDA
































@Copyright 1998
by
Chhandomay Mandal

































To
My Parents,
and My Brother


















ACKNOWLEDGMENTS


I am deeply grateful to my advisor Dr. Baba C. Vemuri, whose guidance, en-

couragement and support made the completion of this work possible. I thank him

for his academic as well as non-academic advice and assistance during my study at

the University of Florida, and for his impressive example of integrity and academic

responsibility. I wish to express sincere appreciation to my co-advisor, Dr. Hong Qin,

for introducing me to the wonders of computer graphics, for patiently guiding me all

the way till the end, and for assisting me in various ways. I would also like to thank

Dr. Sartaj Sahni, Dr. Paul A. Fishwick and Dr. Douglas A. Cenzer for serving in

my supervisory committee and advising me on various aspects of this dissertation.

My dissertation research was supported in part by the NSF grant ECS-9210648

and the NIH grant RO1-LM05944 of Dr. Vemuri, and by the NSF CAREER award

CCR-9702103 and DMI-9700129 of Dr. Qin. I wish to acknowledge Dr. Tim McIn-

erney, Dr. Gregoire Malandain, Dr. Hughes Hoppe, Dr. Kari Pulli, Dr. Dimitry

Goldgof, Dr. Christina Leonard and Dr. Shang-Hong Lai for helping me at various

stages of the work by providing various data sets, software and figures.

I would like to thank the faculty, staff and students of the computer and in-

formation science and engineering department, who were always there to help me.









In my years with the computational vision, graphics and medical imaging group, I

have enjoyed working with a set of bright students, including Yanlin Guo, Yi Cao,

Li Chen, Li Wang, Shuangying Huang, Chinar Kapoor, Fengting Chen, Arun Srini-

vasan, Jundong Liu and Jun Ye. Special thanks goes to my long time roommates,

Raja Chatterjee and Kingshuk Majumdar, for helping me in various ways during my

stay at Gainesville.

Finally, I would like to thank my parents and my brother for their constant en-

couragement and support, for cheering me up during the hard days and for providing

an ever increasing incentive to finish my graduate study. I would have been nowhere

near my goals without them.














TABLE OF CONTENTS




ACKNOWLEDGMENTS ........... .................. iv

LIST OF FIGURES .. ..... ............... .. .... ix

ABSTRACT ............................... ...... xiii

CHAPTERS

1 INTRODUCTION ..................... ........ 1

1.1 Problem Statement ............................ 2
1.2 Proposed Solution ............ ..... ........... 3
1.3 Contributions ............................... 5
1.4 Outline of Dissertation .......................... 8

2 BACKGROUND ................................ 11

2.1 Subdivision Surfaces ........................... 11
2.2 Physics-based Deformable Surface Models ............... 17
2.3 Shape Modeling Using Physics-based Subdivision-surface Model 20
2.4 Shape and Motion Estimation Using Physics-based Subdivision-surface
Model ........... ....................... 22
2.5 Unified Framework for Shape Recovery and Shape Modeling ..... 24

3 DYNAMIC CATMULL-CLARK SUBDIVISION SURFACES ....... 27

3.1 Overview of the Catmull-Clark Subdivision Scheme ......... 27
3.2 Formulation ................. ..... ......... 30
3.2.1 Assigning Patches to Regular Faces .............. 31
3.2.2 Assigning Patches to Irregular Faces .............. 33
3.2.3 Kinematics ............. ................ 36
3.2.4 Dynamics ................... ...... 44
3.2.5 Multilevel Dynamics ............ ........ .... 46
3.3 Finite Element Implementation ... .... .. 49
3.3.1 Data Structures .................. .. ..... .. 49
3.3.2 Normal Elements .................. .... .. 51
3.3.3 Special Elements ........................ .. 51









3.3.4 Force Application .................. ....... ..52
3.3.5 Discrete Dynamic Equation .................. ..54
3.3.6 Model Subdivision .......... ........ 55
3.4 Generalization of the Approach ...... .. 56

4 DYNAMIC BUTTERFLY SUBDIVISION SURFACES ... 57

4.1 Overview of the (Modified) Butterfly Subdivision Scheme ...... 57
4.2 Formulation ............. 61
4.2.1 Local Parameterization ... ..62
4.2.2 Dynamics ................ 72
4.2.3 Multilevel Dynamics ...................... .. 74
4.3 Finite Element Implementation ..... ... 76
4.3.1 Elemental Mass and Damping matrices .... 77
4.3.2 Elemental Stiffness Matrix ..... .. 79
4.3.3 Other Implementation Issues ... ..... 81
4.4 Generalization of the Approach ... ..... 82

5 UNIFIED DYNAMIC FRAMEWORK FOR SUBDIVISION SURFACES 83

5.1 Overview of the Unified Approach .. .. 83
5.2 Unified Dynamic Framework for Catmull-Clark Subdivision Surfaces 85
5.2.1 Local Parameterization ... ... 88
5.2.2 Dynamics ........ . 92
5.2.3 Finite Element Implementation ..... 92
5.3 Unified Dynamic Framework for Loop Subdivision Surfaces 95
5.3.1 Local Parameterization ............ .... .. .. 97
5.3.2 Dynamics ............. ..... .... ......... 100
5.3.3 Finite Element Implementation .... 100
5.4 A General Outline of the Framework for Approximating Subdivision
Schemes ............... .................. 101
5.5 A General Outline of the Framework for Interpolatory Subdivision
Schemes ............ ............ 102

6 MULTIRESOLUTION DYNAMICS .................. .. ..103

6.1 Overview of Multiresolution Analysis and Wavelets .... 109
6.2 Multiresolution Analysis for Surfaces of Arbitrary Topology 114
6.2.1 Nested Spaces through Subdivision ... 116
6.2.2 Inner Product ................ ........ 118
6.2.3 Biorthogonal Surface Wavelets on Arbitrary Manifold 119
6.3 Multiresolution Representation ... ... .. 125
6.4 Dynamics ................................ .. 128
6.5 Implementation Details .................. ....... .. 129









7 SYSTEM ARCHITECTURE ........... ............... 133

7.1 Topological Information Processing Module ........... 135
7.2 Subdivision Module .................. .. ....... .. 135
7.3 Finite Element Analysis Module ... 135
7.4 Force Synthesis Module .... 136
7.5 Update Engine ................. 136
7.6 Display Module ............ ...... 137

8 APPLICATIONS ............... 138

8.1 Geometric Modeling .... . ..... 138
8.2 Shape Recovery from Range and Volume Data ... 142
8.3 Non-rigid Motion Estimation. ....... 148
8.4 Multiresolution Editing .................. ..... .. 150
8.5 Natural Terrain Modeling .................. ... ..151

9 CONCLUSIONS AND FUTURE WORK. ...... ..... 155

9.1 Conclusions ............... ....... 155
9.2 Future Directions .......... . 156
9.2.1 Automatic Change of Topology .... 156
9.2.2 Local Refinement ........ ...... 157
9.2.3 Automatic Model Parameter Selection .... 157
9.2.4 Constraint Imposition ..................... .. 157
9.2.5 Recovery of Sharp Features ... ..... 158
9.2.6 Automatic Model Initialization 158
9.2.7 Improved Synthesized Force Fields .. 158

REFERENCES .................. ... 159

BIOGRAPHICAL SKETCH .................. .......... .. 168














LIST OF FIGURES


2.1 Refinement of an initial control mesh to obtain the limit surface.
(Courtesy : H. Hoppe) .................. ....... 12

3.1 Catmull-Clark subdivision : (a) initial mesh, (b) mesh obtained after
one step of Catmull-Clark subdivision, and (c) mesh obtained after
another subdivision step. .................. ...... ..28

3.2 A rectangular mesh and its limit surface consisting of 4 bicubic surface
patches. .................................. 32

3.3 A mesh with an extraordinary point of degree 3 and its limit surface. 34

3.4 Local subdivision around the extraordinary point and the limit surface. 35

3.5 Local subdivision around the extraordinary point and the correspond-
ing patches in the limit surface from different levels of subdivision. .. 43

3.6 The data structure used for dynamic subdivision surface implementation. 50

3.7 Model subdivision to increase the degrees of freedom : (a) evolution
of the initial mesh, and (b) the corresponding limit surface evolution
perceived by the user. ...... ........ 55

4.1 (a) The control polygon with triangular faces; (b) mesh obtained after
one subdivision step using butterfly subdivision rules. ... 58

4.2 (a) The contributing vertices in the j-th level for the vertex in the j+1-
th level corresponding to the edge between v and v2; (b) the weighing
factors for different vertices. .. ...... 58

4.3 (a) The weighing factors of contributing vertex positions for an edge
connecting two vertices of degree 6; (b) the corresponding case when
one vertex is of degree n and the other is of degree 6. ... 60

4.4 The smoothing effect of the subdivision process on the triangles of the
initial mesh. ................... ............... 62









4.5 Tracking a point x through various levels of subdivision : (a) initial
mesh, (b) the selected section (enclosed by dotted lines) of the mesh
in (a), after one subdivision step, (c) the selected section of the mesh
in (b), after another subdivision step. 64

4.6 Different set of control points at a subdivided level obtained by apply-
ing different subdivision matrices on a given set of control points in a
coarser mesh. ............. .... 67 68

4.7 (a) An initial mesh, and (b) the corresponding limit surface. The do-
mains of the shaded elements in the limit surface are the corresponding
triangular faces in the initial mesh. The encircled vertices in (a) are
the degrees of freedom for the corresponding element. ... 75

5.1 A control mesh with an extraordinary vertex of degree 5 and the corre-
sponding limit surface : (a) using the concepts developed in Chapter 3,
where the limit surface consists of quadrilateral normal elements and
a pentagonal special element; (b) using the unified concept developed
in this chapter, where the limit surface consists of one single type of
quadrilateral finite element. ................... 86

5.2 In Catmull-Clark subdivision, each non-boundary quadrilateral face in
the control mesh has a corresponding quadrilateral patch in the limit
surface : (a) control mesh, (b) limit surface. ... 87

5.3 (a) The marked 16 control vertices define the shaded quadrilateral
patch associated with the shaded regular face in the control mesh; (b)
the marked 14 control vertices define the shaded quadrilateral patch
associated with the shaded irregular face in the control mesh.. 91

5.4 (a) The control polygon with triangular faces; (b) mesh obtained after
one subdivision step using Loop's subdivision rules. ... 96

5.5 An initial mesh and the corresponding limit surface obtained using
Loop's subdivision rules. The domains of the shaded triangular patches
in the limit surface are the corresponding triangular faces in the initial
mesh. The encircled vertices are the control vertices for the corre-
sponding triangular patch in the limit surface. .... 98

5.6 Each triangular patch in the limit surface can be associated with a
non-boundary triangular face in the initial mesh, which in turn can be
parameterized over a triangle with vertices at (0, 0), (1, 0) and (0, 1). 99

6.1 Schematic block diagram of the multilevel dynamics approach. 104









6.2 Multilevel dynamics vs. multiresolution dynamics. ... 106

6.3 Representation of the degrees of freedom in multilevel dynamics and
multiresolution dynamics approach. ..... 108

6.4 The filter bank. .............. 112

6.5 Subdivision refinement of a tetrahedron. .. ...... 119

6.6 Wavelet construction. ................ ...... ..123

7.1 System architecture. ......................... .. 134

8.1 (a), (b), (c) and (d) : Initial shapes (obtained applying Catmull-Clark
subdivision rules on control meshes); (e), (f), (g) and (h) : the corre-
sponding modified shapes after interactive force application. 139

8.2 (a), (b), (c) and (d) : Initial shapes (obtained applying modified but-
terfly subdivision rules on control meshes); (e), (f), (g) and (h) : the
corresponding modified shapes after interactive sculpting via force ap-
plication ................... .. .............. 141

8.3 (a) Range data of a bulb along with the initialized model, (b) an
intermediate stage of evolution, and (c) the fitted dynamic Catmull-
Clark subdivision surface model. ..143

8.4 (a) Range data of an anvil along with the initialized model, (b) an
intermediate stage of evolution, and (c) the fitted dynamic Catmull-
Clark subdivision surface model. ..... 144

8.5 (a) Range data of a head along with the initialized model, (d) the
fitted dynamic butterfly subdivision model, and (c) visualization of
the shape from another view point. . 144

8.6 (a) A slice from a brain MRI, (b) initialized model inside the region
of interest superimposed on the slice, (c) the fitted model superim-
posed on the slice, and (d) a 3D view of the dynamic Catmull-Clark
subdivision surface model fitted to the cerebellum. .. ... 146

8.7 (a) Data points identifying the boundary of the caudate nucleus on
a MRI slice of human brain, (b) data points (from all slices) in 3D
along with the initialized model, and (c) the fitted dynamic butterfly
subdivision model. .......................... .. 146









8.8 Snapshots from the animation of canine heart motion over a cardiac
cycle using the dynamic butterfly subdivision model. ... 149

8.9 (a) The initialized model along with the data set; (b) the fitted model
with two subdivisions on the initial mesh, along with attached springs
for editing. The model after editing (c) at lower resolution, (d) at the
same resolution of the fitted model, and (e) at higher resolution. 150

8.10 (a) Discrete elevation data set (4096 points), (b) fitted dynamic but-
terfly subdivision surface model with 841 vertices (without noise ad-
dition), and (c) fitted dynamic subdivision surface model with 841
vertices when Gaussian noise is added. 153

8.11 Synthesized natural terrain of different roughness using the dynamic
butterfly subdivision surface model with 841 vertices from a data set
of 3099 elevation values. ........ 154













Abstract of Dissertation
Presented to the Graduate School of the University of Florida
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy



A DYNAMIC FRAMEWORK FOR SUBDIVISION SURFACES

By

Chhandomay Mandal

December 1998


Chairman: Dr. Baba C. Vemuri
Cochairman : Dr. Hong Qin
Major Department: Computer and Information Science and Engineering


Subdivision surfaces are extensively used to model smooth shapes of arbitrary

topology. Recursive subdivision on a user-defined initial control mesh generates a

visually pleasing smooth surface in the limit. However, users have to carefully select

the initial mesh and/or manipulate the control vertex positions at different levels

of subdivision hierarchy to satisfy the functional and aesthetic requirements in the

smooth limit surface. This modeling drawback results from the lack of direct manip-

ulation tools for the limit surface. In this dissertation, techniques from physics-based

modeling are integrated with geometric subdivision methodology, and a dynamic

framework is presented for direct manipulation of the smooth limit surface generated

by the subdivision schemes using physics-based "force" tools.









In a typical subdivision scheme, the user starts with an initial control mesh

which is refined recursively using a fixed set of subdivision rules, and a smooth sur-

face is produced in the limit. Most often this limit surface does not have an analytic

expression, and hence poses challenging problems in incorporating mass and damp-

ing distribution functions, internal deformation energy, forces, and other physical

quantities required to develop a physics-based subdivision surface model. In this

dissertation, local parameterization techniques suitable for embedding the geometric

subdivision surface model in a physics-based modeling framework have been devel-

oped. Specific local parameterization techniques have been fully developed for the

Catmull-Clark, modified butterfly and the Loop subdivision schemes. Techniques

for assigning material properties to geometric subdivision surfaces are derived, and

a motion equation for the dynamic model has been formulated using Lagrangian

dynamics. Furthermore, advantages of the physics-based deformable models are in-

corporated into the conventional subdivision schemes, and a dynamic hierarchical

control of this model is introduced. Finally, a multiresolution representation of the

control mesh is developed and a unified approach for deriving subdivision surface-

based finite elements is presented.

The proposed dynamic framework enhances the applicability of the subdivision

surfaces in modeling applications. It is also useful for hierarchical shape recovery

from large range and volume data sets, as well as for non-rigid motion tracking from

a temporal sequence of data sets. Multiresolution representation of the initial mesh

controlling the smooth limit surface enables global and local editing of the shape as









desired by the modeler. This dynamic framework has also been used for synthesizing

natural terrain models from sparse elevation data.

















CHAPTER 1
INTRODUCTION



Generating smooth surfaces of arbitrary topology poses a grand challenge for

the computer graphics, computer-aided geometric design and scientific visualization

researchers. The existing techniques for modeling smooth complex shapes can be

broadly classified into two distinct categories namely, (1) explicit patching using

parametric surfaces and (2) subdivision surfaces.

The explicit patching technique involves partitioning the model surface into a

collection of individual parametric surface patches. Adjacent surface patches are

then explicitly stitched together using continuity constraints. This explicit stitching

process is very complicated for modeling smooth surfaces of arbitrary topology due

to the continuity constraints which need to be satisfied along the patch boundaries.

Subdivision surfaces are simple procedural models which offer an alternate rep-

resentation for the smooth surfaces of arbitrary topology. A typical recursive subdi-

vision scheme produces a visually pleasing smooth surface in the limit by repeated

application of a fixed set of refinement rules on an user-defined initial control mesh.

The initial control mesh is a simple polygonal mesh of the same topological type

as that of the smooth surface to be modeled. At each step of subdivision, a finer

polygonal mesh with more vertices and faces is constructed from the previous one via









the refinement process, and the smooth surface is obtained in the limit. However, a

few subdivision steps on the initial mesh generally suffice to approximate the smooth

surface for all practical purposes. Various sets of rules lead to different subdivision

schemes which mainly differ in the smoothness property of the resulting limit surface

and/or the type of initial mesh (i.e. triangular, quadrilateral etc.) chosen.

1.1 Problem Statement

To model a smooth surface of arbitrary topology using subdivision surfaces, first

a polygonal mesh of same topology needs to be chosen as the initial mesh. This ini-

tial mesh, also known as the control mesh, is refined via recursive subdivision using

a fixed set of rules, and a smooth surface of the desired topology is obtained in the

limit. However, users have to carefully select the initial mesh and/or manipulate

the control vertex positions at different levels of subdivision hierarchy in order to

satisfy the functional and aesthetic requirements on the smooth limit surface. For

example, to obtain a desired effect on the smooth limit surface, it might be necessary

to reposition a handful of vertices in the mesh obtained after one subdivision step,

or a large number of vertices might need to be moved in the mesh produced after

three subdivision steps! This process is not intuitive and at best laborious. Despite

the presence of a variety of subdivision schemes in the computer graphics and geo-

metric modeling literature, there is no direct and natural way of manipulating the

limit surface. The current state-of-the-art only permits the modeler to interactively

obtain the desired effects on the smooth surface by kinematically manipulating the

vertex positions at various levels of the subdivision hierarchy. In this dissertation, the









challenging problem of directly manipulating the smooth limit surface at arbitrary

locations/areas is addressed and a novel solution is presented where the modeler can

not only manipulate the smooth limit surface directly but also control the extent of

the manipulation effect (i.e., global or local manipulation) on the limit surface.

Subdivision surfaces have also been used for recovering shapes from a given set of

points in 3D. However, most of the existing subdivision surface-based shape recovery

techniques resort to complex algorithms to derive a mesh for the underlying shape,

and then mesh optimization techniques are used to obtain a compact representation

of the same. Nevertheless, this process yields a control mesh of the smooth subdi-

vision surface representing the underlying shape that typically uses a large number

of degrees of freedom (control vertices) for representation. In this dissertation, an

efficient hierarchical method of recovering shapes from range and volume data sets is

proposed where the control mesh of the smooth limit surface will use very few degrees

of freedom for representation.

1.2 Proposed Solution

Physics-based modeling techniques offer a potential solution to the problem

of directly manipulating the smooth limit surface generated by the recursive sub-

division procedure. In the physics-based modeling paradigm, a deformable model

is derived by assigning mass, damping, stiffness and other physical properties to a

surface model. The model is deformed by applying synthesized forces and this de-

formation is governed by physical laws. Now, if the purely geometric subdivision

surface models can be embedded in a physics-based modeling paradigm, then the









smooth limit surface can be deformed directly by using physics-based "force" tools.

However, this procedure-based surface model obtained through the subdivision pro-

cess does not have a closed-form analytic formulation in general, and hence poses

challenging problems to incorporate mass and damping distribution functions, inter-

nal deformation energy, forces, and other physical quantities required to develop a

physics-based model. Techniques of locally parameterizing the smooth limit surface

generated by various subdivision schemes are proposed in this dissertation. Once a

suitable local parameterization scheme is developed for a specific subdivision scheme,

a dynamic framework is provided for the corresponding subdivision scheme where

the modeler can directly manipulate the smooth limit surface by using synthesized

forces. At the same time, a dynamic framework for the wavelets derived using subdi-

vision schemes will assist in adopting a multiresolution representation of the control

mesh defining the smooth limit surface, and the modeler can control the extent of

manipulation effect by choosing a desired level of editing. Synthesized force applica-

tion at a lower resolution will yield a global effect whereas manipulations at a finer

resolution will have localized effects on the limit surface. The motion of this physics-

based deformable subdivision surface model is governed by a second-order differential

equation, which is solved numerically using the finite element method. New types of

finite elements for the chosen subdivision scheme are also presented for representing

the smooth limit surface.

The proposed dynamic framework for the subdivision surfaces provides an effi-

cient solution to the shape recovery problem as well. A simple subdivision surface









model with very few vertices in the control mesh can be initialized (positioned) fully

inside a given set of points in 3D. The initialized model will be deformed by applying

forces synthesized from the given data points. When an equilibrium is obtained, the

number of vertices in the control mesh can be increased via a subdivision step on the

current control mesh thereby increasing the degrees of freedom for model representa-

tion, and a new equilibrium with a better fit to the given data set can be obtained.

This process can be repeated till a prescribed error in fit is achieved. Similar approach

can be taken for shape recovery from volume data sets as well where a different type

of synthesized force needs to be specified. The hierarchical shape recovery process

ensures a compact representation of the recovered smooth limit surface using very

few degrees of freedom.

1.3 Contributions

In this dissertation, techniques from physics-based modeling are integrated with

geometric subdivision methodology to present a scheme for directly manipulating

the smooth limit surface generated by the subdivision process. As a result, unlike

the existing geometric solutions that only allow the operations on control vertices,

the proposed methodology and algorithms permit the user to physically modify the

shape of subdivision surfaces at desired locations via application of forces. This gives

the user a "virtual" clay/play-dough modeling environment. The proposed model

can be edited directly in a hierarchical fashion using synthesized forces. Also, this

physics-based subdivision surface model efficiently recovers shapes as well as non-

rigid motions from large range and volumetric data sets. Note that this dissertation









neither proposes a new subdivision technique nor provides a different interpretation

of any existing subdivision technique, but integrates the advantages of subdivision

surface-based and physics-based modeling techniques to solve important theoretical

and practical problems. Although the principles of physics-based modeling are well

understood by the graphics experts and modeling researchers, this dissertation will

greatly advance the state of the art in physics-based shape modeling due to the

contributions listed below.


Local parameterization techniques for the smooth limit surface generated by

various subdivision schemes are systematically derived in a hierarchical frame-

work, and subsequently the initial control polyhedron can be viewed as the

parametric domain of the physics-based smooth limit surface.


The smooth dynamic subdivision surface in the limit is treated as a "real"

physical object represented by a set of novel finite elements. The basis (shape)

functions of these new variety of finite elements are derived using the subdivision

schemes. The proposed finite element methods will prove to be useful not only

in the areas of computer graphics and geometric design, but also in engineering

analysis as well.


The subdivision techniques are used to create a surface model that incorporates

mass and damping distribution functions, internal deformation energy, forces,

and other physical quantities. The motion equations are also systematically

derived for this dynamic subdivision surface model whose degrees of freedom are

the collection of initial user-specified control vertices. Therefore, the advantages









of both the physics-based modeling philosophy and the geometric subdivision

schemes are incorporated within a single unified framework.


Users will be able to manipulate this physics-based model in an arbitrary region,

and the model will respond naturally (just like a real-world object would) to

this force application. This shape deformation is quantitatively characterized

by the time-varying displacement of control points that uniquely define the

geometry of the limit surface.


The dynamic framework for wavelets derived using subdivision schemes en-

ables a multiresolution representation of the control mesh defining the smooth

limit surface. This provides additional flexibility to the physics-based modeling

framework since the modeler is free to choose the desired editing level. If a

lower resolution editing level is chosen, the synthesized force application on the

smooth limit surface will have global effects, whereas editing with force-based

tools at a finer resolution will yield a localized effect.


Algorithms and procedures are developed which approximate the proposed new

finite elements using a collection of linear and/or bilinear finite elements subject

to the implicit geometric constraints enforced by the subdivision rules. This

hierarchically-structured approximation can satisfy any user-specified error tol-

erance.


The proposed dynamic framework enhances the applicability of subdivision sur-

face models in various application areas. It provides a direct and intuitive way of









manipulating shapes in a hierarchical fashion for geometric modeling applications.

It has also been successfully used for efficient and hierarchical shape recovery from

range and volume data sets as well as for tracking a shape of interest from a time

sequence of range or volume data sets. Finally, the dynamic framework of subdivision

surfaces is combined with a variant of a fractal surface synthesis technique to present

a novel natural terrain modeling method.

1.4 Outline of Dissertation

The rest of the dissertation is organized as follows :

Chapter 2 contains an overview of the subdivision surfaces along with a review of

the related literature. The motivation for embedding the subdivision surface models

in a physics-based modeling framework is discussed, and the proposed model is com-

pared with the existing physics-based models. The advantages of a unified framework

for shape modeling and shape recovery are also pointed out in this chapter.

Chapter 3 provides a dynamic framework for the Catmull-Clark subdivision

scheme, which is one of the most popular subdivision techniques for modeling compli-

cated objects of arbitrary genus. Analytic formulation of the limit surface generated

by the Catmull-Clark subdivision scheme is derived and the "physical" quantities

required to develop the dynamic model are introduced [57, 60, 76]. The governing

dynamic differential equation is derived using Lagrangian mechanics and is imple-

mented using a finite element technique.









In Chapter 4, a dynamic framework is provided for the butterfly subdivision

scheme, another popular subdivision technique for modeling smooth surfaces of arbi-

trary topology. A local parameterization scheme for the butterfly subdivision surface

is derived in a hierarchical style. The physics-based butterfly subdivision model is

formulated as a set of novel finite elements which are optimally approximated by

a collection of standard finite elements subjected to implicit geometric constraints

[58, 61].

Chapter 5 presents an unified approach for providing a dynamic framework for

subdivision surfaces in general. In particular, it has been shown that the limit surface

obtained using any subdivision scheme can be viewed as a collection of either quadri-

lateral or triangular finite elements whose basis (shape) functions can be derived

using the chosen subdivision scheme [59].

In Chapter 6, a dynamic framework is presented for the subdivision surface-

based wavelet schemes. This dynamic framework enables a multiresolution represen-

tation of the control mesh defining the smooth limit surface and consequently, the

modeler can control the extent of the effect of manipulation on the limit surface by

choosing a proper editing level.

In Chapter 7, a system that integrates the implementation of the concepts pro-

posed in earlier chapters is presented. Several modules comprise the entire system,

and the functionality of these modules are discussed in this chapter.





10


The proposed modeling framework is used in various application areas and the

results are presented in Chapter 8. The proposed dynamic model has been success-

fully used in geometric modeling, shape recovery and non-rigid motion estimation

applications. A novel technique for natural terrain modeling is also presented where

the dynamic framework is combined with a variant of a fractal surface synthesis

technique.

Finally, conclusions are drawn in Chapter 9 where future directions of research

are also pointed out.
















CHAPTER 2
BACKGROUND


In this chapter, a brief overview of the subdivision surfaces is presented followed

by a review of the previous work done in the area of subdivision surfaces. This is

followed by an overview of the physics-based models along with a review of the re-

lated literature. The motivating factors for embedding geometric subdivision surface

models in a physics-based modeling paradigm are presented in the next section. The

advantages of the proposed dynamic models over the existing physics-based models

for shape recovery are discussed and finally advantages of a unified framework for

shape modeling and shape recovery are also presented.

2.1 Subdivision Surfaces

The input to any subdivision scheme is an initial mesh (also known as control

mesh), MO = (Vo, Fo), which is a collection of vertices V and a collection of faces

Fo. The subdivision surface S(MO) associated with the initial mesh MO is defined

as the limit of the recursive application of the refinement R as shown in Fig.2.1.

The refinement R, when applied on a mesh Mk = (VI, Fk), creates a refined mesh

Mk+ l (Vk+l, Fk+1) where the vertices in Vk+' are computed as affine combinations

of the vertices in Vk and the faces in Fk+1 are obtained by splitting each face in Fk

into a fixed number of sub-faces.




























(b) Meh M = R(M)


(c) Mesh M = R(R(M)) (d) Limit snlrfax S(M) = R(M)

Figure 2.1. Refinement of an initial control mesh to obtain the limit surface. (Cour-
tesy : H. Hoppe)









It may be noted that the vertices introduced through the subdivision process

have a fixed degree in general (4 in case of quadrilateral meshes and 6 in case of tri-

angular meshes). The vertices which do not have this degree are called extraordinary

vertices. The number of extraordinary vertices are very few as these vertices are not

introduced via subdivision refinement in general. For example, in the Catmull-Clark

subdivision scheme (defined on quadrilateral meshes) the number of extraordinary

vertices does not change after the first subdivision on the initial mesh, whereas in the

modified butterfly subdivision scheme (defined on triangular meshes) the subdivision

process never introduces an extraordinary vertex. The limit surface defined by the

subdivision process is most often C (first derivative continuous) or C' (second deriva-

tive continuous) depending on the subdivision rules expect at very few extraordinary

points corresponding to the extraordinary vertices in the mesh. Specific subdivision

schemes along with the properties of their limit surfaces will be discussed in later

chapters.

A lot of research have been carried out on subdivision surfaces, which are mainly

focused either on the development of a new subdivision technique or on analyzing

the smoothness properties of the limit surface generated by a subdivision scheme.

However, in this dissertation the focus is entirely different, namely, embedding the

subdivision surfaces in a physics-based modeling paradigm to enhance the applicabil-

ity of these surface models. In the rest of this section, these "traditional" techniques

of subdivision surfaces are reviewed.









Chaikin [14] first introduced the concept of subdivision to the computer graphics

community for generating a smooth curve from a given control polygon. During the

last two decades, a wide variety of subdivision schemes for modeling smooth surfaces

of arbitrary topology have been derived following Chaikin's pioneering work on curve

generation. In general, these subdivision schemes can be categorized into two distinct

classes namely, (1) approximating subdivision techniques and (2) interpolatory sub-

division techniques. The limit surface obtained using an approximating subdivision

technique only approximates the initial mesh, whereas the limit surface interpolates

the initial mesh in case of interpolatory subdivision techniques.

Among the approximating schemes, the techniques of Doo and Sabin [27, 83]

generalize the idea of obtaining uniform biquadratic patches from a rectangular con-

trol mesh. This scheme is an example of vertex subdivision scheme where a vertex

surrounded by k faces is split into k sub-vertices, one for each face. Doo and Sabin's

subdivision technique can be applied to any mesh of arbitrary topology, and the re-

sulting smooth limit surface would be biquadratic B-splines, except at extraordinary

points corresponding to the extraordinary vertices in the control mesh.

Catmull and Clark [10] developed a method for recursively generating a smooth

surface from a polyhedral mesh of arbitrary topology. The Catmull-Clark subdivision

surface, defined by an arbitrary initial mesh, can be reduced to a set of standard

bicubic B-spline patches except at a finite number of degenerate points. This scheme

is an example of face subdivision scheme, where a k-sided face is split into k sub-

faces. The details of this subdivision scheme are discussed in Section 3.1. Halstead









et al. [38] proposed an algorithm to construct a Catmull-Clark subdivision surface

that interpolates the vertices of a mesh of arbitrary topology.

Loop [51] presented a subdivision scheme based on the generalization of quartic

triangular B-splines for triangular meshes. Loop's subdivision rules are presented in

Section 5.3. Hoppe et al. [40] extended Loop's work to produce piecewise smooth

surfaces with discontinuities in selected areas of the limit surface. They introduced

new subdivision rules that allow for sharp features such as creases and corners in the

limit surface.

Peters and Reif [73] proposed a simple subdivision scheme for smoothing poly-

hedra. Their refinement rules yield a C1 surface that has a piecewise quadratic pa-

rameterization except at isolated extraordinary points. Most recently, non-uniform

Doo-Sabin and Catmull-Clark surfaces that generalize non-uniform tensor product

B-spline surfaces to arbitrary topologies were introduced by Sederberg et al. [88]. All

the schemes mentioned above generalize recursive subdivision schemes for generating

limit surfaces with a known parameterization. Various issues involved with character

animation using these approximating subdivision schemes were discussed at length

by DeRose et al. [25].

The most well-known interpolation-based subdivision scheme is the "butterfly"

algorithm proposed by Dyn et al. [30]. Butterfly subdivision method, like other sub-

division schemes, makes use of a small number of neighboring vertices for subdivision.

It requires simple data structures and is extremely easy to implement. However, it

needs a topologically regular setting of the initial (control) mesh in order to obtain









a smooth C' limit surface. A variant of this scheme with better smoothness prop-

erties can be found in Dyn et al. [29]. Zorin et al. [109] developed an improved

interpolatory subdivision scheme (which we call the modified butterfly scheme) that

retains the simplicity of the butterfly scheme and results in much smoother surfaces

even from irregular initial meshes. These interpolatory subdivision schemes have

extensive applications in wavelets on manifolds, multiresolution decomposition of

polyhedral surfaces and multiresolution editing.

A variational approach for interpolatory refinement has been proposed by Kobbelt

[43, 44] and by Kobbelt and Schrider [46]. In this approach, the vertex positions in

the refined mesh at each subdivision step are obtained by solving an optimization

problem. Therefore, these schemes are global, i.e., every new vertex position depends

on all the vertex positions of the coarser level mesh. The local refinement property

which makes the subdivision schemes attractive for implementation in the computer

graphics applications is not retained in the variational approach.

The derivation of various mathematical properties of the smooth limit surface

generated by the subdivision algorithms is rather complex. Doo and Sabin [28] first

analyzed the smoothness behavior of the limit surface using the Fourier transform and

an eigen-analysis of the subdivision matrix. Ball and Storry [3, 4] and Reif [80] further

extended Doo and Sabin's prior work on continuity properties of subdivision surfaces

by deriving various necessary and sufficient conditions on smoothness for different

subdivision schemes. Specific subdivision schemes were analyzed by Schweitzer [87],

Habib and Warren [37], Peters and Reif [74] and Zorin [108]. Most recently, Stam









[90] presented a method for exact evaluation of Catmull-Clark subdivision surfaces

at arbitrary parameter values. It may be noted that the focus of this dissertation

is not on deriving a new subdivision technique or analyzing a subdivision technique,

but on deriving a deformable surface model using these subdivision schemes.

2.2 Physics-based Deformable Surface Models

Shape models can be broadly categorized into two types lumped parameter

models and distributed parameter models. A lumped parameter model uses small

number of parameters to represent model geometry, whereas a distributed parameter

model uses large number of degrees of freedom for model representation. An example

of a lumped parameter surface model is a superquadric. Spline surfaces serve as a

good example for distributed parameter model.

Deformable surfaces are distributed parameter models which use large number

of degrees of freedom for representing the model geometry. A large variety of shapes

can be modeled using this type of model, but handling a large number of degrees of

freedom can be cumbersome. However, the large degrees of freedom of a deformable

model are embedded in a physics-based framework to allow only a "physically mean-

ingful" model behavior. Various types of energies are assigned to the model using the

degrees of freedom, and the model "deforms" to find an equilibrium state with min-

imal energy. The motion equation is derived using Lagrangian dynamics [36], where

various energies associated with model gives rise to internal and external forces. The

equilibrium state of the model is a model position where the internal deformation

force becomes equal to the externally applied force. These physics-based deformable








models are useful for modeling where the modeler can deform a surface by applying

synthesized forces, and for data fitting where external forces are synthesized from a

given data set such that the model approximates the given data at equilibrium.

A deformable surface model is typically parameterized over the domain [0, 1]2.

Let s(u, v, p) be a deformable surface model, where 0 < u, v < 1 and p is a collection

of the degrees of freedom p,, i = 1, 2,..., n associated with the model (assuming the

model has n degrees of freedom). The degrees of freedom pi(t) is a set of generalized

coordinates which are functions of time and are assembled into the vector p. Let

fi(t) be the generalized force represented by the vector fp and acting on p,. Let T be

the kinetic energy, F be the dissipation energy and U be the potential energy of the

deformable surface model. Then, the Lagrangian equation of motion for the model

can be expressed as
d T 8T OF 9U
T + + f (2.1)
dt( p, op, 8p, op,

Let p(u, v) be the mass density function of the surface. Then the kinetic energy

of the surface is



T = J s sdudv = ipMp, (2.2)



where M is called the mass matrix. Similarly, let y(u, v) be the damping density

function of the surface. Then the dissipation energy is



F = Tsi f dudv = TDp, (2.3)
2I 2









where D is called the damping matrix. The potential energy of the model can be

defined using the thin-plate-under-tension energy model [96], and is given by


T f OT Ns s 9s'T s 82sT 2's
U = ( + a22 + i1u1
82sT 02s _2sT 2s 1
+12 a + 22 )dudv -pTKp, (2.4)
u9uv Yaby a2v 82v 2


where aii(u, v) and Aj(u, v) are elasticity functions which control tension and rigidity

respectively of the deformable surface model, and K is known as the stiffness matrix.

The discretized equation of motion can be derived using the expressions of the kinetic,

dissipation and potential energy in Eqn.2.1, which is given by



Mp + Dp + Kp = fp, (2.5)



where fp is the generalized force vector.

The free-form deformable models discussed above were first introduced to com-

puter graphics and visualization in Terzopoulos et al. [99] and further developed by

Terzopoulos and Fleischer [97], Pentland and Williams [72], Metaxas and Terzopou-

los [65] and Vemuri and Radisavljevic [104]. Celniker and Gossard [11] developed a

system for interactive free-form design based on the finite element optimization of en-

ergy functionals proposed in Terzopoulos and Fleischer [97]. Bloor and Wilson [8, 9],

Celniker and Welch [12] and Welch and Witkin [105] proposed deformable B-spline

curves and surfaces which can be designed by imposing the shape criteria via the

minimization of the energy functionals subject to hard or soft geometric constraints









through Lagrange multipliers or penalty methods. Qin and Terzopoulos [77, 78, 100]

developed dynamic NURBS (D-NURBS) which are very sophisticated models suit-

able for representing a wide variety of free-form as well as standard analytic shapes.

The D-NURBS have the advantage of interactive and direct manipulation of NURBS

curves and surfaces, resulting in physically meaningful hence intuitively predictable

motion and shape variation.

Deformable models are also widely used for shape recovery, segmentation, mo-

tion tracking and other computer vision and medical imaging applications. A detailed

survey of deformable models used in these techniques can be found in McInerney and

Terzopoulos [64] and the references therein. Some specific existing deformable models

used for shape recovery and non-rigid motion estimation will be reviewed in Section

2.4.

2.3 Shape Modeling Using Physics-based Subdivision-surface Model

Recursive subdivision surfaces are powerful for representing smooth geometric

shapes of arbitrary topology. However, they constitute a purely geometric represen-

tation, and furthermore, conventional geometric modeling with subdivision surfaces

may be difficult for representing highly complicated objects. For example, modelers

are faced with the tedium of indirect shape modification and refinement through time-

consuming operations on a large number of (most often irregular) control vertices

when using typical subdivision surface-based modeling schemes. Despite the advent

of advanced 3D graphics interaction tools, these indirect geometric operations remain

non-intuitive and laborious in general. In addition, it may not be enough to obtain









the most "fair" surface that interpolates a set of (ordered or unorganized) data points.

A certain number of local features such as bulges or inflections may be strongly de-

sired while requiring geometric objects to satisfy global smoothness constraints in

geometric modeling and computer graphics applications. In contrast, physics-based

modeling provides a superior approach to shape modeling that can overcome most of

the limitations associated with traditional geometric modeling approaches. Free-form

deformable models governed by the laws of continuum mechanics are of particular

interest in this context. These dynamic models respond to externally applied forces

in a very intuitive manner. The dynamic formulation marries the model geometry

with time, mass, damping and constraints via a force balance equation. Dynamic

models produce smooth, natural motions which are easy to control. In addition, they

facilitate interaction especially direct manipulation of complex geometries. Fur-

thermore, the equilibrium state of the model is characterized by a minimum of the

deformation energy of the model subject to the imposed constraints. The deformation

energy functionals can be formulated to satisfy local and global modeling criteria, and

geometric constraints relevant to shape design can also be imposed. The dynamic

approach subsumes all of the aforementioned modeling capabilities in a formulation

which grounds everything in real-world physical behavior.

A severe limitation of the existing deformable models, including D-NURBS, is

that they are defined on a rectangular parametric domain. Hence, it can be very

difficult to model surfaces of arbitrary genus using these models. In a recent work,









DeRose et al. [25] assigned material properties to control meshes at various subdivi-

sion levels in order to simulate cloth dynamics using subdivision surfaces. Note that,

they assign physical properties on the control meshes at various levels of subdivision

and not on the limit surface itself, and hence can not solve the modeling goal we

are trying to achieve. In this dissertation, a dynamic framework is presented for

subdivision surfaces which combines the benefits of subdivision surfaces for modeling

arbitrary topology as well as that of dynamic splines for interactive shape manip-

ulation by applying synthesized forces. The proposed dynamic framework presents

the modeler a formal mechanism of direct and intuitive manipulation of the smooth

limit surface, as if they were seamlessly sculpting a piece of real-world "clay." A

dynamic framework for subdivision surface-based wavelets adds flexibility in the pro-

posed modeling paradigm. It enables a multiresolution representation of the evolving

control mesh defining the smooth limit surface, and the modeler can control the

extent of manipulation effect by choosing a proper editing level. The formulation

and implementation details of these dynamic frameworks are discussed in subsequent

chapters.

2.4 Shape and Motion Estimation Using Physics-based Subdivision-surface Model

The dynamic subdivision surface model has been developed primarily for mod-

eling arbitrary (known) topology where modelers can directly manipulate the limit

surface by applying synthesized forces in an intuitive fashion. However, another im-

portant application of the dynamic subdivision surfaces is in non-rigid shape and









motion reconstruction/recovery. Accurate shape recovery requires distributed pa-

rameter models which typically possess a large number of degrees of freedom. On

the other hand, efficient shape representation imposes the requirement of geometry

compression, i.e., models with fewer degrees of freedom. These requirements are

conflicting and numerous researchers have been seeking to strike a balance between

these contradicting requirements [5, 21, 47, 49, 62, 89, 100, 104]. Another important

criterion in the design of shape models is that the initialization of the model during

the shape recovery process should not be restricted to globally parameterized input

meshes since it may be infeasible to globally parameterize shapes of arbitrary topol-

ogy. A physics-based model best satisfying the above mentioned criteria is an ideal

candidate for a solution to the shape recovery problem for obvious reasons.

Deformable models which come in many varieties, have been used to solve this

problem in the physics-based modeling paradigm. These models involve the use of

either fixed size [21, 66, 71, 98, 104] or adaptive size [15, 41, 47, 62, 95, 101] grids.

The models with fixed grid size generally use less number of degrees of freedom for

representation at the cost of accuracy of the recovered shape. On the other hand,

models using adaptive grids generally need large number of degrees of freedom to

recover the shapes accurately. Note that the shapes being recovered from the image

data are smooth in most of the medical applications, i.e. the anatomical structures

even with considerable amount of details have more or less a C1 surface. Under

these circumstances, the finite element approaches as in Cohen and Cohen [21] and

McInerney and Terzopoulos [62] need a large number of degrees of freedom for deriving









a smooth and accurate representation. In addition, they can not represent shapes

with known arbitrary topology. Moreover, almost all of these schemes require a

globally parameterized mesh as their input which may be infeasible at times.

The proposed model solves the shape recovery problem very effectively as it

can recover shapes from large range and volume data sets using very few degrees

of freedom (control vertices) for its representation and can cope with any arbitrary

input mesh, not necessarily parameterized. The initialized model deforms under the

influence of synthesized forces to fit the data set by minimizing its energy. Once the

approximate shape is recovered, the model is further subdivided automatically and a

better approximation to the input data set is achieved using more degrees of freedom.

The process of subdivision after achieving an approximate fit is continued till a pre-

scribed error criteria for fitting the data points is achieved. The proposed dynamic

subdivision surface models have also been successfully used in motion tracking and

visualization of a dynamically deforming shape from a time sequence of volumetric

data sets.

2.5 Unified Framework for Shape Recovery and Shape Modeling

Currently, shape recovery and shape modeling are viewed as two distinct areas in

computer graphics and geometric modeling literature. However, there are potential

benefits if these two can be combined in an unified framework. For example, the

modeler starts from scratch to build a specific model in a typical geometric modeling

scenario. First, a rough shape is modeled and then it is fine-tuned by manipulating

control vertex positions to obtain the desired effects. This turns out a cumbersome









process in general. On the other hand, shape recovery using state-of-the-art methods

yield large polygonal meshes which are very difficult to manipulate, especially for

global changes in shape. Most often the resulting meshes from a shape recovery

application are not directly amenable to multiresolution analysis. Computationally

expensive re-meshing techniques are needed to convert these meshes into a specific

type of meshes on which multiresolution analysis can be performed. This specific type

of mesh is known as mesh with "subdivision-connectivity", implying a topologically

equivalent mesh with the same connectivity as of the given mesh can be obtained by

recursive subdivision of a very simple known initial mesh.

The proposed dynamic framework combines shape recovery and shape modeling

in an unified framework where the modeler can scan in 3D points of a prototype

model, recover the shape using the proposed dynamic subdivision surface model, and

edit at any desired resolution using physics-based force tools. Thus, the modeler is re-

lieved of the burden of both building an initial model and editing a cumbersome huge

polygonal mesh. The shape recovery process starts with a smooth subdivision surface

model which has a simple initial mesh. This model is deformed using synthesized

forces from the given data points and is automatically refined using some pre-defined

error in fit criteria. The initial mesh of the final recovered smooth shape has subdivi-

sion connectivity in-built, as it is obtained by subdivision refinement, and therefore

no re-meshing is needed for multiresolution analysis. The dynamic framework for

subdivision surface-based wavelets makes multiresolution editing using physics-based





26


force tools easy to perform. These advantages of shape modeling and shape recovery

in a unified framework will be further illustrated in later chapters.

















CHAPTER 3
DYNAMIC CATMULL-CLARK SUBDIVISION SURFACES



Subdivision surfaces, as mentioned earlier, can be broadly classified into two

categories approximating and interpolatory subdivision schemes. The approxi-

mating subdivision schemes typically generalize recursive subdivision techniques for

generating limit surfaces with known parameterizations. In this chapter, a dynamic

framework will be presented for Catmull-Clark subdivision scheme, a popular approx-

imating subdivision technique. The Catmull-Clark subdivision scheme generalizes the

idea of obtaining a bicubic B-spline surface by recursive subdivision of a rectangu-

lar initial mesh. Before discussing the dynamic framework, first the Catmull-Clark

subdivision scheme is briefly reviewed.

3.1 Overview of the Catmull-Clark Subdivision Scheme

Catmull-Clark subdivision scheme, like any other subdivision scheme, starts

with a user-defined mesh of arbitrary topology. It refines the initial mesh by adding

new vertices, edges and faces with each step of subdivision following a fixed set of

subdivision rules. In the limit, a recursively refined polyhedral mesh will converge to

a smooth surface. The subdivision rules are as follows:





















(b) (c)


Figure 3.1. Catmull-Clark subdivision : (a) initial mesh, (b) mesh obtained after one
step of Catmull-Clark subdivision, and (c) mesh obtained after another subdivision
step.


For each face, a new face vertex is introduced which is the average of all the

old vertices defining the face.


For each (non-boundary) edge, a new edge vertex is introduced which is the

average of the following four points: two old vertices defining the edge and two

new face vertices of the faces adjacent to the edge.


For each (non-boundary) vertex V, a new vertex is introduced whose position

is E + 2E + n-3, where F is the average of the new face vertices of all faces

adjacent to the old vertex V, E is the average of the midpoints of all the edges

incident on the old vertex V and n is the number of the edges incident on the

vertex.


New edges are formed by connecting each new face vertex to the new edge

vertices of the edges defining the old face and by connecting each new vertex









corresponding to an old vertex to the new edge vertices of all old edges incident

on the old vertex.


New faces are defined as faces enclosed by new edges.


An example of Catmull-Clark subdivision on an initial mesh is shown in Fig.3.1.

The most important property of the Catmull-Clark subdivision surfaces is that a

smooth surface can be generated from any control mesh of arbitrary topology. There-

fore, this subdivision scheme is extremely valuable for modeling various complicated

geometric objects of arbitrary topology. Catmull-Clark subdivision surfaces include

standard bicubic B-spline surfaces as their special case (i.e., the limit surface is a

bicubic B-spline surface for a rectangular mesh with all non-boundary vertices of

degree 4). In addition, the aforementioned subdivision rules generalize the recur-

sive bicubic B-spline patch subdivision algorithm. For non-rectangular meshes, the

limit surface converges to a bicubic B-spline surface except at a finite number of ex-

traordinary points. These extraordinary points correspond to extraordinary vertices

(vertices whose degree is not equal to 4) in the mesh. Note that, after the first sub-

division, all faces are quadrilaterals, hence all the new vertices created subsequently

will have four incident edges. The number of extraordinary points on the limit sur-

face is a constant, and is equal to the number of extraordinary vertices in the refined

mesh obtained after applying one step of the Catmull-Clark subdivision on the initial

mesh. The limit surface is curvature-continuous everywhere except at extraordinary

vertices, where only tangent plane continuity is achieved. In spite of the popularity

of Catmull-Clark subdivision surfaces for representing complex geometric shapes of









arbitrary topology, these subdivision surfaces may not be easily parameterizable and

deriving a closed-form analytic formulation of the limit surface can be very difficult.

These deficiencies preclude their immediate pointwise manipulation and hence may

restrain the applicability of these schemes. In the rest of the chapter, a dynamic

framework is developed for the Catmull-Clark subdivision surfaces which offers a

closed-form analytic formulation and allows users to manipulate the model directly

and intuitively. It may be noted that recently a scheme was proposed by Stam [90]

for directly evaluating the Catmull-Clark subdivision surfaces at arbitrary parameter

values, but the work presented here on dynamic Catmull-Clark subdivision surfaces

[57, 60, 76] was published prior to the publication of the above-mentioned work.

3.2 Formulation

In this section, a systematic formulation of the new dynamic model based on

Catmull-Clark subdivision surfaces is presented. This dynamic model treats the

smooth limit surface as a function of its initial mesh. However, the control vertex po-

sitions need to be updated continually at every level of subdivision in order to develop

the dynamic framework. Note that, all the vertices introduced through subdivision

are obtained as an affine combination of control vertex positions in the initial mesh.

Therefore, the dynamic behavior of the limit surface can be controlled by formulat-

ing the dynamic model on the initial mesh itself, the only exception being the case

when the initial mesh has non-rectangular faces. This problem can be circumvented

by taking the mesh obtained through one step of subdivision as the initial mesh.

It may be noted that an additional subdivision step may be needed in some cases









to isolate the extraordinary points, and the resulting mesh is treated as the initial

mesh. A typical example of the above mentioned scenario is when the initial mesh is

a tetrahedron where two subdivision steps are needed, and the dynamic framework

can be developed by treating the mesh obtained after two subdivision steps as the

initial mesh. To define the limit surface using the vertices of the initial mesh, the

enumeration of the bicubic patches in the limit surface is necessary. In the next two

sections, various schemes are presented for assigning the bicubic patches to various

faces of the initial mesh.

3.2.1 Assigning Patches to Regular Faces

In Fig.3.2, a rectangular control mesh is shown along with the bicubic B-spline

surface (4 patches) in the limit after an infinite number of subdivision steps. Note

that, each of the bicubic patches in the limit surface is defined by a rectangular face

with each vertex of degree four, thereby accounting for 16 control points (from its

8 connected neighborhood) needed to define a bicubic surface patch in the limit.

Therefore, for each rectangular face in the initial mesh with a degree of 4 at each

vertex, the corresponding bicubic surface patch can be assigned to it in a straight

forward way. In Fig.3.2, the surface patches S1,S2, S3 and S4 are assigned to face

F1, F2, F3 and F4 respectively. The 16 control points for the patch S1, corresponding

to face F1, are highlighted in Fig.3.2. Note that, the initial control mesh can be viewed

as the parametric domain of the limit surface. Therefore, face Fi can be thought of

as the portion of the parametric domain over which the patch Si is defined, i.e., has

non-zero values. Nevertheless, each rectangular face (e.g. Fi) can be parametrically





32

















U-J J - -













[-- -- "-- ~- _




Figure 3.2. A rectangular mesh and its limit surface consisting of 4 bicubic surface
patches.









defined over [0, 1]2, and hence, all bicubic B-spline patches defined by 16 control

points are locally parameterized over [0, 1]2.

3.2.2 Assigning Patches to Irregular Faces

In Fig.3.3, a mesh containing an extraordinary point of degree 3 and its limit sur-

face are shown. The faces Fo, FI, ..., F are assigned to bicubic patches So, Si,..., Ss

respectively (as they all have vertices of degree 4) following the aforementioned

scheme. However, the central smooth surface enclosed by the patches So, S1,..., Ss

consists of infinite number of bicubic patches converging to a point in the limit. A re-

cursive way of enumerating these bicubic patches and assigning them to various faces

at different levels need to be developed in order to formulate the dynamic framework

for Catmull-Clark subdivision surface model.

The idea of enumerating the bicubic patches corresponding to faces having an

extraordinary vertex is shown in Fig.3.4 where a local subdivision of the mesh con-

sisting of faces Fo, F,..., F. PF P P2 (and not the other boundary faces) of Fig.3.3

is carried out. Topologically, the resulting local subdivision mesh (shown as dotted

mesh) is exactly the same as the mesh in Fig.3.3 and hence exactly the same number

of bicubic patches can be assigned to its faces with vertices of degree 4 as is evident

from Fig.3.4 (the new faces and the corresponding patches are marked by "p" and

"n" respectively). This process of local subdivision and assignment of bicubic patches

around an extraordinary point can be carried out recursively and in the limit, the

enclosed patch corresponding to faces sharing the extraordinary point will converge

to a point. However, there is no need to carry out an infinite number of subdivision





34

































S 4
So S S S
-----t ----- --- p-i




F 3n






Figure 3.3. A mesh with an extraordinary point of degree 3 and its limit surface.














selected mesh for
S local subdivision
- -


\7


Figure 3.4. Local subdivision around the extraordinary point and the limit surface.









steps. This description is for formulation purposes only and the exact implementation

will be detailed in a later section.

3.2.3 Kinematics

In this section, the mathematics for the kinematics of the limit surface is devel-

oped via illustrative examples and then the generalized formulas are presented. It

may be noted that a single bicubic B-spline patch is obtained as the limiting process

of the Catmull-Clark subdivision algorithm applied to an initial 4 by 4 rectangular

control mesh. Let sp(u, v), where (u, v) E [0, 1]2, denote this bicubic B-spline patch

which can be expressed analytically as


3 3
sp(u, v) = (x(u, v), y(u, v), z(u, v)) = E E di,j B,4(u)Bj,4(v), (3.1)
i=0 j=0


where dj represents a 3-dimensional position vector at the (i,j)th control point

location and Bi,4(u) and B,4 (v) are the cubic B-spline basis functions. The subscript

p on s denotes the patch under consideration. Expressing Eqn.3.1 in a generalized

coordinate system we have

sp = Jq,, (3.2)


where Jp is a transformation matrix storing the basis functions of a bicubic B-spline

patch, and is of size (3, 48). Vector q, is the concatenation of all control points

defining a B-spline patch in 3D. Note that in the concatenation of the control points,

each control point has an (x, y, z) component. For example, the (x, y, z) components

of the control point (i,j) correspond to positions 3k,3k + 1,3k + 2 where, k =









4i + j respectively in the vector qp. We can express the entries of Jp explicitly

in the following way: Jp(0, k) = Jp(1, k + 1) = Jp(2, k + 2) = Bi,4(u)Bj,4(v) and

J,(0, k + 1) = J,(0, k + 2) = Jp(1, k) = Jp(1, k + 2) = Jp(2, k) = Jp(2, k + 1) = 0.

Limit surface with many bicubic patches from a rectangular initial mesh

Now a limit surface consisting of many bicubic surface patches obtained after

applying an infinite number of subdivision steps to a rectangular initial mesh is

considered. For example, let the limit surface of Fig.3.2 be sm which can be written

as


1 1 1 1
sm(u, v) = sm,(2u, 2v)+s,.(2(u- ),2v)+s, (2(u--),2(v-2))+sa(2u,2(v- )),

(3.3)

where s, (2u, 2v) sm(u, v) for 0 < u, v < and 0 otherwise. Similarly, s2,,Sm3

and s., are also equal to sm(u, v) for an appropriate range of values of u, v and 0

outside. It may be noted that SMn,, ,2,, sm,, sn correspond to patches S1, S, S3, S4

respectively in Fig.3.2. Eqn.3.3 in generalized coordinates becomes


4
Sm = Jlq + J2q2 + J3q3 + J44 = Z Jjqi, (3.4)
i=1


where Jis are the transformation matrices of size (3, 48) and qis are the (x,y,z) com-

ponent concatenation of a subset of the control points of sm defining sm,, i = 1, 2,3









and 4. A more general expression for Sm is


4
sm = JAq + J2A2qm + JaA3qm + J4A4qm = JiAfqm = Jmqm. (3.5)
i=l


Where, q, is the 75-component vector of 3D positions of the 25 vertex control mesh

defining the limit surface si. Matrices Ai, 1 < i < 4, are of size (48, 75) each

row consisting of a single nonzero entry (= 1) and the (3,75)-sized matrix Jm =

E=1 JA,.

Limit surface with many bicubic patches from an arbitrary initial mesh

The stage is now set to define the limit surface s using the vertices of initial

mesh M for any arbitrary topology, assuming all faces are quadrilateral and no face

contains more than one extraordinary point as its vertex (i.e., extraordinary points

are isolated). As mentioned earlier, if these assumptions are not satisfied, one or two

steps of global subdivision may be required and the resulting mesh can be treated

as the initial mesh. Let the number of vertices in the initial mesh M be a, and let

I of these be the extraordinary vertices. The number of faces in the initial mesh are

assumed to be b, and k of these faces have vertices with degree 4 (henceforth termed a

"normal face") and each of the remaining (b- k) faces have one of the I extraordinary

vertices (henceforth termed a "special face"). Let p be the 3a = N dimensional vector

containing the control vertex positions in 3D. Using the formulations in Section 3.2.1









and Section 3.2.2, the smooth limit surface can be expressed as


k
s = En,+ s, (3.6)
i=1 j=1


where n, is a single bicubic patch assigned to each of the normal faces and sj is a

collection of infinite number of bicubic patches corresponding to each of the extraor-

dinary points. As s is a surface defined over the faces of the initial control mesh, each

n, can be viewed as a bicubic patch defined over the corresponding regular face in the

initial control mesh. Similarly, each sj can be defined over the corresponding irregu-

lar faces in the initial mesh (refer to Fig.3.2 Fig.3.4 for the detailed description on

parametric domains of subdivision surfaces). For the simplicity of the following math-

ematical derivations on the dynamic model, from now on the parametric domain will

not be explicitly provided in the formulation. Employing the same approach taken

before to derive Eqn.3.5, it can be shown that


k k k
Sni = E(n')(P" ) = ( n J.)("A,))p = ("J)p, (3.7)
I=1 i=1 i=l


where "Ji,n pi and "nA are the equivalent of Ji, p, in Eqn.3.4 and A, in Eqn.3.5 respec-

tively. The pre-superscript n is used to indicate that these mathematical quantities

describe bicubic patch in the limit surface corresponding to normal faces.

The following notational convention will be used for describing various math-

ematical quantities used in the derivation of the expression for a collection of infinite

number of bicubic patches around an extraordinary vertex. The pre-superscript s









is used to represent a collection of bicubic patches around an extraordinary vertex,

the subscript j is used to indicate the j-th extraordinary point, the post-superscript

represents the exponent of a mathematical quantity and the level indicator (to rep-

resent various levels of subdivision in the local control mesh around an extraordinary

vertex) is depicted via a subscript on the curly braces.

The expression for sj is derived using the recursive nature of local subdivision

around an extraordinary vertex as shown in Section 3.2.2. First, sj can be expressed

as

si = {'Jj}l{Spi + {s,}A, (3.8)


where the first term of Eqn.3.8 is the generalized coordinate representation of the

bicubic B-spline patches corresponding to the normal faces of the new local subdivi-

sion mesh obtained after one subdivision step on the local control mesh (similar to

those patches marked n in Fig.3.4). {sj}, represents the rest of the infinite bicubic

B-spline patches surrounding the extraordinary point (similar to the central patch

enclosed by patches marked n in Fig.3.4). The vertices in the newly obtained lo-

cal subdivision mesh {'pjl} can be expressed as a linear combination of a subset

of the vertices of the initial mesh M (which will contribute to the local subdivi-

sion) following the subdivision rules. We can name this subset of initial control

vertices {'pj}o. Furthermore, there exists a matrix {'Bi,} of size (3c,3d), such that

({Bj}l('pj}o = {'Pj}i where {'pj}1 and {'pj}o are vectors of dimension 3c and 3d

respectively. Applying the idea of recursive local subdivision again on {sj}l, sj can








be further expanded as



s = {*Jj3}1{'B3j},'pj} + {SJj03 'Bj}2{%jp} + {si}1. (3.9)


In the above derivation, {'fj}1 is a vector of dimension 3d, comprising of a

subset of the vertices defining the 3c dimensional vector {'pj}i. Note that, {'fij}

has the same structure as {%pj}o, therefore, there exists a (3d,3d) matrix {*Cj}1

such that {'Cj),{'pj}, = {ij},. Each subdivision of a local mesh with d vertices

creates a new local mesh with c vertices which contributes a fixed number of bicubic

B-spline patches. Proceeding one step further, it can be shown that



sj = {*Jl}1{'Bj,}'p}o+ {'Jj,}2{Bj}2{"'C}1{'pj}o

+{'Jj}3{Bj}3{'Cj}1{'pj}o + {sj}3. (3.10)


Because of the intrinsic property of the local recursive subdivision around the ex-

traordinary point, we have {'Jj} { SJj}2 = ... = {S'Jj} = ... = {'Jj}o In

addition, the subdivision rules remain the same throughout the refinement process,

we also have {1Bj}1 = {'Bj}2 = ... = {'Bj} = ... = {'Bj}m. So, the above

equations can be further simplified leading to


s = {SJj}1{Bj}{'pj+}o+ {'Jj}{'SBj}] {C,}{'pj}o

+{"'Jj}1{'Bj}( {'Cj}){'pj}o +.

= {*JJ}{iB,} ,( {SC)})P"{p "o, (3.11)
-=0









Vector sj can be written as



sj = ('J, i'p,1, (3.12)



where 'Jj = {'Jj,}l{'Bj}l(= 0o {'Cj,}) and 'pj = {'pj}0. The ideaoflocal recursive

subdivision around an extraordinary point is illustrated in Fig.3.5. Note that, each

vertex position in the subdivided mesh is obtained by an affine combination of some

vertices in the previous level and hence any row of ({Cj}1 sums to 1. The largest

eigenvalue of such a matrix is 1 and it can be shown that the corresponding infinite

series is convergent following a similar approach as in Halstead et al. [38]. The rest

of the derivation leading to an expression for s is relatively straight forward. Using

the same approach used to derive the Eqn.3.7, it can be shown that



sj = E(,'j)(p)) = (-('Jj)('Aj))p = ('J)p. (3.13)
j=1 j=1 j=1


From Eqn.3.6, Eqn.3.7 and Eqn.3.13,



s = ("J)p + ('J)p. (3.14)



Let J = ("J) + ('J), and hence


s = Jp.


(3.15)




















{SB}









i J1 0



Figure 3.5. Local subdivision around the extraordinary point and the corresponding
patches in the limit surface from different levels of subdivision.
f<'>;?B^S>^ (Sj.

WU ^^^ ^l
^ps''^^ *







patches in the limit surface from different levels of subdivision.








3.2.4 Dynamics

In the previous section, an expression is derived for the smooth limit surface

obtained via infinite subdivision steps. This expression treats the smooth limit surface

as a function of the control vertex positions in the initial mesh. Now the vertex

positions in the initial mesh defining the smooth limit surface s are treated as a

function of time in order to embed the Catmull-Clark subdivision model in a dynamic

framework. The velocity of this surface model can be expressed as



s(u,v,p)= Jp, (3.16)



where an overstruck dot denotes a time derivative. The physics of the dynamic

subdivision surface model is based on the work-energy version of Lagrangian dynamics

[36] and is formulated in an analogous way to the dynamic framework presented in

Section 2.2.

Let ju(u, v) be the mass density function of the surface. Then the kinetic energy

of the surface is



T = Jf sdudv = p Mp, (3.17)



where (using Eqn.3.16)


M= f pJJdudv


(3.18)








is the mass matrix of size (N, N). Similarly, let y(u, v) be the damping density

function of the surface. Then the dissipation energy is


F = J fy sdudv = j1 Dp, (3.19)


where

D= f / 7JTjdudv (3.20)


is the damping matrix of size (N, N).

The potential energy of the smooth limit surface can be expressed as


U = pTKp, (3.21)


where K is the stiffness matrix of size (N, N). A thin-plate-under-tension energy

model [96] is used to compute the elastic potential energy of the dynamic subdivision

surface. The corresponding expression for the stiffness matrix K is


K = J(anJ J + a22J J +011 .TJ. + J12J J. + 223J J,)dudv. (3.22)


where the subscripts on J denote the parametric partial derivatives. The aii(u, v)

and i., ,,. v)s are elasticity functions controlling local tension and rigidity in the two

parametric coordinate directions. Note that the thin-plate energy expression might

diverge at extraordinary points on the limit surface for Catmull-Clark subdivision

scheme as shown in Halstead et al. [38]. Several methods have been suggested in









Halstead et al. [38] to overcome the problem of the divergent series. However, in this

dissertation, the problem is circumvented by setting the corresponding rigidity coef-

ficients to be zero at these points. Therefore, the thin-plate energy at extraordinary

points is zero. The effect is negligible to the overall thin-plate energy as very few

extraordinary points are present in the smooth limit surface.

Using the expressions for the kinetic, dissipation and potential energy in Eqn.2.1,

the same motion equation as in Eqn.2.5 can be obtained. The generalized force vector

fp, which can be obtained through the principle of virtual work [36], is expressed as



f, = f Jf(u, v, t)dudv (3.23)



Different types of forces can be applied on the smooth limit surface, and the limit

surface would evolve over time according to Eqn.2.5 to obtain an equilibrium position

characterized by a minimum of the total model energy.

3.2.5 Multilevel Dynamics

At this point, a dynamic framework is developed for the Catmull-Clark subdi-

vision scheme where the smooth limit surface evolves over time in response to the

applied forces. The entire process can be described as follows. Given an initial mesh,

a smooth surface is obtained in the limit. Users can directly apply synthesized forces

to this smooth limit surface to enforce various functional and aesthetic constraints.

This direct manipulation is then transferred back as virtual forces acting on the initial

mesh through a transformation matrix (Eqn.3.23), and the initial mesh (as well as









the underlying smooth limit surface) deforms continuously over time until an equilib-

rium position is obtained. The deformation of the surface in response to the applied

forces is governed by the motion equation (Eqn.2.5). Within the physics-based mod-

eling framework, the limit surface evolves as a consequence of the evolution of the

initial mesh. One can apply various types of forces on the limit surface to obtain

a desired effect, but the possible level of details appearing in a shape that can be

obtained through evolution is constrained by the number of control vertices in the

initial mesh. It might be necessary to increase the number of control vertices in the

initial mesh in order to obtain a detailed shape through this evolution process.

The number of control vertices defining the same smooth limit surface can be

increased by simply replacing the initial mesh with a mesh obtained after one subdi-

vision step applied to the initial mesh. This new mesh has more number of vertices

but defines the same limit surface. Therefore, the dynamic Catmull-Clark surface

model can be subdivided globally to increase the number of vertices (control points)

of the model. For example, after one step of global subdivision, the initial degrees of

freedom p (refer to Eqn.3.15 and Eqn.3.16) in the dynamic system will be replaced

by a larger number of degrees of freedom q, where q = Ap. A is a global subdivi-

sion matrix of size (M, N) whose entries are uniquely determined by Catmull-Clark

subdivision rules (see Section 3.1 for the details about the rules). Thus, p, expressed

as a function of q, can be written as


p = (ATA)-IATq = Atq,


(3.24)









where At = (ATA)-'AT. Therefore, Eqn.3.15 and Eqn.3.16 can be rewritten as



s = (JAt)q, (3.25)



and

s(u, v, q) = (JAt)q, (3.26)


respectively. Now the equation of motion needs to be derived for this new subdivided

model involving a larger number of control vertices namely q. The mass, damping

and stiffness matrices need to be recomputed for this "finer" level. The structure of

the motion equation as given by Eqn.2.5 remains unchanged, but the dimensionality

and the entries of M, D, K, p and fp change correspondingly in this newly obtained

subdivided level. In particular the motion equation, explicitly expressed as a function

of q, can be written as

Mq + Dq + Kq = fq, (3.27)


where Mq f f p(At)TJTJAtdudv and the derivation of D,, K, and f, follow suit.

It may be noted that further subdivision, if necessary, can be carried out in a

similar fashion. Therefore, multilevel dynamics is achieved through recursive subdi-

vision on the initial set of control vertices. Users can interactively choose the level of

detail representation of the dynamic model as appropriate for their modeling and de-

sign requirements. Alternatively, the system can automatically determine the level of

subdivision most suitable for an application depending on some application-specific

criteria.









3.3 Finite Element Implementation

The evolution of the generalized coordinates for our new dynamic surface model

can be determined by the second-order differential equation as given by Eqn.2.5. An

analytical solution of the governing differential equation can not be derived in general.

However, an efficient numerical implementation can be obtained using the finite ele-

ment method [42]. The limit surface of the dynamic Catmull-Clark subdivision model

is a collection of bicubic patch elements. We use two types of finite elements for this

purpose normal elements (bicubic patches assigned to the normal faces of the initial

mesh) and special elements (collection of infinite number of bicubic patches assigned

to each extraordinary vertex of the initial mesh). The shape functions for the nor-

mal element are the basis functions of the bicubic surface patch, whereas the shape

functions for the special element are the collection of basis functions corresponding

to the bicubic patches in the special element. In the current implementation, the

M, D and K matrices for each individual normal and special elements are calculated

and they can be assembled into the global M, D and K matrices that appear in the

corresponding discrete equation of motion. In practice, we never assemble the global

matrices explicitly in the interest of time performance. The detailed implementation

is explained in the following sections.

3.3.1 Data Structures

The limit surface of the dynamic Catmull-Clark subdivision model is a collection

of bicubic patches, and hence requires us to keep track of the control polygons defining

such patches. We can get the control polygons for the normal elements in the limit










Slist of faces
list of edges
Subdiv object
S e list of vertices
list of normal elements


S subdivision aro




list of faces list of facdges
lubd of d ---- ofbdi ob
Jim ofnorml elemCls i st of ononml element




Figure 3.6. The data structure used for dynamic subdivision surface implementation.



surface from the initial control mesh itself. However, we need to locally subdivide the

initial control mesh around the extraordinary vertices to obtain the control polygons

of the bicubic patches in the special elements on the limit surface.

A subdivision surface defined by a control mesh at any level is designed as a

class which has a pointer to its parent mesh, a set of pointers to its offspring meshes

(arising out of local subdivision around the extraordinary vertices at that level), a list

of faces, edges, vertices and normal elements Face, edge, vertex and normal elements

are, in turn, classes which store all the connectivity and other information needed to

either enumerate all the patches or locally subdivide around an extraordinary vertex

in that level. The implementation takes the initial mesh as the base subdivision

surface object (with its parent pointer set to NULL) and locally subdivides the initial

mesh upto a user-defined maximum level around each extraordinary vertex to create









offspring objects at different levels (Fig.3.6). At this point, let's take a closer look at

the normal and special element data structures and computation of the corresponding

local M, D and K matrices.

3.3.2 Normal Elements

Each normal element is a bicubic surface patch and hence is defined by 16

vertices (from the 8-connected neighborhood of the corresponding normal face) in the

initial control mesh. Each normal element keeps a set of pointers to those vertices

of the initial mesh which act as control points for the given element. For a normal

element, the mass, damping and stiffness matrices are of size (16,16) and can be

computed exactly by carrying out the necessary integration analytically. The matrix

J in Eqn.3.18, 3.20 and 3.22 needs to be replaced by J, (of Eqn.3.2) for computation

of the local M, D and K matrices respectively of the corresponding normal element.

3.3.3 Special Elements

Each special element consists of an infinite number of bicubic patches in the

limit. We have already described a recursive enumeration of the bicubic patches of

a special element in Section 3.2.2. Let us now consider an arbitrary bicubic patch of

the special element in some level j. The mass matrix M, of this patch can be written

as

M, = n' Mps (3.28)


where Mp is the normal element mass matrix (scaled by a factor of y to take into

account of the area shrinkage in bicubic patches at higher level of subdivision) and









f, is the transformation matrix of the control points of that arbitrary patch from the

corresponding control points in the initial mesh. The damping and stiffness matrices

for the given bicubic patch can be derived in a similar fashion. These mass, damping

and stiffness matrices from various levels of (local) subdivision can then be assembled

to form the mass, damping and stiffness matrices of the special element. It may be

noted that the stiffness energy due to plate terms diverges at the extraordinary points

on the limit surface. We solve the problem using a slightly different approach than

the one used in Halstead et al. [38]. When the area of the bicubic patch obtained via

local subdivision of the initial mesh around an extraordinary vertex becomes smaller

than the display resolution, the contribution from such a bicubic patch is ignored in

computing the physical matrices of the corresponding special element. The number

of extraordinary points in the limit surface is very few, and the above mentioned

approximation is found to work well in practice.

3.3.4 Force Application

The force f(u, v, t) in Eqn.3.23 represents the net effect of all externally applied

forces. The current implementation supports spring, inflation as well as image-based

forces. However, other types of forces like repulsion forces, gravitational forces etc.

can easily be implemented.

To apply spring forces, a spring of stiffness k can be connected from a point do

to a point xo on the limit surface (or, to the j-th level approximation mesh), the net








applied spring force being


f(x, t) = s k(d0 s(x, t))6(x xo)dx, (3.29)



where 6 is the unit impulse function implying f(xo, t) = kido s(xo, t) and vanishes

elsewhere on the surface. However, the 6 function can be replaced with a smooth

kernel to spread the force over a greater portion on the surface. The spring forces can

be applied interactively using the computer mouse or the points from which forces

need to be applied can be read in from a file.

To recover shapes from 3D image data, we synthesize image-based forces. A 3D

edge detection is performed on a volume data set using the 3D Monga-Deriche(MD)

operator [69] to produce a 3D potential field P(x, y, z), which we use as an external

potential for the model. The force distribution is then computed as



f(x,y,z)= A VP(x, z) (3.30)
1 VP(lx, y, z) II'


where A controls the strength of the force. The applied force on each vertex at the j-

th approximation level is computed by trilinear interpolation for evaluating Eqn.3.23

in Cartesian coordinates. More sophisticated image-based forces which incorporate

region-based information such as gradients of a thresholded fuzzy voxel classification

can also be used to yield better and more accurate shape recovery. It may be noted

that we can apply spring forces in addition with the image-based forces by placing









points on the boundary of the region of interest in each slices of the 3D volume (MR,

CT etc.) image data.

3.3.5 Discrete Dynamic Equation

The differential equation given by Eqn.2.5 is integrated through time by dis-

cretizing the time derivative of p over time steps At. The state of the dynamic

subdivision surface at time t + At is integrated using prior states at time t and

t At. An implicit time integration method is used in the current implementation

where discrete derivatives of p are calculated using


( + At) p(t + At) 2p(t) + p(t At)
p(t+ At) =At2 (3.31)



and
S p(t + At) p(t At)
p(t + At) = (3.32)
2At

The elemental mass, damping and stiffness matrices can be assembled to get the

global mass, damping and stiffness matrix for the smooth subdivision surface model.

However, we do not assemble these global sparse matrices explicitly for efficiency

reasons. For the time varying stiffness matrix, we recompute K at each time step.

Using Eqn.2.5, Eqn.3.31 and Eqn.3.32, the discrete equation of motion is obtained as



(2M + DAt + 2At2K)p(t + At) = 2At2(t + At) + (DAt -2M)p(t At) + 4Mp(t).

(3.33)









A
evoslton evLion



subdivson same limit surface
at equilium with more patches



eoion evolio



(a) (b)

Figure 3.7. Model subdivision to increase the degrees of freedom : (a) evolution of
the initial mesh, and (b) the corresponding limit surface evolution perceived by the
user.


This linear system of equations is solved iteratively between each time step using the

conjugate gradient method [34, 75].

3.3.6 Model Subdivision

The initialized model grows dynamically according to the equation of motion

(Eqn.2.5). The degrees of freedom of the initialized model is equal to the number

of control vertices in the initial mesh as mentioned earlier. When an equilibrium is

achieved for the model, the number of control vertices can be increased by replacing

the original initial mesh by a new initial mesh obtained through one step of Catmull-

Clark subdivision. This increases the number of degrees of freedom to represent the

same (deformed) smooth limit surface and a new equilibrium position for the model

can be obtained. This process is depicted schematically in Fig.3.7. Model subdivision

might be needed to obtain a very localized effect on a smooth limit surface. For a









shape recovery application, one may start with a very simple initial model, and when

an approximate shape is recovered, the degrees of freedom can be increased to obtain

a new equilibrium position for the model with a better fit to the given data set. The

error of fit criteria for the discrete data is based on distance between the data points

and the points on the limit surface where the corresponding springs are attached.

In the context of image-based forces, if the model energy does not change between

successive iterations indicating an equilibrium for the given resolution, the degrees of

freedom for the model can be increased by the above-mentioned replacement scheme

until the model energy is sufficiently small and the change in model energy between

successive iterations becomes less than a pre-specified tolerance.

3.4 Generalization of the Approach

The proposed approach can be generalized for other approximating subdivision

schemes. However, a more general approach is presented in Chapter 5 for deriving the

dynamic framework. Dynamic Catmull-Clark subdivision surface model is reformu-

lated in Section 5.2 using this general approach. A dynamic framework for another

popular approximating subdivision scheme namely, Loop's subdivision scheme, is also

discussed in Chapter 5.
















CHAPTER 4
DYNAMIC BUTTERFLY SUBDIVISION SURFACES


In the previous chapter, a dynamic framework has been presented for an approx-

imating subdivision scheme namely, Catmull-Clark subdivision scheme. In this chap-

ter, a dynamic framework is developed for (modified) butterfly subdivision scheme,

which is an interpolatory subdivision technique. First, a brief overview of the (mod-

ified) butterfly subdivision scheme is presented. Next, a local parameterization tech-

nique for the limit surface obtained via (modified) butterfly subdivision is discussed.

This parameterization scheme is then used to derive the dynamic model. The imple-

mentation details are also discussed.

4.1 Overview of the (Modified) Butterfly Subdivision Scheme

The butterfly subdivision scheme [30], like many other subdivision schemes used

in geometric design literature/applications, starts with an initial triangular mesh

which is also known as the control mesh. The vertices of the control mesh are known

as the control points. In each step of subdivision, the initial (control) mesh is refined

through the transformation of each triangular face into a patch with four smaller

triangular faces. After one step of refinement, the new mesh in the finer level retains

the vertices of each triangular face in the previous level and hence, interpolates the

coarser mesh in the previous level. In addition, every edge in each triangular face





























(a)


Figure 4.1. (a) The control polygon with triangular faces;
one subdivision step using butterfly subdivision rules.




6 v -w
V6 V V
1 5




4 12 3




Vj V2j VI -w


(b)


(b) mesh obtained after


-w


Figure 4.2. (a) The contributing vertices in the j-th level for the vertex in the j+1-
th level corresponding to the edge between v{ and v2; (b) the weighing factors for
different vertices.


I






I



"

:. .-









is spilt by adding a new vertex whose position is obtained by an affine combination

of the neighboring vertex positions in the coarser level. For instance, the mesh in

Fig.4.1(b) is obtained by subdividing the initial mesh shown in Fig.4.1(a) once. It

may be noted that all the newly introduced vertices corresponding to the edges in the

original mesh have degree 6, whereas the position and degree of the original vertices

do not change in the subdivided mesh.

In the original butterfly scheme, the new vertices corresponding to the edges in

the previous level are obtained using an eight-point stencil as shown in Fig.4.2(a). To

show how this stencil is used, a vertex to be introduced and the contributing vertices

from the stencil are highlighted in Fig.4.1(a). The name of the scheme originated

from the "butterfly"-like configuration of the contributing vertices. The weighing

factors for different contributing vertex positions are shown in Fig.4.2(b). The vertex

e~1 in the j + 1-th level of subdivision, corresponding to the edge connecting vertices

v{ and v2 at level j, is obtained by



ejf' = 0.5(vi + vN) + 2w(v~ + vJ) w(v3 + v) + v' + vj), (4.1)



where 0 < w < 1, and vj denotes the position of the i-th vertex at the j-th level.

The butterfly subdivision scheme produces a smooth C' surface in the limit

except at the extraordinary points corresponding to the extraordinary vertices (ver-

tices with degree not equal to 6) in the initial mesh [109]. All the vertices introduced

through subdivision have degree 6, and therefore, the number of extraordinary points

in the smooth limit surface equals the number of extraordinary vertices in the initial








































(a) (b)

Figure 4.3. (a) The weighing factors of contributing vertex positions for an edge
connecting two vertices of degree 6; (b) the corresponding case when one vertex is of
degree n and the other is of degree 6.









mesh. Recently, this scheme has been modified by Zorin et al. [109] to obtain better

smoothness properties at the extraordinary points. This modified scheme is known

as modified butterfly subdivision technique. In Zorin et al. [109], all the edges have

been categorized into three classes : (i) edges connecting two vertices of degree 6 (a

10 point stencil, as shown in Fig.4.3(a), is used to obtain the new vertex positions

corresponding to these edges), (ii) edges connecting a vertex of degree 6 and a vertex

of degree n 0 6 (the corresponding stencil to obtain new vertex position is shown in

Fig.4.3(b), where q = .75 is the weight associated with the vertex of degree n 0 6,

and si = (0.25 + cos(2ri/n) + 0.5cos(47ri/n))/n, i = 0, 1,..., n 1, are the weights

associated with the vertices of degree 6), and (iii) edges connecting two vertices of

degree n 5 6. The last case can not occur except in the initial mesh as the newly

introduced vertices are of degree 6, and the new vertex position in this last case is

obtained by averaging the positions obtained through the use of stencil shown in

Fig.4.3(b) at each of those two extraordinary vertices.

4.2 Formulation

In this section, a systematic formulation of the dynamic framework for the modi-

fied butterfly subdivision scheme is presented. Unlike the approximating schemes, the

limit surface obtained via modified butterfly subdivision does not have any analytic

expression even in a regular setting. Therefore, derivation of a local parameterization

suitable for developing the dynamic framework for the limit surface is the key issue

here which is discussed next.

















Figure 4.4. The smoothing effect of the subdivision process on the triangles of the
initial mesh.


4.2.1 Local Parameterization

The smooth limit surface defined by the modified butterfly subdivision technique

is of arbitrary topology where a global parameterization is impossible. Nevertheless,

the limit surface can be locally parameterized over the domain defined by the initial

mesh following a similar approach described in Lounsbery et al. [53]. The idea

is to track any arbitrary point on the initial mesh across the meshes obtained via

the subdivision process (see Fig.4.4 and Fig.4.5), so that a correspondence can be

established between the point being tracked in the initial mesh and its image on the

limit surface.

The modified butterfly subdivision scheme starts with an initial mesh consisting

of a set of triangular faces. The recursive application of the subdivision rules smoothes

out each triangular face, and in the limit, a smooth surface is obtained which can also

be considered as a collection of smooth triangular patches. The subdivision process

and the triangular decomposition of the limit surface is depicted in Fig.4.4. Note

that, the limit surface can be represented by the same number of smooth triangular

patches as that of the triangular faces in the initial mesh. Therefore, the limit surface


- 0 0









s can be expressed as

s = Sk, (4.2)
k=1

where n is the number of triangular faces in the initial mesh and sk is the smooth

triangular patch in the limit surface corresponding to the k-th triangular face in the

initial mesh.

The stage is now set for describing the parameterization of the limit surface over

the initial mesh. The process is best explained through the following example. A

simple planar mesh shown in Fig.4.5(a) is chosen as the initial mesh. An arbitrary

point x inside the triangular face abc is tracked over the meshes obtained through

subdivision. The vertices in the initial mesh are darkly shaded in Fig.4.5. After

one step of subdivision, the initial mesh is refined by addition of new vertices which

are lightly shaded. Another subdivision step on this refined mesh leads to a finer

mesh with introduction of new vertices which are unshaded. It may be noted that

any point inside the smooth triangular patch in the limit surface corresponding to the

face abc in the initial mesh depends only on the vertices in the initial mesh which

are within the 2-neighborhood of the vertices a, b and c due to the local nature of

the subdivision process. For example, the vertex d, introduced after first subdivision

step, can be obtained using the 10 point stencil shown in Fig.4.3(a) on the edge ab.

All the contributing vertices in the initial mesh are within the 1-neighborhood of the

vertices a and b. A 10 point stencil can be used again in the next subdivision step

on the edge db to obtain the vertex g. Some of the contributing vertices at this level

of subdivision, for example, the (lightly shaded) 1-neighbors of the vertex b (except



























(a) (b)


Figure 4.5. Tracking a point x through various levels of subdivision : (a) initial
mesh, (b) the selected section (enclosed by dotted lines) of the mesh in (a), after one
subdivision step, (c) the selected section of the mesh in (b), after another subdivision
step.









d and e) in Fig.4.5(b), depend on some vertices in the initial mesh which are within

the 2-neighborhood of the vertices a, b and c in the initial mesh.

In the rest of the formulation, superscripts are used to indicate the subdivision

level. For example, v-7 denotes the collection of vertices at level j which control

the smooth patch in the limit surface corresponding to the triangular face uvw at

the j-th level of subdivision. Let vob be the collection of vertices in the initial mesh

which are within the 2-neighborhood of the vertices a, b and c (marked black in

Fig.4.5(a)). Let the number of such vertices be r. Then, the vector v0bc, which is the

concatenation of the (x, y, z) positions for all the r vertices, is of dimension 3r. These

r vertices control the smooth triangular patch in the limit surface corresponding to

the triangular face abc in the initial mesh. Now, there exists four subdivision matrices

(Aa,),, (Anac)t, (Aa.e), and (Ab), of dimension (3r, 3r) such that



v = (abc),rVac

Vbd = (Aoa)iVabc,

vfe = (Aoate),ac,

v, = (Aabc)mva (4.3)


where the subscripts t, 1, r and m denote top, left, right and middle triangle positions

respectively (indicating the relative position of the new triangle with respect to the

original triangle), and Vid, V vi and vYef are the concatenation of the (x, y, z)

positions for the vertices in the 2-neighborhood of the corresponding triangle in the








newly obtained subdivided mesh. The new vertices in this level of subdivision are

lightly shaded in Fig.4.5(b). The 2-neighborhood configuration of the vertices in the

newly obtained triangles is exactly the same as that of the original triangle, hence

local subdivision matrices are square and the vector dimensions on both sides of

Eqn.4.3 are the same. This concept is further illustrated in Fig.4.6.

Carrying out one more level of subdivision, a new set of vertices which are

unshaded in Fig.4.5(c) are obtained along with the old vertices. Adopting a similar

approach as in the derivation of Eqn.4.3, it can be shown that



vd = (Abed)tVed

Vbg = (Abed)ived

6vh = (Abed),vted

vjh, = (Abed)mVed (4.4)


The relative position of the triangular face dgi in Fig.4.5(c) with respect to the

triangular face bed is topologically the same as of the triangular face adf in Fig.4.5(b)

with respect to the triangular face abc. Therefore, one can write (Abed), = (Ab~)t,.

Using similar reasoning, Eqn.4.4 can be rewritten as



v~, = (AMd)tved = (Aabc)(tVed

vbh = (Abed)Ved = (Aab)rVld

v2 = (Abed)rVd = (Aabc)rVed
Vethbe





67






---------------.?'!' ,b ,






















(b)
4 v




















































Figure 4.6. Different set of control points at a subdivided level obtained by applying
different subdivision matrices on a given set of control points in a coarser mesh.









vg2 = (Ad)mVV = (Aak)mV d. (4.5)


Combining Eqn.4.3 and Eqn.4.5, it can be shown that



vd = (Ask)t(Aok),vo,

= (Aa~)l(A.kc)iVO,

vh = (Ao.),(Aai),v ,

vh = (A.)m(Aa ),v.bc (4.6)


Let x be a point with barycentric coordinates (ac fac, Y7 ) inside the trian-

gular face abc. When the initial mesh is subdivided, x becomes a point inside the

triangular face bed with barycentric coordinates (-i,. r.,.,y*). Another level of

subdivision causes x to be included in the triangular face dgi with barycentric coor-

dinates (admi, },"y/,i). Let sj denote the j-th level approximation of the smooth

triangular patch Sab in the limit surface corresponding to the triangular face abc in

the initial mesh. Now v0 can be written as

r r r
c = C ., b, ,..., b,,c,,..


where the subscripts x, y and z indicate the x, y and z coordinates respectively of the

corresponding vertex position. The expressions for v d and vi can also be written








in a similar manner. Next, the matrix Bb can be constructed as follows:

r r r
"oa, b, 7, 0,... 0,0,...,0,0,..., 0
r r r
Bob(x) = 0- .,,
0,1 ..., 0, bc r bc 7el ,0, .,0,0,..., 0
r r
0F. 0,0,0.. a5b- 0 7l 0,.. ,


The matrices B1d and B' can also be constructed in a similar fashion. Now sobc(X),

s9c(x), and sla(x) can be written as


Sabc(X) B0
Sk(X) = Bn (x)vab,
Sabc(x) = B' (x)v = Bd(x)(A )
= Bbed(X)Vbed = Bbed(X)(Aabc)ljVab,
s(x) = B' (x)v, = B'd(X)(Aabc)tVd g B (X)(Aabc)t(Aasc)iVac.

(4.7)



Proceeding in a similar way, the expression for s8c(x), j-th level approximation

of sa(x), is given by



sI(x) = B,,,(x)(Abc)... (Aak)t(Aabc) Ve c

= Bu,(x)(Alb-)vbc

= Bjc(x)vObc, (4.8)


where x is inside the triangular face uvw at level j (with an assumption that uvw is the

triangular face in the middle with respect to its coarser level original triangular face in









the previous level), (Ab) = (Ab,,)m... (A(Aa,,)(Aa,) and Bi (x) = Bj,(x)(A- ).

It may be noted that the sequence of applying (Aabs)t, (Aask), (Abc), and (Aabc)m

depends on the triangle inside which the tracked point x falls after each subdivision

step. Finally, the local parameterization process can be completed by writing



Sab(x) = (lim BI.(x))v = Back(x)v (4.9)


In the above equation, Bae is the collection of basis functions at the vertices of

v b. It may also be noted that the modified butterfly subdivision scheme is a sta-

tionary subdivision process, and hence new vertex positions are obtained by affine

combinations of nearby vertices. This guarantees that each row of the matrices

(Anb),, (Aabe), (Aab), and (Aatc)m sums to one. The largest eigenvalue of such ma-

trices is 1 and therefore the limit in Eqn.4.9 exists. Now, assuming the triangular

face abc is the k-th face in the initial mesh, Eqn.4.9 can be rewritten as



Sk(x) = Bk(x)v0 = Bk(x)Akp, (4.10)


where p is the concatenation of the (x,y,z) positions of all the vertices in the initial

mesh and the matrix Ak, when post-multiplied by p, selects the vertices vA controlling

the k-th smooth triangular patch in the limit surface. If there are t vertices in the

initial mesh and r of them control the k-th patch, then p is a vector of dimension 3t,

Ak is a matrix of dimension (3r, 3t), and Bk(x) is a matrix of dimension (3,3r).









Combining Eqn.4.2 and Eqn.4.10, it can be shown that



s(x) =( Bk(x)Ak)p = J(x)p, (4.11)
k=1


where J, a matrix of dimension (3,3t), is the collection of basis functions for the

corresponding vertices in the initial mesh. The vector p is also known as the degrees

of freedom vector of the smooth limit surface s.

4.2.2 Dynamics

An expression of the limit surface obtained via modified butterfly subdivision

process is derived in the last section. Now a dynamic framework for the limit surface

can be derived using this expression. The derivation of the dynamic model from this

point is exactly similar in spirit as that of the dynamic Catmull-Clark subdivision

model presented in Section 3.2.4. The vertex positions in the initial mesh defining

the smooth limit surface s are treated as a function of time in order to embed the

modified butterfly subdivision model in a dynamic framework. The velocity of this

surface model can be expressed as



9(x,p) = J(x)p, (4.12)


where an overstruck dot denotes a time derivative and x e So, So being the domain

defined by the initial mesh. As in the case of the dynamic Catmull-Clark subdivision

surfaces, the Lagrangian equation of motion (Eqn.2.1) is used to derive the dynamic

modified butterfly subdivision surface model.








Let p be the mass density function of the surface. Then the kinetic energy of

the surface is



T= j ( p(x)s(x)x) = p Mp, (4.13)


where (using Eqn.4.12) M = fx5so p(x)JT(x)J(x)dx is the mass matrix of dimen-

sion (3t, 3t). Similarly, let 7 be the damping density function of the surface. The

dissipation energy is


F = 1 y(x)s(x)s(x)dx -= TDp, (4.14)
2 JeSO 2


where D= JxES y(x)JT(x)J(x)dx is the damping matrix of dimension (3t, 3t). The

potential energy of the smooth limit surface can be expressed as



U pTKp, (4.15)


where K is the stiffness matrix of dimension (3t, 3t) obtained by assigning various

internal energies to a discretized approximation of the limit surface and is detailed

in Section 4.3.

Using the expressions for the kinetic, dissipation and potential energy in Eqn.2.1,

the same motion equation as in Eqn.2.5 can be obtained. The generalized force vector









f, can be obtained is a similar fashion described in Section 3.2.4 and is given by



f,= feso J'(x)f(x, t)dx. (4.16)


4.2.3 Multilevel Dynamics

The initial mesh of the dynamic modified butterfly subdivision surface model

can be subdivided to increase the degrees of freedom for model representation. The

situation is exactly similar to that of the dynamic Catmull-Clark subdivision surface

model as described in Section 3.2.5 and a similar example is used for illustration.

After one step of modified butterfly subdivision, the initial degrees of freedom p (refer

to Eqn.4.11 and Eqn.4.12) in the dynamic system is replaced by a larger number of

degrees of freedom q, where q = Bp. B is a global subdivision matrix of size (3s, 3t)

whose entries are uniquely determined by the weights used in the modified butterfly

subdivision scheme (see Section 4.1 for the weights). Thus, p, expressed as a function

of q, can be written as

p = (BTB)-IBTq = Btq, (4.17)


where Bt = (BTB)-'BT. Therefore, Eqn.4.11 and Eqn.4.12 is modified as



s(x) = (J(x)Bt)q, (4.18)


and


s(x,q) = (J(x)Bt)il,


(4.19)

































(b)

Figure 4.7. (a) An initial mesh, and (b) the corresponding limit surface. The domains
of the shaded elements in the limit surface are the corresponding triangular faces in
the initial mesh. The encircled vertices in (a) are the degrees of freedom for the
corresponding element.


respectively. The motion equation, explicitly expressed as a function of q, can be

written as


Mq + Dq + Kq = fq,


(4.20)


where Mq = fxEsl p(x)(Bt)TJT(x)J(x)Btdx, S' being the domain defined by the

newly obtained subdivided mesh. The derivations of Dq, K, and f, are similar.


t ----









4.3 Finite Element Implementation

In the previous section, the smooth limit surface obtained via modified butterfly

subdivision process is expressed as a function of the control vertex positions in its

initial mesh; mass and damping distribution, internal deformation energy and forces

are assigned to the limit surface in order to develop the corresponding physical model.

In this section, the implementation of this physical model is detailed. It may be noted

that even though the formulation of the dynamic framework for modified butterfly

subdivision scheme is similar to that of Catmull-Clark subdivision scheme once the

local parameterization is derived, the finite element implementation is totally different

as the limit surface does not have any analytic expression in case of the modified

butterfly subdivision scheme.

In Section 4.2 it was pointed out that the smooth limit surface obtained by the

recursive application of the modified butterfly subdivision rules can be represented by

a set of smooth triangular patches, each of which is represented by a finite element.

The shape (basis) function of this finite element is obtained by smoothing a hat

function through repeated application of the modified butterfly subdivision rules.

The number of finite elements in the smooth limit surface is equal to the number

of triangular faces in the initial mesh as mentioned earlier (refer Fig.4.4 and 4.7).

A detailed discussion on how to derive the mass, damping and stiffness matrices

for these elements is presented next. These elemental matrices can be assembled

to obtain the global physical matrices M, D and K, and a numerical solution to

the governing second-order differential equation as given by Eqn.2.5 can be obtained









using finite element analysis techniques [42]. The same example as in Section 4.2

(refer Fig.4.5) is used to develop the related concepts. The concept of elements along

with the control vertices and their corresponding domains is further illustrated in

Fig.4.7.

Now it will be shown how to derive the mass, damping and stiffness matrices

for the element corresponding to the triangular face abc in Fig.4.5. The derivations

hold for any element in general.

4.3.1 Elemental Mass and Damping matrices

The mass matrix for the element given by Sbc (Eqn.4.9) can be written as



Mb = f p(x) {B.(x)}T{Ba(x)}dx. (4.21)



However, from Eqn.4.9 it is known that the basis functions corresponding to the

vertices in vk which are stored as entries in Ba, are obtained as a limiting process.

These basis functions do not have any analytic form, hence computing the inner

product of such basis functions as needed in Eqn.4.21 is a challenging problem. In

Lounsbery et al. [53], an outline is provided on the computation of these inner

products without performing any integration. In this dissertation, a different yet

simpler approach is developed to solve this problem. The smooth triangular patch in

the limit surface corresponding to the face abc in the initial mesh is approximated by

a triangular mesh with 43 faces obtained after j levels of subdivision of the original

triangular face abc (each subdivision step splits one triangular face into 4 triangular








faces). Then the mass matrix can be expressed as



Mbc= f L (x){BJi(x)}T{B,(x)}dx. (4.22)


The j-th level approximation of the corresponding basis functions can be explicitly

evaluated (refer Eqn.4.8 for an expression of B1 ). An important point to note is

that Eqn.4.8 involves several matrix multiplications and hence can be very expensive

to evaluate. However, the matrix (AJ()(= (Aab) ... (Aoac)t(Aa.)t) in Eqn.4.8

encodes how vertices in the 2-neighborhood of the triangular face uvw at level j are

related to the vertices in the 2-neighborhood of the triangular face abc in the initial

mesh. In the implementation, how a new vertex is obtained from the contributing

vertices in its immediate predecessor level is tracked. The information stored in (Abc)

can be obtained by tracking the contributions from level j to level 0 recursively. Thus

the entries of (Abc) can be determined without forming any local subdivision matrices

and thereby avoiding subsequent matrix multiplications.

By choosing a sufficiently high value of j, a reasonably good approximation of

the elemental mass matrices is achieved. The computations involved in evaluating

the integrals in Eqn.4.22 are eliminated by using discrete mass density function which

has non-zero values only at the vertex positions of the j-th subdivision level mesh.

Therefore, the approximation of the mass matrix for the element can be written as



Ma,. = k ),(V, T B v) B )}, (4.23)
i=1









where k is the number of vertices in the triangular mesh with 43 faces. This approxi-

mation has been found to be very effective and efficient for implementation purposes.

The elemental damping matrix can be obtained in a similar fashion.

4.3.2 Elemental Stiffness Matrix

The internal energy is assigned to each element in the limit surface, thereby

defining the internal energy of the smooth subdivision surface model. A similar

approach as in the derivation of the elemental mass and damping matrix is taken

and the internal energy is assigned to a j-th level approximation of the element.

Three types of internal energy are used tension, stiffness and spring energy.

For the examples used throughout the paper, this energy at the j-th level of approx-

imation can be written as



Ekac Ej- = (Ek), + (Ek),, + (Ef )p,, (4.24)



where the subscripts t, st and sp denote tension, stiffness and spring respectively.

The expression for the tension energy, which is essentially equivalent to the first

order strain (membrane) energy [96], is



(E )t = ktEvi vI2
2 0
= 1ktv .T (K )Jt(Vk}, (4.25)



where kt is a constant, vf and vi, the 1-th and m-th vertex in the j-th level mesh,

are in the 1-neighborhood of each other, Q is the domain defined by all such vertex








pairs, and vb is the concatenation of the (x,y,z) positions of all the vertices in the

j-th subdivision level of the triangular face abc in the initial mesh.

Similarly, the expression for stiffness energy, which is equivalent to the second

order strain (thin plate) energy [96], can be written as



(E ) k., Iv -2v + +vi2
n
1= k{v (L v }, (4.26)



where vJ, vi and v- are the three vertices of a triangular face. The summation in-

volves three terms corresponding to each triangular face, and is over all the triangular

faces in the mesh at the j-th level of subdivision.

Finally, a spring energy term is added which is given by


1 (kim) (I vi v e1)
(E bc)sp 2 i -- vm (vj VIm)

S2{v } (Ki ),,{v }, (4.27)



where (km),, is the spring constant, Q is the domain as in Eqn.4.25 and ,im is the

natural length of the spring connected between vi and vi. It may be noted that

the entries in (Kb ), depend on the distance between the connected vertices and

hence, (Kib),, unlike other elemental matrices, is a function of time which needs to

be recomputed in each time step.








Combining the expressions for tension, stiffness and spring energy, it can be

shown that


Ekb ab c} {kt(Ka + ket(K c),+(K ),j}{v }

= {v }(K abc) tva
T1 j T j j
2{(A f){v }} (K bc)(A bc){V }

I T KT j )j (4.28)
= {v } {(AIse) (K )(A ac)} vac}, (4.28)


where (Abe) and vb are same as in Eqn.4.8. Therefore, the expression for the

elemental stiffness matrix is given by


Ka- = (Ab)T(Ka)(Ab). (4.29)


It may be noted that the matrix multiplications for constructing Kabc are avoided in

the implementation by following the same technique described in Section 4.3.1.

4.3.3 Other Implementation Issues

The techniques of applying synthesized forces on the smooth limit surface ob-

tained via modified butterfly subdivision process are similar to those described in

Section 3.3.4 in the context of dynamic Catmull-Clark subdivision surface model.

The techniques of model subdivision and discrete dynamic equation derivation are

also the same.





82


4.4 Generalization of the Approach

The proposed approach of deriving the dynamic framework for modified butter-

fly subdivision scheme is generalized for other interpolatory subdivision schemes in

the next chapter, Chapter 5.

















CHAPTER 5
UNIFIED DYNAMIC FRAMEWORK FOR SUBDIVISION SURFACES



Dynamic framework for specific subdivision schemes was presented in the last

two chapters. In Chapter 3, Catmull-Clark subdivision scheme, a representative of

the approximating subdivision schemes, is embedded in a physics-based modeling en-

vironment. In Chapter 4, a dynamic framework is provided for the modified butterfly

subdivision scheme, a popular interpolatory subdivision technique. In this chapter,

a unified approach is presented for embedding any subdivision scheme in a dynamic

framework.

5.1 Overview of the Unified Approach

The key concept of the unified approach is the ability to express the smooth

limit surface generated by a subdivision scheme as a collection of a single type of

finite elements. The type of this finite element depends on the subdivision scheme

involved. For example, the limit surface generated by the Catmull-Clark subdivision

scheme can be expressed as a collection a quadrilateral finite element patches. This

approach differs form the one taken in Chapter 3 where the limit surface was a

collection of normal elements (corresponding to the quadrilateral bicubic B-spline

patches away from the extraordinary points) and special elements (corresponding to









the n-sided patches near extraordinary vertices of degree n). This concept will be

further elaborated in the next few sections.

The unified approach adopts two different strategies depending on the subdivi-

sion scheme involved. One strategy is for the approximating schemes, whose limit

surface is some type of B-spline surface away from the extraordinary points, and

the other is for the interpolatory subdivision schemes, where the limit surface does

not have any analytical expression. It may be recalled that recursive subdivision is

a method for smoothing polyhedra, and hence corresponding to each n-sided face

in the control mesh, there would be a smooth n-sided patch in the limit surface

in general. For example, the control mesh (after at most one subdivision step) of

the Catmull-Clark subdivision scheme consists of quadrilateral faces, which when

smoothed through subdivision lead to quadrilateral patches in the limit surface (see

Fig.5.1(b)). In the case of approximating subdivision schemes, the faces which do

not contain any extraordinary vertex lead to spline patches with analytic expressions.

On the other hand, faces which have one extraordinary vertex, lead to patches whose

analytic expressions are difficult to determine. Continuing the example with the

Catmull-Clark subdivision scheme, the faces without extraordinary vertices lead to

quadrilateral bicubic B-spline patches in the limit, whereas the faces with one extraor-

dinary vertex lead to quadrilateral patches whose analytic expressions are difficult

to obtain. In Chapter 3, n such adjacent patches corresponding to an extraordinary

vertex of degree n are effectively combined together to obtain an expression for the re-

sulting n-sided patch (see Fig.5.1(a)). Recently, Stam [90, 91] has introduced schemes









for exactly evaluating the limit surface, even near the extraordinary points, for ap-

proximating subdivision schemes. Analytic expressions for the patches resulting from

smoothing faces with an extraordinary vertex can be obtained using these schemes.

Therefore, the limit surface can be expressed as a collection of a single patch type,

and consequently a single type of finite element can be used to provide the dynamic

framework for the approximating subdivision schemes. This unified framework for

approximating subdivision schemes is fully worked out in this chapter for Catmull-

Clark and Loop subdivision schemes and a general outline is also presented on how

to carry out the same strategy for other approximating subdivision schemes.

The limit surface resulting from an interpolatory subdivision scheme does not

have an analytic expression in general, and hence a strategy similar to the one adopted

in Chapter 4 needs to be used. An outline of the strategy to be used for interpolatory

subdivision schemes is also presented in this chapter.

5.2 Unified Dynamic Framework for Catmull-Clark Subdivision Surfaces

A systematic formulation along with the implementation details of the unified

dynamic framework for Catmull-Clark subdivision surfaces is presented in this sec-

tion. The key difference between the framework developed in Chapter 3 and the one

presented here is the representation of the limit surface, which leads to different types

of finite elements used for developing the dynamic framework. This is illustrated with

a schematic diagram in Fig.5.1.




Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EZKAYFT03_SDXRFY INGEST_TIME 2013-01-23T14:39:46Z PACKAGE AA00012984_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES