A Novel FEMBased Dynamic Framework For Subdivision Surfaces
Chhandomay Mandal* Hong Qint Baba C. Vemuri*
*Department of Computer and Information Science and Engineering
University of Florida
tDepartment of Computer Science
State University of New York at Stony Brook
Abstract
Subdivision surfaces have been extensively used to model smooth
shapes of arbitrary topology. Recursive subdivision on an user
defined initial control mesh generates a visually pleasing smooth
surface in the limit. However, users have to carefully specify
the initial mesh and/or painstakingly manipulate the control ver
tices at different levels of subdivision hierarchy to satisfy various
functional and aesthetic requirements in the limit surface. This
modeling drawback results from the lack of direct manipulation
tools for the limit surface. In this paper, we integrate physics
based modeling techniques with geometric subdivision method
ology and present an unified approach for arbitrary subdivision
schemes. Our dynamic framework permits users to directly manip
ulate the limit subdivision surface via physicsbased "force" tools.
The key contribution of this unified approach is to formulate the
smooth limit surface of any subdivision scheme as a single type
of novel finite elements. The geometric and dynamic features of
our subdivisionbased finite elements depend on the subdivision
scheme involved. We present our finite element method (FEM) and
formulation for the modified butterfly and CatmullClark subdivi
sion schemes, and further generalize our dynamic framework for
any subdivision scheme. Our FEMbased approach significantly
advances the stateoftheart of physicsbased geometric modeling
because (1) our dynamic framework provides a universal physics
based solution to any subdivision scheme beyond frequentlyused
and popular splinelike subdivision techniques; (2) we systemat
ically devise a natural mechanism that allows users to intuitively
deform any subdivision surface; (3) we represent the smooth limit
surface of any subdivision scheme using a single type of novel
subdivisionbased finite elements. Our experiments demonstrate
that the new unified FEMbased framework promises a greater po
tential of subdivision techniques for solid modeling, finite element
analysis, and engineering design.
1 INTRODUCTION
Efficiently modeling and manipulating smooth surfaces of arbitrary
topology is a grand challenge to scientists and engineers in solid
modeling, computeraided design, and interactive graphics. The
recursive subdivision scheme, which produces a visually pleasing
smooth surface in the limit by repeated application of a fixed set of
refinement rules on an userspecified initial control mesh, is well
suited for this task. Despite the prevalence of diverse subdivision
schemes in the computer graphics and geometric modeling litera
ture, it is almost impossible to manipulate the limit surface (ob
tained through procedurebased subdivision) in a direct and natural
way. The current stateoftheart only permits modelers to interac
tively obtain the desired effects on the smooth surface by kinemati
cally manipulating the control vertices at various levels of subdivi
sion hierarchy. In this paper, we address the challenging problem of
directly manipulating the limit subdivision surface at arbitrary loca
tions/areas, and offer a novel solution to this problem by embedding
purely geometric subdivision schemes in a physicsbased modeling
framework. Unlike the existing geometric solutions that only al
low operations on control vertices, our methodology and algorithms
permit users to physically modify the shape of subdivision surfaces
at desired locations via application of forces. Consequently, this
gives the user an intuitive and natural feeling that is uniquely pro
duced while modeling with real clay/playdough. Additionally, we
will demonstrate that the proposed model can efficiently recover
shapes from a cloud of 3D points. However, prior to the technical
details of our novel scheme, we shall briefly review the previous
work on subdivision surfaces.
1.1 Background
In [4], Chaikin first introduced the concept of subdivision to the
modeling community for generating a smooth curve from an arbi
trary control polygon. During the last two decades, a wide vari
ety of subdivision schemes for modeling smooth surfaces of arbi
trary topology have been derived following Chaikin's pioneering
work on curve generation. The existing subdivision schemes can
be broadly categorized into two distinct classes namely, (1) approx
imating subdivision techniques, and (2) interpolating subdivision
techniques.
Among the approximating schemes, the techniques of Doo and
Sabin [6] and Catmull and Clark [3] generalize the idea of obtain
ing uniform biquadratic and bicubic Bspline patches, respectively,
from a rectangular control mesh. In [3], Catmull and Clark devel
oped an algorithm for recursively generating a smooth surface from
a polyhedral mesh of arbitrary topology. The CatmullClark subdi
vision surface, defined by an arbitrary initial mesh, can be reduced
to a set of standard Bspline patches except at a finite number of
degenerate points. In [12], Loop presented a similar subdivision
scheme based on the generalization of quartic triangular Bsplines
for triangular meshes. Hoppe et al. [11] further extended Loop's
work to produce piecewise smooth surfaces with selected disconti
nuities. Halstead et al. [10] proposed an algorithm to construct a
CatmullClark subdivision surface that interpolates the vertices of
a mesh of arbitrary topology. Peters and Reif [16] proposed a sim
ple subdivision scheme for smoothing polyhedra. Most recently,
nonuniform DooSabin and CatmullClark surfaces that generalize
nonuniform tensorproduct Bspline surfaces to arbitrary topolo
gies were introduced by Sederberg et al. [22]. All the aforemen
tioned schemes generalize recursive subdivision schemes for gener
ating limit surfaces with a known parameterization. Various issues
involved with the use of these approximating subdivision schemes
for character animation were discussed at length by DeRose et al.
[5].
The most wellknown interpolationbased subdivision scheme is
the "butterfly" algorithm proposed by Dyn et al. [8]. Butterfly
method, like other subdivision schemes, makes use of a small num
ber of neighboring vertices for subdivision. It requires simple data
structures and is rather straightforward to implement. Neverthe
less, it needs a topologically regular setting of the initial (control)
mesh in order to obtain a smooth C1 limit surface. Zorin et al.
[25] has developed an improved interpolatory subdivision scheme
(which we call the modified butterfly scheme) that retains the sim
plicity of the butterfly scheme and results in much smoother sur
faces even from irregular initial meshes. These interpolatory sub
division schemes have extensive applications in wavelets on man
ifolds, multiresolution decomposition of polyhedral surfaces, and
multiresolution editing.
The derivation of various mathematical properties of the smooth
limit surface generated by the subdivision algorithms is rather com
plex. Doo and Sabin [7] first analyzed the smoothness behavior of
the limit surface using the Fourier transform and an eigenanalysis
of the subdivision matrix. Ball and Storry [1, 2] and Reif [20] fur
ther 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. Spe
cific subdivision schemes were analyzed by Schweitzer [21], Habib
and Warren [9], Peters and Reif [17] and Zorin [26]. Most recently,
Stam [23] presented a method for exact evaluation of CatmullClark
subdivision surfaces at arbitrary parameter values.
1.2 Motivation
Although recursive subdivision surfaces are extremely powerful for
representing smooth geometric shapes of arbitrary topology, they
constitute a purely geometric representation, and furthermore, con
ventional geometric modeling with subdivision surfaces may be dif
ficult for effectively representing and deforming highly complicated
objects. For example, modelers are faced with the tedium of in
direct shape modification and refinement through timeconsuming
operations on a large number of (oftentimes irregular) control ver
tices when utilizing typical subdivisionbased modeling techniques.
Despite the advent of many modern 3D graphics interaction tools,
these indirect geometric operations remain nonintuitive and labo
rious in general. In addition, oftentimes 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 desired while requiring ge
ometric objects to satisfy global smoothness criteria in solid mod
eling and/or interactive graphics applications. In contrast, physics
based modeling provides a superior approach to shape modeling
that can overcome most of the limitations associated with tradi
tional geometric modeling approaches. Freeform deformable mod
els governed by the laws of continuum mechanics are of particular
relevance in this context. Dynamic models respond to externally
applied forces in a very intuitive manner. The dynamic formula
tion 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 com
plex geometries. Furthermore, 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 model
ing criteria, and geometric constraints relevant to shape design can
also be properly imposed. The dynamic approach subsumes all of
the aforementioned modeling capabilities in a formulation which
grounds everything in realworld physical behavior.
Freeform deformable models were first introduced to the mod
eling community by Terzopoulos et al. [24], and were improved by
a number of researchers over the years. Qin and Terzopoulos [18]
developed DNURBS which are very sophisticated physicsbased
models suitable for representing a wide variety of freeform as well
as standard analytic shapes. The DNURBS have the advantage of
interactive and direct manipulation of NURBS curves and surfaces,
resulting in physically meaningful thus intuitively predictable mo
tion and shape variation. However, a severe limitation of the exist
ing deformable models, including DNURBS, is that they are de
fined on a rectangular parametric domain. Therefore, it can be very
difficult to model surfaces of arbitrary genus using these models.
Subdivision schemes, in contrast, can model complex surfaces of
arbitrary topology, and hence are a good candidate for developing a
novel physicsbased model where the modeler can directly manip
ulate the (complicated) limit surface in an intuitive way.
Previously we had introduced dynamic CatmullClark subdivi
sion surfaces [14, 19] where the smooth limit surface generated by
the CatmullClark subdivision scheme was embedded in a physics
based modeling framework. The current research differs signifi
cantly from our prior work because the approach taken in this pa
per is much more general. It aims to develop a systematic and
universal mechanism with which any subdivision scheme can be
formulated within the physicsbased framework. The critical math
ematical technique we resort to is finite element analysis. We shall
first formulate a dynamic representation and equation for an inter
polatory subdivision scheme the modified butterfly subdivision
method where the limit surface, unlike other generalized spline
based subdivision schemes, does not have any closedform analytic
formulation. Moreover, we shall reformulate the dynamic Catmull
Clark subdivision surface model using this novel methodology, and
describe how to develop an unified dynamic framework for any sub
division scheme. The key contribution of this unified approach is to
represent the smooth limit surface of any subdivision scheme using
a single type of novel finite elements. The geometric and physi
cal features of our subdivisionbased finite elements depend only
on the subdivision scheme involved. Our FEMbased approach sig
nificantly advances the stateoftheart of physicsbased geometric
modeling because (1) it provides a universally physicsbased solu
tion to any subdivision schemes beyond prevalent splinelike sub
division techniques; (2) a natural mechanism that allows users to
intuitively deform any subdivision surface has been systematically
devised; (3) the limit surface of any subdivision schemes has been
represented using a single type of novel subdivisionbased finite
elements; and (4) our subdivisionbased finite elements are poten
tially of great interest to FEM communities.
1.3 Overview
The rest of the paper is organized as follows. A dynamic frame
work for the interpolatory (modified) butterfly subdivision scheme
is presented in Section 2. We reformulate the dynamic framework
for the approximating CatmullClark subdivision scheme using the
proposed approach in Section 3. Section 4 presents a solution on
how to develop a dynamic framework for any subdivision scheme.
Experimental applications are presented in Section 5. Finally, we
conclude the paper in Section 6.
2 Dynamic Butterfly Subdivision Sur
faces
This section discusses a dynamic framework for an interpolatory
subdivision scheme namely, the (modified) butterfly subdivision
technique. First, a brief overview of the (modified) butterfly sub
division scheme is presented. Next, a local geometric parameter
ization technique for the limit surface of the (modified) butterfly
subdivision is detailed. Our parameterization method is then used
to derive the new triangular finite element model for iii,,, i i .,.,. I
subdivision. Finally, the implementation details are described. Note
that, we will further generalize our physicsbased formulation for
other interpolatory subdivision schemes in Section 4.
2.1 Geometry of The (Modified) Butterfly Subdivi
sion
(a) (b)
Figure 1: (a) The control polygon with triangular faces. (b) The re
fined mesh obtained after one subdivision step using butterfly sub
division rules.
Figure 2: (a) The weighing factors of contributing vertex positions
for an edge connecting two vertices of degree 6; (b) the correspond
ing case when one vertex is of degree n and the other is of degree
6.
The butterfly subdivision scheme [8], like other subdivision tech
niques used in geometric modeling and graphics, starts with an ini
tial triangular mesh (a.k.a. the control mesh) defined by a set of
control vertices. 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 is split by adding a new vertex whose position
is obtained by an affine combination of the neighboring vertex po
sitions in the coarser level. For instance, the mesh in Fig.l(b) is
obtained by subdividing the initial mesh shown in Fig.l(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 posi
tion and degree of all original vertices do not change in the refined
mesh.
In the original butterfly scheme, the new vertices corresponding
to the edges in the previous level are obtained using an eightpoint
stencil. It produces a smooth C' surface in the limit except at the
extraordinary points corresponding to the extraordinary vertices
(vertices with degree not equal to 6) in the initial mesh [25]. Since
all the vertices introduced through subdivision have degree 6, the
number of extraordinary points in the smooth limit surface equals to
the number of extraordinary vertices in the initial mesh. Recently,
the original butterfly scheme has been modified by Zorin et al. [25]
to obtain better smoothness properties at the extraordinary points.
In this modified butterfly subdivision technique, all the edges had
been categorized into three classes: (i) edges connecting two ver
tices of degree 6 (a 10 point stencil, as shown in Fig.2(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 / 6
(the corresponding stencil to obtain new vertex position is shown in
Fig.2(b), where q = .75 is the weight associated with the vertex of
degree n $ 6, and si = (0.25+cos(27ri/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 $ 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.2(b) at each of those two
extraordinary vertices.
2.2 Formulation
In this section, we systematically formulate the dynamic framework
for the modified butterfly subdivision scheme. Unlike the approx
imating schemes, the geometry of the limit surface obtained via
modified butterfly subdivision does not have any closedform ana
lytic expression even for a regular mesh. Therefore, the key issue
is to define an appropriate parametric domain and derive a local pa
rameterization for iii,,, i. ... / subdivision. These relevant geo
metric components are critical to the development of our physics
based finite element model for the limit surface of butterfly scheme.
The smooth limit surface defined by the modified butterfly sub
division technique is of arbitrary topology where a global parame
terization is impossible. Nevertheless, the limit surface can be lo
cally parameterized over the geometric domain defined by the ini
tial mesh. The idea is to track an arbitrary point on the initial mesh
across the mesh hierarchy obtained via the subdivision process (see
Fig.3 and Fig.4), so that a correspondence can be established be
tween the point being tracked in the initial mesh and its image on
the limit surface.
Figure 3: The smoothing effect of the subdivision process on the
triangles of the initial mesh.
The modified butterfly subdivision scheme starts with an initial
mesh consisting of a set of triangular faces. The recursive appli
cation of the subdivision rules smoothes out each triangular face,
and in the limit, we obtain a smooth surface consists of a collec
tion of smooth triangular patches. The subdivision process and the
triangular decomposition of the limit surface is depicted in Fig.3.
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 s can be expressed as
S= ISk,
k=l
.X Xf
X
(a) (b)
Figure 4: 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.
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 kth triangular face in the initial mesh.
After the above geometric decomposition, we now describe the
parameterization of the limit surface over the initial mesh. The pro
cedure can be best explained through the following example. A
simple planar mesh shown in Fig.4(a) is chosen as the initial mesh.
An arbitrary point x inside the triangular face abe is tracked over
the meshes obtained through subdivision. The vertices in the initial
mesh are darkly shaded in Fig.4. 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 abe in the initial mesh
depends only on the vertices in the initial mesh which are within
the 2neighborhood of the vertices a, b and c due to the local na
ture of the subdivision process (the kneighborhood of a vertex in
cludes all the vertices that can be reached following at most k edges
from the given vertex). For example, the vertex d, introduced af
ter first subdivision step, can be obtained using the 10 point stencil
shown in Fig.2(a) on the edge ab. All the contributing vertices in
the initial mesh are within the 1neighborhood 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 contribut
ing vertices at this level of subdivision, for example, the (lightly
shaded) 1neighbors of the vertex b (except d and e) in Fig.4(b),
depend on some vertices in the initial mesh which are within the
2neighborhood 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,,,, 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 jth level
of subdivision. Let v ,b be the collection of vertices in the initial
mesh that are within the 2neighborhood of the vertices a, b and c
(marked black in Fig.4(a)). Let the number of such vertices be r.
Then, the vector v 0, which is the concatenation of the (x, y, z)
positions for all the r vertices, is of dimension 3r. Based on the
above observation of the 2neighborhood property, the geometry of
the smooth triangular patch in the limit surface corresponding to
the triangular face abe in the initial mesh is uniquely determined
by these r vertices. Because of the recursive characteristic, there
now exists four subdivision matrices (Aab c), (Aabc)j, (Aabc),
and (Aa,,c) of dimension (3r, 3r) such that
1 0
vadf = (Aabc),V)bc,
Vbd = (Aabc)lvabc,
Vcfe = (Aabc),vc.c,
Vdef = (Aabc),Vbc, (2)
where the subscripts t, 1, r and m denote top, left, right and mid
dle triangle positions, respectively (indicating the relative posi
tion of the new triangle with respect to the original triangle), and
vadf, V d, Vcfe and vd f are the concatenation of the (x, y, z) po
sitions for the vertices in the 2neighborhood of the corresponding
triangle within the newly obtained refined mesh after one subdi
vision. Note that, the new vertices in this level of subdivision are
lightly shaded in Fig.4(b). The 2neighborhood 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.2 are the same.
Carrying out one more level of subdivision, a new set of vertices
which are unshaded in Fig.4(c) are obtained along with the old ver
tices. Adopting a similar approach as in the derivation of Eqn.2, it
can be shown that
vdg = (Abed),Vbed
\ = (Abed) Vbed
vih = (Abed),vbed
2
Vghi = (Abed),Vbed (3)
The relative position and geometric structure of the triangu
lar face dgi in Fig.4(c) with respect to the triangular face bed is
topologically the same as of the triangular face adf in Fig.4(b)
with respect to the triangular face abc. Therefore, we can obtain
(Abed), = (Aabc)i. Based on the similar reasoning, Eqn.3 can be
rewritten as
2 1 1
"dgj = (Abed)1Vbed = (Aabc)1Vbed
\ = (Abed);Vbed = (Aabc);Vbed
Ve ih = (Abed),vbed = (Aabc),vbed
2 1
Vhi = (Abed),Vt) d = (Aabc),vCd. (4)
Combining Eqn.2 and Eqn.4, it can be shown that
vdgi = (Aabc)t(Aabc)jVabc,
S = (Aabc)j(Aabc)V bc,
2 0
Veih = (Aabc),r(Aabc)j Vabc)
Vghi = (Aabc),(Aabc)l Vbc (5)
Let x be a point with barycentric coordinates (0bc, bc'
inside the triangular face abc. When the initial mesh is refined, x
becomes a point inside the triangular face bed with barycentric co
ordinates (bed,. bd). Another level of subdivision causes x
to be included in the triangular face dgi with barycentric coordi
nates (a i, It) 2 ). Let sj denote the jth level approxima
nates (gidgi abc
tion of the smooth triangular patch Sabc in the limit surface corre
sponding to the triangular face abc in the initial mesh. Now v0bc
can be written as
r r r
vabc = [ax, bx, cx, ...,y, by,, c...,0a, b,, c. ...]
where the subscripts x, y and z indicate the x, y and z coordinates
of the corresponding vertex position,respectively. The expressions
for vbed and v di can also be written in a similar manner. Next, the
matrix Bae can be constructed as follows:
r r r
a ,bc, P bc, bc,0, . ., 0, ...,0, 0,...
r r
B c(X) o
0,.... 0 ,, abc ,bc bc,, 0, 0, O, ,. ., O
r r r
0,..., 0, 00, .., 0, a ,bc, abc, Yabc, , 0,..., O0
The matrices Bd and B ,i can also be constructed in a similar
fashion. Now sOb(x), Sb(x), and s b(x) can be written as
aabc (X) can be written as
S bc
Sbc (x)
Sabc (x)
Sabc (x)
B' I\'I d
B' i.A
= B1 \'IA
' (A bc) abc,
(Aabc)t, (Aabc)j, (Aabc), and (Aabc)m sums to one. The largest
eigenvalue of such matrices is 1 and therefore the mathematical
limit in Eqn.8 exists. Now, assuming the triangular face abc is the
kth face in the initial mesh, Eqn.8 can be rewritten as
Sk(x) = Bk((x)v = Bk(x)Akp, (9)
where p is the concatenation of the (x,y,z) positions of all the ver
tices in the initial mesh and the matrix Ak, when postmultiplied by
p, only selects the vertices vk defining the kth smooth triangular
patch in the limit surface. If there are t vertices in the initial mesh
and r of them control the kth 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.1 and Eqn.9, it can be shown that
s(x) = (Y B/k(x)Ak)p = J(x)p,
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.
We now treat the vertex positions in the initial mesh defining the
smooth limit surface s as time variables in order to develop the new
dynamic butterfly subdivision model. The velocity of the surface
model can be expressed as s(x, p) = J(x)p, where an overstruck
dot denotes a time derivative and x C So, So being the domain
defined by the initial mesh. Note that, So is the parametric domain
of the limit surface, each triangle of the initial control mesh serves
as a local parametric domain for its corresponding triangular patch.
2.3 Finite Element Procedure
Vabc,
' bed
Proceeding in a similar way, the expression for s,b(x), jth
level approximation of sab, (x), is given by
,() = (x)(A c)... (AbC), (A.bc); Vbc
S B (,,(x)(A, b) Vbc
= B (' \) ',c, (7)
where x is inside the triangular face uvw at level j (with an assump
tion that uvw is the triangular face in the middle with respect to its
coarserlevel original triangular face in the previous level), (Abc)
= (Abc)m ... (Abc),(Abc), and Bbc(x) = BUw(x)(Ab,).
It may be noted that the sequence of applying (Abc),, (Aabc)l,
(Aabc) and (Aab,)m depends on the triangle inside which the
tracked point x falls after each subdivision step. Finally, the local
geometric parameterization procedure can be completed by writing
Sabc(x) = (lim B (x))vbc = Babc(x)v bc. (8)
Note that, Bab, is the collection of basis functions at the ver
tices of vbc. It may also be noted that the modified butter
fly subdivision scheme is a stationary subdivision process, and
hence new vertex positions are obtained by affine combinations
of nearby vertices. This guarantees that each row of the matrices
(b)
Figure 5: (a) An initial mesh, and (b) the corresponding limit sur
face. 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 el
ement.
In Section 2.2 we have demonstrated that the smooth limit sur
face of butterfly subdivision can be represented by a collection of
smooth triangular patches. In our dynamic framework, we now con
sider each patch of the limit surface as a finite element. The number
of such patches is equal to the number of triangular faces in the ini
tial mesh as mentioned earlier. The concept of decomposing the
smooth limit surface into a collection of elements is illustrated in
Fig.5. We also show the parametric domain and control vertices for
shaded elements in Fig.5. The governing motion equation of this
subdivisionbased FEM model is given by
Mp + Dp + Kp = fp, (11)
where fp is the generalized force vector, and M, D, and K are the
mass, damping and stiffness matrices of the physical model. We
provide an outline on how to derive the mass, damping and stiffness
matrices for these finite elements so that a numerical solution to
the governing secondorder differential equation can be obtained
using popular finite element analysis techniques. We use the same
example as in Section 2.2 (refer to Fig.4) to introduce the relevant
concepts and derive our FEM model.
The mass matrix for the element Sab, corresponding to the tri
angular face abc, can be expressed as
fESobc
I(x)B,,C (x)BB,,G (x)dx.
However, the basis functions (stored as entries in Babc) do not have
any analytic form, hence computing this integral is a difficult propo
sition. We solve this problem by approximating the smooth trian
gular patch in the limit surface corresponding to the face abc in
the initial mesh by a triangular mesh with 42 faces obtained after j
levels of subdivision of the original triangular face abc (each sub
division step splits one triangular face into 4 triangular faces). In
addition, we choose a discretized form of mass distribution function
which has nonzero values only at the vertex positions of the jth
subdivision level mesh. Then the mass matrix can be approximated
k
_l =a Iza(IB( )} {Bc( )},
where k is the number of vertices in the triangular mesh with 43
faces. This approximation has been found to be very effective and
efficient for the implementation of FEM procedure. The computa
tion of elemental damping matrix follows suit.
Physicsbased models have both kinetic and potential energies.
We now define the internal (e.g., elastic) energy of the subdivision
based dynamic model by assigning deformation energy to each ele
ment. We take a similar approach as shown above and consider the
jth level approximation of the element. Note that, a wide range of
functional formulations can be employed to describe various ma
terial and physical behaviors such as linear elastic deformation and
nonlinear plastic deformation. Throughout this paper, in particular,
we assign springlike energy to the approximated model because of
its simplicity and efficient computation. The energy at the jth level
of approximation can be defined as
Eabc z Eabc
1 k;m i(vi ,,l )
2 [2 vF (v1
1 { }T (K*){vy}.
where ki, is the springcontrolling variable, v and v,, the 1th
and mth vertex in the jth level mesh, are in the 1neighborhood
of each other, Q is the domain defined by all such vertex pairs, fl,
is the natural length of the spring connected between v, and v,,
and v, is the concatenation of the (x,y,z) positions of all the ver
tices in the jth subdivision level of the triangular face abc in the
(a) (b)
Figure 6: CatmullClark subdivision: (a) initial mesh, (b) mesh
obtained after one step of CatmullClark subdivision, and (c) mesh
obtained after another subdivision step.
initial mesh. The vertex positions in v_, are obtained by a lin
ear combination of the vertex positions in v0bC, and hence we can
write v3,b = (A3,C)vC where (Aa,.c) is the transformation (sub
division) matrix. Therefore, the expression for the elemental stiff
ness matrix is given by Kabc = (Ajb,,T) (Kjac)(AbC). It may
be noted that this approach is applicable for modeling isotropic as
well as anisotropic phenomena because kl,, the springcontrolling
variable, can be a timedependent function in general.
3 Dynamic CatmullClark Subdivision
Surfaces
In this section, we consider a new FEM model based on an ap
proximating subdivision scheme, namely, CatmullClark subdivi
sion technique. Please note, the dynamic formulation of Catmull
Clark subdivision previously proposed in [14, 19] could not be gen
eralized for other approximating subdivision schemes. The frame
work developed in this section can be generalized to other approx
imating subdivision schemes as shown in Section 4. In fact, a dy
namic framework for Loop's technique (another popular approxi
mating subdivision scheme) has been presented in [13] using the
algorithm proposed here. We limit our discussion to CatmullClark
subdivision surfaces only in this paper due to space restrictions.
We first outline the CatmullClark subdivision scheme. Next, we
present the dynamic formulation. In particular, we address the dif
ference between the current work and prior results [14, 19]. Finally,
we discuss the finite element implementation of the physicsbased
model.
3.1 CatmullClark Subdivision Scheme
CatmullClark subdivision scheme, like any other subdivision
scheme, starts with an userdefined 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 sequence of recursively refined polyhedral meshes
will converge to a smooth surface. The subdivision rules are as fol
lows:
(1) For each face, a new face point is introduced which is the aver
age of all the old vertices defining the face.
(2) For each (nonboundary) edge, a new edge point is introduced
which is the average of the following four points: two old vertices
defining the edge and two new face points of the faces adjacent to
the edge.
(3) For each (nonboundary) vertex V, new vertex is introduced
whose position is + l+ + 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 edges incident on the old vertex V
and n is the number of the edges incident on the vertex.
(4) New edges are formed by connecting each new face point to the
new edge points of the edges defining the old face and by connect
ing each new vertex point to the new edge points of all old edges
incident on the old vertex point.
(5) New faces are defined as faces enclosed by new edges.
An example of CatmullClark subdivision on an initial mesh is
shown in Fig.6. The most important property of the CatmullClark
subdivision surfaces is that a smooth surface can be generated from
any control mesh of arbitrary topology. CatmullClark subdivision
surfaces include standard bicubic Bspline surfaces as their special
case (i.e., the limit surface is a bicubic Bspline surface for a rect
angular mesh with all nonboundary vertices of degree 4). In addi
tion, the aforementioned subdivision rules generalize the recursive
bicubic Bspline patch subdivision algorithm. For nonrectangular
meshes, the limit surface converges to a bicubic Bspline surface
except at a finite number of extraordinary points. These extraor
dinary 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 new vertices created
subsequently will have four incident edges. The number of extraor
dinary points on the limit surface is a constant, and is equal to the
number of extraordinary vertices in the refined mesh obtained after
applying one step of the CatmullClark subdivision on the initial
mesh. The limit surface is curvaturecontinuous everywhere ex
cept at extraordinary vertices, where only tangent plane continuity
is achieved.
3.2 Formulation
A systematic formulation of the newly proposed dynamic frame
work for CatmullClark subdivision surfaces is presented in this
section. The key difference between the dynamic model devel
oped in [14, 19] and the one presented here is the representation
of the limit surface. The previously proposed approach leads to di
verse types of finite elements, whereas the present approach leads
to single type of finite elements. This is illustrated with a schematic
diagram in Fig.7.
Following the concepts developed in [14, 19], the limit surface
of the control mesh shown in Fig.7, consists of quadrilateral bicubic
Bspline patches corresponding to the faces marked 'n' (faces with
no extraordinary points), and a pentagonal patch corresponding to
the faces marked 's' (faces having one extraordinary vertex of de
gree 5) (Fig.7(a)). However, in this section, it has been shown that
the entire limit surface can be expressed as a collection of quadri
lateral patches as shown in Fig.7(b) using the algorithm proposed in
[23]. We next discuss a local parameterization of the limit surface
which is critical to embed the limit surface in a dynamic framework.
As mentioned earlier, the control mesh (after at most one subdi
vision step) for the CatmullClark subdivision scheme consists of
quadrilateral faces which lead to quadrilateral patches in the limit
surface. For the sake of formulation simplicity, it has been assumed
that each face has at most one extraordinary vertex. If this assump
tion is not valid, then one more subdivision step needs to be per
formed on the current control mesh in order to obtain a new control
mesh on which the following analysis can be carried out. The num
ber of quadrilateral patches in the limit surface is equal to the num
ber of nonboundary quadrilateral faces in the control mesh (Fig.8).
Therefore, the smooth limit surface s can be expressed as
s= s1, (14)
=
where n is the number of nonboundary faces in the control mesh
(a) (b)
Figure 7: A control mesh with an extraordinary vertex of degree 5
and the corresponding limit surface: (a) using the concepts devel
oped in [14, 19], where the limit surface consists of quadrilateral
normal elements and a pentagonal special element; (b) using the
unified approach developed in this paper, where the limit surface
consists of one single type of quadrilateral finite element.
and sl is the smooth quadrilateral patch corresponding to the lth
nonboundary quadrilateral face in the control mesh. Each of these
quadrilateral patches can be parameterized over the correspond
ing nonboundary quadrilateral face in the control mesh. How
ever, since a quadrilateral face can easily be reparameterized over
a [0, 1]2 domain, each quadrilateral patch is locally parameterized
over [0, 1]2.
The nonboundary quadrilateral faces are of two types : (a)
faces having no extraordinary vertices (dubbed as "regular" faces
in [14, 19], marked as K in Fig.8(a)) and (b) faces with one extraor
dinary vertex (dubbed as "irregular" faces in [14, 19], marked as (
in Fig.8(a)). If there are m regular and n m irregular faces, then
Eqn. 14 can be rewritten as
s= s + s sj, (15)
i=1 j=1
where si is the quadrilateral patch corresponding to the ith regu
lar face and sj is the quadrilateral patch corresponding to the jth
irregular face.
The quadrilateral patch in the limit surface corresponding to
each regular face is a bicubic Bspline patch, which is defined over
[0, 1]2. The set of control vertices defining this bicubic Bspline
patch can be obtained using the adjacent face information. There
fore, the quadrilateral patches in the smooth limit surface corre
sponding to the regular faces in the control mesh can be easily ex
pressed analytically, which are essentially bicubic Bspline patches
defined by 16 control vertices over a [0, 1]2 domain. The analytic
expression for the quadrilateral patch corresponding to the regular
face i is given by
si = Jb(u,v)pi
S(Jb(u, v)Ai)p
= Ji(u,v)p, (16)
(b)
Figure 8: In CatmullClark subdivision, each nonboundary quadri
lateral face in the control mesh has a corresponding quadrilateral
patch in the limit surface : (a) control mesh, (b) limit surface.
where 0 < u, v < 1, Jb(u, v) is the collection of the bicubic B
spline basis functions, pi is the concatenation of the 16 control
vertex positions defining the bicubic Bspline patch, Ai is a se
lection matrix which when multiplied with p, the concatenation of
all the control vertex positions defining the smooth limit surface,
selects the corresponding set of control vertices, and Ji(u, v)
Jb(u, v))Ai.
By contrast, the analytic expression of the quadrilateral patches
corresponding to the irregular faces in the control mesh was difficult
to derive, and hence an alternative approach was taken in [14, 19].
However, very recently an efficient scheme for evaluating Catmull
Clark subdivision surfaces at arbitrary parameter values has been
proposed by Stam [23]. The proposed approach, involving eigen
analysis of the subdivision matrix, leads to an analytic expression of
the quadrilateral patches which are parameterized over an irregular
face in the control mesh, and hence over [0, 1]2 after reparameteri
zation. Following the scheme developed by Stam [23], the quadri
lateral patch corresponding to the irregular face j is given by
Si = Jd, (u, v)pj
= (Jd,(u,v)Aj)p
= Jj(u, )p, (17)
where 0 < u, v < 1 as before. Jd, (u, v) is the collection of basis
functions for the corresponding quadrilateral patch in the smooth
limit surface. The subscript dk is used to denote the fact that the
irregular face has an extraordinary vertex of degree k. The de
tailed derivation and the analytic expressions of these basis func
tions involving the eigenvalues and eigenvectors of the subdivision
matrix can be found in [23]. The other symbols used in Eqn.17
have the usual meaning: pj is the concatenation of the 2k + 8 con
trol vertices defining the quadrilateral patch in the limit surface, p
is the concatenation of all the control vertex positions defining the
smooth limit surface, Aj is a selection matrix which when multi
plied with p selects the corresponding set of control vertices, and
Jy(u, v) =Jd (u, v)Aj.
(a) (b)
Figure 9: (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.
It may be noted that the number of control vertices in the ini
tial mesh defining a quadrilateral patch in the smooth limit surface
is 2k + 8, where k = 4 in case the associated quadrilateral face
in the control mesh is regular, or k = degree of the extraordinary
vertex if the associated quadrilateral face is irregular. For example,
the shaded quadrilateral patch is associated with the shaded regu
lar face in Fig.2(a), and the 16 control vertices defining this patch
(which is actually a bicubic Bspline patch) are marked. Similarly,
the shaded quadrilateral patch is associated with the shaded irregu
lar face in Fig.2(b), and the 14 control vertices defining this patch
are highlighted. Now an expression of the smooth limit surface can
be formulated. Using Eqn.15, 16 and 17, it can be shown that
s Jp + JjP
i=1 j=1
= (ZJ + j)p
i=1 j=l
= Jp, (18)
where J = ( J + =L Jj). Note that even though the
initial mesh serves as the parametric domain of the smooth limit
surface, each quadrilateral face in the initial mesh and consequently
the smooth limit surface can be defined over a [0, 1]2 domain.
Once an analytic expression of the smooth limit surface of
CatmullClark subdivision is derived, we then develop the dynamic
model by considering the control vertex positions as timevarying
variables. The velocity of the surface model can be expressed as
s(x, u, v) = J(x, ', i where an overstruck dot denotes a time
derivative and x C S", ." being the domain defined by the initial
mesh.
3.3 Finite Element Implementation
The smooth limit surface of CatmullClark subdivision comprises
a collection of quadrilateral patches. Each quadrilateral patch is
considered as a finite element. Therefore, within the unified frame
work the limit surface can be decomposed into one single type of
finite elements rather than two different types as in [14, 19]. Our
new FEM technique significantly simplifies the data structure and
system architecture. Consequently, more efficient algorithms for
finiteelement assembly, dynamic simulation, etc. can be devised
using this unified approach. The motion equation of the dynamic
model is same as that of the dynamic model of butterflybased sub
division:
Mp + Dp + Kp = fp, (19)
where fp is the generalized force vector and M, D, and K are the
mass, damping and stiffness matrices of the model. The expressions
of the mass, damping and stiffness matrices for a quadrilateral ele
ment (which is a bicubic Bspline) can be written as
x1 = ~pJ Jbdudv, (20)
De = 7JBJbdudv, (21)
0 0
and
Ke, = / (a{(Jb) T{ (Jb))} + 22b{(Jb)}T{(J b)}
+ 311{(Jb)uI (Jb)u.} + {b(Jb)uu. T(Jb)uu}
+3)22{(Jb,)UU} {(Jb),u})dudv (22)
respectively, where Jb is the bicubic Bspline basis matrix, P(u, v)
is the mass density, y(u, v) is the damping density, aii(u, v) and
3ij (u, v) are the tension and rigidity functions respectively. The
subscript u and v denote partial derivatives with respect to u and
v respectively. The subscript e is used to indicate elemental matri
ces which are of size (16, 16). Note that, the mass, damping and
stiffness matrices for these elements can be evaluated analytically,
provided the material properties (e.g., mass, damping, rigidity and
bending distributions) have analytic expressions. In some cases,
these distribution functions can be assumed to be constant to sim
plify the matter.
The mass, damping and stiffness matrices for the quadrilateral
elements which are not bicubic Bsplines (corresponding to the ir
regular faces) can also be expressed analytically by simply replac
ing the matrix Jb in Eqn.20, 21 and 22 with the matrix Jd, (refer
to Eqn.17), where k denotes the degree of the extraordinary vertex
associated with the corresponding irregular face. These elemental
matrices are of size (2k + 8, 2k + 8). The generalized force vector
for these elements can also be determined in a similar fashion. It
may be noted that the limits of integration need to be chosen care
fully for elemental stiffness matrices as the second derivative di
verges near the extraordinary points for CatmullClark subdivision
surfaces.
Even though an analytical expression for a nonBspline quadri
lateral element in the limit surface exists, it is cumbersome to ac
tually evaluate the elemental matrix expressions. Numerical inte
gration using Gaussian quadrature may be used to obtain approx
imations of these elemental matrices. However, in this paper, an
approach similar to the FEM procedure presented in Section 2 is
utilized because of its simplicity and effectiveness. An approxima
tion of the smooth limit surface is obtained by refining the initial
control mesh j times, and a springmass system is developed on
this jth approximation level in a similar fashion as in Section 2.3.
The physical matrices of this system is then used as an approxima
tion to the actual physical matrices. This approximation has been
found to be very efficient for implementation purposes.
4 Unified Approach For Any Subdivision
Scheme
The dynamic framework for modified butterfly and CatmullClark
subdivision scheme can be generalized to any subdivision scheme.
The key observation is that the smooth limit surface can be viewed
as a collection of a single type finite element. Because of the
nature of recursive refinement, any subdivisionbased scheme es
sentially defines a "natural" correspondence which leads to a lo
cal parameterization of the smooth limit surface. The unique type
of the associated finite element results from the local parameteri
zation scheme. This is evident from the triangular finite element
patches developed for the modified butterfly subdivision scheme
and from the quadrilateral finite element patches developed for
CatmullClark subdivision scheme. We shall present a general out
line on how to provide a dynamic framework for interpolatory and
approximating subdivision schemes.
4.1 Interpolatory subdivision schemes
Most of the interpolatory subdivision schemes are obtained by mod
ifying the butterfly subdivision scheme [8]. Therefore, the frame
work for the modified butterfly subdivision scheme in Section 2
and its principles can be applied to other interpolatory subdivision
schemes. The only difference is that the basis functions as well
as the set of control vertices of arbitrary patch in the limit surface
depend on the chosen interpolatory subdivision rules. It may also
be noted that unlike the approximating schemes, the physical ma
trices can not be obtained analytically as the basis functions cor
responding to interpolatory subdivision schemes do not have any
analytic expressions in general. Even though these matrices can
be obtained via numerical integration, the pointmass system con
nected by springs as developed in Section 2 is more preferable for
implementation purposes because of efficiency reasons.
4.2 Approximating subdivision schemes
The unified approach for a dynamic model of CatmullClark sub
division can be generalized for other approximating subdivision
schemes as well. This generalized approach involves three steps:
(a) The limit surface obtained via an approximating subdivision
scheme can be expressed as a collection of smooth patches which
can be locally parameterized over a corresponding face in the con
trol mesh. Each patch is nsided if it is locally parameterized over a
nsided face. Analytic expressions for each of these patches can be
derived even in the presence of extraordinary vertices in the control
mesh, and hence an expression of the limit surface can be obtained.
(b) Once an expression of the limit surface is obtained, the dynamic
framework can be developed by considering control vertex posi
tions as a function of time. The corresponding motion equation can
be derived.
(c) Each patch in the limit surface is treated as a finite element in im
plementation. The elemental mass, damping and stiffness matrices
along with the generalized force vector can be obtained by either
analytic or numerical integration. Alternatively, the control mesh
can be subdivided j times to obtain an approximation of the smooth
limit surface, and a springmass system can be developed on this
approximation mesh. The physical matrices of this system provide
an approximation to the original physical matrices and works well
in practice.
5 Solid Modeling Applications
The proposed FEMbased dynamic subdivision models can be used
to represent a wide variety of smooth shapes with arbitrary genus.
The smooth limit object can be sculpted by applying synthesized
forces in a direct and intuitive way in shape design applications
for solid modeling. The underlying shape from a cloud of 3D
points can also be recovered hierarchically using our FEM mod
els. For data fitting applications, springs are attached to the initial
ized model from the data points in 3D, and the initialized model
evolves dynamically according to the equation of motion subject to
the applied spring forces and various geometric constraints. When
an optimal fit to the given data set is achieved, the number of con
trol vertices can be increased by replacing the original initial mesh
by a new initial mesh obtained by applying a single subdivision
step. This increases the number of degrees of freedom to repre
sent the same limit surface and a new equilibrium position for the
model with a much better fit to the given data set can be achieved.
The fittingerror criteria for the discrete data can be computed ac
cording to distance between the data points and the points on the
limit surface where the corresponding springs are attached. We now
demonstrate modeling and data fitting examples using our dynamic
FEM model.
In a shape modeling application, the user can specify any mesh
as the initial (control) mesh, and the corresponding limit surface
can be sculpted directly and interactively by applying synthesized
forces in realtime. We show several initial surfaces obtained from
different control meshes and the corresponding deformed surfaces
after interactive sculpting on the limit surface in Fig.10 (see last
page). To change the shape of an initial surface, the user can attach
springs from different points in 3D to the nearest point on the limit
surface such that the limit surface deforms towards these locations
to generate the desired shape. The limit surface here consists of
a single type of smooth triangular finite element patches, irrespec
tive of the number of extraordinary vertices in the control mesh.
The initial mesh of the smooth surface shown in Fig.10(a) has 125
faces and 76 vertices (degrees of freedom), which is deformed to
the smooth shape shown in Fig.10(b) by interactive spring force ap
plication. The initial mesh of the closed solid shape in Fig.10(c) has
24 faces and 14 vertices. This solid shape is deformed to the shape
shown in Fig.10(d). The one hole torus in Fig.10(e) and the corre
sponding modified shape in Fig.10(f) have initial meshes with 64
faces and 32 vertices. A two hole torus with a control mesh of 272
faces and 134 vertices, shown in Fig.10(g), is dynamically sculpted
to the shape shown in Fig.10(h).
We have also performed several experiments testing the applica
bility of our model to recover the underlying shapes from a cloud of
points in 3D. In all the experiments, the initialized dynamic model
has a control mesh comprising of 24 triangular faces and 14 ver
tices whereas the control mesh of the fitted model has 384 triangu
lar faces and 194 vertices. It may be noted that once an optimal
shape defined by a fixed number of control vertices (determined by
subdivision levels) is recovered, the limit smooth model is capable
of refining itself in accordance with the datafitting criteria, thereby
increasing the degrees of freedom of the recovered shape only when
necessary. For the fittingerror (defined as the maximum distance
between a data point and the nearest point on the limit surface ex
pressed as a percentage of the diameter of the smallest sphere en
closing the object) of approximately 3%, the initialized model is
refined twice. The datafitting examples are shown in Fig.11 (see
last page). In the first data fitting experiment, range data acquired
from multiple views of a light bulb is used and the model was ini
tialized inside the 1000 data points (Fig.ll(a)). The fitted dynamic
model is shown in Fig.1 l(b). In the next experiment, the shape of a
mechanical part is recovered from a range dataset containing 2031
data points (Fig.11(c) and (d)). We also recover the shape of a hu
man head from the data set as shown in Fig.1 l(e). The head data set
has 1779 3D points. It may be noted that the final shape with a very
low error tolerance is recovered using very few number of control
points in comparison to the large number of data points present in
the original range data set.
6 Conclusions
In this paper, we have presented a new FEMbased dynamic frame
work where a single type of subdivisionbased finite elements are
used to represent the smooth limit surface generated by any sub
division scheme. The primary objective is to integrate physics
based modeling techniques with geometric subdivision methodol
ogy for the interactive sculpting and direct manipulation of the
limit surface of prevalent subdivision schemes. We have proposed
an unified approach and demonstrated how to transform any sub
division scheme into our dynamic modeling framework. Model
ers can physically sculpt virtual objects defined through arbitrary
procedurebased subdivision techniques in a natural and intuitive
manner within the proposed framework. Users can also directly
enforce various functional and aesthetic requirements in the limit
surface without the need to explicitly handle control vertices. Fur
thermore, this dynamic framework permits physicsbased models
to be refined adaptively in a hierarchical fashion which is an intrin
sic feature of subdivision geometry. Our experiments have demon
strated the applicability of the new unified FEMbased framework
in solid modeling and data fitting applications. This unified method
will offer a greater potential for popular subdivision techniques in
solid and geometric modeling, interactive graphics, finite element
analysis, and engineering design applications.
7 Acknowledgments
This research was supported in part by the NSF grant ECS9210648
and the NIH grant RO1LM05944 to B.C. Vemuri; the NSF CA
REER award CCR9702103, the NSF grant DMI9896170, and a
research grant from Ford Motor Company to H. Qin. We wish to
acknowledge Dr. Hughes Hoppe and Dr. Kari Pulli for the data
sets.
References
[1] A.A. Ball and D.J.T. Storry, "Conditions for tangent plane
continuity over recursively generated Bspline surfaces,"
ACM Transactions on Graphics, vol. 7, no. 2, pp. 83 102,
1988.
[2] A.A. Ball and D.J.T. Storry, "An investigation of curvature
variations over recursively generated Bspline surfaces," ACM
Transactions on Graphics, vol. 9, no. 4, pp. 424 437, 1990.
[3] E. Catmull and J. Clark, "Recursively generated Bspline sur
faces on arbitrary topological meshes," Computer AidedDe
sign, vol. 10, no. 6, pp. 350 355, 1978.
[4] G.M. Chaikin, "An algorithm for high speed curve genera
tion," Computer Vision, Graphics and Image Processing, vol.
3, no. 4, pp. 346 349, 1974.
[5] T. DeRose, M. Kass and T. Truong, "Subdivision surfaces
in character animation," in Computer Graphics Proceedings,
ACM SIGGRAPH, Annual Conference Series, pp. 85 94,
July, 1998.
[6] D. Doo, "A subdivision algorithm for smoothing down ir
regularly shaped polyhedrons," in Proceedings on Interactive
Techniques in Computer AidedDesign, pp. 157 165, 1978.
[7] D. Doo and M. Sabin, "Analysis of the behavior of recursive
division surfaces near extraordinary points," Computer Aided
Design, vol. 10, no. 6, pp. 356 360, 1978.
[8] N. Dyn, D. Levin and J.A. Gregory, "A butterfly subdivision
scheme for surface interpolation with tension control," ACM
Transactions on Graphics, vol. 9, no. 2, pp. 160 169, April,
1990.
[9] A. Habib and J. Warren, "Edge and vertex insertion for a
class of C1 subdivision surfaces," Computer Aided Geometric
Design, to appear.
[10] M. Halstead, M. Kass and T. DeRose, "Efficient, fair interpo
lation using CatmullClark surfaces," in Computer Graphics
Proceedings, ACM SIGGRAPH, Annual Conference Series,
pp. 35 44, August, 1993.
[11] H. Hoppe, T. DeRose, T. Duchamp, M. Halstead, H. Jin,
J. McDonald, J. Schweitzer and W. Stuetzle, "Piecewise
smooth surface reconstruction," in Computer Graphics Pro
ceedings, ACM SIGGRAPH, Annual Conference Series, pp.
295 302, July, 1994.
[12] C. Loop, Smooth subdivision surfaces based on triangles,
M.S. thesis, University of Utah, Department of Mathematics,
1987.
[13] C. Mandal, A dynamic framework for subdivision surfaces,
Ph.D. thesis, University of Florida, Gainesville, 1998.
[14] C. Mandal, H. Qin and B.C. Vemuri, "Dynamic smooth sub
division surfaces for data visualization," in IEEE Visualiza
tion '97 Conference Proceedings, Phoenix,AZ, pp. 371 377,
October, 1997.
[15] C. Mandal, H. Qin and B.C. Vemuri. "Direct manipula
tion of butterfly subdivision surfaces : A physicsbased ap
proach," Technical Report CISETR98009, University of
Florida, 1998.
[16] J. Peters and U. Reif, "The simplest subdivision scheme for
smoothing polyhedra," ACM Transactions on Graphics, vol.
16, no. 4, pp. 420 431, October, 1997.
[17] J. Peters and U. Reif, "Analysis of generalized
Bspline subdivision algorithms," SIAM Journal
of Numerical Analysis, to appear, available at
ftp://ftp.cs.purdue.edu/pub/jorg/9697agb.ps.Z.
[18] H. Qin and D. Terzopoulos, "DNURBS : A physicsbased
framework for geometric design," IEEE Transactions on Vi
sualization and Computer Graphics, vol. 2, no. 1, pp. 85 96,
January March, 1996.
[19] H. Qin, C. Mandal and B.C. Vemuri, "Dynamic Catmull
Clark subdivision surfaces," IEEE Transactions on Visual
ization and Computer Graphics, vol. 4, no. 3, pp. 215 229,
July September, 1998.
[20] U. Reif, "A unified approach to subdivision algorithms near
extraordinary points," Computer Aided Geometric Design,
vol. 12 no. 2, pp. 153 174, 1995.
[21] J.E. Schweitzer, Analysis andApplication of Subdivision Sur
faces, Ph.D. thesis, University of Washington, Seattle, 1996.
[22] T.W. Sederberg, J. Zheng, D. Sewell and M. Sabin, "Non
uniform recursive subdivision surfaces," in Computer Graph
ics Proceedings, ACM SIGGRAPH, Annual Conference Se
ries, pp. 387 394, July, 1998.
[23] J. Stam, "Exact evaluation of CatmullClark subdivision sur
faces at arbitrary parameter values," in Computer Graphics
Proceedings, ACM SIGGRAPH, Annual Conference Series,
pp. 395 404, July, 1998.
[24] D. Terzopoulos, J. Platt, A. Barr and K. Fleischer, "Elasti
cally deformable models," in Computer Graphics Proceed
ings, ACM SIGGRAPH, Annual Conference Series, pp. 205
214, 1987.
[25] D. Zorin, P. Schr6der and W. Sweldens, "Interpolating sub
division for meshes with arbitrary topology," in Computer
Graphics Proceedings, ACM SIGGRAPH, Annual Confer
ence Series, pp. 189 192, August, 1996.
[26] D. Zorin, "Smoothness of stationary subdivision on irregular
meshes," Constructive Approximation, submitted, available as
Stanford Computer Science Lab Tech. Rep. CSLTR98752,
1998.
(a) (c) (e) (g)
D) (0)) (n)
Figure 10: (a), (c), (e) and (g) : Initial shapes; (b), (d), (f) and (h) : the corresponding modified shapes after interactive sculpting via force
application.
(c) (e)
(U) (I)
Figure 11: (a), (c) and (e) : Collection of points in 3D along with the initialized model; (b), (d) and (f) : the corresponding fitted dynamic
subdivision surface model.
