Dynamic CatmullClark Subdivision
Surfaces
Hong Qin Chhandomay Mandal Baba C. Vemuri
CISE Technical Report # TR97020
Department of Computer and Information Science and Engineering
P.O. Box 116120
University of Florida, Gainesville, Florida 32611.
Email: {qin Icmandal I vemuri}'i., .ufl.edu
Tel: (352) 3921482, (352) 3921220 (fax).
10th November, 1997
November 10, 1997
DRAFT
Abstract
Recursive subdivision schemes have been extensively used in computer graphics, computer
aided geometric design and scientific visualization for modeling smooth surfaces of arbitrary
'I" 'I . Recursive subdivision generates a visually pleasing smooth surface in the limit from an
initial userspecified i" ..1 i ,1 mesh through the repeated application of a fixed set of subdivision
rules. In this paper, we present a new dynamic surface model based on the Catmull(' I !:
subdivision scheme, which is a very popular method to model complicated objects of arbitrary
genus because of 11ii i of its nice properties. Our new dynamic surface model inherits the
attractive properties of the Catmull( l: subdivision scheme as well as that of the I'! ', 
based models. This new model provides a direct and intuitive means of manipulating geometric
shapes, a fast, robust, and hierarchical approach for recovering complex geometric shapes from
range and volume data sets using very few degrees of freedom (control vertices). We provide
an :lii 11 I, formulation and introduce the i'1! i. 1 quantities required to develop the dynamic
subdivision surface model which can be interactively deformed by ;,1, 1 ii, synthesized forces in
real time. The governing dynamic differential equation is derived using Lagrangian mechanics
and a finite element discretization. Our experiments demonstrate that this new dynamic model
has a promising future in computer graphics, geometric shape design and scientific visualization.
Keywords
Computer Graphics, CAGD, Visualization, Subdivision Surfaces, Deformable Models, Dy
namics, I !ii!I 1.I ii! I Interactive Techniques.
I. INTRODUCTION
Generating smooth surfaces of arbitrary I. 1, .1 . is a grand challenge in geometric modeling, computer
graphics and visualization. The recursive subdivision scheme first introduced by ('I ,ii:ii [1] is very well
suited for this purpose. During the past two decades, a wide i of subdivision schemes for modeling
smooth surfaces of arbitrary , .1.  have been derived in geometric modeling after ( '!I i pioneering
work on the curve generation. A recursive subdivision algorithm I i' 11 generates a smooth surface
November 10, 1997
DRAFT
which is the limit of a sequence of recursively refined polyhedral surfaces based on a userdefined initial
control mesh. At each step of the subdivision, a finer polyhedral surface with more vertices and faces will
be constructed from the previous one via a refinement process (also called i'ipp'!i < .! 'i ). In general,
subdivision schemes can be categorized into two distinct classes namely, (1) approximating subdivision
methods and (2) interpolating subdivision methods.
A. Background
Among the approximating schemes, the techniques of Doo and Sabin [2], [3], [4] and Catmull and ('I! :
[5] generalize the idea of obtaining biquadratic and bicubic Bspline patches from rectangular control
meshes. In [5], Catmull and ('I !: developed a method for recursively generating a smooth surface from
a polyhedral mesh of arbitrary 1 1 .1,_ . The Catmull('I !: subdivision surface, defined by an arbitrary
nonrectangular mesh, can be reduced to a set of standard Bspline patches except at a finite number
of extraordinary points, where the indegree of the vertex in the mesh is not equal to four. Doo and
Sabin [3] further :i!I .1 .1 the smoothness behavior of the limit surface near extraordinary points using
Fourier transforms and an eigenvalue :;i i1 i of the subdivision matrix. Ball and 1I ,*i [6], [7] and Reif
[8] further extended the prior work on !il li properties of subdivision surfaces by deriving various
necessary and itL!!i, !0 conditions on smoothness for different subdivision schemes. In [9], Loop presented
a similar subdivision scheme based on the generalization of quartic triangular Bsplines for triangular
meshes. Halstead, Kass and Derose [10] proposed an algorithm to construct a Catmull('I! !: subdivision
surface that interpolates the vertices of a mesh of arbitrary 1 I 1 _1 . In [11], Taubin developed a signal
processingbased approach to fair polyhedral surfaces of arbitrary l .1 1.._ .
The most wellknown interpolationbased subdivision scheme is the "1. 111. il algorithm proposed by
Dyn, Gregory and Levin [12]. Butterfly subdivision method 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 for the initial polygonal meshes in order to obtain a smooth limit
surface. A variant of this scheme was proposed by Dyn, Hed and Levin [13]. Recently, Zorin, Schroder
and Sweldens [14] further developed an improved interpolatory subdivision scheme that can retain the
November 10, 1997
DRAFT
ii! 1.ii '1i of the butterfly scheme and result in much smoother surfaces even from initial polygonal meshes
that are irregular.
B. Motivation
Although recursive subdivision surfaces are extremely powerful to represent smooth geometric shapes of
arbitrary I .1 .,_ , they constitute a purely geometric representation, and furthermore, conventional geo
metric modeling with subdivision surfaces 11i be infeasible 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  1, 1
splinebased modeling schemes. In addition, it 11 , 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. o' .i!.  ') ii , be strongly desired while making geometric objects satisfy global
smoothness requirements in geometric modeling and graphics applications. In contrast, 1,1i i. based
modeling provides a superior approach to shape modeling that can overcome most of the limitations
associated with traditional geometric modeling approaches. Freeform deformable models governed by
1I1! ., 1 laws are of particular interest in this context. These models respond dynamically to applied
forces in a very intuitive manner. The equilibrium state of the model is characterized by a minimum of
the potential i. of the model ,L,.j. I to imposed constraints. The potential ii. i functionals can
be formulated to satisfy local and global modeling criteria and impose geometric constraints relevant to
shape design.
Freeform deformable models were first introduced to computer graphics and visualization in Terzopou
los et al. [15] and further developed by Terzopoulos and Fleischer [16], Pentland and \\ ill !! [17],
Metaxas and Terzopoulos [18] and Vemuri and I ,l1 li. 1i [19]. Celniker and Gossard [20] developed
a * , i i for interactive freeform design based on the finite element optimization of. i I function
als proposed in [16]. Bloor and \\ !I. ., [21], [22], Celniker and Welch [23] and Welch and \\ i1!: i [24]
proposed deformable Bspline curves and surfaces which can be designed by imposing the shape crite
ria via the minimization of the i. functionals 11lj. i I to hard or soft geometric constraints through
November 10, 1997
DRAFT
Lagrange multipliers or i,. ii 11 methods. Recently, Qin and Terzopoulos '2], [26], [27] have developed
dynamic NURBS (DNURBS) which are very sophisticated models suitable for representing a wide i 1I
of freeform as well as standard :I ,! I1 I shapes. The DNURBS have the advantage of interactive and
direct manipulation of NURBS curves and surfaces, resulting in 1,1i 11 1! 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 parametric domain. Hence, it is almost impossible to model surfaces of arbitrary genus using these
models. In this paper, we develop a dynamic generalization of recursive subdivision schemes based on
Catmull( 'I i I: subdivision surfaces. Our new dynamic model combines the benefits of subdivision surfaces
for modeling arbitrary l ., .1 _ as well as the dynamic splines for direct and interactive manipulation of
shapes by ,ii, 1i i1 simulated forces. Note that, the derivation of our dynamic subdivision surface poses a
significant technical challenge because of the fact that no closedform parameterization of the limit surface
exists near the extraordinary points. We present the details of our formulation in a later section.
The dynamic Catmull('I i!: subdivision surface has been developed primarily for modeling arbitrary
,l '.1,1. However, another important application of the developed model is in shape recovery. In a
S I,'. ,1 shape reconstruction application, we need to recover shapes of arbitrary l. I 1_ from large data
sets. P! i, based models are often used for this purpose. However, the model used for fitting should
be able to recover the shape accurately. At the same time the number of degrees of freedom for model
representation should be kept low. Another important criterion is that the model initialization should
not be restricted to parameterized input meshes since it is infeasible to parameterize shapes of arbitrary
, ,! 1,,_. A 1.1! 1 based model ;rif' ;_ the aforementioned criteria is a good candidate for a solution
to the shape recovery problem.
Pr! i, based deformable models used to solve shape recovery problem involve either fixed size [19], [28],
[29], [30], [31] or adaptive size [32], [33], [34], ,;], [36], [37] grids. The models with fixed grid size generally
use less number of degrees of freedom for representation, but the accuracy of the recovered shape is lacking
in ii i cases. On the other hand, the number of degrees of freedom used for shape representation of
the model is generally very high and computationally expensive ad hoc schemes are used in models with
November 10, 1997
DRAFT
adaptive grid size methods. The recovered shape is however satisfactory in the context of accuracy. The
hierarchical shape representation using locally adaptive finite elements discussed in [34] can ti1. i0 I;
represent the shape of an object of genus zero with a small number of nodal points. However, this scheme
can not be easily extended to cope with arbitrary shapes. The balloon model for describing the shape
of complex objects [32] also adapts the mesh surface to local surface shapes and is purely driven by an
applied inflation force towards the object surface from the interior of the object. This scheme involves
a large number of nodal points for representing complex shapes. Moreover, all the existing models using
either a fixed or an adaptive grid size require a parameterized mesh as their input.
The proposed model solves the shape recovery problem very tin. wi as it can recover shapes from
large range and volume data sets using very few degrees of freedom (control vertices) for its representation
and can cope with ;~,i arbitrary input mesh, not necessarily parameterized, with an arbitrary number
of extraordinary points. The initialized model deforms under the influence of synthesized forces to fit
the data set by minimizing its ii. i Once the approximate shape is recovered, the model is further
subdivided automatically and a better approximation to the input data set is achieved using more degrees
of freedom. The process of subdivision after achieving an approximate fit is continued till a prescribed
error criteria for fitting the data points is achieved.
In a nutshell, the dynamic Catmull('I !!: subdivision surface model has been motivated by its I '111
to model arbitrary 1..11.._** where modelers can directly manipulate the smooth limit surface in an
intuitive fashion and by its applicability to the shape recovery problem.
C. Overview
The rest of the paper is organized as follows: Section II presents the detailed formulation of the dynamic
Catmull( 'I i I: subdivision Surfaces. The implementation details are provided in Section III. Experimental
results can be found in IV. Ii II ,11 !, we make concluding remarks and point out future directions of research
in Section V.
November 10, 1997
DRAFT
II. FOI:. II NATION
In this section we present a I in, i i formulation of our new dynamic model based on Catmull( 'I i!
subdivisions. 1 i i1 we ',i fi review the Catmull('I! I: subdivision scheme. Then, we demonstrate how
to assign a bicubic patch in the limit surface to a nonboundary face in a rectangular setting. We further
generalize this idea to assign the infinite number of bicubic patches in the limit surface to faces that are in
the i! ii of an extraordinary point/vertex. Next, we formulate a closed form ;,!i , 1 i, 1 representation
of the limit smooth surface which can be viewed as a function of its (initial) polyhedral control vertices.
I i i!! ,11 we introduce ,1i 1. i1 quantities into our dynamic model in order to derive its motion equation.
A. Catmull(l I subdivision .... f....
Catmull('I !: subdivision scheme, like ;,i, other subdivision scheme, starts with a userdefined mesh
of arbitrary l, ,1, .1, _ 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 follows:
For each face, introduce a new face point which is the average of all the old vertices defining the face.
For each (nonboundary) edge, introduce a new edge point which is the average of the following four
points: two old vertices defining the edge and two new face points of the faces .'li 1 ,!I to the edge.
For each (nonboundary) vertex, introduce a new face point obtained from the average F+ 2iE+ (i3
where F is the average of the new face points of all faces ,li I !i to the old vertex point, E is the
average of the midpoints of all edges incident on the old vertex and n is the number of the edges
incident on the vertex.
Form new edges by connecting each new face point to the new edge points of the edges defining the
old face and by connecting each new vertex point to the new edge points of all old edges incident on
the old vertex point.
Define new faces as those enclosed by new edges.
The most important property of Catmull(' Il: subdivision surfaces is that the smooth surface can
be generated from control meshes of arbitrary I ,1i .1. . *. Therefore, this subdivision scheme is extremely
November 10, 1997
DRAFT
valuable for modeling various complicated geometric objects of arbitrary l. 1. 1,_ . Catmull('I I: subdi
vision surfaces include standard bicubic Bspline surfaces as their special case (i.e., the limit surface is a
tensorproduct Bspline surface for a rectangular control point mesh). In addition, the aforementioned sub
division 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. Note that, after the first subdivision, all faces are quadrilaterals, hence all new vertices created
subsequently will have four incident edges. The number of extraordinary points on the surfaces remains a
constant which is determined by the refined meshes after one subdivision. The limit surface is curvature
continuous everywhere except at extraordinary vertices, where only tangent plane .i0ii,,i I is achieved.
In spite of the popularity of Catmull( 'I i!: subdivision surfaces for representing complex geometric shapes
of arbitrary I. 1 i .1. , these subdivision surfaces are not parameterizable and lack closedform i,! i1 i for
mulations. These deficiencies preclude their immediate pointwise manipulation and hence 11i i restrain
the applicability of these schemes. We develop a new dynamic model based on Catmull( 'I !: subdivision
surfaces which offer modelers a closedform I, 1 i,1 111 formulation and allows users to manipulate the model
directly and intuitively.
To develop the dynamic model which treats the limit smooth surface as a function of its control mesh in
a hierarchical fashion, we need to update control vertex positions continually at ;,i! given level. However,
all the vertices introduced through subdivision are obtained as an ;.!!i. combination of control vertex
positions of the initial mesh. Therefore, we can control the dynamic behavior of the limit surface by
formulating the dynamic model on the initial mesh itself, the only exception being the case when the
initial mesh has nonrectangular faces. This problem can be circumvented by taking the mesh obtained
through one step of subdivision as the initial mesh. To define the limit surface using the vertices of the
initial mesh, the enumeration of the bicubic patches in the limit surface is necessary. In the next two
subsections, we present a scheme of assigning the bicubic patches to various faces of the initial mesh. It
11i be noted that one additional subdivision step 11i be needed in some cases to isolate the extraordinary
points and treat the obtained mesh as the initial mesh (one I 1. ,1 example is when the initial mesh is a
tetrahedron).
November 10, 1997
DRAFT
B. Assigning patches to regular faces

(i
Fig. 1. A rectangular mesh and its limit surface consisting of 4 bicubic surface patches.
In Fig.1, a rectangular control mesh is shown along with the bicubic Bspline surface (4 patches) in
the limit after an infinite number of subdivision steps. Note that, each of the bicubic patches in the
limit surface is defined by a rectangular face with each vertex of degree four, thereby accounting for
16 control points (from its 8 connected neighborhood) needed to define a bicubic surface patch in the
limit. Therefore, for each rectangular face in the initial mesh with a valence of 4 at each vertex, the
corresponding bicubic surface patch can be assigned to it in a straight forward way. In Fig.1, the surface
patches S1, S2, S3 and S4 are assigned to face F1, F2, F3 and F4 respectively. The 16 control points for
the patch SI, corresponding to face F1, are highlighted in Fig.l.
November 10, 1997
DRAFT
S
2 \
A
SS 6
  I  
S' i S4 i
__  _  
Fig. 2. A mesh with an extraordinary point of valence 3 and its limit surface.
C. Assigning patches to irregular faces
In Fig.2, a mesh containing an extraordinary point of valence 3 and its limit surface are shown. The
faces Fo, FI, ..., F8 are assigned to bicubic patches So, S1,..., S8 respectively (as they all have vertices
of valence 4) following the aforementioned scheme. However, the central smooth surface enclosed by the
patches So, S,..., S8 consists of infinite number of bicubic patches converging to a point in the limit. We
need to develop a recursive way of enumerating these bicubic patches and assigning them to various faces
at !t.11 i. Ii levels in order to develop the dynamic subdivision surface model.
The idea of enumerating the bicubic patches corresponding to faces having an extraordinary vertex
is shown in 1 i ; where a local subdivision of the mesh consisting of faces F0, F, ..., F8, P, P1, P2 (and
not the other boundary faces) of Fig.2 is carried out. Topologically, the resulting local subdivision mesh
(shown as dotted mesh) is exactly the same as the mesh in Fig.2 and hence exactly the same number of
November 10, 1997
DRAFT
selected mesh for
local subdivision
P S
nn
n n
HR
A\
\\
n n
n n
\ \
i,''
        
nf
Fig. 3. Local subdivision around the extraordinary point and the limit surface.
bicubic patches can be assigned to its faces with vertices of valence 4 as is evident from 1 !, ; (the new faces
and the corresponding patches are marked by "p" and "n" respectively). This process of local subdivision
and assignment of bicubic patches around an extraordinary point can be carried out recursively and in
the limit, the enclosed patch corresponding to faces sharing the extraordinary point will converge to a
point. However, there is no need to carry out an infinite number of subdivision steps. This description is
for formulation purposes only and the exact implementation will be detailed in a later section.
D. Kinematics of the limit .... f..,
In this section we develop the mathematics for the kinematics of the limit surface via illustrative
examples and then present the generalized formulas. We start the illustration with a single bicubic B
November 10, 1997
DRAFT
spline patch which is obtained as the limiting process of the Catmull( 'I 1i subdivision algorithm applied
to an initial 4 by 4 rectangular control mesh. Let sp(u, v), where (u,v) E [0,1]2, denote this bicubic
Bspline patch which can be expressed;! 1 11. as
3 3
s p(u, v) = (x(u, v), y(u, v), z(u, v)) = di,jBi,4(u)Bj,4(v) (1)
i=0 j=0
where dij represents a 3dimensional position vector at the (i, j)th control point location and Bi,4 (u) and
Bj,4(v) are the cubic Bspline basis functions. The subscript p on s denotes the patch under consideration.
Expressing Eqn.1 in a generalized coordinate  . ii we have
Sp = Jpqp (2)
where Jp is the standard Jacobian matrix of a bicubic Bspline patch, and is of size (3, 48). Vector qp is the
concatenation of all control points defining a Bspline patch in 3D. Note that in the concatenation of the
control points, each control point has an (x, y, z) component. For example, the (x, y, z) components of the
control point (i, j) correspond to positions 3k, 3k + 1, 3k + 2 where, k = 4i +j respectively in the vector
qp. We can express the entries of Jp explicitly in the following way: Jp(0, k) = Jp(1, k+1) = Jp(2, k+2) =
Bi,4(u)Bj,4(v) and Jp(0, k + 1) = Jp(0, k + 2) = Jp(1, k) = Jp(1, k + 2) = Jp(2, k) = Jp(2, k + 1) = 0.
D.1 Limit surface with i ii i bicubic patches from a rectangular initial mesh
Now let's consider a limit surface consisting of i ,ii~ bicubic surface patches obtained after .'1,1 i,!
an infinite number of subdivision steps to a rectangular initial mesh. For example, let the limit surface of
Fig.1 be s, which can be written as
1 1 1 1
s,, (u, v) =s, (2,, 2v) + s, (2(u ), 2v) +s,, (2(u ), 2(v )) +s ( 2(v )) (3)
where s,, (.,, 2v) = s,,(u, v) for 0 < u, v < and 0 otherwise. Similarly, s,, S,I3 and s,4 are also equal
to s, (u, v) for an appropriate range of values of u, v and 0 outside. It i1 be noted that sM, s,,, Ms,, s, 4
correspond to patches S1, S2, S3, S4 respectively in Fig.1. Rewriting Eqn.3 in generalized coordinates we
have
4
s, = Jiqi + J2q2 +J3q3 + J4q4 = Jii (4)
i=1
November 10, 1997
DRAFT
where Jis are the Jacobian matrices of size (3,48) and qis are the (x,y,z) component concatenation of a
subset of the control points of s, defining s,,, i = 1, 2, 3 and 4. A more general expression for s,, is
4
sm = JiAlqm + J2A2qm + J3A3qm + J4A4qm = JiAiqm = Jrqm. (5)
i=1
Where, q, is the 75component vector of 3D positions of the 25 vertex control mesh defining the limit
surface s,.n Matrices Ai, 1 < i < 4, are of size (48,75) each row consisting of a single nonzero entry
(= 1) and the (3, 75)sized matrix J, = i=l JiA.
D.2 Limit surface with i! i , bicubic patches from an arbitrary initial mesh
The stage is now set to define the limit surface s using the vertices of initial mesh 4M for i~ arbitrary
1..**I' ,!_, assuming all faces are rectangular and no face contains more than one extraordinary point as
its vertex (i.e., extraordinary points are isolated). As mentioned earlier, if these assumptions are not
satisfied, one or two steps of global subdivision 11i be required and the resulting mesh can be treated
as the initial mesh. Let the number of vertices in the initial mesh 4M be a, and let I of these be the
extraordinary vertices. Let us assume that the number of faces in the initial mesh are b, and that k of
these have vertices with valence 4 (henceforth termed a "ii.. i I I f ) and each of the remaining (b k)
faces have one of the I extraordinary vertices (henceforth termed a "1 i '! i ). Let p be the 3a = N
dimensional vector containing the control vertex positions in 3D. Using the formulations in subsections
IIB and IIC, the smooth limit surface can be expressed as
k I
s= n + S (6)
i=1 j=l
where ni is a single bicubic patch assigned to each of the normal faces and sj is a collection of infinite num
ber of bicubic patches corresponding to each of the extraordinary points. "',i '1~,!'! ; the same approach
taken before to derive Eqn.5, it can be shown that
k k k
>ni = E("Ji)("pi)= ( '("Ji)("Ai))p = ("J)p (7)
i=1 i=1 i=1
where "Ji,"pi and "Ai are the equivalent of Ji,pi in Eqn.4 and Ai in Eqn.5 respectively. The pre
superscript n is used to indicate that these mathematical quantities describe bicubic patch in the limit
surface corresponding to normal faces.
We will use the following notational convention for describing various mathematical quantities used
November 10, 1997
DRAFT
in the derivation of the expression for a collection of infinite number of bicubic patches around an ex
traordinary vertex. The presuperscript s is used to represent a collection of bicubic patches around an
extraordinary vertex, the subscript j is used to indicate the jth extraordinary point, the postsuperscript
represents the exponent of a mathematical !I i0 I and the level indicator (to represent various levels of
subdivision in the local control mesh around an extraordinary vertex) is depicted via subscripts on the
curly braces.
The expression for sj is derived using the recursive nature of local subdivision around an extraordinary
vertex as shown in subsection IIC. i 1, sj can be expressed as
sj = {'Jj}1{SP}I + {sj} (8)
where the first term of Eqn.8 is the generalized coordinate representation of the bicubic Bspline patches
corresponding to the normal faces of the new local subdivision mesh obtained after one subdivision step
on the local control mesh (similar to those patches marked n in Fig.3). {sjy} represents the rest of the
infinite bicubic Bspline patches surrounding the extraordinary point (similar to the central patch enclosed
by patches marked n in Fig.3). The vertices in the newly obtained local subdivision mesh {"pj } can be
expressed as a linear combination of a subset of the vertices of the initial mesh M (which will contribute to
the local subdivision) following the subdivision rules. We can name this subset of initial control vertices
{ PJ}0. Furthermore, there exists a matrix {"Bj}1 of size (3c,3d), such that {8Bj5}{5pj,} = {5pj}i
where {Spj}l and {fpj}o are vectors of dimension 3c and 3d respectively. AI.I1. i,_ the idea of recursive
local subdivision again on {sy},, sy can be further expanded as
Si = {SJj}l{SBj}j{spj}o + {'SJj}2{SB}2 }1 + {sj}2 (9)
In the above derivation, {'y I}, is a vector of dimension 3d, comprising of a subset of the vertices defining
the 3c dimensional vector {5py}1. Note that, {fpj}j has the same structure as {pjjY}, therefore, there
exists a (3d, 3d) matrix {"Cjy} such that {'Cjy 1} pjy = {SP} 1l. Each subdivision of a local mesh with
d vertices creates a new local mesh with c vertices which contributes a fixed number of bicubic Bspline
patches. So, if we proceed one step further, we obtain
s e = {SJth iSBt}i spj}o + fJj2 t 2e cl r 'J B 2 fCu p s 0 + 3' arJ B Cd t3 e oI {s} (10)j
Because of the intrinsic property of the local recursive subdivision around the extraordinary point, we
November 10, 1997
DRAFT
have {SJj = {,= {J = }2 = {Jj} = ... = {'Jj} In addition, the subdivision rules remain the same
throughout the refinement process, we also have {"Bj}1 = {"Bj}2 ... = {"B},2 = ...= { B3j}. So,
we can further simplify the above equations leading to
sj = {SJ ll SIBll JSpj}o + SJll S }Jll SCJl Sp o + SJll S } JSC l p }o ...
Oo
i= 0
0 ,,
{scj}1 {s 0
i I fpo
][ 12
s
fScjI f}2 s?}
{ ][3}3
Fig. 4. Local subdivision around the extraordinary point and the corresponding patches in the limit
surface from I,!!. i. 1 levels of subdivision.
November 10, 1997
DRAFT
We can rewrite sj as
si = (SJj)(SPj) (12)
where SJj = {SJj }{SB }I,(' {sCj/}) and "pj = {SP/ }0* The idea of local recursive subdivision
around an extraordinary point is illustrated in Fig.4. Note that, each vertex position in the subdivided
mesh is obtained by an ;~!hil combination of some vertices in the previous level and hence ;~i row of
{"Cj } sums to 1. The largest eigenvalue of such a matrix is 1 and it can be shown that the corresponding
infinite series is convergent following a similar approach as in [10]. The rest of the derivation leading to
an expression for s is relatively straight forward. Using the same approach used to derive the Eqn.7, it
can be shown that
sj = (Jj)(Spj) = ( (SJj)(SAj))p = (SJ)p (13)
j=1 j=1 j=1
From Eqn.6,7 and 13,
s = ("J)p + (SJ)p (14)
Let J = ("J) + ("J), hence
s = Jp (15)
E. D ,,.. .. ',.
We now treat the control point positions (alternatively, the vertex positions in the initial mesh) defining
the limit surface s as a function of time in order to develop our new dynamic model. The velocity of the
surface model can be expressed as
s(u, v,p) = Jp (16)
where an overstruck dot denotes a time derivative. The 1 11! i. of the dynamic subdivision surface model
is based on the 1 i _ version of Lagrangian dynamics [38] and is formulated in an analogous way
to that in [27].
In an abstract 1,1i 1  in 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. The Lagrangian equation of motion can then be expressed as
Mpi + Dp + Kp = fp (17)
November 10, 1997
DRAFT
Let pt(u, v) be the mass 1. i!i function of the surface. Then
M= ff IjJdudv (18)
is an N x N mass matrix. Similarly the expression for damping matrix is
D= JJJ Jdudv (19)
where (u, v) is the damping 1L i!.iI
A thinplateundertension i model [39] is used to compute the elastic potential !, i of the
dynamic subdivision surface. The corresponding expression for the stiffness matrix K is
K = (a nJJu + 22 v + 11J JJu +12 JTu 22J vJ)ddudv (20)
where the subscripts on J denote the parametric partial derivatives. The aii(u,v) and ij (u, v)s are
S1 I i. i functions controlling local tension and i_ii I in the two parametric coordinate directions. The
generalized force vector fp can be obtained through the principle of virtual work [38] done by the applied
force distribution f(u, v, t) and can be expressed as
fp= Tf(u, v, t)dudv (21)
E.1 Multilevel Dynamics
Our dynamic Catmull( 'I I: surface model can be subdivided globally to increase the number of vertices
(control points) of the model. For example, after one step of global subdivision, the initial degrees of
freedom p (refer to Eqn.15 and Eqn.16) in the dynamic  I. I, will be replaced by a larger number of
degrees of freedom q, where q = Ap. A is a global subdivision matrix of size (M, N) whose entries are
uniquely determined by Catmull('I I: subdivision rules (see Section IIA for the details about the rules).
Thus, p, expressed as a function of q, can be written as
p = (ATA) ATq = Bq (22)
where B = (ATA) 'AT. Therefore, we can rewrite Eqn.15 and Eqn.16 as
s = (JB)q (23)
and
(u, v, q) = (JB) (24)
November 10, 1997
DRAFT
respectively. Now we need to 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 l i1 ,. 1!i..i il 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
Mqi + Dq + Kqq = fq (25)
where Mq = f /pBTJTJBdudv and the derivation of Dq, Kq and f, follow suit.
It ii , be noted that further subdivision, if necessary, can be carried out in a similar fashion. Therefore,
multilevel dynamics is achieved through recursive subdivision on the initial set of control vertices. Users
can interactively choose the level of detail representation of the dynamic model as appropriate for their
modeling and design requirements. Alternatively, the  I. i. can automatically determine the level of
subdivision most suitable for an application depending on some applicationspecific criteria.
III. FINITE 1.1.1.. I.NT I. IP1.1.. I.NTATION
The evolution of the generalized coordinates for our new dynamic surface model can be determined
by the secondorder differential equation as given by Eqn.17. An iii1 1 I, G1 solution of the governing
differential equation can not be obtained in general. However, an. mtn. ii numerical implementation can
be obtained using finite element ; !! ,1 1 techniques [40]. For the dynamic subdivision surface model, two
I of finite elements are considered normal elements (bicubic patches assigned to the normal faces of
the initial mesh) and special elements (collection of infinite number of bicubic patches assigned to each
extraordinary vertex of the initial mesh). In the current implementation, the M, D and K matrices for
each individual normal and special elements are calculated and they can be assembled into the global M,
D and K matrices that appear in the corresponding discrete equation of motion. In practice, we never
assemble the global matrices explicitly in the interest of time performance. The detailed implementation
is explained in the following subsections.
November 10, 1997
DRAFT
A. Data S, '. ,.,. .*
A subdivision surface defined by a control mesh at ;,i level is designed as a class which has a pointer
to its parent mesh, a set of pointers to its offspring meshes (arising out of local subdivision around the
extraordinary vertices at that level), a list of faces, edges, vertices and normal elements Face, edge,
vertex and normal elements are, in turn, classes which store all the connectivity and other information
needed to either enumerate all the patches or locally subdivide around an extraordinary vertex in that
level. The implementation takes the initial mesh as the base subdivision surface object (with its parent
pointer set to NULL) and locally subdivides the initial mesh upto a userdefined maximum level around
each extraordinary vertex to create offspring objects at .l !. i. !!I levels. At this point, let's take a closer
look at the normal and special element data structures and computation of the corresponding local M, D
and K matrices.
A.1 Normal i.. !,. !i 
Each normal element is a bicubic surface patch and is hence defined by 16 vertices (from the 8connected
neighborhood of the corresponding normal face). Each normal element keeps a set of pointers to those
vertices of the initial mesh which act as control points for the given element. For a normal element, the
mass, damping and stiffness matrices are of size (16,16) and can be computed exactly by i!! i, out
the necessary integration .1 ,! II i.11 The matrix J in Eqn.18,19 and 20 needs to be replaced by Jp
(of Eqn.2) for computation of the local M, D and K matrices respectively of the corresponding normal
element.
A.2 Special i. !i 
Each special element consists of an infinite number of bicubic patches in the limit. We have already
described a recursive enumeration of the bicubic patches of a special element in Section IIC. Let us now
consider an arbitrary bicubic patch of the special element in some level j. The mass matrix M, of this
patch can be written as
M = Nf I !0 (26)
November 10, 1997
DRAFT
where Mp is the normal element mass matrix (scaled by a factor of to take into account of the area
shrinkage in bicubic patches at higher level of subdivision) and QS is the transformation matrix of the
control points of that arbitrary patch from the corresponding control points in the initial mesh. The
damping and stiffness matrices for the given bicubic patch can be derived in an exactly similar fashion.
Now, these mass, damping and stiffness matrices can be assembled to form the mass, damping and
stiffness matrices of the special element. As mentioned in Section IID.2, the infinite series summation is
convergent. However, it has been found that the contribution from bicubic patches of a special element
at a higher level of subdivision to the mass, damping and stiffness matrices becomes negligible and in the
implementation, the local subdivision is carried out until the contribution is small enough to be ignored.
B. Force Application
The force f(u, v, t) in Eqn.21 represents the net effect of all applied forces. The current implementation
supports spring, inflation as well as imagebased forces. However, other 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 (uo, Vo) on
the limit surface, the net applied spring force being
f(u, v, t) = k(do s(u, v, t))(u ,o,v vo)dudv (27)
where 6 is the unit impulse function i1 i11 f(uo, v, t) = k(do s(uo, v0, t)) and vanishes elsewhere in
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 a mouse button or
the points from which forces need to be applied can be read in from a t!.
To recover shapes from 3D image data, we synthesize imagebased forces. A 3D edge detection is
performed on a Gaussian smoothed volume data set using the 3D MongaD(, i !,. .II)) operator [41] 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
VP(x, y, z)
f(x,y,z)= k P (28)
H VP( x, y, z) z
where k controls the strength of the force. The applied force on each element is computed using Gaussian
November 10, 1997
DRAFT
quadrature for evaluating Eqn.21 in Cartesian coordinates. It !i be noted that we can apply spring
forces in addition with the imagebased forces by placing points near the region of interest in the slices of
the 3D image data.
C. Discrete D., ...... Equation
The 1!l, i i!!ii 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
where discrete derivatives of p are calculated using
p(t + At) 2p(t) + p(t At) (29)
j(t+ At) = A2()
At2
and
(t At) = p(t + At) p(t At) (30)
p^(+ At) = (30)
2At
Using Eqn.17,29 and 30, the discrete equation of motion is obtained as
(2M + DAt + 2At2K)p(t + At) = 2At2fp(t + At) + (DAt 2M)p(t At) + 4Mp(t) (31)
This linear I. i 1 of equations is solved iteratively between each time step using the ..I il il. gradient
method. For a first order  1. i, with no mass, the above equation reduces to
(D + 2AtK)p(t + At) = 2Atfp(t + At) + Dp(t At) (32)
which gives a faster convergence.
D. Model Subdivision
The initialized model grows dynamically according to the equation of motion (Eqn.17) and when an
equilibrium is achieved at a given level of subdivision, the model can be subdivided, if necessary, according
to the Catmull('I I: subdivision rules to increase the number of vertices (control points) and a better fit
to the data can be achieved. Currently the error of fit criteria is based on distance between the data points
and the points on the limit surface where the corresponding springs are attached. However, other I i" 
of error criterion can also be defined and used in this context. For example, in the context of imagebased
forces, if the model. i does not change between successive iterations indicating an equilibrium for the
given resolution, the model can be subdivided further until the model !!. i _ is it!l, !. !,I I1 small and the
November 10, 1997
DRAFT
change in i i between successive iterations becomes less than a prespecified tolerance.
IV. RESULTS
The proposed dynamic subdivision surface can be used to represent a wide ,i. of shapes with
arbitrary genus. In this section we demonstrate the power of our modeling scheme via model fitting
examples to a i i. of data sets of ii, degree of,! !.ii1 .i In all the experiments, normal elements
are shaded yellow, while special elements are colored green.
In i " .) an open limit surface defined by an initial mesh of 61 vertices and 45 faces is shown. The
mesh has one extraordinary point of valence 5. The limit surface is acted upon by spring forces as shown
in I 'i, .(b). The evolving model and its control mesh is shown in Fig.5(c) and(d). The final fitted model
is depicted in Fig.5(e) and (f). It 1i i be noted that the model controlled by the initial mesh reached local
minimum without fitting the points exactly. In order to obtain an exact fit ( i, '(f)), the control mesh
is subdivided once thereby increasing the degrees of freedom (control vertices) of the ,11 L!. 1 i! model.
Thus the dynamics can be applied in a hierarchical fashion. The developed model can be used to obtain a
very fast approximate fitting with fewer number of vertices and an exact fit after more subdivision steps
as needed.
In the next experiment, we show the fitting process using spring forces with a closed surface of genus
two(Fig.6). The smooth surface is controlled by an initial mesh of 544 faces and 542 vertices, 8 of them
being extraordinary points of valence 5. In this experiment, the model has it!, !. !!, degrees of freedom
and fitted the data points exactly without needing further subdivision of its control mesh.
In all the experiments to follow, the initialized model had 96 faces and 98 vertices, 8 of them being
extraordinary vertices of valence 3. The final fitted model, obtained through one step of subdivision, has
a control polygon of 384 faces with 386 vertices. The tolerance level of the error in fit, which is defined
as the maximum distance between a data point and the nearest point on surface as a percentage of the
object diameter, was set to be 1 .
In Fig.7, we demonstrate the model fitting algorithm applied to laser range data acquired from multiple
views of a light bulb. Prior to ,,I'1 i,, our algorithm, the data were transformed into a single reference
November 10, 1997
DRAFT
coordinate  . i!, The model was initialized inside the 1000 range data points on the surface of the bulb.
In the next experiment, the shape of a human head is recovered from a range data set as shown in
1i, ', The range data set has 1779 points in 3D. It ii 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 data
points present in the original range data set. I II!!, example with an anvil data set is shown in Fig.9.
The anvil data set has 2031 data points.
We show the application of our model to anatomical shape recovery from 3D volumetric 1, li1 data in
the last two experiments. 1 i we fit the model to a cerebellum (a cortical structure in brain) given an
input of 30 sagittal slices from a I: 1 brain scan. Fig.10(a) depicts a slice from this 1: i1 scan and the
model initialization is shown in Fig.10(b). Continuous image based forces are applied to the model and
the model deforms under the influence of these forces until maximum conformation to the boundaries of
the desired cerebellum shape. Fig.10(c) depicts an intermediate stage of the model evolution during the
fitting process and the final fitted model is shown in Fig.10(d). Arbitrary 3D views of the fitted model
from different viewing angles are depicted in Fig.10(e) and (f).
In the last experiment, we present the shape extraction of a caudate nucleus (another cortical structure
in human brain) from 64 I l1 I slices, each of size (256, 256). Fig.11(a) depicts a slice from this I I I scan
along with the points placed by an expert neuroscientist on the boundary of the shape of interest. Fig.11(b)
depicts the data points (placed in each of the slices depicting the boundary of the shape of interest) in 3D.
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. Fig.11(c) depicts the initialized
model and the data points. 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 maximum conformation to the
boundaries of the desired caudate shape. Fig.11(d) depicts an intermediate stage of the model evolution
during the fitting process and two arbitrary views of the final fitted model in 3D is shown in Fig.11(e)
and (f).
November 10, 1997
DRAFT
V. CONCLUSIONS
In this paper, a dynamic generalization of the Catmull('I l: subdivision surfaces is presented which
has numerous applications in geometric modeling, computer graphics and visualization. Apart from
providing a direct and intuitive way of manipulating shapes, it facilitates the modeling and shape ;i! I1 i
of objects contained in range and volume data sets using very few degrees of freedom. We have presented
an ;li 1 11, formulation of the subdivision scheme, incorporated the advantages of freeform deformable
models in subdivision scheme, introduced hierarchical dynamic control and shown the advantages of our
model via experiments. However, the current scheme can not recover very sharp edges in the data. Also,
the initialization is interactive; ideally, initialization should be done automatically on the basis of the
input data set. Our future efforts will be focused toward addressing these issues.
VI. ACKNOWLEDG(I.. !1.NTS
This research was supported in part by the NSF grant ECS9210648 and the NIH grant RO1LM05944
to BCV, the NSF ('.\il:.i.i award CCR9702103 and DMI9700129 to HQ. We also wish to acknowledge
Dr. H. Hoppe and Dr. K. Pulli for the range data and Dr. C.M. Leonard for the brain 11I data.
REFERENCES
[1] G.M. ('.ii.! iii "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] 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.
[4] M. Sabin, !' use of piecewise forms for the numerical representation of shape, Ph.D. thesis, Hungarian
Academy of Sciences, Budapest, 1976.
[5] E. Catmull and J. ('l.! I .... 1 I generated Bspline surfaces on arbitrary topological meshes," Computer
Aided Design, vol. 10, no. 6, pp. 350 355, 1978.
[6] A.A. Ball and D.J.T. '..ii", "Conditions for tangent plane continuity over recursively generated Bspline
surfaces," ACM I ... .. on Graphics, vol. 7, no. 2, pp. 83 102, 1988.
November 10, 1997
DRAFT
[7] 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.
[8] U. Reif, "A unified approach to subdivision algorithms near extraordinary points," Computer Aided Geometric
Design, vol. 12 no. 2, pp. 153 174, 1995.
[9] C. Loop, !ii ...i subdivision surfaces based on triangles," M.S. thesis, University of Utah, Department of
Mathematics, 1987.
[10] M. Halstead, M. Kass and T. DeRose, I.1 t. ii fair interpolation using Catmull('!.Cl surfaces," in Computer
Graphics Proceedings. ACM SIGGRAPH, August 1993, Annual Conference Series, pp. 35 44.
[11] G. Taubin, "A signal processing approach to fair surface design," in Computer Graphics Proceedings. ACM
SIGGRAPH, August 1995, Annual Conference Series, pp. 351 358.
[12] 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, .Apl 1990.
[13] N. Dyn, S. Hed and D. Levin, 'I I.I i..hi schemes for surface interpolation," in Workshop in Computational
Geometry, A. C. et al., Ed., pp. 97 118. 1993.
[14] D. Zorin, P. Schr6der and W. Sweldens, i!i. ii,..!.ii ii subdivision for meshes with arbitrary topology," in
Computer Graphics Proceedings. ACM SIGGRAPH, August 1996, Annual Conference Series.
[15] D.Terzopoulos, J.Platt, A.Barr and K.Fleischer, !.1. l...11r deformable models," in Computer Graphics
Proceedings. ACM SIGGRAPH, 1987, Annual Conference Series, pp. 205 214.
[16] D. Terzopoulos and K. Fleischer, "Deformable models," '. Visual Computer, vol. 4, no. 6, pp. 306 331,
1988.
[17] A. Pentland and J. Williams, "Good l.II Mt..I : Modal dynamics for graphics and animation," in Computer
Graphics Proceedings. ACM SIGGRAPH, 1989, Annual Conference Series, pp. 215 222.
[18] 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.
[19] B. C. Vemuri and A. Radisavljevic, I. lii! i i ..i. ..i stochastic hybrid shape models with fractal priors,"
ACM I ... .. on Graphics, vol. 13, no. 2, pp. 177 207, Apil 1994.
[20] 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.
[21] M.I.J. Bloor and M.J. \\ !..!i, I,. p. , m.ii i PDE surfaces in terms of Bsplines," Computer Aided Design,
vol. 22, no. 6, pp. 324 331, 1990.
[22] M.I.J. Bloor and M.J. \\ !,..i i I ,1 partial
Aided Design, vol. 22, no. 4, pp. 202 212, 1990.
November 10, 1997
DRAFT
[23] G. ('. i!. ii and W. \\. !1. i "Linear constraints for deformable Bspline surfaces," in Proceedings of the
Symposium on Interactive 3D Graphics, pp. 165 170. ACM, .'.. York, 1992.
[24] W. \\. !I. i and A. Witkin, \ ..i.I ..i.!l surface modeling," in Computer Graphics Proceedings. ACM SIG
GRAPH, July 1992, Annual Conference Series, pp. 157 166.
[25] H. Qin and D. Terzopoulos, "Dynamic NURBS swung surfaces for physicsbased shape design," Computer
Aided Design, vol. 27, no. 2, pp. 111127, 1995.
[26] 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.
[27] D. Terzopoulos and H. Qin, "Dynamic NURBS with geometric constraints for interactive sculpting," ACM
S.. .. .* on Graphics, vol. 13, no. 2, pp. 103 136, Apil 1994.
[28] L.D. Cohen and I. Cohen, I ini. !i. 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.
[29] D. Metaxas and D. Terzopoulos, "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.
[30] D. Metaxas and D. Terzopoulos, !i.,i., and nonrigid motion estimation through physicsbased synthesis,"
IEEE ... .. on Pattern Analysis and Machine Intelligence, vol. 15, no. 6, pp. 580 591, June 1993.
[31] 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.
[32] Y. ('!I. 1 and G. Medioni, 1i! !... description of complex objects from multiple range images," in Proceedings
of the Conference on Computer Vision and Pattern Recognition, Seattle,NVA, June 1994, pp. 153 158.
[33] 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. I !lp.! _i, IL, June
1992, pp. 833 835.
[34] E. Koh, D. Metaxas and N. Badler, 1L. i.h I!, ..1 shape representation using locally adaptive finite elements,"
in Lecture ... in computer Science, Computer Vision ECCV 1994, J.O..! !i.EIIil, Ed., Berlin, Germany,
1994, vol. 800, pp. 441 446, SpringerVerlag.
[35] T. McInerney and D. Terzopoulos, "A finite element model for 3d shape reconstruction and nonrigid motion
tracking," in Proceedings of the International Conference on Computer Vision, Berlin, Germany, May 1993,
pp. 518 523.
[36] H. Tanaka and F. Kishino, "Adaptive mesh generation for surface reconstruction : Parallel hierarchical
November 10, 1997
DRAFT
triangulation without discontinuities," in Proceedings of the Conference on Computer Vision and Pattern
Recognition, .'.. York City, NY, June 1993, pp. 88 94.
[37] M1. N .! . 11 and D. Terzopoulos, "Adaptive meshes and shells : Irregular triangulation, discontinuities and
hierarchical subdivision," in Proceedings of the Conference on Computer Vision and Pattern Recognition,
UrbanaC'!h..ii,..ij,!_ IL, June 1992, pp. 829 832.
[38] B.R. Gossick, Hamilton's Principle and Physical Systems, Academic Press,.'.. York, 1967.
[39] D. Terzopoulos, I i.i. i '.ii .1. of inverse visual problems involving discontinuities," IEEE .. .. on
Pattern Analysis and Machine Intelligence, vol. 8, no. 4, pp. 413 424, 1986.
[40] H. Kardestuncer, '. L..'.  Handbook, McGrawHill, .'. York, 1987.
[41] 0. Monga and R. Deriche, "3d edge detection using recursive ti I. i ii. in Proceedings of IEEE Conference
on Computer Vision and Pattern Recognition, June 1989, pp. 28 35.
November 10, 1997
DRAFT
(a) (b)
(c) (d)
(e) (f)
Fig. 5. 1 !I ii, the dynamic open surface model to discrete points in 3D : (a) model initialization depicting
the associated control mesh, (b) model and data points, (c) & (d) intermediate stages of the fitting,
(e) final fit depicting the control mesh and (f) fitted model without control mesh.
November 10, 1997
DRAFT
(a) (b)
(c) (d)
(e) (f)
Fig. 6. 1 !,II1, the dynamic closed surface model to discrete points in 3D : (a) model initialization
depicting the associated control mesh, (b) model and data points, (c) & (d) intermediate stages of
the fitting, (e) final fit depicting the control mesh and (f) fitted model without control mesh.
November 10, 1997
DRAFT
(a) (b)
(c) (d)
(e) (f)
Fig. 7. (a) Range data of a bulb, (b)initialized model and the associated control polygon, (c) data and
superimposed initialized model, (d) intermediate stage of evolution, (e) fitted model and (f) fitted
model with its control polygon.
November 10, 1997
DRAFT
(a) (b)
(c) (d)
(e) (f)
Fig. 8. (a) Range data of a head, (b)initialized model and the associated control i.'.1 _.,i (c) data and
initialized model, (d) intermediate stage of evolution, (e) fitted model and (f) fitted model with its
control I'*1 **. i
November 10, 1997
DRAFT
(a) (b)
(c) (d)
(e) (f)
Fig. 9. (a) Range data of an anvil, (b)initialized model and the associated control i" .1 _ *i (c) data and
initialized model, (d) intermediate stage of evolution, (e) fitted model and (f) fitted model with its
control I,1, _ i
November 10, 1997
DRAFT
(c) (d)
(e) (f)
Fig. 10. (a) A slice from a brain 11i 1i (b) initialized model inside the region of interest superimposed on
the slice, (c) intermediate stage of the evolving model, (d) fitted model, (e) & (f) arbitrary 3D views
of the model fitted to the cerebellum.
November 10, 1997
DRAFT
(a) (b)
(C) (d)
(e) (f)
Fig. 11. (a) Data points ;.1. i li ii the boundary of the region of interest (a caudate nucleus) on a liI
slice of human brain, (b) data points (from all the slices) in 3D, (c) data and initialized model, (d)
intermediate stage of evolution, (e) fitted model and (f) another view of the fitted model.
November 10, 1997
DRAFT
