Dynamic Catmull-Clark Subdivision Surfaces.
AbstractRecursive subdivision schemes have been extensively used in computer graphics, computer-aided geometric design, and scientific visualization for modeling smooth surfaces of arbitrary topology. Recursive subdivision generates a visually pleasing smooth surface in the limit from an initial user-specified polygonal 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-Clark subdivision scheme, a popular technique for modeling complicated objects of arbitrary genus. Our new dynamic surface model inherits the attractive properties of the Catmull-Clark subdivision scheme, as well as those of the physics-based models. This new model provides a direct and intuitive means of manipulating geometric shapes, and an efficient hierarchical approach for recovering complex shapes from large range and volume data sets using very few degrees of freedom (control vertices). We provide an analytic formulation and introduce the "physical" quantities required to develop the dynamic subdivision surface model which can be interactively deformed by applying synthesized forces. The governing dynamic differential equation is derived using Lagrangian mechanics and the finite element method. Our experiments demonstrate that this new dynamic model has a promising future in computer graphics, geometric shape design, and scientific visualization.
INTRODUCTION
Generating smooth surfaces of arbitrary topology is a grand challenge in geometric modeling, computer
graphics and visualization. The recursive subdivision scheme first introduced by Chaikin [1] is very well
suited for this purpose. During the past two decades, a wide variety of subdivision schemes for modeling
smooth surfaces of arbitrary topology have been derived in geometric modeling after Chaikin's pioneering
work on the curve generation. A recursive subdivision algorithm typically generates a smooth surface
November 10, 1997
which is the limit of a sequence of recursively refined polyhedral surfaces based on a user-defined 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 "chopping corners"). 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 Clark
[5] generalize the idea of obtaining biquadratic and bicubic B-spline patches from rectangular control
meshes. In [5], Catmull and Clark developed a method for recursively generating a smooth surface from
a polyhedral mesh of arbitrary topology. The Catmull-Clark subdivision surface, defined by an arbitrary
non-rectangular mesh, can be reduced to a set of standard B-spline patches except at a finite number
of extraordinary points, where the in-degree of the vertex in the mesh is not equal to four. Doo and
Sabin [3] further analyzed the smoothness behavior of the limit surface near extraordinary points using
Fourier transforms and an eigenvalue analysis of the subdivision matrix. Ball and Storry [6], [7] and Reif
[8] further extended the prior work on continuity properties of subdivision surfaces by deriving various
necessary and sufficient conditions on smoothness for different subdivision schemes. In [9], Loop presented
a similar subdivision scheme based on the generalization of quartic triangular B-splines for triangular
meshes. Halstead, Kass and Derose [10] proposed an algorithm to construct a Catmull-Clark subdivision
surface that interpolates the vertices of a mesh of arbitrary topology. In [11], Taubin developed a signal
processing-based approach to fair polyhedral surfaces of arbitrary topology.
The most well-known interpolation-based subdivision scheme is the "butterfly" 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
simplicity 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 topology, they constitute a purely geometric representation, and furthermore, conventional geometric
modeling with subdivision surfaces may be infeasible for representing highly complicated objects.
For example, modelers are faced with the tedium of indirect shape modification and refinement through
time-consuming operations on a large number of (most often irregular) control vertices when using typical
spline-based modeling schemes. In addition, it may not be enough to obtain the most "fair" surface that
interpolates a set of (ordered or unorganized) data points. A certain number of local features such as
bulges or inflections ("roughness") may be strongly desired while making geometric objects satisfy global
smoothness requirements in geometric modeling and graphics applications. In contrast, physics-based
modeling provides a superior approach to shape modeling that can overcome most of the limitations
associated with traditional geometric modeling approaches. Free-form deformable models governed by
physical 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 energy of the model subject to imposed constraints. The potential energy functionals can
be formulated to satisfy local and global modeling criteria and impose geometric constraints relevant to
shape design.
Free-form deformable models were first introduced to computer graphics and visualization in Terzopoulos
et al. [15] and further developed by Terzopoulos and Fleischer [16], Pentland and Williams [17],
Metaxas and Terzopoulos [18] and Vemuri and Radisavljevic [19]. Celniker and Gossard [20] developed
a system for interactive free-form design based on the finite element optimization of energy functionals
proposed in [16]. Bloor and Wilson [21], [22], Celniker and Welch [23] and Welch and Witkin [24]
proposed deformable B-spline curves and surfaces which can be designed by imposing the shape criteria
via the minimization of the energy functionals subject to hard or soft geometric constraints through
November 10, 1997
Lagrange multipliers or penalty methods. Recently, Qin and Terzopoulos [25], [26], [27] have developed
dynamic NURBS (D-NURBS) which are very sophisticated models suitable for representing a wide variety
of free-form as well as standard analytic shapes. The D-NURBS have the advantage of interactive and
direct manipulation of NURBS curves and surfaces, resulting in physically meaningful hence intuitively
predictable motion and shape variation.
A severe limitation of the existing deformable models, including D-NURBS, 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-Clark subdivision surfaces. Our new dynamic model combines the benefits of subdivision surfaces
for modeling arbitrary topology as well as the dynamic splines for direct and interactive manipulation of
shapes by applying simulated forces. Note that, the derivation of our dynamic subdivision surface poses a
significant technical challenge because of the fact that no closed-form parameterization of the limit surface
exists near the extraordinary points. We present the details of our formulation in a later section.
The dynamic Catmull-Clark subdivision surface has been developed primarily for modeling arbitrary
topology. However, another important application of the developed model is in shape recovery. In a
typical shape reconstruction application, we need to recover shapes of arbitrary topology from large data
sets. Physics-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
topology. A physics-based model satisfying the aforementioned criteria is a good candidate for a solution
to the shape recovery problem.
Physics-based deformable models used to solve shape recovery problem involve either fixed size [19], [28],
[29], [30], [31] or adaptive size [32], [33], [34], [35], [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 many 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
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 efficiently
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 efficiently as it can recover shapes from
large range and volume data sets using very few degrees of freedom (control vertices) for its representation
and can cope with any arbitrary input mesh, not necessarily parameterized, 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 energy. Once the approximate shape is recovered, the model is further
subdivided automatically and a better approximation to the input data set is achieved using more degrees
of freedom. The process of subdivision after achieving an approximate fit is continued till a prescribed
error criteria for fitting the data points is achieved.
In a nutshell, the dynamic Catmull-Clark subdivision surface model has been motivated by its capability
to model arbitrary topology 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-Clark subdivision Surfaces. The implementation details are provided in Section III. Experimental
results can be found in IV. Finally, we make concluding remarks and point out future directions of research
in Section V.
II. FORMULATION
In this section we present a systematic formulation of our new dynamic model based on Catmull-Clark
subdivisions. First, we briefly review the Catmull-Clark subdivision scheme. Then, we demonstrate how
to assign a bicubic patch in the limit surface to a non-boundary 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 vicinity of an extraordinary point/vertex. Next, we formulate a closed form analytical representation
of the limit smooth surface which can be viewed as a function of its (initial) polyhedral control vertices.
Finally, we introduce physical quantities into our dynamic model in order to derive its motion equation.
A. Catmull-Clark subdivision surfaces
Catmull-Clark subdivision scheme, like any other subdivision scheme, starts with a user-defined mesh
of arbitrary topology. It refines the initial mesh by adding new vertices, edges and faces with each step
of subdivision following a fixed set of subdivision rules. In the limit, a sequence of recursively refined
polyhedral meshes will converge to a smooth surface. The subdivision rules are as follows:
ffl For each face, introduce a new face point which is the average of all the old vertices defining the face.
ffl For each (non-boundary) 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 adjacent to the edge.
ffl For each (non-boundary) vertex, introduce a new face point obtained from the average F
where F is the average of the new face points of all faces adjacent 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.
ffl 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.
ffl Define new faces as those enclosed by new edges.
The most important property of Catmull-Clark subdivision surfaces is that the smooth surface can
be generated from control meshes of arbitrary topology. Therefore, this subdivision scheme is extremely
November 10, 1997
valuable for modeling various complicated geometric objects of arbitrary topology. Catmull-Clark subdivision
surfaces include standard bicubic B-spline surfaces as their special case (i.e., the limit surface is a
tensor-product B-spline surface for a rectangular control point mesh). In addition, the aforementioned sub-division
rules generalize the recursive bicubic B-spline patch subdivision algorithm. For non-rectangular
meshes, the limit surface converges to a bicubic B-spline surface except at a finite number of 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 continuity is achieved.
In spite of the popularity of Catmull-Clark subdivision surfaces for representing complex geometric shapes
of arbitrary topology, these subdivision surfaces are not parameterizable and lack closed-form analytic for-
mulations. These deficiencies preclude their immediate pointwise manipulation and hence may restrain
the applicability of these schemes. We develop a new dynamic model based on Catmull-Clark subdivision
surfaces which offer modelers a closed-form analytic 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 any given level. However,
all the vertices introduced through subdivision are obtained as an affine 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 non-rectangular 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
may be noted that one additional subdivision step may be needed in some cases to isolate the extraordinary
points and treat the obtained mesh as the initial mesh (one typical example is when the initial mesh is a
tetrahedron).
B. Assigning patches to regular faces
F
F
3Fig. 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 B-spline surface (4 patches) in
the limit after an infinite number of subdivision steps. Note that, each of the bicubic patches in the
limit surface is defined by a rectangular face with each vertex of degree four, thereby accounting for
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 are assigned to face F 1 respectively. The 16 control points for
the patch S 1 , corresponding to face F 1 , are highlighted in Fig.1.
F
F
F
F
F
F
F
F
F
S1024681357Fig. 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
are assigned to bicubic patches S 0 respectively (as they all have vertices
of valence 4) following the aforementioned scheme. However, the central smooth surface enclosed by the
patches 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 different 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 Fig.3 where a local subdivision of the mesh consisting of faces F
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
selected mesh for
local subdivision
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 Fig.3 (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 surface
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-spline
2patch which is obtained as the limiting process of the Catmull-Clark subdivision algorithm applied
to an initial 4 by 4 rectangular control mesh. Let s p (u; v), where (u; v) 2 [0;
B-spline patch which can be expressed analytically as
where d i;j represents a 3-dimensional position vector at the (i; j)th control point location and B i;4 (u) and
are the cubic B-spline basis functions. The subscript p on s denotes the patch under consideration.
Expressing Eqn.1 in a generalized coordinate system we have
where J p is the standard Jacobian matrix of a bicubic B-spline patch, and is of size (3; 48). Vector q p is the
concatenation of all control points defining a B-spline patch in 3D. Note that in the concatenation of the
control points, each control point has an (x; component. For example, the (x; components of the
control point (i; correspond to positions 3k; 3k +1; 3k +2 - where, respectively in the vector
. We can express the entries of J p explicitly in the following way: J p (0;
D.1 Limit surface with many bicubic patches from a rectangular initial mesh
Now let's consider a limit surface consisting of many bicubic surface patches obtained after applying
an infinite number of subdivision steps to a rectangular initial mesh. For example, let the limit surface of
Fig.1 be s m which can be written as
where s m1 (2u;
2 , and 0 otherwise. Similarly, s m2 ; s m3 and s m4 are also equal
to s m (u; v) for an appropriate range of values of u; v and 0 outside. It may be noted that s
correspond to patches S 1 respectively in Fig.1. Rewriting Eqn.3 in generalized coordinates we
have
are the Jacobian matrices of size (3; 48) and q i s are the (x,y,z) component concatenation of a
subset of the control points of s m defining s m i
4. A more general expression for s m is
Where, qm is the 75-component vector of 3D positions of the 25 vertex control mesh defining the limit
surface s m . Matrices A are of size (48; 75) , each row consisting of a single nonzero entry
(= 1) and the (3; 75)-sized matrix
D.2 Limit surface with many bicubic patches from an arbitrary initial mesh
The stage is now set to define the limit surface s using the vertices of initial mesh M for any arbitrary
topology, assuming all faces are 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 may be required and the resulting mesh can be treated
as the initial mesh. Let the number of vertices in the initial mesh M be a, and let l 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 "normal face") and each of the remaining (b \Gamma
faces have one of the l extraordinary vertices (henceforth termed a "special face"). Let p be the
dimensional vector containing the control vertex positions in 3D. Using the formulations in subsections
II-B and II-C, the smooth limit surface can be expressed as
l
where n i is a single bicubic patch assigned to each of the normal faces and s j is a collection of infinite number
of bicubic patches corresponding to each of the extraordinary points. Employing the same approach
taken before to derive Eqn.5, it can be shown that
are the equivalent of J in Eqn.4 and A i 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
in the derivation of the expression for a collection of infinite number of bicubic patches around an extraordinary
vertex. The pre-superscript s is used to represent a collection of bicubic patches around an
extraordinary vertex, the subscript j is used to indicate the j-th extraordinary point, the post-superscript
represents the exponent of a mathematical quantity and the level indicator (to 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 s j is derived using the recursive nature of local subdivision around an extraordinary
vertex as shown in subsection II-C. First, s j can be expressed as
where the first term of Eqn.8 is the generalized coordinate representation of the bicubic B-spline 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). fs j g 1
represents the rest of the
infinite bicubic B-spline patches surrounding the extraordinary point (similar to the central patch enclosed
by patches marked n in Fig.3). The vertices in the newly obtained local subdivision mesh f s
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
. Furthermore, there exists a matrix f s
of size (3c; 3d), such that f s
are vectors of dimension 3c and 3d respectively. Applying the idea of recursive
local subdivision again on fs j g 1 , s j can be further expanded as
In the above derivation, f s ~ p j g 1 is a vector of dimension 3d, comprising of a subset of the vertices defining
the 3c dimensional vector f s
. Note that, f s ~
has the same structure as f s
, therefore, there
exists a (3d; 3d) matrix f s C
such that f s C
. Each subdivision of a local mesh with
d vertices creates a new local mesh with c vertices which contributes a fixed number of bicubic B-spline
patches. So, if we proceed one step further, we obtain
Because of the intrinsic property of the local recursive subdivision around the extraordinary point, we
November 10, 1997
have f s J
In addition, the subdivision rules remain the same
throughout the refinement process, we also have f s
we can further simplify the above equations leading to
{
{
s
s
{
{ j
{ j
Fig. 4. Local subdivision around the extraordinary point and the corresponding patches in the limit
surface from different levels of subdivision.
We can rewrite s j as
. 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 affine combination of some vertices in the previous level and hence any row of
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
l
l
l
From Eqn.6,7 and 13,
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
where an overstruck dot denotes a time derivative. The physics of the dynamic subdivision surface model
is based on the work-energy version of Lagrangian dynamics [38] and is formulated in an analogous way
to that in [27].
In an abstract physical system, let p i (t) be a set of generalized coordinates which are functions of time
and are assembled into the vector p. Let f i (t) be the generalized force assembled into the vector f p and
acting on p i . The Lagrangian equation of motion can then be expressed as
Let -(u; v) be the mass density function of the surface. Then
Z Z
is an N \Theta N mass matrix. Similarly the expression for damping matrix is
Z Z
where fl(u; v) is the damping density.
A thin-plate-under-tension energy model [39] is used to compute the elastic potential energy of the
dynamic subdivision surface. The corresponding expression for the stiffness matrix K is
Z Z
where the subscripts on J denote the parametric partial derivatives. The ff ii (u; v) and fi ij (u; v)s are
elasticity functions controlling local tension and rigidity in the two parametric coordinate directions. The
generalized force vector f p 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
Z Z
Multilevel Dynamics
Our dynamic Catmull-Clark surface model can be subdivided globally to increase the number of vertices
(control points) of the model. For example, after one step of global subdivision, the initial degrees of
(refer to Eqn.15 and Eqn.16) in the dynamic system will be replaced by a larger number of
degrees of freedom q, where Ap. A is a global subdivision matrix of size (M; N) whose entries are
uniquely determined by Catmull-Clark subdivision rules (see Section II-A for the details about the rules).
Thus, p, expressed as a function of q, can be written as
A T
Therefore, we can rewrite Eqn.15 and Eqn.16 as
and
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 "finer" level. The structure of the motion equation as given by Eqn.17 remains unchanged, but
the dimensionality and the entries of M;D;K;p and f p change correspondingly in this newly obtained
subdivided level. In particular the motion equation, explicitly expressed as a function of q, can be written
as
R R -B T J T JBdudv and the derivation of D q , K q and f q follow suit.
It may be noted that further subdivision, if necessary, can be carried out in a similar fashion. Therefore,
multilevel dynamics is achieved through recursive 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 system can automatically determine the level of
subdivision most suitable for an application depending on some application-specific criteria.
III.
The evolution of the generalized coordinates for our new dynamic surface model can be determined
by the second-order differential equation as given by Eqn.17. An analytical solution of the governing
differential equation can not be obtained in general. However, an efficient numerical implementation can
be obtained using finite element analysis techniques [40]. For the dynamic subdivision surface model, two
types 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.
A. Data Structures
A subdivision surface defined by a control mesh at any level is designed as a class which has a pointer
to its parent mesh, a set of pointers to its offspring meshes (arising out of local subdivision around the
extraordinary vertices at that level), a list of faces, edges, vertices and normal elements . Face, edge,
vertex and normal elements are, in turn, classes which store all the connectivity and other information
needed to either enumerate all the patches or locally subdivide around an extraordinary vertex in that
level. The implementation takes the initial mesh as the base subdivision surface object (with its parent
pointer set to NULL) and locally subdivides the initial mesh upto a user-defined maximum level around
each extraordinary vertex to create offspring objects at different levels. 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 Elements
Each normal element is a bicubic surface patch and is hence defined by 16 vertices (from the 8-connected
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 carrying out
the necessary integrations analytically. The matrix J in Eqn.18,19 and 20 needs to be replaced by J p
(of Eqn.2) for computation of the local M;D and K matrices respectively of the corresponding normal
element.
A.2 Special Elements
Each special element consists of an infinite number of bicubic patches in the limit. We have already
described a recursive enumeration of the bicubic patches of a special element in Section II-C. Let us now
consider an arbitrary bicubic patch of the special element in some level j. The mass matrix M s of this
patch can be written as
where M p is the normal element mass matrix (scaled by a factor of 1
to take into account of the area
shrinkage in bicubic patches at higher level of subdivision)
and\Omega s 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.
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 II-D.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 image-based forces. However, other types of forces like repulsion
forces, gravitational forces etc. can easily be implemented.
To apply spring forces, a spring of stiffness k can be connected from a point d 0 to a point
the limit surface, the net applied spring force being
Z Z
where ffi is the unit impulse function implying f vanishes elsewhere in
the surface. However, the ffi 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 file.
To recover shapes from 3D image data, we synthesize image-based forces. A 3D edge detection is
performed on a Gaussian smoothed volume data set using the 3D Monga-Deriche(MD) operator [41] to
produce a 3D potential field P (x; which we use as an external potential for the model. The force
distribution is then computed as
where k controls the strength of the force. The applied force on each element is computed using Gaussian
November 10, 1997
quadrature for evaluating Eqn.21 in Cartesian coordinates. It may be noted that we can apply spring
forces in addition with the image-based forces by placing points near the region of interest in the slices of
the 3D image data.
C. Discrete Dynamic Equation
The differential equation given by Eqn.17 is integrated through time by discretizing the time derivative
of p over time steps \Deltat. The state of the dynamic subdivision surface at time t + \Deltat is integrated using
prior states at time t and t\Gamma\Deltat. An implicit time integration method is used in the current implementation
where discrete derivatives of p are calculated using
\Deltat 2
and
2\Deltat
Using Eqn.17,29 and 30, the discrete equation of motion is obtained as
This linear system of equations is solved iteratively between each time step using the conjugate gradient
method. For a first order system with no mass, the above equation reduces to
(D
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-Clark 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 types
of error criterion can also be defined and used in this context. For example, in the context of image-based
forces, if the model energy does not change between successive iterations indicating an equilibrium for the
given resolution, the model can be subdivided further until the model energy is sufficiently small and the
November 10, 1997
change in energy between successive iterations becomes less than a pre-specified tolerance.
IV. RESULTS
The proposed dynamic subdivision surface can be used to represent a wide variety of shapes with
arbitrary genus. In this section we demonstrate the power of our modeling scheme via model fitting
examples to a variety of data sets of varying degree of complexity. In all the experiments, normal elements
are shaded yellow, while special elements are colored green.
In Fig.5(a) 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 Fig.5(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 may 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 (Fig.5(f)), the control mesh
is subdivided once thereby increasing the degrees of freedom (control vertices) of the underlying 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 sufficient 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 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 applying our algorithm, the data were transformed into a single reference
November 10, 1997
coordinate system. 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
Fig.8. The range data set has 1779 points in 3D. It may be noted that the final shape with a very low
error tolerance is recovered using very few number of control points in comparison to the number of data
points present in the original range data set. Fitting 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 MRI data in
the last two experiments. First, we fit the model to a cerebellum (a cortical structure in brain) given an
input of sagittal slices from a MR brain scan. Fig.10(a) depicts a slice from this MRI 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 MRI slices, each of size (256; 256). Fig.11(a) depicts a slice from this MRI 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).
V. CONCLUSIONS
In this paper, a dynamic generalization of the Catmull-Clark 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 analysis
of objects contained in range and volume data sets using very few degrees of freedom. We have presented
an analytic formulation of the subdivision scheme, incorporated the advantages of free-form 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 focussed toward addressing these issues.
VI.

ACKNOWLEDGEMENTS

This research was supported in part by the NSF grant ECS-9210648 and the NIH grant RO1-LM05944
to BCV, the NSF CAREER award CCR-9702103 and DMI-9700129 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 MRI data.



--R

"An algorithm for high speed curve generation,"
"A subdivision algorithm for smoothing down irregularly shaped polyhedrons,"
"Analysis of the behavior of recursive division surfaces near extraordinary points,"
The use of piecewise forms for the numerical representation of shape
"Recursively generated B-spline surfaces on arbitrary topological meshes,"
"Conditions for tangent plane continuity over recursively generated B-spline surfaces,"
"An investigation of curvature variations over recursively generated B-spline surfaces,"
"A unified approach to subdivision algorithms near extraordinary points,"
"Smooth subdivision surfaces based on triangles,"
"Efficient, fair interpolation using Catmull-Clark surfaces,"
"A signal processing approach to fair surface design,"
"A butterfly subdivision scheme for surface interpolation with tension control,"
"Subdivision schemes for surface interpolation,"
"Interpolating subdivision for meshes with arbitrary topology,"
"Elastically deformable models,"
"Deformable models,"
"Good vibrations : Modal dynamics for graphics and animation,"
"Dynamic deformation of solid primitives with constraints,"
"Multiresolution stochastic hybrid shape models with fractal priors,"
"Deformable curve and surface finite elements for free-form shape design,"
"Representing PDE surfaces in terms of B-splines,"
"Using partial differential equations to generate free-form surfaces,"
"Linear constraints for deformable B-spline surfaces,"
"Variational surface modeling,"
"Dynamic NURBS swung surfaces for physics-based shape design,"
physics-based framework for geometric design,"
"Dynamic NURBS with geometric constraints for interactive sculpting,"
"Finite-element methods for active contour models and balloons for 2-d and 3-d images,"
"Dynamic 3d models with local and global deformations: Deformable superquadrics,"
"Shape and nonrigid motion estimation through physics-based synthesis,"
"Recovery of non-rigid motion and structure,"
"Surface description of complex objects from multiple range images,"
"Adaptive-size physically-based models for nonrigid motion analysis,"
"Hierarchical shape representation using locally adaptive finite elements,"
"A finite element model for 3d shape reconstruction and nonrigid motion tracking,"

"Adaptive meshes and shells : Irregular triangulation, discontinuities and hierarchical subdivision,"
Hamilton's Principle and Physical Systems
"Regularization of inverse visual problems involving discontinuities,"
The Finite Element Handbook
"3d edge detection using recursive filtering,"
--TR

--CTR
Kaihuai Qin , Zhengyi Chang , Huawei Wang , Denggao Li, Physics-based loop surface modeling, Journal of Computer Science and Technology, v.17 n.6, p.851-858, November 2002
Frank Dachille, IX , Hong Qin , Arie Kaufman , Jihad El-Sana, Haptic sculpting of dynamic surfaces, Proceedings of the 1999 symposium on Interactive 3D graphics, p.103-110, April 26-29, 1999, Atlanta, Georgia, United States
Ye Duan , Hong Qin, Intelligent balloon: a subdivision-based deformable model for surface reconstruction of arbitrary topology, Proceedings of the sixth ACM symposium on Solid modeling and applications, p.47-58, May 2001, Ann Arbor, Michigan, United States
Han , Raffaele De Amicis , Giuseppe Conti, Interactive spline-driven deformation for free-form surface styling, Proceedings of the 2006 ACM symposium on Solid and physical modeling, June 06-08, 2006, Cardiff, Wales, United Kingdom
Seth Green , George Turkiyyah , Duane Storti, Subdivision-based multilevel methods for large scale engineering simulation of thin shells, Proceedings of the seventh ACM symposium on Solid modeling and applications, June 17-21, 2002, Saarbrcken, Germany
Kevin T. McDonnell , Hong Qin , Robert A. Wlodarczyk, Virtual clay: a real-time sculpting system with haptic toolkits, Proceedings of the 2001 symposium on Interactive 3D graphics, p.179-190, March 2001
Guiqing , Li Hua, Blending parametric patches with subdivision surfaces, Journal of Computer Science and Technology, v.17 n.4, p.498-506, July 2002
Gerold Wesche , Hans-Peter Seidel, FreeDrawer: a free-form sketching system on the responsive workbench, Proceedings of the ACM symposium on Virtual reality software and technology, November 15-17, 2001, Baniff, Alberta, Canada
Chhandomay Mandal , Hong Qin , Baba C. Vemuri, A novel FEM-based dynamic framework for subdivision surfaces, Proceedings of the fifth ACM symposium on Solid modeling and applications, p.191-202, June 08-11, 1999, Ann Arbor, Michigan, United States
Ye Duan , Hong Qin, A subdivision-based deformable model for surface reconstruction of unknown topology, Graphical Models, v.66 n.4, p.181-202, July 2004
Gerold Wesche , Marc Droske, Conceptual free-form styling on the responsive workbench, Proceedings of the ACM symposium on Virtual reality software and technology, October 22-25, 2000, Seoul, Korea
Steve Capell , Seth Green , Brian Curless , Tom Duchamp , Zoran Popovi, A multiresolution framework for dynamic deformations, Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation, July 21-22, 2002, San Antonio, Texas
Zo J. Wood , Peter Schrder , David Breen , Mathieu Desbrun, Semi-regular mesh extraction from volumes, Proceedings of the conference on Visualization '00, p.275-282, October 2000, Salt Lake City, Utah, United States
Chhandomay Mandal , Hong Qin , Baba C. Vemuri, Dynamic Modeling of Butterfly Subdivision Surfaces, IEEE Transactions on Visualization and Computer Graphics, v.6 n.3, p.265-287, July 2000
Denis Zorin, Modeling with multiresolution subdivision surfaces, ACM SIGGRAPH 2006 Courses, July 30-August 03, 2006, Boston, Massachusetts
finite elements;subdivision surfaces;interactive techniques;computer graphics;visualization;dynamics;CAGD;deformable models