Direct Manipulation of Butterfly Subdivision
Surfaces: A Physicsbased Approach
Chhandomay Mandal*, Hong Qint and Baba C. Vemuri*
*Authors' address: Department of Computer and Information Science and Engineering, University of Florida,
Gainesville, Florida, 32611. Email:{cmandal  vemuri}lcise.ufl.edu.
tAuthor's address : Department of Computer Science, I.. University of .'". York at '..I Brook, ',..I
Brook, .' York, 11794. i. .il 1 e!.  Ii. lI, I
August 20, 1998
DRAFT
Abstract
Subdivision surfaces are extensively used to model smooth shapes of arbitrary l I 1 !_ . Re
cursive subdivision on an userdefined 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 manipulation tools for the limit surface. In this paper, we integrate techniques
from ,1!i i, based modeling with geometric subdivision i, 1 .1.i.1, ._, and present a scheme for
direct manipulation of the smooth limit surface generated by the (modified) butterfly scheme
using 1i1! .i based '!,.!,, tools. This procedurebased surface model obtained through but
terfly subdivision does not have a closedform ,i!i ,1 I formulation (unlike other well known
splinebased models), and hence poses challenging problems to incorporate mass and damping
distribution functions, internal deformation ii. i forces, and other i1!. . i1 quantities required
to develop a 1'1! i. based model. Our primary contributions to computer graphics and geomet
ric modeling include : (1) a new hierarchical formulation for locally parameterizing the butterfly
subdivision surface over its initial control polyhedron, (2) formulation of 1,1! i. based butterfly
subdivision surface as a set of novel finite elements, and (3) optimal approximation of this new
I "' of finite elements by a collection of existing finite elements ,. I 1, ,i 1 to implicit geometric
constraints. Our new dynamic model can be sculpted directly by ;,1l', i~l~ synthesized forces,
and its equilibrium is characterized by a minimum of the deformation ii. i L .ii. I to the im
posed constraints. We demonstrate that this novel dynamic framework not only provides a direct
and natural means of manipulating geometric shapes, but also facilitates hierarchical shape and
nonrigid motion estimation from large range and volumetric data sets using very few degrees
of freedom (control vertices that define the initial polyhedron). This new ,1!i ;. based model
promises a greater potential of subdivision schemes in computer graphics, geometric modeling,
and virtual environments.
August 20, 1998
DRAFT
Keywords
Computer Graphics, P! ; based Modeling, Geometric Modeling, CAGD, Scientific Visu
alization, Subdivision Surfaces, Deformable Models, I !!!!i~ 1. I ii I Interactive Techniques.
I. INTRODUCTION
The concept of subdivision is ubiquitous in computer science. In visual computing areas, subdivision
surfaces are extensively used to model smooth shapes of arbitrary l. ,',.1 _ for computer graphics, an
imation, and geometric design applications. A I 1i.; 1 recursive subdivision 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. Despite the presence of a ! i;. of subdivision schemes in the computer
graphics and geometric modeling literature, there is no direct and natural way of manipulating the limit
surface. The current stateoftheart only permits the modeler to interactively obtain the desired effects
on the smooth surface by kinematically manipulating the vertex positions at various levels of subdivision
hierarchy. In this paper, we tackle the challenging problem of direct manipulation of the limit subdivision
surface at arbitrary locations/areas and offer a novel solution to this problem by embedding the modified
butterfly subdivision scheme in a 1!i ;1 based modeling framework. As a result, unlike the existing
geometric solutions that only allow the operations on control vertices, our i. !. ,1.l .1._ and algorithms
permit the user to '11! "; 1! modify the shape of subdivision surfaces at desired locations via application
of forces. This gives the user an intuitive and natural feeling that is produced while modeling with real
clay/' I '1.. il1! We will also demonstrate that the proposed model t I !. II, recovers shapes as well as
nonrigid motions from large range and volumetric data sets. Note that, this paper neither proposes a
new subdivision technique nor provides a different interpretation of ;!!, existing subdivision technique,
but integrates the advantages of subdivision surfacebased and 1i1! ;i based modeling techniques to solve
important theoretical and practical problems. Although the principles of 1'1! i based modeling are well
understood by the graphics experts and modeling researchers, this paper will greatly advance the state
of the art in i1. i; based shape modeling due to the following contributions : (1) a hierarchical local
August 20, 1998
DRAFT
parameterization scheme for (modified) butterfly subdivision surfaces is derived; (2) material properties
are assigned to the smooth limit surface in order to embed this geometric model in a '11! i. based mod
eling framework; (3) the smooth limit surface is decomposed into a set of novel finite elements; (4) the
motion equations that govern the behavior of (modified) butterfly subdivision surfaces are derived; (5)
users can sculpt this 1i'1! i. based model at desired locations and the model responds naturally (as a
realworld object would); and (6) algorithms and procedures are developed to best approximate this novel
finite element using a collection of existing lower order finite elements. WVe will address the ii.. 1 of
the proposed model, detail our contributions, and compare this new model with our previous research
results on dynamic subdivision surfaces in later sections. in I we will ',i f[ review the previous work
on subdivision surfaces.
A. Background
In [1], ('!i ii!:!n first introduced the concept of subdivision to the computer graphics .!!' .iil,! for
generating a smooth curve from a given control i..1 _.i During the last two decades, a wide i I, 
of subdivision schemes for modeling smooth surfaces of arbitrary 1. ,..1. _ have been derived following
('I! ii,!:, pioneering work on curve generation. In general, these subdivision schemes can be catego
rized into two distinct classes namely, (1) approximating subdivision techniques and (2) interpolating
subdivision techniques.
Among the approximating schemes, the techniques of Doo and Sabin [2], [3] and Catmull and ('! !:
[4] generalize the idea of obtaining uniform biquadratic and bicubic Bspline patches respectively from
a rectangular control mesh. In [4], Catmull and ('I i !: developed a method for recursively generating a
smooth surface from a polyhedral mesh of arbitrary Il. .1. . . The Catmull('I! I : subdivision 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 [5], Loop presented a similar subdivision scheme based on the
generalization of quartic triangular Bsplines for triangular meshes. Hoppe et al. [6] extended his work
to produce piecewise smooth surfaces with selected discontinuities. Halstead et al. [7] proposed an
algorithm to construct a Catmull( '! !: subdivision surface that interpolates the vertices of a mesh of
August 20, 1998
DRAFT
arbitrary I. 1I,. .1 Peters and Reif [8] proposed a simple subdivision scheme for smoothing i" .1 1!. 1
Most recently, nonuniform DooSabin and Catmull('!I !: surfaces that generalize nonuniform tensor
product Bspline surfaces to arbitrary topologies were introduced by Sederberg et al. [9]. 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. [10].
The most wellknown interpolationbased subdivision scheme is the "1.111i. i!l algorithm proposed by
Dyn et al. [11]. Butterfly subdivision method, like other subdivision 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 C1 limit surface. A variant of this scheme with better smoothness properties can
be found in [12]. Zorin et al. [13] has developed an improved interpolatory subdivision scheme (which
we call the !i... ..r' 1 butterfly scheme) that retains the ii.wli, '1 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 [14], [15] and by Kobbelt and Schr6der [16]. 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 sub
division algorithms is rather complex. Doo and Sabin [17] first i!i! 1 ,1 the smoothness behavior of
the limit surface using the Fourier transform and an ._. i_ 1, ,1 i of the subdivision matrix. Ball and
1i..i, [18], [19] and Reif [20] further extended Doo and Sabin's prior work on .. .ii,,,il I properties of
subdivision surfaces by deriving various necessary and it!!li wI conditions on smoothness for .[l!!. i. II
subdivision schemes. Specific subdivision schemes were :i! ~1 .1 by Schweitzer [21], Habib and Warren
August 20, 1998
DRAFT
[22], Peters and Reif [23] and Zorin [24]. Most recently, i, i 1i [25] presented a method for exact evaluation
of Catmull('I! I : subdivision surfaces at arbitrary parameter values.
B. s.., modeling using the ;,,,' based subdivisionsurface model
Although recursive subdivision surfaces are powerful for representing smooth geometric shapes of ar
bitrary l.l, ,., they constitute a purely geometric representation, and furthermore, conventional geo
metric modeling with subdivision surfaces 11i be [t! t111 for representing highly complicated objects.
For example, modelers are faced with the tedium of indirect shape modification and refinement through
timeconsuming operations on a large number of (most often irregular) control vertices when using I ,i
cal subdivision surfacebased modeling schemes. Despite the advent of advanced 3D graphics interaction
tools, these indirect geometric operations remain nonintuitive and laborious in general. In addition, it
ii not be enough to obtain the most ''I !! surface that interpolates a set of (ordered or unorganized)
data points. A certain number of local features such as bulges or inflections i! be strongly desired while
requiring geometric objects to satisfy global smoothness constraints in geometric modeling and computer
graphics applications. In contrast, i'1! 1, based modeling provides a superior approach to shape model
ing that can overcome most of the limitations associated with traditional geometric modeling approaches.
Freeform 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. Furthermore,
the equilibrium state of the model is characterized by a minimum of the deformation ii. i of the model
 ,i, 'r i to the imposed constraints. The deformation i i 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  i li!i_ in realworld i,1! ii ,1 behavior.
Freeform deformable models were first introduced to computer graphics and visualization in Terzopou
August 20, 1998
DRAFT
los et al. [26] and further developed by Terzopoulos and Fleischer [27], Pentland and \\ ill ,!!! [28],
Metaxas and Terzopoulos [29] and Vemuri and I ,1l li. ; [30]. Celniker and Gossard [31] developed
a I, i, for interactive freeform design based on the finite element optimization of i functionals
proposed in [27]. Bloor and \\ i! . [32], [33], Celniker and Welch [34] and Welch and \\ i !:,i [,;.] proposed
deformable Bspline curves and surfaces which can be designed by imposing the shape criteria via the
minimization of the ii. i functionals ,i. lr I to hard or soft geometric constraints through Lagrange
multipliers or i" ii ,l methods. Qin and Terzopoulos [36], [37], [38] developed dynamic NURBS (D
NURBS) which are very sophisticated models suitable for representing a wide i. of freeform as well
as standard !! 1 1 i shapes. The DNURBS have the advantage of interactive and direct manipulation of
NURBS curves and surfaces, resulting in ,1!. ii11 meaningful hence intuitively predictable motion and
shape variation.
A severe limitation of the existing deformable models, including DNURBS, is that they are defined on
a rectangular parametric domain. Hence, it can be very lIi 1iII to model surfaces of arbitrary genus using
these models. In [10], DeRose et al. assigned material properties to control meshes at various subdivision
levels in order to simulate cloth dynamics using subdivision surfaces. Note that, they assign i1 i, ,1
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 I i 11 to achieve. Previously we introduced a dynamic
Catmull('C! I subdivision surface model [39], [40] which combined the benefits of subdivision surfaces
for modeling arbitrary l,,.I.1._ as well as that of dynamic splines for interactive shape manipulation
by 1',,1, i 1 synthesized forces. The dynamic (modified) butterfly subdivision surface model formulated
and developed in this paper aims to achieve the same longterm objective, i.e., a formal mechanism of
allowing the modeler to directly and intuitively manipulate the smooth limit surface of arbitrary I .1. 
as if they were seamlessly sculpting a piece of realworld " !I '. However, this new model is superior to
our previously reported research in several significant aspects which will be discussed in detail later in
Section ID. In this paper, we derive a novel technique for locally parameterizing the smooth limit surface
generated by the modified butterfly subdivision surface algorithm which embeds the proposed model in
a dynamic framework in a straightforward manner. The model can be initialized interactively by an
August 20, 1998
DRAFT
userdefined control mesh and is amenable to further sculpting via direct application of !i l i ,1 forces
to ,i i region of object surface. The formulation and implementation details are discussed in subsequent
sections.
C. !ji,1, and motion estimation using the ;/,', .based subdivisionsurface model
The dynamic subdivision surface model has been developed primarily for modeling arbitrary (known)
l I'1" ,! ,_ where modelers can directly manipulate the limit surface by ;,.1,, l  ,l. I .1 ,1 forces in an
intuitive fashion. However, as we have shown in our earlier work [41], another important application of the
dynamic subdivision surfaces is in nonrigid shape and motion reconstruction/recovery. Accurate shape
recovery requires distributed parameter models which I i. ,11 possess a large number of degrees of free
dom. On the other hand, Iti ! ii 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 [30], [38], [41], [42], [43],
[44], [45], l11], [47]. 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 in be infeasible to globally parameterize shapes of arbitrary 1i i.1 ._ . A 1i.1 . based
model best 11r i the above mentioned criteria is an ideal candidate for a solution to the shape recovery
problem for obvious reasons.
Deformable models which come in in, i, varieties, have been used to solve the problem in the I'1! ' 
based modeling paradigm. These models involve the use of either fixed size [30], [43], [48], [49], [50]
or adaptive size [44], i], [51], [52],, ;], [54] 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 [43], [46] need a
large number of degrees of freedom for deriving a smooth and accurate representation. In addition, they
August 20, 1998
DRAFT
can not represent shapes with known arbitrary I. 1. 1, _. Moreover, almost all of these schemes require a
globally parameterized mesh as their input which 11i be infeasible at times.
Our previous dynamic subdivision surface model [39], [40], [41] offered an elegant solution to the above
mentioned problem as it could recover complex shapes in a hierarchical fashion using very few degrees
of freedom without requiring parameterized input mesh. However, the model proposed in this paper
outperforms the previous one in the compactness of the model representation. We will show experimental
results in support of this claim. We will also demonstrate the potential of this model in motion tracking
and visualization of a dynamically deforming shape from a time sequence of volumetric data sets. Like
the previous model, the dynamic modified butterfly subdivision surface model also deforms under the
influence of ,il i 1 forces to fit the i1 il i,1 i!_ shape in the given range or volumetric data set via
the principle of i. minimization.
D. Contributions
In this section, we summarize the contributions of the present work along with the advantages of
the proposed dynamic modified butterfly subdivision surface model over the dynamic Catmull('I! !:
subdivision surface model [39], [40]. The primary contributions of the present work are as follows.
We  l. 1 1, ,1 i11 derive a local parameterization scheme for the modified butterfly subdivision sur
faces in a hierarchical  1. and subsequently the initial control polyhedron can be viewed as the
parametric domain of the 1i'1! i based smooth limit surface.
We treat the dynamic subdivision surface in the limit as a "'. i!" 1'1 1,1 1 object and represent the
smooth limit surface as a set of novel finite elements whose shape (basis) functions are derived using
the modified butterfly subdivision scheme. We envision that this new finite element method will prove
to be useful not only in the areas of computer graphics and geometric design, but also in engineering
We use modified butterfly subdivision techniques to create a surface model that incorporates mass and
damping distribution functions, internal deformation i!. forces, and other 1i1!1 11 i1 quantities. We
also I. I ,11 formulate the motion equations for this (modified) butterfly subdivision surface
August 20, 1998
DRAFT
whose degrees of freedom are the collection of initial userspecified control vertices. Therefore, the
advantages of both the '1! i. based modeling I'!'i1.... i! and the geometric subdivision schemes are
presented within a single unified framework.
Users will be able to manipulate this i,1! i1 based model in an arbitrary region, and the model
responds naturally (just like the realworld object) to this force application. This shape deformation
is quantitatively characterized by the timevarying displacement of control points that uniquely define
the geometry of the limit surface.
We develop algorithms and procedures which approximate our novel finite element using a collection
of linear and/or bilinear finite elements il. i I to the implicit geometric constraints enforced by the
butterfly subdivision rules. This hierarchicallystructured approximation is optimal in the sense that
it can fall within :,,i userspecified error tolerance.
Although the longterm goals inherent in the lo' ;* 1 1 1, i. .1 dynamic Catmull('I! !: model [39],
[40] are the same as our current endeavor of deriving dynamic modified butterfly subdivision scheme, the
research presented in this paper achieves them in a much more elegant fashion. I i! I of all, the finite
element implementation of the dynamic Catmull(' : subdivision surface is specific to the subdivision
technique involved where a diversity of complicated finite elements must be i il1 .1 in order to account
for the special cases and can not be readily generalized to other approximating subdivision schemes in
a straightforward way. However, the finite element techniques developed in this paper can be easily
generalized to other interpolatory (as well as I I '! i !!1 i1 i subdivision schemes involving triangular (or
nsided) meshes. Second, for some specific cases the plate i i of Catmull('I! !: subdivision surface
diverges as shown in [7], hence various case :iI 1  need to be performed to derive the internal i. of
the dynamic Catmull('I ! subdivision surface model. The internal deformation ii. i of the dynamic
butterfly scheme can be derived in an unified fashion for a spectrum of scenarios. Third, for both models,
we need to derive the closest point on the limit surface from a given point in 3D for force applications.
The calculation overhead involved in this process is significantly less for the butterfly case as it is an
interpolatory scheme where all vertices at various levels of subdivision lie on the limit surface, and the
search space can be reduced rapidly in a hierarchical fashion. The situation is quite 1!I! i, II for Catmull
August 20, 1998
DRAFT
('I! iil subdivision scheme since it is an approximating scheme, and the technique used for finding the
closest model point [48] is computationally expensive. The force application is vital to ; i' I'1! '. 
based model, and hence the adopted computationally inexpensive method to find the closest point for
the proposed scheme in this paper has very significant advantages over the previous model. Lastly, it
has been empirically found that the recovered shape is more compact (less number of degrees of freedom)
when using the proposed model in comparison with our earlier reported model using the same data set
and modelfitting criteria.
E. Overview
The rest of the paper is organized as follows: Section II presents the detailed formulation of the dynamic
modified butterfly subdivision surfaces. The implementation details are presented in Section III followed
by experimental results in Section IV. I ii I1I we make the concluding remarks in Section V.
II. FORMULATION
In this section, we provide a I. i!,, formulation of the dynamic subdivision surface model. I i, I
we ', i. f review the modified butterfly subdivision scheme. Next we introduce a local parameterization
scheme which will facilitate the formulation of the smooth limit surface as a function of the control point
positions defining the initial mesh. This parameterization scheme is then used to derive the dynamic
model. It ii , be noted that these techniques can be generalized to define and construct a generic
dynamic framework for other trianglebased subdivision surface schemes as well, but we will focus only
on the modified butterfly subdivision technique in this paper.
A. I !, ( ...... "1', ,I) butterfly subdivision scheme
The butterfly subdivision scheme [11], like 11ii i~ 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
August 20, 1998
DRAFT
\ \ I \
triangular face in the previous level and hence, interpolates the coarser mesh in the previous level. In
S\ \1 \ I \ / \
I II" I
(a) (b)
Fig. 1. (a) The control position an with triangular faces; (b) Mvesh obtained after one subdivision step.
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 spilt by adding a new vertex whose position is obtained by
an ;I !!!r, combination of the neighboring vertex positions 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 c' be noted that all the
newly introduced vertices (marked in blue) corresponding to the edges in the original mesh have valence
(degree) 6, whereas the position and valence of the original vertices (marked in red) 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 eightpoint stencil as shown in i 2, , i, The name of the scheme originated from the
" 1.i I. i II 'like configuration of the contributing vertices. The weighing factors for I h!! I. i!. contributing
vertex positions are shown in Fig.2(b). The vertex e l in the j + 1th level of subdivision, corresponding
to the edge connecting vertices vi and v2 at level j, is obtained by
ei+1 = 0.5(v + v?) + D. \1 + v ) w(v 3 + v + v + v3), (1)
S< w < and v denotes the position of the h ver at the jh level.
where 0 < w < 1, and vj denotes the position of the ith vertex at the jth level.
August 20, 1998
DRAFT
j w w
5 0.5
7 2 8 w w
(a) (b)
Fig. 2. (a) The contributing vertices in the jth level for the vertex in the j+lth level corresponding to
the edge between v1 and v2; (b) the weighing factors for different vertices.
The butterfly subdivision scheme produces a smooth C1 surface in the limit except at the extraordinary
points corresponding to the extraordinary vertices (vertices with valence not equal to 6) in the initial
mesh [13]. All the vertices introduced through subdivision have valence 6, and therefore, the number of
extraordinary points in the smooth limit surface equals the number of extraordinary vertices in the initial
mesh. Recently, this scheme has been modified by Zorin et al. [13] to obtain better smoothness properties
at the extraordinary points. In [13], all the edges have been categorized into three classes : (i) edges
connecting two vertices of valence 6 (a 10 point stencil, as shown in I i, ;, i is used to obtain the new
vertex positions corresponding to these edges), (ii) edges connecting a vertex of valence 6 and a vertex of
valence n r 6 (the corresponding stencil to obtain new vertex position is shown in Fig.3(b), where q = .75
is the weight associated with the vertex of valence n r 6, and si = (0.25 + cos(27i/n) + 0.5cos(47i/n))/n,
i = 0, 1,..., n 1, are the weights associated with the vertices of valence 6), and (iii) edges connecting two
vertices of valence n r 6. The last case can not occur except in the initial mesh as the newly introduced
vertices are of valence 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.3(b) at each of those two extraordinary vertices.
August 20, 1998
DRAFT
(a) (b)
Fig. 3. (a) The weighing factors of contributing vertex positions for an edge connecting two vertices of
valence 6; (b) the corresponding case when one vertex is of valence n and the other is of valence 6.
B. Local parameterization of the limit '.. f f,..I
Oftentimes, the smooth limit surface defined by the recursive subdivision process is of;' ,1 'i 1, 1. I. 1 1_ *
where a global parameterization in , not be possible. Nevertheless, we can locally parameterize the limit
surface over the domain defined by the initial mesh following a similar approach described in [55]. The
idea is to track ; i arbitrary point on the initial mesh across the meshes obtained via the subdivision
process (see Fig.4 and 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. Note that, the limit surface can be represented by the same number of smooth triangular patches
August 20, 1998
DRAFT
x X * *
Fig. 4. The smoothing effect of the subdivision process on the triangles of the initial mesh.
as that of the triangular faces in the initial mesh. Therefore, we can express the limit surface s as
a
s = k, (2)
k=l
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.
We are now ready to describe the parameterization of the limit surface over the initial mesh. The
process is best explained through the following example. We choose a simple planar mesh shown in
I i, ,, 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 drawn in black in Fig.5. After
one step of subdivision, the initial mesh is refined by addition of new vertices which are colored green.
Another subdivision step on this refined mesh leads to a finer mesh with introduction of magenta colored
new vertices. It !1 i be noted that ., ,. point inside the smooth triangular patch in the limit r ,,tII..
corresponding to the face abc in the initial mesh depends (..l', on the vertices in the initial mesh which
are within the 2neighborhood 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 I i, ;,: 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 contributing vertices at this level of subdivision, for
example, the (green colored) 1neighbors of the vertex b (except d and e) in I i, .(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,
August 20, 1998
DRAFT
(a) (b)
Fig. 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.
August 20, 1998
DRAFT
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 vobc be the collection of
vertices in the initial mesh which are within the 2neighborhood of the vertices a, b and c (marked black
in i i, ,, ,). Let the number of such vertices be r. Then, the vector vObc, 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 (3r x 3r) subdivision matrices (Aabc)t, (Aabc)l, (Aabc)r and (Aabc)m such that
vadf = (Aabc)tVabc,
Vbed = (Aabc) Vabc,
IVfe = (Aabc),V bc
VCf( abc,
V f = (Aabc),, b (3)
where the subscripts t, I, r and m denote top, left, right and middle triangle positions respectively (indi
cating the relative position of the new triangle with respect to the original triangle), andVadfi, v d, v1fel
and Vef are the concatenation of the (x, y, z) positions for the vertices in the 2neighborhood of the cor
responding triangle in the newly obtained subdivided mesh. The new vertices in this level of subdivision
are marked green in Fig.5(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.3 are the same. This concept is further illustrated in
Fig.6.
C i11 out one more level of subdivision, along with the old vertices, we get a new set of vertices
which are marked in magenta in I i, .(c). Adopting a similar approach as in the derivation of Eqn.3, we
obtain
Vdgi = (Abed)tVbed
2
V hg = (Abed)1V ed
V eih = (Abed),vred
V hi = (Abed),,Ved (4)
August 20, 1998
DRAFT
August 20, 1998

DRAFT
Fig. 6. Different set of control points at a subdivided level obtained by ;1:''1;! ,ltl' i' !ii subdivision
matrices on a given set of control points in a coarser mesh.
August 20, 1998
DRAFT
The relative position of the triangular face ~I.,' in Fig.5(c) with respect to the triangular face bed is
topologically the same as of the triangular face adf in i, .(b) with respect to the triangular face abc.
Therefore, we can write (Abed)f = (Aabc)1. Using similar reasoning, Eqn.4 can be rewritten as
vdgi = (Abed)tVled = (Aabc)tled
Vhg = (Abed)lVed = (Aabc),vbed
h = (Abed)rVed = (Aabc)rVLed
Vghi = (Abed),nVbed = (Aabc),nVed. (5)
Combining Eqn.3 and Eqn.5, we get
Vgi = (Aabc)t(Aabc)lv b,
bg= (Aabc)l(Aabc)lVabc>
Veih = (Aabc)r(Aabc)lVabc,
Vghi = (Aabc),(Aabc)lVabc. (6)
Let x be a point with '. "1 i. i il,.; coordinates (abc, ab, Yabc) inside the triangular face abc. When the
initial mesh is subdivided, x becomes a point inside the triangular face bed with '.i id i i. coordinates
(ab~d,>Pled,' .). Another level of subdivision causes x to be included in the triangular face 1, with
i I. !i coordinates (a di,dgi, T gi). Let sabc denote the jth level approximation of the smooth
triangular patch Sabc in the limit surface corresponding to the triangular face abc in the initial mesh. Now
bc can be written as
r r r
vObc bx ,c,....,ay,by,cy, ..,a,, b ,c ,, .. ]T
where the subscripts x, y and z indicate the x, y and z coordinates respectively of the corresponding
vertex position. The expressions for v ed and vdg, can also be written in a similar manner. Next, we
construct the matrix Bbc as follows:
bcc, ,bc, T bc, 0,...,0, 0,...,0,0,...,0
B. 0abc ab, 0,bc, b 0, ...0, 00, ..,0
0, 0. ,0 ... ab, O bc bc ,..., 0
August 20, 1998
DRAFT
The matrices Bed and B g. can also be constructed in a similar fashion. We can now write sbc(X),
abc(x), and S bc(x) as
sbc(x) = B0 bc(X)V
bc(x) = Bd(x)vbed = B (x)(Albc)l AVb, (7)
S 2.. B 0
Sbc(x) = B2. x = B xi(Aabc)Ved = B .x(Aabc)(Aabc)lVabc.
Proceeding in a similar way, the expression for s bc(x), jth level approximation of Sabc(x), is given by
3
Sbc(X) = B,.(x) (Aabc), ...(Aabc)t(Aabc)V bV
= B4 ,\xi(Abc)VOb
= bc(X)Vab, (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), (Abc) =
(Aabc),n ... (Aabc)t(Aabc)1 and BIbc(x) = B1 ,fi(Ab). It ii, be noted that the sequence of i.I 'I ,
(Aabc)t, (Aabc)1, (Aabc)r and (Aabc) depends on the triangle inside which the tracked point x falls after
each subdivision step. I ii!! we can complete the local parameterization process by writing
0 +
Sabc(x) = ( lim Bb,(X))Vbc = Babc(X)Vbc. (9)
In the above equation, Bab, is the collection of basis functions at the vertices of vbc. It II also be
noted that the modified butterfly subdivision scheme is a stationary subdivision process, and hence new
vertex positions are obtained by ,t !!ii. combinations of nearby vertices. This guarantees that each row of
the matrices (Aabc)t, (Aabc)1, (Aabc)r and (Aabc), sums to one. The largest eigenvalue of such matrices
is 1 and therefore the limit in Eqn.9 exists. Now, if we assume that the triangular face abc is the kth face
in the initial mesh, then Eqn.9 can be rewritten as
sk(x) = Bk(x = Bk(x)Akp, (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 postmultiplied by p, selects the vertices vA controlling 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 (3r x 3t) matrix and Bk (x) is a (3 x 3r) matrix.
August 20, 1998
DRAFT
Combining Eqn.2 and Eqn.10, we get
s(x) = ( B (x)Ak)p = J(x)p, (11)
k=l
where the (3 x 3t) matrix J 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.
C D ..... ', ,
In the previous section, we have derived an expression of the smooth limit surface obtained via infinite
subdivision steps as a function of the control vertex positions in the initial mesh. We now treat the vertex
positions in the initial mesh defining the smooth limit surface s 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
s(x,p) =J(x)p, (12)
where an overstruck dot denotes a time derivative and x E So, So being the domain defined by the initial
mesh. The 1i1! '. of the dynamic subdivision surface model is based on the i. i!:i i relationship of
Lagrangian dynamics [,1,] and is formulated in an analogous way to that in [39], [40].
In an abstract 1i1! 1  .. i let pi(t) be a set of generalized coordinates which are functions of time
and are assembled into the vector p. Let fi(t) be the generalized force assembled into the vector fp and
acting on pi. Then, the Lagrangian equation of motion can be expressed as
d OT OT OF OU
dti p+ + = fi, (13)
dt dpi dpi dpi dpi
where T, F and U are the kinetic, dissipation and potential. . i _ respectively.
Let p be the mass l. i,i function of the surface. Then the kinetic i. i, of the surface is
T = 1 p(x)sT(x)s(x)dx = IMp, (14)
2 xeSO 2
where (using Eqn.12) M = Jxo P(x)JT(x)J(x)dx is the (3t x 3t) mass matrix. Similarly, let 7 be the
damping 1I. i! i function of the surface. The dissipation !i. i is
F = 1 f (x)T (x)s(x)dx = IT D, (15)
2 xxxx e 2
where D = fS, o 7(x)JT(x)J(x)dx is the (3t x 3t) damping matrix. The potential i i of the smooth
August 20, 1998
DRAFT
limit surface can be expressed as
U pTKp, (16)
2
where the (3t x 3t) stiffness matrix K is obtained by assigning various internal energies to a discretized
approximation of the limit surface and is detailed in Section III.
Using the expressions for the kinetic, dissipation and potential i i in Eqn.13, we get the motion
equation given by
Mpi + Dp + Kp = fp. (17)
The generalized force vector fp, which can be obtained through the principle of virtual work 'i,]j, is
expressed as
f= J(x)f(x, t)dx. (18)
We can apply various I I," of forces on the smooth limit surface, and the limit surface would evolve over
time according to Eqn.17 to obtain an equilibrium position characterized by a minimum of the total model
 I I 1
C.1 Multilevel Dynamics
We have developed a dynamic framework 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.18), and the initial
mesh (as well as the 1l. i1 i!i smooth limit surface) deforms continuously over time until an equilibrium
position is obtained. The deformation of the surface in response to the applied forces is governed by
the motion equation (Eqn.17). \\ !iili our 1i1! i, based modeling framework, the limit surface evolves
as a consequence of the evolution of the initial mesh. One can apply various I I," 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
August 20, 1998
DRAFT
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 subdivision step applied to the initial mesh.
This new mesh has more number of vertices but defines the same limit surface. For example, after one
step of modified butterfly subdivision, the initial degrees of freedom p (refer to Eqn.11 and Eqn.12) in
the dynamic 1. I will be replaced by a larger number of degrees of freedom q, where q = Ap. A is
a global subdivision matrix of size (3s x 3t) whose entries are uniquely determined by the weights used
in the modified butterfly subdivision scheme (see Section IIA for the weights). Thus, p, expressed as a
function of q, can be written as
p = (ATA) ATq = Aq, (19)
where At = (ATA) AAT. Therefore, we can rewrite Eqn.11 and Eqn.12 as
s(x) = (J(x)At)q, (20)
and
s(x,q) = (J(x)At)q, (21)
respectively. Now we derive the equation of motion for this new subdivided model involving a larger
number of control vertices namely q. We need to recompute the mass, damping and stiffness matrices
for this " i level. The structure of the motion equation as given by Eqn.17 remains unchanged, but
the .1i1 ,. 1!..i! ,,1 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
Mql + Dq Kqq = fq, (22)
where Mq = fxC, p(x)(At) JT(x)J(x)Atdx, S1 being the domain defined by the newly obtained sub
divided mesh. The derivation of D,, Kq and fq follow suit.
It ii be noted that further increase in the number of control vertices, if necessary, can be obtained via
another level of subdivision. Therefore, multilevel dynamics is achieved through recursive subdivision on
the initial set of control vertices. Users can interactively choose ;,i subdivided mesh as the control mesh
for the dynamic model depending on their needs. Alternatively, the  1. i can automatically determine
August 20, 1998
DRAFT
the most suitable control mesh for certain applications based on an applicationspecific criteria.
III. FINITE ELEMENT IMPLEMENTATION
(b)
Fig. 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.
In the previous section, we have expressed the smooth limit surface as a function of the control vertex
positions in the initial mesh, and have assigned mass and damping distribution, internal deformation
. i. and forces to the limit surface to develop the corresponding ,1!. i. 1 model. In this section, we
describe the implementation of this 1'1! i, ,1 model using the finite element method.
In Section II we 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
August 20, 1998
DRAFT
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 and 7). We now provide a detailed discussion on how to derive the
mass, damping and stiffness matrices for these elements. These elemental matrices can be assembled to
obtain the global '1!i i1 ,1 matrices M, D and K, and a numerical solution to the governing secondorder
differential equation as given by Eqn.17 can be obtained using finite element ; ,! I1 i techniques [57].
We use the same example as in Section II (refer I i, ) to develop the related concepts. The concept of
elements along with the control vertices and their corresponding domains is further illustrated in Fig.7.
We will now show how to derive the mass, damping and stiffness matrices for the element corresponding
to the triangular face abc in I i, ' The derivations hold for ;,i element in general.
A. ./ ./. !,,/I mass and damping matrices
The mass matrix for the element given by Sabc (Eqn.9) can be written as
I = s p(x){Babc(x)} {Babc(x)}dx. (23)
JXESabc
However, from Eqn.9 we know that the basis functions corresponding to the vertices in v b which are
stored as entries in Bab, are obtained as a limiting process. These basis functions do not have :i 111 1 iG
form, hence computing the inner product of such basis functions as needed in Eqn.23 is a challenging
problem. In [55], an outline is provided on the computation of these inner products without performing ;, i
integration. In this paper, we develop a different yet simpler approach 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 4J 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
43
NI ,i = / (x){BJbc(x)}T{Bbc(x)}dx. (24)
The jth level approximation of the corresponding basis functions can be explicitly evaluated (refer Eqn.8
for an expression of BIbc). An important point to note is that Eqn.8 involves several matrix multiplications
and hence can be very expensive to evaluate. However, the matrix (AJbc)(= (Aabc), .. (Aabc)f(Aabc)1)
in Eqn.8 encodes how vertices in the 2neighborhood of the triangular face uvw at level j are related to
August 20, 1998
DRAFT
the vertices in the 2neighborhood of the triangular face abc in the initial mesh. In the implementation,
we keep track of how a new vertex is obtained from the contributing vertices in its immediate predecessor
level. If we move up from level j to level 0, we get the information stored in (AJbc) without forming iii
local subdivision matrices and thus avoiding subsequent matrix multiplications.
By choosing a , L!i .I Ii high value of j, we achieve a reasonably good approximation of the elemental
mass matrices. We eliminate the computations involved in evaluating the integrals in 51!! 21 by using
discrete mass 1. i function which has nonzero values only at the vertex positions of the jth subdivision
level mesh. Therefore, the approximation of the mass matrix for the element can be written as
k
NI = (v ){B b(vJIT ) {BJbc(vj)}, (25)
i=l
where k is the number of vertices in the triangular mesh with 4J faces. This approximation has been
found to be very effective and. tin. ii for implementation purposes. The elemental damping matrix can
be obtained in a similar fashion.
B. I./ .....I! 'i stiffness matrix
We assign the internal. ii. to each element in the limit surface, thereby defining the internal. . i
of the smooth subdivision surface model. We take a similar approach as in the derivation of the elemental
mass and damping matrix and assign the internal. ,. i i to a jth level approximation of the element.
In this paper, we use three I ,. of internal i.  tension, stiffness and spring ii i _. For the
examples used throughout the paper, this i at the jth level of approximation can be written as
E = (E b t + (E +bct + (E b9p (26)
where the subscripts t, st and sp denote tension, stiffness and spring respectively.
The expression for the tension !. i which is essentially equivalent to the first order strain (membrane)
. [58], is
1 2
(E bc t 2 kt vI v
1 kt1{v bc T(KT b9 Vbc}, (27)
where kt is a constant, v3 and v the 1th and mth vertex in the jth level mesh, are in the 1
neighborhood of each other, Q is the domain defined by all such vertex pairs, and vibc is the concatenation
August 20, 1998
DRAFT
of the (x,y,z) positions of all the vertices in the jth subdivision level of the triangular face abc in the
initial mesh.
Similarly, the expression for stiffness. . i which is equivalent to the second order strain (thin plate)
* i [58], can be written as
1 2
(Ebct 2kt v 2v + vI
Q
= I bc tVJ}, kI (28)
2n ab, 1 abc)J V ab,
where v3, v3 and v3 are the three vertices of a triangular face. The summation involves three terms
corresponding to each triangular face, and is over all the triangular faces in the mesh at the jth level of
subdivision.
I i!i !! we add a spring. !i. term which is given by
1 (k,n),,( Iv v , (v)
2
= iv bc (K c sp bc}, (29)
where (k,im),p is the spring constant, Q is the domain as in Eqn.27 and l is the natural length of the
spring connected between v and vJ. It 11 i be noted that the entries in (K c)sp depend on the distance
between the connected vertices and hence, (Kbc)s 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. i we get
jb1 kT ab)st (I(b {
E3 bc = 2Vbc} {kt(Kbc)t + kst(Kbc)st + (Kbc), p}{Vab}
= {v ab ba(Keb{v b
= {(A bc) IV T }} (K b)(A b){v
= { abc} {(ab) (Kab )(abc)} Vab,, (30)
where (Abc) and v are same as in Eqn.8. Therefore, the expression for the elemental stiffness matrix
is given by
Kabc = (Aibc) (Kibc)(Abc). (31)
It i i be noted that the matrix multiplications for constructing Kabc are avoided in the implementation
by following the same technique described in Section IIIA.
August 20, 1998
DRAFT
C. Force Application
The force f(x,t) in Eqn.18 represents the net effect of all externally applied forces. The current im
plementation supports spring, inflation as well as imagebased forces. However, other I i" 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 jth level approximation mesh), the net applied spring force being
f(x, t) = k(do s(x, t))(x xo)dx, (32)
where 6 is the unit impulse function i!i'l, il! f(xo,t) = k do 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 t!l
To recover shapes from 3D image data, we synthesize imagebased forces. A 3D edge detection is
performed on a volume data set using the 3D MongaDT i i. II)) operator [59] 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 P(x,,) (33)
I vP(x,y, z) /
where A controls the strength of the force. The applied force on each vertex at the jth approximation level
is computed by trilinear interpolation for evaluating Eqn.18 in Cartesian coordinates. More sophisticated
imagebased forces which incorporate regionbased information such as gradients of a thresholded fuzzy
voxel classification can also be used to yield better and more accurate shape recovery. It i!i be noted
that we can apply spring forces in addition with the imagebased forces by placing points on the boundary
of the region of interest in each slices of the 3D volume (MR, CT etc.) image data.
D. Discrete D,. ..... '', Equation
The t1! 1. i i I equation given by Eqn.17 is integrated through time by discretizing 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
August 20, 1998
DRAFT
where discrete derivatives of p are calculated using
p(t + At) 2p(t) + p(t At) (34)
pj(t + At) = ) 2(34)
At2
and
(t At) (t + At) p(t At)
p(t+Af)= (. )
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 tin i reasons. For the time i stiffness matrix, we recompute K at
each time step. Using Eqn.17, Eqn.34 and Eqn.35, the discrete equation of motion is obtained as
(2M + DAt + 2At2K)p(t + At) = 2At2fp(t + At) + (DAt 2M)p(t At) + 4Mp(t). (36)
This linear ,. in of equations is solved iteratively between each time step using the '.. il_ . gradient
method [60], [61].
E. Model Subdivision
evolution evolution
subdiv sion same limit surface
at equili rium with more patches
evol tion evol tion
(a) (b)
Fig. 8. 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.
The initialized model grows dynamically according to the equation of motion (Eqn.17). 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
August 20, 1998
DRAFT
replacing the original initial mesh by a new initial mesh obtained through one step of butterfly 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
1 i, Model subdivision might be needed to obtain a very localized effect on a smooth limit surface. For
a shape recovery application, one 111 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 imagebased forces, if the model !! i 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 abovementioned replacement scheme until the model !! i is itLil, ! small and the
change in model < i between successive iterations becomes less than a prespecified tolerance.
IV. RESULTS
The proposed dynamic (modified) butterfly subdivision surface model can be used to represent a wide
,i. I of smooth shapes with arbitrary genus. The smooth limit surface can be sculpted by ;1''1 i!_
synthesized forces in a direct and intuitive way in shape design applications. The [1!Il. i! i,_ shape in
a range or volume data set can also be recovered hierarchically using the proposed dynamic (modified)
butterfly subdivision surface model. Before illustrating the application results of the proposed model, we
would like to point out the advantages of our model over the existing techniques. It is known that shapes
of arbitrary lI ,.1. _ can be modeled either by explicit patching or by subdivision surfaces. C. .iJ ii[il
constraints across patches need to be imposed with care while modeling arbitrary I, 1,. .1 _ using explicit
patching techniques. Subdivision surfaces do not face this problem, but one has to manipulate control
vertex positions at various levels of subdivision to obtain a desired smooth shape of arbitrary l 1 I ,1_ * us
ing subdivision surfaces. Our dynamic subdivision surface model overcomes this problem by allowing the
user to obtain the desired effect on the limit surface by direct sculpting via the application of synthesized
forces. To recover shapes from a given set of points in 3D, the existing subdivision surface based tech
August 20, 1998
DRAFT
(b) (c)
(f) (g)
Fig. 9. (a), (b), (c) and (d) : Initial shapes; (e), (f), (g) and (h) : the corresponding modified shapes
after interactive sculpting via force application.
niques resort to complex techniques to derive a mesh for the '1i.1 l i1 il shape, and then I i.,, 11 mesh
optimization techniques are used to obtain a compact representation of the same. Our model recovers the
shape from a set of points in an t. !. ii hierarchical way ,i, simple mesh of the same I, 1. . can be
used as an initial mesh which will evolve over time to fit the given data and depending on the error of
fit achieved, it will automatically refine itself until prescribed error of fit is obtained. In the rest of this
section, we elaborate on these points via examples.
In a shape modeling application, the user can specify :i!i mesh as the initial (control) mesh, and
the corresponding limit surface can be sculpted interactively by ;,1,1,1  ,lli .11 forces. In I i, '
August 20, 1998
DRAFT
we show several initial surfaces obtained from different control meshes and the corresponding modified
surfaces after interactive sculpting. To change the shape of an initial surface, the user can attach springs
from h !!. I. ,I points in 3D to the nearest points on the limit surface such that the limit surface deforms
towards these points to generate the desired shape. The initial mesh of Fig.9(a) has 125 faces and 76
vertices (degrees of freedom), the initial mesh of the closed cubelike shape in Fig.9(b) has 24 faces and
14 vertices. The one hole torus in Fig.9(c) and the corresponding modified shape in Fig.9(g) 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.9(d), is dynamically sculpted to the shape shown in Fig.9(h).
We have performed several experiments testing the applicability of our model to recover the I 1111. i l
shapes in range and volume data sets. In all the experiments, the initialized model has a control mesh
comprising of 24 triangular faces and 14 vertices whereas the control mesh of the fitted model has 384
triangular faces and 194 vertices. It i!! be noted that once an approximate shape is recovered, the model
is refined depending on the datafitting criteria, thereby increasing the degrees of freedom of the recovered
shape only when necessary. For an error in fit (defined as the maximum distance between a data point
and the nearest point on the limit surface expressed as a percentage of the diameter of the smallest sphere
enclosing the object) of approximately 3%, the initialized model is refined twice following the technique
described in Section IIIE. Also, the limit surface of i,! control mesh (of the desired genus) can be used
as the initialized model. However, an initial mesh with few degrees of freedom usually performs better in
terms of recovering a compact representation of the [L.II. ! i ,_ shape.
In the first shape recovery experiment, we depict the laser range data acquired from multiple views
of a light bulb in Fig.10(a). Prior to ;,1I1' ii, our algorithm, the data were transformed into a single
reference coordinate I. ii, The model was initialized inside the 1000 range data points on the surface
of the bulb. The fitted dynamic (modified) butterfly subdivision surface model is shown in Fig.10(b) and
(c). In the next experiment, the shape of a mechanical part is recovered from a range data set having
2031 data points (Fig.11). We also recover the shape of a human head from a range data set as shown in
Fig.12. The head range data set has 1779 points in 3D. It i!! 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 number of
August 20, 1998
DRAFT
(b) (c)
Fig. 10. (a) Range data of a bulb along with the initialized model, (d) the fitted dynamic butterfly
subdivision model, and (c) visualization of the shape from another view point.
(a) (b) (c)
Fig. 11. (a) Range data of a mechanical part along with the initialized model, (d) the fitted dynamic
butterfly subdivision model, and (c) visualization of the shape from another view point.
(a) (b) (c)
Fig. 12. (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.
August 20, 1998
DRAFT
(b) (c)
Fig. 13. (a) Data points i.1. iiflf i _ the boundary of the caudate nucleus on a : 11: 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.
data points present in the original range data set.
WVe present the recovery of the caudate nucleus (a cortical structure in human brain) from 64 1: Il
slices, each of size 256 x 256 in our next experiment. Fig.13(a) depicts a slice from this 11: 1 scan along
with a sparse set of points placed by an expert neuroscientist on the boundary of the shape of interest.
Fig.13(b) depicts the sparse data points (placed in each of the slices depicting the boundary of the shape
of interest) in 3D along with the initialized model. Note that points had to be placed on the boundary of
the caudate nucleus due to lack of image gradients delineating the caudate from the surrounding tissue in
parts of the image. Continuous image based forces as well as spring forces are applied to the model and
the model deforms under the influence of these forces until maximal conformation to the data is achieved.
The final fitted model is shown in Fig.13(c). We like to point out the fact that the recovered shape in [41]
using our previous dynamic subdivision surface model for the same data set has 386 degrees of freedom
and therefore, we achieve a factor of 2 improvement in the number of degrees of freedom required to
represent the model in this particular example.
In the last experiment, we animate the motion of the leftventricular chamber of a canine heart over
a complete cardiac cycle. The data set comprised of 16 3D CT images, with each volume image having
118 slices of 128 x 128 pixels. in I, we have recovered the shape from one data set using imagebased
(gradient) as well as pointbased forces. Once the shape is recovered from one data set, this fitted model
August 20, 1998
DRAFT
(b) (c)
(f) (g)
(j) (k)
(I') (o)
Fig. 14. Snapshots from the animation of canine heart motion over a cardiac cycle using the dynamic
butterfly subdivision model.
August 20, 1998
DRAFT
is used as an initialization for the next data set to track the shape of interest. The snapshots from motion
tracking are shown in Fig.14 for the 16 volume data sets. It 11 i' be noted that the control mesh describing
the smooth surfaces shown in Fig.14 has only 384 triangular faces with a total of 194 vertices as mentioned
earlier. This experiment clearly shows that our model can be used to track a shape of interest from a set
of time dependent volume data sets in an tin i. manner. Note that no other existing purely geometric
subdivision surface technique can be used with (time i i,_i continuous data sets.
V. CONCLUSIONS
In this paper, we have presented a novel finite element method to derive and :,i ,1 the new dynamic
model based on the modified butterfly subdivision surface scheme. The new ,1! based surface model
provides a direct and intuitive way of manipulating smooth shapes of arbitrary I. ..I i,1 and is very useful
for directly extracting and visualizing the ,I1ii,1 I1 i shapes in large range and volume data sets. The
proposed model has also been used successfully for nonrigid motion tracking from a temporal sequence of
volume data sets. We have developed a hierarchical local parameterization of the subdivision scheme which
is critical to the formulation of our dynamic model. We have combined material properties with geometric
entities, formulated the motion equations for our dynamic model, and incorporated the advantages of free
form deformable models into conventional subdivision scheme. Moreover, we have introduced an tin. III
hierarchical dynamic control for various applications. Our experiments indicate a promising future of the
proposed model in computer graphics, geometric modeling and scientific visualization. Furthermore, the
finite element techniques proposed in this paper should be of great interest to the engineering design and
, i i 1 . I ,~,~! , ii as w ell.
VI. ACKNOWLEDGMENTS
This research was supported in part by the NSF grant ECS9210648 and the NIH grant RO1LM05944
to B.C. Vemuri, the NSF ('.C\I:I..I award CCR9702103 and DMI9700129 to H. Qin. We wish to
acknowledge Dr. Tim Mi Ii. i. for his help, Dr. Gregoire Malandain for the edge detection software,
Dr. Hughes Hoppe and Dr. Kari Pulli for the range data, Dr. Dmitry Goldgof for the 4D CT images and
Dr. ('I! ,1 i,, Leonard for brain I, 11:1 images.
August 20, 1998
DRAFT
REFERENCES
[I] G.M. ('h.i.! ii "An algorithm for high speed curve generation," Computer Vision, Graphics and Image
Processing, vol. 3, no. 4, pp. 346349, 1974.
[2] D. Doo, "A subdivision algorithm for smoothing down irregularly shaped polyhedrons," in Proceedings on
Interactive Techniques in Computer Aided Design, 1978, pp. 157 165.
[3] M. Sabin, !' use of piecewise forms for the numerical representation of shape, Ph.D. thesis, Hungarian
Academy of Sciences, Budapest, 1976.
[4] E. Catmull and J. ('!l. I i1... ,i  generated Bspline surfaces on arbitrary topological meshes," Computer
Aided Design, vol. 10, no. 6, pp. 350 355, 1978.
[5] C. Loop, Smooth subdivision surfaces based on triangles, M.S. thesis, University of Utah, Department of
Mathematics, 1987.
[6] H. Hoppe, T. DeRose, T. Duchamp, 11. Halstead, H. Jin, J. McDonald, J. Schweitzer and WV. Ii. ,.
I'l... . smooth surface reconstruction," in Computer Graphics Proceedings, ACM SIGGRAPH, July
1994, Annual Conference Series, pp. 295 302.
[7] M. Halstead, M. Kass and T. DeRose, I.!1. i. !, fair interpolation using Catmull('!.l. I surfaces," in Computer
Graphics Proceedings, ACM SIGGRAPH, August 1993, Annual Conference Series, pp. 35 44.
[8] J. Peters and U. Reif, "The simplest subdivision scheme for smoothing polyhedra," ACM .... .. on
Graphics, vol. 16, no. 4, pp. 420 431, October 1997.
[9] T.W. Sederberg, J. Zheng, D. Sewell and M. Sabin, .'... i!!..i ii. i iri e subdivision surfaces," in Computer
Graphics Proceedings, ACM SIGGRAPH, July 1998, Annual Conference Series, pp. 387 394.
[10] T. DeRose, M. Kass and T. 1!i.._ i 411,..1 i ..ii surfaces in character animation," in Computer Graphics
Proceedings, ACM SIGGRAPH, July 1998, Annual Conference Series, pp. 85 94.
[11] N. Dyn, D. Levin and J.A. Gregory, "A butterfly subdivision scheme for surface interpolation with tension
control," ACM I ... .. on Graphics, vol. 9, no. 2, pp. 160 169, ApIl 1990.
[12] N. Dyn, S. Hed and D. Levin, 'iI. it ..!i schemes for surface interpolation," in Workshop in Computational
Geometry, A. C. et al., Ed., pp. 97 118. 1993.
[13] D. Zorin, P. Schrider and W. Sweldens, IiwI.. Ip'.l.!ii. subdivision for meshes with arbitrary topology," in
Computer Graphics Proceedings, ACM SIGGRAPH, August 1996, Annual Conference Series, pp. 189 192.
[14] L. Kobbelt, I iI. I 1..I..!I.. refinement by variational methods," in Approximation .. VIII, C. ('Ciii and
L. Schumaker, Eds., vol. 2 of Wavelets and Multilevel Approximation, \\..il!I Scientific Publishing Co., 1995,
pp. 217 224.
August 20, 1998
DRAFT
[15] L. Kobbelt, "A variational approach to subdivision," ComputerAided Geometric Design, vol. 13, pp. 743 
761, 1996.
[16] L. Kobbelt and P. Schr6der, "Constructing variationally optimal curves through subdivision," Tech. Rep.
CS 1. '.711. California Institute of Technology Computer Science Department Technical Report, 1997.
[17] 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.
[18] A.A. Ball and D.J.T. ',..i;, "Conditions for tangent plane continuity over recursively generated Bspline
surfaces," ACM I ... .. .. on Graphics, vol. 7, no. 2, pp. 83 102, 1988.
[19] A.A. Ball and D.J.T. '..i i, "An investigation of curvature variations over recursively generated Bspline
surfaces," ACM ... .. .. on Graphics, vol. 9, no. 4, pp. 424 437, 1990.
[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 and Application of Subdivision Surfaces, Ph.D. thesis, University of \\.dhi,! ..I
Seattle, 1996.
[22] A. il1.,il., and J. \\.! 1h, "Edge and vertex insertion for a class of C1 subdivision surfaces," Computer Aided
Geometric Design, to appear.
[23] J. Peters and U. Reif, "Analysis of generalized Bspline subdivision algorithms," SIAM Journal of .
Analysis, to appear, I, ..il..1I. at ftp://ftp.cs.purdue.edu/pub/jorg/9697agb.ps.Z.
[24] D. Zorin, ii....i l. !! of stationary subdivision on irregular meshes," Constructive Approximation, submit
ted, available as i .!1.. .I Computer Science Lab Tech. Rep. CSLI1 i' ._'
[25] J. '.In1 I.... I evaluation of Catmull('!., I subdivision surfaces at arbitrary parameter values," in Computer
Graphics Proceedings, ACM SIGGRAPH, July 1998, Annual Conference Series, pp. 395 404.
[26] D. Terzopoulos, J. Platt, A. Barr and K. Fleischer, I.!.1.1 ..17 deformable models," in Computer Graphics
Proceedings, ACM SIGGRAPH, 1987, Annual Conference Series, pp. 205 214.
[27] D. Terzopoulos and K. Fleischer, "Deformable models," '. Visual Computer, vol. 4, no. 6, pp. 306 331,
1988.
[28] A. Pentland and J. Williams, "Good l..iI It.i : Modal dynamics for graphics and animation," in Computer
Graphics Proceedings, ACM SIGGRAPH, 1989, Annual Conference Series, pp. 215 222.
[29] D. Metaxas and D. Terzopoulos, "Dynamic deformation of solid primitives with constraints," in Computer
Graphics Proceedings, ACM SIGGRAPH, July 1992, Annual Conference Series, pp. 309 312.
[30] B.C. Vemuri and A. Radisavljevic, l .liii i..,i !.. stochastic hybrid shape models with fractal priors,"
ACM I ... .. on Graphics, vol. 13, no. 2, pp. 177 207, Apil 1994.
August 20, 1998
DRAFT
[31] G. Celniker and D. Gossard, "Deformable curve and surface finite elements for freeform shape design," in
Computer Graphics Proceedings, ACM SIGGRAPH, July 1991, Annual Conference Series, pp. 257 266.
[32] M.I.J. Bloor and M.J. \\ i,..i, i,. p' ii ..; iiPDE surfaces in terms of Bsplines," Computer Aided Design,
vol. 22, no. 6, pp. 324 331, 1990.
[33] M.I.J. Bloor and M.J. \\ !,..,1 I i_, partial
Aided Design, vol. 22, no. 4, pp. 202 212, 1990.
[34] G. ('. ii,! i. and W. \\ !. i, "Linear constraints for deformable Bspline surfaces," in Proceedings of the
Symposium on Interactive 3D Graphics, pp. 165 170. ACM, .'.. York, 1992.
[35] W. \\. !. i and A. Witkin, \ .I1 .Ii ...l surface modeling," in Computer Graphics Proceedings, ACM SIG
GRAPH, July 1992, Annual Conference Series, pp. 157 166.
[36] H. Qin and D. Terzopoulos, "Dynamic NURBS swung surfaces for physicsbased shape design," Computer
Aided Design, vol. 27, no. 2, pp. 111127, 1995.
[37] H. Qin and D. Terzopoulos, "DNURBS : A physicsbased framework for geometric design," IEEE .
tions on Visualization and Computer Graphics, vol. 2, no. 1, pp. 8596, March 1996.
[38] D. Terzopoulos and H. Qin, "Dynamic NURBS with geometric constraints for interactive sculpting," ACM
I .. .. ...* on Graphics, vol. 13, no. 2, pp. 103 136, Apil 1994.
[39] C. Mandal, H. Qin and B.C. Vemuri, "Dynamic smooth subdivision surfaces for data visualization," in IEEE
Visualization'97 Conference Proceedings, Phoenix,AZ, October 1997, pp. 371 377.
[40] H. Qin, C. Mandal and B.C. Vemuri, "Dynamic Catmull('!.,! I subdivision surfaces," IEEE ... .. on
Visualization and Computer Graphics, vol. 4, no. 3, September 1998, to appear.
[41] C. Mandal, B.C. Vemuri and H. Qin, 'I.ip. recovery using dynamic subdivision surfaces," in Proceedings
of the International Conference on Computer Vision, Bombay, India, January 1998, pp. 805 810.
[42] E. Bardinet, L.D. Cohen and N. Ayache, '1q,. !i.p. Il! and freeform deformations : A global model to fit
and track 3D medical data," in Lecture .,' in Computer Science, CVRMed'95, N. Ayache, Ed., Berlin,
Germany, 1995, vol. 905, pp. 319 326, Springer Verlag.
[43] L.D. Cohen and I. Cohen, I ini. , !w i !, methods for active contour models and balloons for 2D and 3D
images," IEEE .. .. . on Pattern Analysis and Machine Intelligence, vol. 15, no. 11, pp. 1131 1147,
.. . i. i 1993.
[44] E. Koh, D. Metaxas and N. Badler, 1 I i..,. I ..! shape representation using locally adaptive finite elements,"
in Lecture ...'  in computer Science, Computer Vision ECCV 1994, J.O..! .1, i.,I, Ed., Berlin, Germany,
1994, vol. 800, pp. 441 446, SpringerVerlag.
[45] F. Leitner and P. Cinquin, "Complex topology 3D objects segmentation," in Modelbased Vision Development
August 20, 1998
DRAFT
and Tools, ". !.' Proceedings, Bellingham, NA, 1991, SPIE, vol. 1609, pp. 16 26.
[46] T. McInerney and D. Terzopoulos, "A dynamic finite element surface model for segmentation and tracking
in multidimensional medical images with application to cardiac 4d image analysis," Computerized Medical
Imaging and Graphics, vol. 19, no. 1, pp. 69 83, 1995.
[47] L.H. '.iil, and J.S. Duncan, I 1..,ii ..I finding with parametrically deformable models," IEEE ... ..
on Pattern Analysis and Machine Intelligence, vol. 14, no. 11, pp. 1061 1075, 1992.
[48] D. Terzopoulos and D. Metaxas, "Dynamic 3D models with local and global deformations: Deformable
superquadrics," IEEE ... .. on Pattern Analysis and Machine Intelligence, vol. 13, no. 7, pp. 703
714, July 1991.
[49] D. Metaxas and D. Terzopoulos, '!..p'. and nonrigid motion estimation through physicsbased synthesis,"
IEEE .... .. on Pattern Analysis and Machine Intelligence, vol. 15, no. 6, pp. 580 591, June 1993.
[50] A. Pentland and B. Horowitz, I...... of nonrigid motion and structure," IEEE ... ... on Pattern
Analysis and Machine Intelligence, vol. 13, no. 7, pp. 730 742, July 1991.
[51] Y. ('!I. i and G. Medioni, '1I .. .. description of complex objects from multiple range images," in Proceedings
of the Conference on Computer Vision and Pattern Recognition, Seattle, A, June 1994, pp. 153 158.
[52] W.C. Huang and D. Goldgof, "Adaptivesize physicallybased models for nonrigid motion analysis," in
Proceedings of the Conference on Computer Vision and Pattern Recognition, Urbana('!i.. nipl,.i, IL, June
1992, pp. 833 835.
[53] H. Tanaka and F. Kishino, "Adaptive mesh generation for surface reconstruction : Parallel hierarchical
triangulation without discontinuities," in Proceedings of the Conference on Computer Vision and Pattern
Recognition, .'.. York City, NY, June 1993, pp. 88 94.
[54] M. N .ii i and D. Terzopoulos, "Adaptive meshes and shells : Irregular triangulation, discontinuities and
hierarchical subdivision," in Proceedings of the Conference on Computer Vision and Pattern Recognition,
Urbana('C!..1I,,1. .,11, IL, June 1992, pp. 829 832.
[55] J.M. Lounsbery, T. DeRose and J. \\.i, l. i, II !,! i.. !ii..,i analysis for surfaces of arbitrary topological
type," ACM i .. ....* on Graphics, vol. 16, no. 1, pp. 34 73, January 1997.
[56] B.R. Gossick, Hamilton's Principle and Physical Systems, Academic Press, .'`. York, 1967.
[57] H. Kardestuncer, !'. /..'.. Handbook, McGrawHill, .'. York, 1987.
[58] D. Terzopoulos, Ih. !..ll. '.oii .'. of inverse visual problems involving discontinuities," IEEE .. .. on
Pattern Analysis and Machine Intelligence, vol. 8, no. 4, pp. 413 424, 1986.
[59] 0. Monga and R. Deriche, "3D edge detection using recursive till. i i. in Proceedings of IEEE Conference
on Computer Vision and Pattern Recognition, June 1989, pp. 28 35.
August 20, 1998
DRAFT
[60] G.H. Golub and V.H. Van Loan, Matrix Computations, The Johns Hopkins University Press, 1989.
[61] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, .. .. Recipes in C, Cambridge Uni
versity Press, 1992.
August 20, 1998
DRAFT
