This item has the following downloads:

1998274 ( PDF )

Direct Manipulation of Butterfly Subdivision



Surfaces: A Physics-based Approach








Chhandomay Mandal*, Hong Qint and Baba C. Vemuri*

















*Authors' address: Department of Computer and Information Science and Engineering, University of Florida,
Gainesville, Florida, 32611. Email:{cmandal | vemuri}lcise.ufl.edu.
tAuthor's address : Department of Computer Science, -I.. University of .'". York at '-..I- Brook, -',..I
Brook, .' York, 11794. i. .il 1 e!. - Ii.- lI, I


August 20, 1998


DRAFT
















Abstract


Subdivision surfaces are extensively used to model smooth shapes of arbitrary l I 1 !_ -. Re-

cursive subdivision on an user-defined initial control mesh generates a visually pleasing smooth

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

the control vertex positions at different levels of subdivision hierarchy to satisfy the functional

and aesthetic requirements in the smooth limit surface. This modeling drawback results from

the lack of direct manipulation tools for the limit surface. In this paper, we integrate techniques

from ,1!i -i, --based modeling with geometric subdivision i, 1 .1.i.1, ._, and present a scheme for

direct manipulation of the smooth limit surface generated by the (modified) butterfly scheme

using 1i1! -.i --based '-!,.!,, tools. This procedure-based surface model obtained through but-

terfly subdivision does not have a closed-form ,i!i ,1- I formulation (unlike other well known

spline-based models), and hence poses challenging problems to incorporate mass and damping

distribution functions, internal deformation ii. i forces, and other i1!. -. i1 quantities required

to develop a 1'1! -i. --based model. Our primary contributions to computer graphics and geomet-

ric modeling include : (1) a new hierarchical formulation for locally parameterizing the butterfly

subdivision surface over its initial control polyhedron, (2) formulation of 1,1! -i. --based butterfly

subdivision surface as a set of novel finite elements, and (3) optimal approximation of this new

I "' of finite elements by a collection of existing finite elements -,. I 1, ,i 1 to implicit geometric

constraints. Our new dynamic model can be sculpted directly by ;,1l', i~l~ synthesized forces,

and its equilibrium is characterized by a minimum of the deformation ii. i -L --.ii. I to the im-

posed constraints. We demonstrate that this novel dynamic framework not only provides a direct

and natural means of manipulating geometric shapes, but also facilitates hierarchical shape and

non-rigid motion estimation from large range and volumetric data sets using very few degrees

of freedom (control vertices that define the initial polyhedron). This new ,1!i -;. --based model

promises a greater potential of subdivision schemes in computer graphics, geometric modeling,

and virtual environments.


August 20, 1998


DRAFT












Keywords


Computer Graphics, P!- -;- --based Modeling, Geometric Modeling, CAGD, Scientific Visu-

alization, Subdivision Surfaces, Deformable Models, I !!!!i~ 1.- I ii I- Interactive Techniques.





I. INTRODUCTION

The concept of subdivision is ubiquitous in computer science. In visual computing areas, subdivision

surfaces are extensively used to model smooth shapes of arbitrary l. ,',.1 _-- for computer graphics, an-

imation, and geometric design applications. A I- 1i.; 1 recursive subdivision scheme produces a visually

pleasing smooth surface in the limit by repeated application of a fixed set of refinement rules on an user-

defined initial control mesh. Despite the presence of a !- i;. of subdivision schemes in the computer

graphics and geometric modeling literature, there is no direct and natural way of manipulating the limit

surface. The current state-of-the-art only permits the modeler to interactively obtain the desired effects

on the smooth surface by kinematically manipulating the vertex positions at various levels of subdivision

hierarchy. In this paper, we tackle the challenging problem of direct manipulation of the limit subdivision

surface at arbitrary locations/areas and offer a novel solution to this problem by embedding the modified

butterfly subdivision scheme in a 1!i -;1 --based modeling framework. As a result, unlike the existing

geometric solutions that only allow the operations on control vertices, our i. !. ,1.l .1._- and algorithms

permit the user to '11! -"; 1-! modify the shape of subdivision surfaces at desired locations via application

of forces. This gives the user an intuitive and natural feeling that is produced while modeling with real

clay/' I -'1.. i-l1! We will also demonstrate that the proposed model t I- !. II, recovers shapes as well as

non-rigid motions from large range and volumetric data sets. Note that, this paper neither proposes a

new subdivision technique nor provides a different interpretation of ;!!, existing subdivision technique,

but integrates the advantages of subdivision surface-based and 1i1! -;i --based modeling techniques to solve

important theoretical and practical problems. Although the principles of 1'1! -i --based modeling are well

understood by the graphics experts and modeling researchers, this paper will greatly advance the state

of the art in i1.- -i; --based shape modeling due to the following contributions : (1) a hierarchical local


August 20, 1998


DRAFT











parameterization scheme for (modified) butterfly subdivision surfaces is derived; (2) material properties

are assigned to the smooth limit surface in order to embed this geometric model in a '11! -i. --based mod-

eling framework; (3) the smooth limit surface is decomposed into a set of novel finite elements; (4) the

motion equations that govern the behavior of (modified) butterfly subdivision surfaces are derived; (5)

users can sculpt this 1i'1! -i. --based model at desired locations and the model responds naturally (as a

real-world object would); and (6) algorithms and procedures are developed to best approximate this novel

finite element using a collection of existing lower order finite elements. WVe will address the ii..- -1 of

the proposed model, detail our contributions, and compare this new model with our previous research

results on dynamic subdivision surfaces in later sections. in -I we will ',i f[- review the previous work

on subdivision surfaces.


A. Background

In [1], ('!i ii!:!n first introduced the concept of subdivision to the computer graphics .!!' .iil,! for

generating a smooth curve from a given control i..1 _.i During the last two decades, a wide i- I, -

of subdivision schemes for modeling smooth surfaces of arbitrary 1. ,..1. _-- have been derived following

('I! ii,!:, pioneering work on curve generation. In general, these subdivision schemes can be catego-

rized into two distinct classes namely, (1) approximating subdivision techniques and (2) interpolating

subdivision techniques.

Among the approximating schemes, the techniques of Doo and Sabin [2], [3] and Catmull and ('! !:

[4] generalize the idea of obtaining uniform biquadratic and bicubic B-spline patches respectively from

a rectangular control mesh. In [4], Catmull and ('I i !: developed a method for recursively generating a

smooth surface from a polyhedral mesh of arbitrary Il. .1. -. -. The Catmull-('I! I : subdivision surface,

defined by an arbitrary initial mesh, can be reduced to a set of standard B-spline patches except at a

finite number of degenerate points. In [5], Loop presented a similar subdivision scheme based on the

generalization of quartic triangular B-splines for triangular meshes. Hoppe et al. [6] extended his work

to produce piecewise smooth surfaces with selected discontinuities. Halstead et al. [7] proposed an

algorithm to construct a Catmull-( '! !: subdivision surface that interpolates the vertices of a mesh of


August 20, 1998


DRAFT












arbitrary I. 1I,. .1 Peters and Reif [8] proposed a simple subdivision scheme for smoothing i" .1 1!. -1

Most recently, non-uniform Doo-Sabin and Catmull-('!I !: surfaces that generalize non-uniform tensor

product B-spline surfaces to arbitrary topologies were introduced by Sederberg et al. [9]. All the schemes

mentioned above generalize recursive subdivision schemes for generating limit surfaces with a known

parameterization. Various issues involved with character animation using these approximating subdivision

schemes were discussed at length by DeRose et al. [10].

The most well-known interpolation-based subdivision scheme is the "1.111i. i!l algorithm proposed by

Dyn et al. [11]. Butterfly subdivision method, like other subdivision schemes, makes use of a small

number of neighboring vertices for subdivision. It requires simple data structures and is extremely easy

to implement. However, it needs a topologically regular setting of the initial (control) mesh in order

to obtain a smooth C1 limit surface. A variant of this scheme with better smoothness properties can

be found in [12]. Zorin et al. [13] has developed an improved interpolatory subdivision scheme (which

we call the !i... ..r' 1 butterfly scheme) that retains the -ii.wli, '1 of the butterfly scheme and results

in much smoother surfaces even from irregular initial meshes. These interpolatory subdivision schemes

have extensive applications in wavelets on manifolds, multiresolution decomposition of polyhedral surfaces

and multiresolution editing. A variational approach for interpolatory refinement has been proposed by

Kobbelt [14], [15] and by Kobbelt and Schr6der [16]. In this approach, the vertex positions in the refined

mesh at each subdivision step are obtained by solving an optimization problem. Therefore, these schemes

are global, i.e., every new vertex position depends on all the vertex positions of the coarser level mesh.

The local refinement property which makes the subdivision schemes attractive for implementation in the

computer graphics applications is not retained in the variational approach.

The derivation of various mathematical properties of the smooth limit surface generated by the sub-

division algorithms is rather complex. Doo and Sabin [17] first i!i! 1- ,1 the smoothness behavior of

the limit surface using the Fourier transform and an ._. i_- 1, ,1- -i- of the subdivision matrix. Ball and

1i..i, [18], [19] and Reif [20] further extended Doo and Sabin's prior work on .. .ii,,,il I properties of

subdivision surfaces by deriving various necessary and -it!!li wI conditions on smoothness for .[l!!. i. II

subdivision schemes. Specific subdivision schemes were :i! ~1- .1 by Schweitzer [21], Habib and Warren


August 20, 1998


DRAFT












[22], Peters and Reif [23] and Zorin [24]. Most recently, i, i 1i [25] presented a method for exact evaluation

of Catmull-('I! I : subdivision surfaces at arbitrary parameter values.


B. s.., modeling using the ;,,,' -based subdivision-surface model


Although recursive subdivision surfaces are powerful for representing smooth geometric shapes of ar-

bitrary l.|l, ,.-, they constitute a purely geometric representation, and furthermore, conventional geo-

metric modeling with subdivision surfaces 11i be [t! t111 for representing highly complicated objects.

For example, modelers are faced with the tedium of indirect shape modification and refinement through

time-consuming operations on a large number of (most often irregular) control vertices when using I- ,i-

cal subdivision surface-based modeling schemes. Despite the advent of advanced 3D graphics interaction

tools, these indirect geometric operations remain non-intuitive and laborious in general. In addition, it

ii not be enough to obtain the most ''I !! surface that interpolates a set of (ordered or unorganized)

data points. A certain number of local features such as bulges or inflections i! be strongly desired while

requiring geometric objects to satisfy global smoothness constraints in geometric modeling and computer

graphics applications. In contrast, i'1! -1, --based modeling provides a superior approach to shape model-

ing that can overcome most of the limitations associated with traditional geometric modeling approaches.

Free-form deformable models governed by the laws of continuum mechanics are of particular interest in

this context. These dynamic models respond to externally applied forces in a very intuitive manner.

The dynamic formulation marries the model geometry with time, mass, damping and constraints via a

force balance equation. Dynamic models produce smooth, natural motions which are easy to control. In

addition, they facilitate interaction especially direct manipulation of complex geometries. Furthermore,

the equilibrium state of the model is characterized by a minimum of the deformation ii. i of the model

- ,i, 'r i to the imposed constraints. The deformation i -i functionals can be formulated to satisfy local

and global modeling criteria, and geometric constraints relevant to shape design can also be imposed.

The dynamic approach subsumes all of the aforementioned modeling capabilities in a formulation which

grounds - i li!i_ in real-world i,1! -ii ,1 behavior.

Free-form deformable models were first introduced to computer graphics and visualization in Terzopou-


August 20, 1998


DRAFT











los et al. [26] and further developed by Terzopoulos and Fleischer [27], Pentland and \\ ill ,!!!- [28],

Metaxas and Terzopoulos [29] and Vemuri and I ,1l- li. -; [30]. Celniker and Gossard [31] developed

a -I, i, for interactive free-form design based on the finite element optimization of i functionals

proposed in [27]. Bloor and \\ i! -. [32], [33], Celniker and Welch [34] and Welch and \\ i !:,i [,;.] proposed

deformable B-spline curves and surfaces which can be designed by imposing the shape criteria via the

minimization of the ii. i functionals -,i. lr I to hard or soft geometric constraints through Lagrange

multipliers or i" ii ,l- methods. Qin and Terzopoulos [36], [37], [38] developed dynamic NURBS (D-

NURBS) which are very sophisticated models suitable for representing a wide i. of free-form as well

as standard !! 1 1 i shapes. The D-NURBS have the advantage of interactive and direct manipulation of

NURBS curves and surfaces, resulting in ,1!. ii-11 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 rectangular parametric domain. Hence, it can be very lIi 1iII to model surfaces of arbitrary genus using

these models. In [10], DeRose et al. assigned material properties to control meshes at various subdivision

levels in order to simulate cloth dynamics using subdivision surfaces. Note that, they assign i1- -i, ,1

properties on the control meshes at various levels of subdivision and not on the limit surface itself, and

hence can not solve the modeling goal we are I i 11 to achieve. Previously we introduced a dynamic

Catmull-('C! I subdivision surface model [39], [40] which combined the benefits of subdivision surfaces

for modeling arbitrary l,,.I.1.--_- as well as that of dynamic splines for interactive shape manipulation

by 1',,1, i 1- synthesized forces. The dynamic (modified) butterfly subdivision surface model formulated

and developed in this paper aims to achieve the same long-term objective, i.e., a formal mechanism of

allowing the modeler to directly and intuitively manipulate the smooth limit surface of arbitrary I .1. -

as if they were seamlessly sculpting a piece of real-world "- !I '. However, this new model is superior to

our previously reported research in several significant aspects which will be discussed in detail later in

Section I-D. In this paper, we derive a novel technique for locally parameterizing the smooth limit surface

generated by the modified butterfly subdivision surface algorithm which embeds the proposed model in

a dynamic framework in a straight-forward manner. The model can be initialized interactively by an


August 20, 1998


DRAFT











user-defined control mesh and is amenable to further sculpting via direct application of !i l -i ,1 forces

to ,i i region of object surface. The formulation and implementation details are discussed in subsequent

sections.


C. !ji,1, and motion estimation using the -;/,', .-based subdivision-surface model

The dynamic subdivision surface model has been developed primarily for modeling arbitrary (known)

l I'1" ,! ,_ where modelers can directly manipulate the limit surface by ;,.1,, l -- ,l. I .1 ,1 forces in an

intuitive fashion. However, as we have shown in our earlier work [41], another important application of the

dynamic subdivision surfaces is in non-rigid shape and motion reconstruction/recovery. Accurate shape

recovery requires distributed parameter models which I- i. ,11- possess a large number of degrees of free-

dom. On the other hand, Iti !- ii shape representation imposes the requirement of geometry compression,

i.e., models with fewer degrees of freedom. These requirements are conflicting and numerous researchers

have been seeking to strike a balance between these contradicting requirements [30], [38], [41], [42], [43],

[44], [45], l11], [47]. Another important criterion in the design of shape models is that the initialization

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

meshes since it in be infeasible to globally parameterize shapes of arbitrary 1i i.1 ._ -. A 1i.1 --. --based

model best 11r- i the above mentioned criteria is an ideal candidate for a solution to the shape recovery

problem for obvious reasons.

Deformable models which come in in, i,- varieties, have been used to solve the problem in the I'1! -' --

based modeling paradigm. These models involve the use of either fixed size [30], [43], [48], [49], [50]

or adaptive size [44], i], [51], [52],, ;], [54] grids. The models with fixed grid size generally use less

number of degrees of freedom for representation at the cost of accuracy of the recovered shape. On the

other hand, models using adaptive grids generally need large number of degrees of freedom to recover the

shapes accurately. Note that the shapes being recovered from the image data are smooth in most of the

medical applications, i.e. the anatomical structures even with considerable amount of details have more

or less a C1 surface. Under these circumstances, the finite element approaches as in [43], [46] need a

large number of degrees of freedom for deriving a smooth and accurate representation. In addition, they


August 20, 1998


DRAFT












can not represent shapes with known arbitrary I. 1-. 1, _-. Moreover, almost all of these schemes require a

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

Our previous dynamic subdivision surface model [39], [40], [41] offered an elegant solution to the above

mentioned problem as it could recover complex shapes in a hierarchical fashion using very few degrees

of freedom without requiring parameterized input mesh. However, the model proposed in this paper

outperforms the previous one in the compactness of the model representation. We will show experimental

results in support of this claim. We will also demonstrate the potential of this model in motion tracking

and visualization of a dynamically deforming shape from a time sequence of volumetric data sets. Like

the previous model, the dynamic modified butterfly subdivision surface model also deforms under the

influence of ,il --i 1 forces to fit the i1 il i,1 i!_- shape in the given range or volumetric data set via

the principle of i. minimization.


D. Contributions

In this section, we summarize the contributions of the present work along with the advantages of

the proposed dynamic modified butterfly subdivision surface model over the dynamic Catmull-('I! !:

subdivision surface model [39], [40]. The primary contributions of the present work are as follows.

We -- -l. 1 1, ,1 i11- derive a local parameterization scheme for the modified butterfly subdivision sur-

faces in a hierarchical -- 1. and subsequently the initial control polyhedron can be viewed as the

parametric domain of the 1i'1! -i --based smooth limit surface.

We treat the dynamic subdivision surface in the limit as a "'. i!" 1'1 -1,1 1 object and represent the

smooth limit surface as a set of novel finite elements whose shape (basis) functions are derived using

the modified butterfly subdivision scheme. We envision that this new finite element method will prove

to be useful not only in the areas of computer graphics and geometric design, but also in engineering



We use modified butterfly subdivision techniques to create a surface model that incorporates mass and

damping distribution functions, internal deformation i!. forces, and other 1i1!1 -11 i1 quantities. We

also -I. I ,11 formulate the motion equations for this (modified) butterfly subdivision surface


August 20, 1998


DRAFT












whose degrees of freedom are the collection of initial user-specified control vertices. Therefore, the

advantages of both the |'1! -i. --based modeling I'!'i1..-.. i! and the geometric subdivision schemes are

presented within a single unified framework.

Users will be able to manipulate this i,1! -i1 --based model in an arbitrary region, and the model

responds naturally (just like the real-world object) to this force application. This shape deformation

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

the geometry of the limit surface.

We develop algorithms and procedures which approximate our novel finite element using a collection

of linear and/or bilinear finite elements il. i I to the implicit geometric constraints enforced by the

butterfly subdivision rules. This hierarchically-structured approximation is optimal in the sense that

it can fall within :,,i user-specified error tolerance.

Although the long-term goals inherent in the lo' ;* -1 -1 1, i. .1 dynamic Catmull-('I! !: model [39],

[40] are the same as our current endeavor of deriving dynamic modified butterfly subdivision scheme, the

research presented in this paper achieves them in a much more elegant fashion. I i! -I of all, the finite

element implementation of the dynamic Catmull-(' : subdivision surface is specific to the subdivision

technique involved where a diversity of complicated finite elements must be i| il1- .1 in order to account

for the special cases and can not be readily generalized to other approximating subdivision schemes in

a straight-forward way. However, the finite element techniques developed in this paper can be easily

generalized to other interpolatory (as well as I I '! i !!1 i1 i subdivision schemes involving triangular (or

n-sided) meshes. Second, for some specific cases the plate i i of Catmull-('I! !: subdivision surface

diverges as shown in [7], hence various case :iI 1- -- need to be performed to derive the internal i. of

the dynamic Catmull-('I ! subdivision surface model. The internal deformation ii. i of the dynamic

butterfly scheme can be derived in an unified fashion for a spectrum of scenarios. Third, for both models,

we need to derive the closest point on the limit surface from a given point in 3D for force applications.

The calculation overhead involved in this process is significantly less for the butterfly case as it is an

interpolatory scheme where all vertices at various levels of subdivision lie on the limit surface, and the

search space can be reduced rapidly in a hierarchical fashion. The situation is quite 1!I! i, II for Catmull-


August 20, 1998


DRAFT












('I! iil subdivision scheme since it is an approximating scheme, and the technique used for finding the

closest model point [48] is computationally expensive. The force application is vital to ; i' I'1! -'. --

based model, and hence the adopted computationally inexpensive method to find the closest point for

the proposed scheme in this paper has very significant advantages over the previous model. Lastly, it

has been empirically found that the recovered shape is more compact (less number of degrees of freedom)

when using the proposed model in comparison with our earlier reported model using the same data set

and model-fitting criteria.


E. Overview

The rest of the paper is organized as follows: Section II presents the detailed formulation of the dynamic

modified butterfly subdivision surfaces. The implementation details are presented in Section III followed

by experimental results in Section IV. I ii I1I we make the concluding remarks in Section V.


II. FORMULATION

In this section, we provide a -I. i!,, formulation of the dynamic subdivision surface model. I i, -I

we ', i. f- review the modified butterfly subdivision scheme. Next we introduce a local parameterization

scheme which will facilitate the formulation of the smooth limit surface as a function of the control point

positions defining the initial mesh. This parameterization scheme is then used to derive the dynamic

model. It ii ,- be noted that these techniques can be generalized to define and construct a generic

dynamic framework for other triangle-based subdivision surface schemes as well, but we will focus only

on the modified butterfly subdivision technique in this paper.


A. I !, ( ...... "1', ,I) butterfly subdivision scheme

The butterfly subdivision scheme [11], like 11ii i~ other subdivision schemes used in geometric design

literature/applications, starts with an initial triangular mesh which is also known as the control mesh.

The vertices of the control mesh are known as the control points. In each step of subdivision, the initial

(control) mesh is refined through the transformation of each triangular face into a patch with four smaller

triangular faces. After one step of refinement, the new mesh in the finer level retains the vertices of each


August 20, 1998


DRAFT











\ \ I \

























triangular face in the previous level and hence, interpolates the coarser mesh in the previous level. In
S\ \1 \ I \ / \

I -II" I










(a) (b)
Fig. 1. (a) The control position an with triangular faces; (b) Mvesh obtained after one subdivision step.


triangular face in the previous level and hence, interpolates the coarser mesh in the previous level. In

addition, every edge in each triangular face is spilt by adding a new vertex whose position is obtained by

an ;I !!!r, combination of the neighboring vertex positions in the coarser level. For instance, the mesh in

Fig.l(b) is obtained by subdividing the initial mesh shown in Fig.l(a) once. It c' be noted that all the

newly introduced vertices (marked in blue) corresponding to the edges in the original mesh have valence

(degree) 6, whereas the position and valence of the original vertices (marked in red) do not change in the

subdivided mesh.

In the original butterfly scheme, the new vertices corresponding to the edges in the previous level are

obtained using an eight-point stencil as shown in i 2, -, i, The name of the scheme originated from the

" 1.i I. i II '-like configuration of the contributing vertices. The weighing factors for I h!! I. i!. contributing

vertex positions are shown in Fig.2(b). The vertex e l in the j + 1-th level of subdivision, corresponding

to the edge connecting vertices vi and v2 at level j, is obtained by

ei+1 = 0.5(v + v?) + D. \1 + v ) w(v 3 + v + v + v3), (1)
S< w < and v denotes the position of the h ver at the jh level.
where 0 < w < 1, and vj denotes the position of the i-th vertex at the j-th level.


August 20, 1998


DRAFT











j -w -w

5 0.5










7 2 8 -w -w


(a) (b)
Fig. 2. (a) The contributing vertices in the j-th level for the vertex in the j+l-th level corresponding to
the edge between v1 and v2; (b) the weighing factors for different vertices.

The butterfly subdivision scheme produces a smooth C1 surface in the limit except at the extraordinary

points corresponding to the extraordinary vertices (vertices with valence not equal to 6) in the initial

mesh [13]. All the vertices introduced through subdivision have valence 6, and therefore, the number of

extraordinary points in the smooth limit surface equals the number of extraordinary vertices in the initial

mesh. Recently, this scheme has been modified by Zorin et al. [13] to obtain better smoothness properties

at the extraordinary points. In [13], all the edges have been categorized into three classes : (i) edges

connecting two vertices of valence 6 (a 10 point stencil, as shown in I i, ;, i is used to obtain the new

vertex positions corresponding to these edges), (ii) edges connecting a vertex of valence 6 and a vertex of

valence n r 6 (the corresponding stencil to obtain new vertex position is shown in Fig.3(b), where q = .75

is the weight associated with the vertex of valence n r 6, and si = (0.25 + cos(27i/n) + 0.5cos(47i/n))/n,

i = 0, 1,..., n 1, are the weights associated with the vertices of valence 6), and (iii) edges connecting two

vertices of valence n r 6. The last case can not occur except in the initial mesh as the newly introduced

vertices are of valence 6, and the new vertex position in this last case is obtained by averaging the positions

obtained through the use of stencil shown in Fig.3(b) at each of those two extraordinary vertices.


August 20, 1998


DRAFT






































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


B. Local parameterization of the limit '.. f f,..I

Oftentimes, the smooth limit surface defined by the recursive subdivision process is of;' ,1 'i 1, 1. I. 1- 1_ *

where a global parameterization in ,- not be possible. Nevertheless, we can locally parameterize the limit

surface over the domain defined by the initial mesh following a similar approach described in [55]. The

idea is to track ; i arbitrary point on the initial mesh across the meshes obtained via the subdivision

process (see Fig.4 and 5), so that a correspondence can be established between the point being tracked in

the initial mesh and its image on the limit surface.

The modified butterfly subdivision scheme starts with an initial mesh consisting of a set of triangular

faces. The recursive application of the subdivision rules smoothes out each triangular face, and in the

limit, a smooth surface is obtained which can also be considered as a collection of smooth triangular

patches. The subdivision process and the triangular decomposition of the limit surface is depicted in

Fig.4. Note that, the limit surface can be represented by the same number of smooth triangular patches


August 20, 1998


DRAFT
















x X -* *-





Fig. 4. The smoothing effect of the subdivision process on the triangles of the initial mesh.


as that of the triangular faces in the initial mesh. Therefore, we can express the limit surface s as
a
s = k, (2)
k=l
where n is the number of triangular faces in the initial mesh and sk is the smooth triangular patch in the

limit surface corresponding to the k-th triangular face in the initial mesh.

We are now ready to describe the parameterization of the limit surface over the initial mesh. The

process is best explained through the following example. We choose a simple planar mesh shown in

I i, ,, as the initial mesh. An arbitrary point x inside the triangular face abc is tracked over the

meshes obtained through subdivision. The vertices in the initial mesh are drawn in black in Fig.5. After

one step of subdivision, the initial mesh is refined by addition of new vertices which are colored green.

Another subdivision step on this refined mesh leads to a finer mesh with introduction of magenta colored

new vertices. It !1 i be noted that ., ,. point inside the smooth triangular patch in the limit r ,,tII..

corresponding to the face abc in the initial mesh depends (..l', on the vertices in the initial mesh which

are within the 2-neighborhood of the vertices a, b and c due to the local nature of the subdivision process.

For example, the vertex d, introduced after first subdivision step, can be obtained using the 10 point

stencil shown in I i, ;,: on the edge ab. All the contributing vertices in the initial mesh are within the

1-neighborhood of the vertices a and b. A 10 point stencil can be used again in the next subdivision step

on the edge db to obtain the vertex g. Some of the contributing vertices at this level of subdivision, for

example, the (green colored) 1-neighbors of the vertex b (except d and e) in I i, -.(b), depend on some

vertices in the initial mesh which are within the 2-neighborhood of the vertices a, b and c in the initial

mesh.

In the rest of the formulation, superscripts are used to indicate the subdivision level. For example,


August 20, 1998


DRAFT


































(a) (b)


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


August 20, 1998


DRAFT











v',, denotes the collection of vertices at level j which control the smooth patch in the limit surface

corresponding to the triangular face uvw at the j-th level of subdivision. Let vobc be the collection of

vertices in the initial mesh which are within the 2-neighborhood of the vertices a, b and c (marked black

in i i, -,, ,). Let the number of such vertices be r. Then, the vector vObc, which is the concatenation

of the (x, y, z) positions for all the r vertices, is of dimension 3r. These r vertices control the smooth

triangular patch in the limit surface corresponding to the triangular face abc in the initial mesh. Now,

there exists four (3r x 3r) subdivision matrices (Aabc)t, (Aabc)l, (Aabc)r and (Aabc)m such that

vadf = (Aabc)tVabc,

Vbed = (Aabc) Vabc,

IVfe = (Aabc),V bc
VCf( abc,

V f = (Aabc),, b (3)

where the subscripts t, I, r and m denote top, left, right and middle triangle positions respectively (indi-

cating the relative position of the new triangle with respect to the original triangle), andVadfi, v d, v1fel

and Vef are the concatenation of the (x, y, z) positions for the vertices in the 2-neighborhood of the cor-

responding triangle in the newly obtained subdivided mesh. The new vertices in this level of subdivision

are marked green in Fig.5(b). The 2-neighborhood configuration of the vertices in the newly obtained

triangles is exactly the same as that of the original triangle, hence local subdivision matrices are square

and the vector dimensions on both sides of Eqn.3 are the same. This concept is further illustrated in

Fig.6.

C i11- out one more level of subdivision, along with the old vertices, we get a new set of vertices

which are marked in magenta in I i, -.(c). Adopting a similar approach as in the derivation of Eqn.3, we

obtain

Vdgi = (Abed)tVbed

2
V hg = (Abed)1V ed

V eih = (Abed),vred

V hi = (Abed),-,Ved (4)


August 20, 1998


DRAFT















































































































August 20, 1998


--------


DRAFT











































































Fig. 6. Different set of control points at a subdivided level obtained by ;1:''1;! ,ltl' i' !ii subdivision
matrices on a given set of control points in a coarser mesh.


August 20, 1998


DRAFT











The relative position of the triangular face ~I.,' in Fig.5(c) with respect to the triangular face bed is

topologically the same as of the triangular face adf in i, -.(b) with respect to the triangular face abc.

Therefore, we can write (Abed)f = (Aabc)1. Using similar reasoning, Eqn.4 can be rewritten as

vdgi = (Abed)tVled = (Aabc)tled

Vhg = (Abed)lVed = (Aabc),vbed

h = (Abed)rVed = (Aabc)rVLed

Vghi = (Abed),nVbed = (Aabc),nVed. (5)

Combining Eqn.3 and Eqn.5, we get

Vgi = (Aabc)t(Aabc)lv b,

bg= (Aabc)l(Aabc)lVabc>

Veih = (Aabc)r(Aabc)lVabc,

Vghi = (Aabc),(Aabc)lVabc. (6)

Let x be a point with '. "1 i. i il,.; coordinates (abc, ab, Yabc) inside the triangular face abc. When the

initial mesh is subdivided, x becomes a point inside the triangular face bed with '.i id i i. coordinates

(ab~d,>Pled,' .). Another level of subdivision causes x to be included in the triangular face 1, with

i I. !i coordinates (a di,dgi, T gi). Let sabc denote the j-th level approximation of the smooth

triangular patch Sabc in the limit surface corresponding to the triangular face abc in the initial mesh. Now

bc can be written as
r r r
vObc bx ,c,....,ay,by,cy, ..,a,, b ,c ,, .. ]T

where the subscripts x, y and z indicate the x, y and z coordinates respectively of the corresponding

vertex position. The expressions for v ed and vdg, can also be written in a similar manner. Next, we

construct the matrix Bbc as follows:

bcc, ,bc, T bc, 0,...,0, 0,...,0,0,...,0

B. 0abc ab, 0,bc, b 0, ...0, 00, ..,0

0, 0. ,0 ... ab, O bc bc ,..., 0


August 20, 1998


DRAFT











The matrices Bed and B g. can also be constructed in a similar fashion. We can now write sbc(X),

abc(x), and S bc(x) as
sbc(x) = B0 bc(X)V

bc(x) = Bd(x)vbed = B (x)(Albc)l AVb, (7)
S 2.. B 0

Sbc(x) = B2. x = B xi(Aabc)Ved = B .x(Aabc)(Aabc)lVabc.


Proceeding in a similar way, the expression for s bc(x), j-th level approximation of Sabc(x), is given by
3
Sbc(X) = B,.(x) (Aabc), ...(Aabc)t(Aabc)V bV

= B4 ,\xi(Abc)VOb

= bc(X)Vab, (8)

where x is inside the triangular face uvw at level j (with an assumption that uvw is the triangular face

in the middle with respect to its coarser level original triangular face in the previous level), (Abc) =

(Aabc),n ... (Aabc)t(Aabc)1 and BIbc(x) = B1 ,fi(Ab). It ii, be noted that the sequence of i.I- 'I ,-

(Aabc)t, (Aabc)1, (Aabc)r and (Aabc) depends on the triangle inside which the tracked point x falls after

each subdivision step. I ii!! we can complete the local parameterization process by writing

0 -+
Sabc(x) = ( lim Bb,(X))Vbc = Babc(X)Vbc. (9)

In the above equation, Bab, is the collection of basis functions at the vertices of vbc. It II also be

noted that the modified butterfly subdivision scheme is a stationary subdivision process, and hence new

vertex positions are obtained by ,t !!ii. combinations of nearby vertices. This guarantees that each row of

the matrices (Aabc)t, (Aabc)1, (Aabc)r and (Aabc), sums to one. The largest eigenvalue of such matrices

is 1 and therefore the limit in Eqn.9 exists. Now, if we assume that the triangular face abc is the k-th face

in the initial mesh, then Eqn.9 can be rewritten as

sk(x) = Bk(x = Bk(x)Akp, (10)

where p is the concatenation of the (x,y,z) positions of all the vertices in the initial mesh and the matrix

Ak, when post-multiplied by p, selects the vertices vA controlling the k-th smooth triangular patch in the

limit surface. If there are t vertices in the initial mesh and r of them control the k-th patch, then p is a

vector of dimension 3t, Ak is a (3r x 3t) matrix and Bk (x) is a (3 x 3r) matrix.


August 20, 1998


DRAFT











Combining Eqn.2 and Eqn.10, we get

s(x) = ( B (x)Ak)p = J(x)p, (11)
k=l
where the (3 x 3t) matrix J is the collection of basis functions for the corresponding vertices in the initial

mesh. The vector p is also known as the degrees of freedom vector of the smooth limit surface s.


C D ..... ', ,

In the previous section, we have derived an expression of the smooth limit surface obtained via infinite

subdivision steps as a function of the control vertex positions in the initial mesh. We now treat the vertex

positions in the initial mesh defining the smooth limit surface s as a function of time in order to embed

the modified butterfly subdivision model in a dynamic framework. The velocity of this surface model can

be expressed as

s(x,p) =J(x)p, (12)

where an overstruck dot denotes a time derivative and x E So, So being the domain defined by the initial

mesh. The 1i1! -'. of the dynamic subdivision surface model is based on the i-. i!:-i i relationship of

Lagrangian dynamics [,1,] and is formulated in an analogous way to that in [39], [40].

In an abstract 1i1! 1 -- .. i let pi(t) be a set of generalized coordinates which are functions of time

and are assembled into the vector p. Let fi(t) be the generalized force assembled into the vector fp and

acting on pi. Then, the Lagrangian equation of motion can be expressed as
d OT OT OF OU
dti p+ + = fi, (13)
dt dpi dpi dpi dpi
where T, F and U are the kinetic, dissipation and potential. -. i _- respectively.

Let p be the mass -l. i,-i- function of the surface. Then the kinetic i. i,- of the surface is

T = 1 p(x)sT(x)s(x)dx = IMp, (14)
2 xeSO 2
where (using Eqn.12) M = Jxo P(x)JT(x)J(x)dx is the (3t x 3t) mass matrix. Similarly, let 7 be the

damping 1I. i! -i function of the surface. The dissipation !i. i is

F = 1 f (x)T (x)s(x)dx = IT D, (15)
2 xxxx e 2
where D = fS, o 7(x)JT(x)J(x)dx is the (3t x 3t) damping matrix. The potential i i of the smooth


August 20, 1998


DRAFT












limit surface can be expressed as

U pTKp, (16)
2
where the (3t x 3t) stiffness matrix K is obtained by assigning various internal energies to a discretized

approximation of the limit surface and is detailed in Section III.

Using the expressions for the kinetic, dissipation and potential i i in Eqn.13, we get the motion

equation given by

Mpi + Dp + Kp = fp. (17)

The generalized force vector fp, which can be obtained through the principle of virtual work 'i,]j, is

expressed as

f= J(x)f(x, t)dx. (18)

We can apply various I I," of forces on the smooth limit surface, and the limit surface would evolve over

time according to Eqn.17 to obtain an equilibrium position characterized by a minimum of the total model

- I I 1


C.1 Multilevel Dynamics


We have developed a dynamic framework where the smooth limit surface evolves over time in response

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

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

to enforce various functional and aesthetic constraints. This direct manipulation is then transferred back

as virtual forces acting on the initial mesh through a transformation matrix (Eqn.18), and the initial

mesh (as well as the 1l. i1- i!i- smooth limit surface) deforms continuously over time until an equilibrium

position is obtained. The deformation of the surface in response to the applied forces is governed by

the motion equation (Eqn.17). \\ !iili our 1i1! -i, --based modeling framework, the limit surface evolves

as a consequence of the evolution of the initial mesh. One can apply various I I," of forces on the

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

obtained through evolution is constrained by the number of control vertices in the initial mesh. It might

be necessary to increase the number of control vertices in the initial mesh in order to obtain a detailed


August 20, 1998


DRAFT











shape through this evolution process.

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

replacing the initial mesh with a mesh obtained after one subdivision step applied to the initial mesh.

This new mesh has more number of vertices but defines the same limit surface. For example, after one

step of modified butterfly subdivision, the initial degrees of freedom p (refer to Eqn.11 and Eqn.12) in

the dynamic -1. I will be replaced by a larger number of degrees of freedom q, where q = Ap. A is

a global subdivision matrix of size (3s x 3t) whose entries are uniquely determined by the weights used

in the modified butterfly subdivision scheme (see Section II-A for the weights). Thus, p, expressed as a

function of q, can be written as

p = (ATA) -ATq = Aq, (19)

where At = (ATA) AAT. Therefore, we can rewrite Eqn.11 and Eqn.12 as

s(x) = (J(x)At)q, (20)

and

s(x,q) = (J(x)At)q, (21)

respectively. Now we derive the equation of motion for this new subdivided model involving a larger

number of control vertices namely q. We need to recompute the mass, damping and stiffness matrices

for this "-- i level. The structure of the motion equation as given by Eqn.17 remains unchanged, but

the .1i1 ,. 1!-..i! ,,1 and the entries of M, D, K, p and fp change correspondingly in this newly obtained

subdivided level. In particular the motion equation, explicitly expressed as a function of q, can be written

as

Mql + Dq Kqq = fq, (22)

where Mq = fxC, p(x)(At) JT(x)J(x)Atdx, S1 being the domain defined by the newly obtained sub-

divided mesh. The derivation of D,, Kq and fq follow suit.

It ii be noted that further increase in the number of control vertices, if necessary, can be obtained via

another level of subdivision. Therefore, multilevel dynamics is achieved through recursive subdivision on

the initial set of control vertices. Users can interactively choose ;,i- subdivided mesh as the control mesh

for the dynamic model depending on their needs. Alternatively, the -- -1. i can automatically determine


August 20, 1998


DRAFT











the most suitable control mesh for certain applications based on an application-specific criteria.


III. FINITE ELEMENT IMPLEMENTATION


(b)
Fig. 7. (a) An initial mesh, and (b) the corresponding limit surface. The domains of the shaded elements
in the limit surface are the corresponding triangular faces in the initial mesh. The encircled vertices
in (a) are the degrees of freedom for the corresponding element.


In the previous section, we have expressed the smooth limit surface as a function of the control vertex

positions in the initial mesh, and have assigned mass and damping distribution, internal deformation

. i. and forces to the limit surface to develop the corresponding ,1!. -i. 1 model. In this section, we

describe the implementation of this 1'1! -i, ,1 model using the finite element method.

In Section II we pointed out that the smooth limit surface obtained by the recursive application of the

modified butterfly subdivision rules can be represented by a set of smooth triangular patches, each of

which is represented by a finite element. The shape (basis) function of this finite element is obtained by

smoothing a hat function through repeated application of the modified butterfly subdivision rules. The


August 20, 1998


DRAFT











number of finite elements in the smooth limit surface is equal to the number of triangular faces in the initial

mesh as mentioned earlier (refer Fig.4 and 7). We now provide a detailed discussion on how to derive the

mass, damping and stiffness matrices for these elements. These elemental matrices can be assembled to

obtain the global '1!i -i1 ,1 matrices M, D and K, and a numerical solution to the governing second-order

differential equation as given by Eqn.17 can be obtained using finite element ; ,! I1 -i- techniques [57].

We use the same example as in Section II (refer I i, -) to develop the related concepts. The concept of

elements along with the control vertices and their corresponding domains is further illustrated in Fig.7.

We will now show how to derive the mass, damping and stiffness matrices for the element corresponding

to the triangular face abc in I i, '- The derivations hold for ;,i- element in general.


A. ./ ./. !,,/I mass and damping matrices

The mass matrix for the element given by Sabc (Eqn.9) can be written as

I = s p(x){Babc(x)} {Babc(x)}dx. (23)
JXE-Sabc
However, from Eqn.9 we know that the basis functions corresponding to the vertices in v b which are

stored as entries in Bab, are obtained as a limiting process. These basis functions do not have :i- 111 1- iG

form, hence computing the inner product of such basis functions as needed in Eqn.23 is a challenging

problem. In [55], an outline is provided on the computation of these inner products without performing ;, i

integration. In this paper, we develop a different yet simpler approach to solve this problem. The smooth

triangular patch in the limit surface corresponding to the face abc in the initial mesh is approximated by

a triangular mesh with 4J faces obtained after j levels of subdivision of the original triangular face abc

(each subdivision step splits one triangular face into 4 triangular faces). Then the mass matrix can be

expressed as
43
NI ,i = / (x){BJbc(x)}T{Bbc(x)}dx. (24)

The j-th level approximation of the corresponding basis functions can be explicitly evaluated (refer Eqn.8

for an expression of BIbc). An important point to note is that Eqn.8 involves several matrix multiplications

and hence can be very expensive to evaluate. However, the matrix (AJbc)(= (Aabc),- .. (Aabc)f(Aabc)1)

in Eqn.8 encodes how vertices in the 2-neighborhood of the triangular face uvw at level j are related to


August 20, 1998


DRAFT












the vertices in the 2-neighborhood of the triangular face abc in the initial mesh. In the implementation,

we keep track of how a new vertex is obtained from the contributing vertices in its immediate predecessor

level. If we move up from level j to level 0, we get the information stored in (AJbc) without forming iii

local subdivision matrices and thus avoiding subsequent matrix multiplications.

By choosing a -, L!i .I Ii high value of j, we achieve a reasonably good approximation of the elemental

mass matrices. We eliminate the computations involved in evaluating the integrals in 51!! 21 by using

discrete mass 1. i function which has non-zero values only at the vertex positions of the j-th subdivision

level mesh. Therefore, the approximation of the mass matrix for the element can be written as
k
NI = (v ){B b(vJIT ) {BJbc(vj)}, (25)
i=l
where k is the number of vertices in the triangular mesh with 4J faces. This approximation has been

found to be very effective and. tin. ii for implementation purposes. The elemental damping matrix can

be obtained in a similar fashion.


B. I./ .....I! 'i stiffness matrix

We assign the internal. ii. to each element in the limit surface, thereby defining the internal. -. i

of the smooth subdivision surface model. We take a similar approach as in the derivation of the elemental

mass and damping matrix and assign the internal. ,. i i to a j-th level approximation of the element.

In this paper, we use three I- ,. of internal i. - tension, stiffness and spring ii i _-. For the

examples used throughout the paper, this i at the j-th level of approximation can be written as

E = (E b t + (E +bct + (E b9p (26)

where the subscripts t, st and sp denote tension, stiffness and spring respectively.

The expression for the tension !. i which is essentially equivalent to the first order strain (membrane)

. [58], is
1 2
(E bc t 2 kt vI v

1 kt1{v bc T(KT b9 Vbc}, (27)


where kt is a constant, v3 and v the 1-th and m-th vertex in the j-th level mesh, are in the 1-

neighborhood of each other, Q is the domain defined by all such vertex pairs, and vibc is the concatenation


August 20, 1998


DRAFT











of the (x,y,z) positions of all the vertices in the j-th subdivision level of the triangular face abc in the

initial mesh.

Similarly, the expression for stiffness. -. i which is equivalent to the second order strain (thin plate)

* i [58], can be written as
1 2
(Ebct 2kt v 2v + vI
Q
= I bc tVJ}, kI (28)
2n ab, 1 abc)J V ab,
where v3, v3 and v3 are the three vertices of a triangular face. The summation involves three terms

corresponding to each triangular face, and is over all the triangular faces in the mesh at the j-th level of

subdivision.

I i!i !! we add a spring. !i. term which is given by
1 (k,-n),,( Iv v ,-- (v)
2

= iv bc (K c sp bc}, (29)

where (k,im),p is the spring constant, Q is the domain as in Eqn.27 and l is the natural length of the

spring connected between v and vJ. It 11 i be noted that the entries in (K c)sp depend on the distance

between the connected vertices and hence, (Kbc)s unlike other elemental matrices, is a function of time

which needs to be recomputed in each time step.

Combining the expressions for tension, stiffness and spring. i we get
jb1 kT ab)st (I(b {
E3 bc = 2Vbc} {kt(Kbc)t + kst(Kbc)st + (Kbc), p}{Vab}

= {v ab ba(Keb{v b

= {(A bc) IV T }} (K b)(A b){v

= { abc} {(ab) (Kab )(abc)} Vab,, (30)

where (Abc) and v are same as in Eqn.8. Therefore, the expression for the elemental stiffness matrix

is given by

Kabc = (Aibc) (Kibc)(Abc). (31)

It i i- be noted that the matrix multiplications for constructing Kabc are avoided in the implementation

by following the same technique described in Section III-A.


August 20, 1998


DRAFT











C. Force Application

The force f(x,t) in Eqn.18 represents the net effect of all externally applied forces. The current im-

plementation supports spring, inflation as well as image-based forces. However, other I- i" of forces like

repulsion forces, gravitational forces etc. can easily be implemented.

To apply spring forces, a spring of stiffness k can be connected from a point do to a point xo on the

limit surface (or, to the j-th level approximation mesh), the net applied spring force being

f(x, t) = k(do s(x, t))(x xo)dx, (32)

where 6 is the unit impulse function i!i'l, il!- f(xo,t) = k do s(xo, t) and vanishes elsewhere on the

surface. However, the 6 function can be replaced with a smooth kernel to spread the force over a greater

portion on the surface. The spring forces can be applied interactively using the computer mouse or the

points from which forces need to be applied can be read in from a t!l

To recover shapes from 3D image data, we synthesize image-based forces. A 3D edge detection is

performed on a volume data set using the 3D Monga-DT i i. II)) operator [59] to produce a 3D potential

field P(x,y,z), which we use as an external potential for the model. The force distribution is then

computed as

f (x, y,z) A P(x,,) (33)
I vP(x,y, z) /
where A controls the strength of the force. The applied force on each vertex at the j-th approximation level

is computed by trilinear interpolation for evaluating Eqn.18 in Cartesian coordinates. More sophisticated

image-based forces which incorporate region-based information such as gradients of a thresholded fuzzy

voxel classification can also be used to yield better and more accurate shape recovery. It i!i be noted

that we can apply spring forces in addition with the image-based forces by placing points on the boundary

of the region of interest in each slices of the 3D volume (MR, CT etc.) image data.


D. Discrete D,. ..... '', Equation

The t1! 1. i i I equation given by Eqn.17 is integrated through time by discretizing the time derivative

of p over time steps At. The state of the dynamic subdivision surface at time t + At is integrated using

prior states at time t and t- At. An implicit time integration method is used in the current implementation


August 20, 1998


DRAFT











where discrete derivatives of p are calculated using

p(t + At) 2p(t) + p(t At) (34)
pj(t + At) = ) 2(34)
At2

and

(t At) (t + At) p(t At)
p(t+Af)= (. -)
2At
The elemental mass, damping and stiffness matrices can be assembled to get the global mass, damping

and stiffness matrix for the smooth subdivision surface model. However, we do not assemble these global

sparse matrices explicitly for tin i reasons. For the time i stiffness matrix, we recompute K at

each time step. Using Eqn.17, Eqn.34 and Eqn.35, the discrete equation of motion is obtained as

(2M + DAt + 2At2K)p(t + At) = 2At2fp(t + At) + (DAt 2M)p(t At) + 4Mp(t). (36)

This linear -,. in of equations is solved iteratively between each time step using the '.. il_ .- gradient

method [60], [61].


E. Model Subdivision



evolution evolution




subdiv sion same limit surface
at equili rium with more patches





evol tion evol tion





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


The initialized model grows dynamically according to the equation of motion (Eqn.17). The degrees of

freedom of the initialized model is equal to the number of control vertices in the initial mesh as mentioned

earlier. When an equilibrium is achieved for the model, the number of control vertices can be increased by


August 20, 1998


DRAFT












replacing the original initial mesh by a new initial mesh obtained through one step of butterfly subdivision.

This increases the number of degrees of freedom to represent the same (deformed) smooth limit surface

and a new equilibrium position for the model can be obtained. This process is depicted schematically in

1 i, Model subdivision might be needed to obtain a very localized effect on a smooth limit surface. For

a shape recovery application, one 111 start with a very simple initial model, and when an approximate

shape is recovered, the degrees of freedom can be increased to obtain a new equilibrium position for the

model with a better fit to the given data set. The error of fit criteria for the discrete data is based on

distance between the data points and the points on the limit surface where the corresponding springs are

attached. In the context of image-based forces, if the model !! i does not change between successive

iterations indicating an equilibrium for the given resolution, the degrees of freedom for the model can be

increased by the above-mentioned replacement scheme until the model !! i is -itLil, ! small and the

change in model < i between successive iterations becomes less than a pre-specified tolerance.


IV. RESULTS

The proposed dynamic (modified) butterfly subdivision surface model can be used to represent a wide

,i. I- of smooth shapes with arbitrary genus. The smooth limit surface can be sculpted by ;1''1- i!_

synthesized forces in a direct and intuitive way in shape design applications. The [1!Il. i! i,_ shape in

a range or volume data set can also be recovered hierarchically using the proposed dynamic (modified)

butterfly subdivision surface model. Before illustrating the application results of the proposed model, we

would like to point out the advantages of our model over the existing techniques. It is known that shapes

of arbitrary lI ,.1. _- can be modeled either by explicit patching or by subdivision surfaces. C. .iJ ii[il-

constraints across patches need to be imposed with care while modeling arbitrary I, 1,. .1 _- using explicit

patching techniques. Subdivision surfaces do not face this problem, but one has to manipulate control

vertex positions at various levels of subdivision to obtain a desired smooth shape of arbitrary l 1 I ,1_ *- us-

ing subdivision surfaces. Our dynamic subdivision surface model overcomes this problem by allowing the

user to obtain the desired effect on the limit surface by direct sculpting via the application of synthesized

forces. To recover shapes from a given set of points in 3D, the existing subdivision surface based tech-


August 20, 1998


DRAFT





























(b) (c)


(f) (g)


Fig. 9. (a), (b), (c) and (d) : Initial shapes; (e), (f), (g) and (h) : the corresponding modified shapes
after interactive sculpting via force application.


niques resort to complex techniques to derive a mesh for the '1i.1 l i1 il shape, and then I- i.,, 11- mesh

optimization techniques are used to obtain a compact representation of the same. Our model recovers the

shape from a set of points in an t. !. ii hierarchical way ,i, simple mesh of the same I, 1. -. can be

used as an initial mesh which will evolve over time to fit the given data and depending on the error of

fit achieved, it will automatically refine itself until prescribed error of fit is obtained. In the rest of this

section, we elaborate on these points via examples.

In a shape modeling application, the user can specify :i!i mesh as the initial (control) mesh, and

the corresponding limit surface can be sculpted interactively by ;,1,1,1- -- -,lli .11 forces. In I i, '


August 20, 1998


DRAFT












we show several initial surfaces obtained from different control meshes and the corresponding modified

surfaces after interactive sculpting. To change the shape of an initial surface, the user can attach springs

from h !!. I. ,I points in 3D to the nearest points on the limit surface such that the limit surface deforms

towards these points to generate the desired shape. The initial mesh of Fig.9(a) has 125 faces and 76

vertices (degrees of freedom), the initial mesh of the closed cube-like shape in Fig.9(b) has 24 faces and

14 vertices. The one hole torus in Fig.9(c) and the corresponding modified shape in Fig.9(g) have initial

meshes with 64 faces and 32 vertices. A two hole torus with a control mesh of 272 faces and 134 vertices,

shown in Fig.9(d), is dynamically sculpted to the shape shown in Fig.9(h).

We have performed several experiments testing the applicability of our model to recover the I 1111. i l

shapes in range and volume data sets. In all the experiments, the initialized model has a control mesh

comprising of 24 triangular faces and 14 vertices whereas the control mesh of the fitted model has 384

triangular faces and 194 vertices. It i!! be noted that once an approximate shape is recovered, the model

is refined depending on the data-fitting criteria, thereby increasing the degrees of freedom of the recovered

shape only when necessary. For an error in fit (defined as the maximum distance between a data point

and the nearest point on the limit surface expressed as a percentage of the diameter of the smallest sphere

enclosing the object) of approximately 3%, the initialized model is refined twice following the technique

described in Section III-E. Also, the limit surface of i,!- control mesh (of the desired genus) can be used

as the initialized model. However, an initial mesh with few degrees of freedom usually performs better in

terms of recovering a compact representation of the [L.II. !- i ,_ shape.

In the first shape recovery experiment, we depict the laser range data acquired from multiple views

of a light bulb in Fig.10(a). Prior to ;,1I1'- ii,- our algorithm, the data were transformed into a single

reference coordinate -I. ii, The model was initialized inside the 1000 range data points on the surface

of the bulb. The fitted dynamic (modified) butterfly subdivision surface model is shown in Fig.10(b) and

(c). In the next experiment, the shape of a mechanical part is recovered from a range data set having

2031 data points (Fig.11). We also recover the shape of a human head from a range data set as shown in

Fig.12. The head range data set has 1779 points in 3D. It i!! be noted that the final shape with a very

low error tolerance is recovered using very few number of control points in comparison to the number of


August 20, 1998


DRAFT
























(b) (c)


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















(a) (b) (c)


Fig. 11. (a) Range data of a mechanical part along with the initialized model, (d) the fitted dynamic
butterfly subdivision model, and (c) visualization of the shape from another view point.















(a) (b) (c)


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


August 20, 1998


DRAFT
























(b) (c)


Fig. 13. (a) Data points i.1. iiflf i _- the boundary of the caudate nucleus on a -: 11: slice of human brain,
(b) data points (from all slices) in 3D along with the initialized model, and (c) the fitted dynamic
butterfly subdivision model.


data points present in the original range data set.

WVe present the recovery of the caudate nucleus (a cortical structure in human brain) from 64 1: Il

slices, each of size 256 x 256 in our next experiment. Fig.13(a) depicts a slice from this 11: 1 scan along

with a sparse set of points placed by an expert neuroscientist on the boundary of the shape of interest.

Fig.13(b) depicts the sparse data points (placed in each of the slices depicting the boundary of the shape

of interest) in 3D along with the initialized model. Note that points had to be placed on the boundary of

the caudate nucleus due to lack of image gradients delineating the caudate from the surrounding tissue in

parts of the image. Continuous image based forces as well as spring forces are applied to the model and

the model deforms under the influence of these forces until maximal conformation to the data is achieved.

The final fitted model is shown in Fig.13(c). We like to point out the fact that the recovered shape in [41]

using our previous dynamic subdivision surface model for the same data set has 386 degrees of freedom

and therefore, we achieve a factor of 2 improvement in the number of degrees of freedom required to

represent the model in this particular example.

In the last experiment, we animate the motion of the left-ventricular chamber of a canine heart over

a complete cardiac cycle. The data set comprised of 16 3D CT images, with each volume image having

118 slices of 128 x 128 pixels. in -I, we have recovered the shape from one data set using image-based

(gradient) as well as point-based forces. Once the shape is recovered from one data set, this fitted model


August 20, 1998


DRAFT


























(b) (c)


(f) (g)


(j) (k)


(I') (o)


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


August 20, 1998


DRAFT












is used as an initialization for the next data set to track the shape of interest. The snapshots from motion

tracking are shown in Fig.14 for the 16 volume data sets. It 11 i' be noted that the control mesh describing

the smooth surfaces shown in Fig.14 has only 384 triangular faces with a total of 194 vertices as mentioned

earlier. This experiment clearly shows that our model can be used to track a shape of interest from a set

of time dependent volume data sets in an tin i. manner. Note that no other existing purely geometric

subdivision surface technique can be used with (time i i,_i continuous data sets.


V. CONCLUSIONS

In this paper, we have presented a novel finite element method to derive and :,i ,1 the new dynamic

model based on the modified butterfly subdivision surface scheme. The new ,1! --based surface model

provides a direct and intuitive way of manipulating smooth shapes of arbitrary I. ..I i,1 and is very useful

for directly extracting and visualizing the ,I1ii,1 I1 i shapes in large range and volume data sets. The

proposed model has also been used successfully for non-rigid motion tracking from a temporal sequence of

volume data sets. We have developed a hierarchical local parameterization of the subdivision scheme which

is critical to the formulation of our dynamic model. We have combined material properties with geometric

entities, formulated the motion equations for our dynamic model, and incorporated the advantages of free-

form deformable models into conventional subdivision scheme. Moreover, we have introduced an tin. III

hierarchical dynamic control for various applications. Our experiments indicate a promising future of the

proposed model in computer graphics, geometric modeling and scientific visualization. Furthermore, the

finite element techniques proposed in this paper should be of great interest to the engineering design and

, i i 1- -. I ,~,~! , ii as w ell.


VI. ACKNOWLEDGMENTS


This research was supported in part by the NSF grant ECS-9210648 and the NIH grant RO1-LM05944

to B.C. Vemuri, the NSF ('.C\I:I..I award CCR-9702103 and DMI-9700129 to H. Qin. We wish to

acknowledge Dr. Tim Mi Ii. i. for his help, Dr. Gregoire Malandain for the edge detection software,

Dr. Hughes Hoppe and Dr. Kari Pulli for the range data, Dr. Dmitry Goldgof for the 4D CT images and

Dr. ('I! ,-1 i,, Leonard for brain I, 11:1 images.


August 20, 1998


DRAFT













REFERENCES

[I] G.M. ('h.i.! ii "An algorithm for high speed curve generation," Computer Vision, Graphics and Image

Processing, vol. 3, no. 4, pp. 346-349, 1974.

[2] D. Doo, "A subdivision algorithm for smoothing down irregularly shaped polyhedrons," in Proceedings on

Interactive Techniques in Computer Aided Design, 1978, pp. 157 -165.

[3] M. Sabin, !' use of piecewise forms for the numerical representation of shape, Ph.D. thesis, Hungarian

Academy of Sciences, Budapest, 1976.

[4] E. Catmull and J. ('!l. I i1... ,i - generated B-spline surfaces on arbitrary topological meshes," Computer

Aided Design, vol. 10, no. 6, pp. 350 -355, 1978.

[5] C. Loop, Smooth subdivision surfaces based on triangles, M.S. thesis, University of Utah, Department of

Mathematics, 1987.

[6] H. Hoppe, T. DeRose, T. Duchamp, 11. Halstead, H. Jin, J. McDonald, J. Schweitzer and WV. -Ii. ,.

I'l... --. smooth surface reconstruction," in Computer Graphics Proceedings, ACM SIGGRAPH, July

1994, Annual Conference Series, pp. 295 -302.

[7] M. Halstead, M. Kass and T. DeRose, I.!1. i. !, fair interpolation using Catmull-('!.l. I surfaces," in Computer

Graphics Proceedings, ACM SIGGRAPH, August 1993, Annual Conference Series, pp. 35 -44.

[8] J. Peters and U. Reif, "The simplest subdivision scheme for smoothing polyhedra," ACM .... .. on

Graphics, vol. 16, no. 4, pp. 420 -431, October 1997.

[9] T.W. Sederberg, J. Zheng, D. Sewell and M. Sabin, .'... -i!!..i ii. i i-ri e subdivision surfaces," in Computer

Graphics Proceedings, ACM SIGGRAPH, July 1998, Annual Conference Series, pp. 387 -394.

[10] T. DeRose, M. Kass and T. 1!i-.._- i 411,..1 i- ..ii surfaces in character animation," in Computer Graphics

Proceedings, ACM SIGGRAPH, July 1998, Annual Conference Series, pp. 85 -94.

[11] N. Dyn, D. Levin and J.A. Gregory, "A butterfly subdivision scheme for surface interpolation with tension

control," ACM I ... .. on Graphics, vol. 9, no. 2, pp. 160 -169, ApIl 1990.

[12] N. Dyn, S. Hed and D. Levin, -'iI.- i-t ..!i schemes for surface interpolation," in Workshop in Computational

Geometry, A. C. et al., Ed., pp. 97 -118. 1993.

[13] D. Zorin, P. Schrider and W. Sweldens, IiwI.. Ip'.l.!ii. subdivision for meshes with arbitrary topology," in

Computer Graphics Proceedings, ACM SIGGRAPH, August 1996, Annual Conference Series, pp. 189 -192.

[14] L. Kobbelt, I iI. I 1..I..!I.. refinement by variational methods," in Approximation .. VIII, C. ('Ciii and

L. Schumaker, Eds., vol. 2 of Wavelets and Multilevel Approximation, \\..il!I Scientific Publishing Co., 1995,

pp. 217 -224.


August 20, 1998


DRAFT













[15] L. Kobbelt, "A variational approach to subdivision," Computer-Aided Geometric Design, vol. 13, pp. 743 -

761, 1996.

[16] L. Kobbelt and P. Schr6der, "Constructing variationally optimal curves through subdivision," Tech. Rep.

CS- 1. -'.7-11. California Institute of Technology Computer Science Department Technical Report, 1997.

[17] D. Doo and M. Sabin, "Analysis of the behavior of recursive division surfaces near extraordinary points,"

Computer Aided Design, vol. 10, no. 6, pp. 356 -360, 1978.

[18] A.A. Ball and D.J.T. -',..i;, "Conditions for tangent plane continuity over recursively generated B-spline

surfaces," ACM I ... .. .. on Graphics, vol. 7, no. 2, pp. 83 -102, 1988.

[19] A.A. Ball and D.J.T. -'..i i, "An investigation of curvature variations over recursively generated B-spline

surfaces," ACM ... .. .. on Graphics, vol. 9, no. 4, pp. 424 -437, 1990.

[20] U. Reif, "A unified approach to subdivision algorithms near extraordinary points," Computer Aided Geometric

Design, vol. 12 no. 2, pp. 153 -174, 1995.

[21] J.E. Schweitzer, Analysis and Application of Subdivision Surfaces, Ph.D. thesis, University of \\.dhi,-! ..I

Seattle, 1996.

[22] A. il1.,il., and J. \\.! 1h, "Edge and vertex insertion for a class of C1 subdivision surfaces," Computer Aided

Geometric Design, to appear.

[23] J. Peters and U. Reif, "Analysis of generalized B-spline subdivision algorithms," SIAM Journal of .

Analysis, to appear, I, ..il..1I. at ftp://ftp.cs.purdue.edu/pub/jorg/9697agb.ps.Z.

[24] D. Zorin, ii....i l. !! of stationary subdivision on irregular meshes," Constructive Approximation, submit-

ted, available as -i .!1.. .I Computer Science Lab Tech. Rep. CSL-I1 i-'-- ._'

[25] J. -'.In1 I.... I evaluation of Catmull-('!., I subdivision surfaces at arbitrary parameter values," in Computer

Graphics Proceedings, ACM SIGGRAPH, July 1998, Annual Conference Series, pp. 395 -404.

[26] D. Terzopoulos, J. Platt, A. Barr and K. Fleischer, I.!.1.-1 ..17 deformable models," in Computer Graphics

Proceedings, ACM SIGGRAPH, 1987, Annual Conference Series, pp. 205 -214.

[27] D. Terzopoulos and K. Fleischer, "Deformable models," '. Visual Computer, vol. 4, no. 6, pp. 306 -331,

1988.

[28] A. Pentland and J. Williams, "Good l..iI It.i- : Modal dynamics for graphics and animation," in Computer

Graphics Proceedings, ACM SIGGRAPH, 1989, Annual Conference Series, pp. 215 -222.

[29] D. Metaxas and D. Terzopoulos, "Dynamic deformation of solid primitives with constraints," in Computer

Graphics Proceedings, ACM SIGGRAPH, July 1992, Annual Conference Series, pp. 309 -312.

[30] B.C. Vemuri and A. Radisavljevic, l .liii i..,i !.. stochastic hybrid shape models with fractal priors,"

ACM I ... .. on Graphics, vol. 13, no. 2, pp. 177 -207, Apil 1994.


August 20, 1998


DRAFT













[31] G. Celniker and D. Gossard, "Deformable curve and surface finite elements for free-form shape design," in

Computer Graphics Proceedings, ACM SIGGRAPH, July 1991, Annual Conference Series, pp. 257 -266.

[32] M.I.J. Bloor and M.J. \\ i,-..i, i,. p' ii ..; iiPDE surfaces in terms of B-splines," Computer Aided Design,

vol. 22, no. 6, pp. 324 -331, 1990.

[33] M.I.J. Bloor and M.J. \\ !,-..,1 I -i_, partial
Aided Design, vol. 22, no. 4, pp. 202 -212, 1990.

[34] G. ('. ii,! i. and W. \\ !. i, "Linear constraints for deformable B-spline surfaces," in Proceedings of the

Symposium on Interactive 3D Graphics, pp. 165 -170. ACM, .'.. York, 1992.

[35] W. \\. !. i and A. Witkin, \ .I1 .Ii ...l surface modeling," in Computer Graphics Proceedings, ACM SIG-

GRAPH, July 1992, Annual Conference Series, pp. 157 -166.

[36] H. Qin and D. Terzopoulos, "Dynamic NURBS swung surfaces for physics-based shape design," Computer-

Aided Design, vol. 27, no. 2, pp. 111-127, 1995.

[37] H. Qin and D. Terzopoulos, "D-NURBS : A physics-based framework for geometric design," IEEE .

tions on Visualization and Computer Graphics, vol. 2, no. 1, pp. 85-96, March 1996.

[38] D. Terzopoulos and H. Qin, "Dynamic NURBS with geometric constraints for interactive sculpting," ACM

I .. .. ...* on Graphics, vol. 13, no. 2, pp. 103 -136, Apil 1994.

[39] C. Mandal, H. Qin and B.C. Vemuri, "Dynamic smooth subdivision surfaces for data visualization," in IEEE

Visualization'97 Conference Proceedings, Phoenix,AZ, October 1997, pp. 371 -377.

[40] H. Qin, C. Mandal and B.C. Vemuri, "Dynamic Catmull-('!.,! I subdivision surfaces," IEEE ... .. on

Visualization and Computer Graphics, vol. 4, no. 3, September 1998, to appear.

[41] C. Mandal, B.C. Vemuri and H. Qin, -'I.ip. recovery using dynamic subdivision surfaces," in Proceedings

of the International Conference on Computer Vision, Bombay, India, January 1998, pp. 805 -810.

[42] E. Bardinet, L.D. Cohen and N. Ayache, -'1q,. !i.p. Il! and free-form deformations : A global model to fit

and track 3D medical data," in Lecture .,-' in Computer Science, CVRMed'95, N. Ayache, Ed., Berlin,

Germany, 1995, vol. 905, pp. 319 -326, Springer Verlag.

[43] L.D. Cohen and I. Cohen, I ini-. -, !w i !, methods for active contour models and balloons for 2D and 3D

images," IEEE .. .. . on Pattern Analysis and Machine Intelligence, vol. 15, no. 11, pp. 1131 -1147,

.. -. i. i 1993.

[44] E. Koh, D. Metaxas and N. Badler, 1 I i..,. I ..! shape representation using locally adaptive finite elements,"

in Lecture ...' -- in computer Science, Computer Vision ECCV 1994, J.O..! .1, i.,I, Ed., Berlin, Germany,

1994, vol. 800, pp. 441 -446, Springer-Verlag.

[45] F. Leitner and P. Cinquin, "Complex topology 3D objects segmentation," in Model-based Vision Development


August 20, 1998


DRAFT













and Tools, ". !.' Proceedings, Bellingham, NA, 1991, SPIE, vol. 1609, pp. 16 26.

[46] T. McInerney and D. Terzopoulos, "A dynamic finite element surface model for segmentation and tracking

in multidimensional medical images with application to cardiac 4d image analysis," Computerized Medical

Imaging and Graphics, vol. 19, no. 1, pp. 69 -83, 1995.

[47] L.H. -'.iil, and J.S. Duncan, I 1..,ii ..I finding with parametrically deformable models," IEEE ... ..

on Pattern Analysis and Machine Intelligence, vol. 14, no. 11, pp. 1061 -1075, 1992.

[48] D. Terzopoulos and D. Metaxas, "Dynamic 3D models with local and global deformations: Deformable

superquadrics," IEEE ... .. on Pattern Analysis and Machine Intelligence, vol. 13, no. 7, pp. 703

714, July 1991.

[49] D. Metaxas and D. Terzopoulos, -'!..p'. and non-rigid motion estimation through physics-based synthesis,"

IEEE .... .. on Pattern Analysis and Machine Intelligence, vol. 15, no. 6, pp. 580 -591, June 1993.

[50] A. Pentland and B. Horowitz, I...... of non-rigid motion and structure," IEEE ... ... on Pattern

Analysis and Machine Intelligence, vol. 13, no. 7, pp. 730 -742, July 1991.

[51] Y. ('!I. i and G. Medioni, -'1I .. .. description of complex objects from multiple range images," in Proceedings

of the Conference on Computer Vision and Pattern Recognition, Seattle, A, June 1994, pp. 153 -158.

[52] W.C. Huang and D. Goldgof, "Adaptive-size physically-based models for non-rigid motion analysis," in

Proceedings of the Conference on Computer Vision and Pattern Recognition, Urbana-('!i.. nipl,.i, IL, June

1992, pp. 833 -835.

[53] H. Tanaka and F. Kishino, "Adaptive mesh generation for surface reconstruction : Parallel hierarchical

triangulation without discontinuities," in Proceedings of the Conference on Computer Vision and Pattern

Recognition, .'.. York City, NY, June 1993, pp. 88 -94.

[54] M. N .i-i i and D. Terzopoulos, "Adaptive meshes and shells : Irregular triangulation, discontinuities and

hierarchical subdivision," in Proceedings of the Conference on Computer Vision and Pattern Recognition,

Urbana-('C!..1I,,1. .,11, IL, June 1992, pp. 829 -832.

[55] J.M. Lounsbery, T. DeRose and J. \\.i, l. i, II !,! i.. !ii..,i analysis for surfaces of arbitrary topological

type," ACM i .. ....* on Graphics, vol. 16, no. 1, pp. 34 -73, January 1997.

[56] B.R. Gossick, Hamilton's Principle and Physical Systems, Academic Press, .'`. York, 1967.

[57] H. Kardestuncer, !'. /..'.. Handbook, McGraw-Hill, .'. York, 1987.

[58] D. Terzopoulos, Ih. !..ll. '.oii .'. of inverse visual problems involving discontinuities," IEEE .. .. on

Pattern Analysis and Machine Intelligence, vol. 8, no. 4, pp. 413 -424, 1986.

[59] 0. Monga and R. Deriche, "3D edge detection using recursive till. i i.- in Proceedings of IEEE Conference

on Computer Vision and Pattern Recognition, June 1989, pp. 28 -35.


August 20, 1998


DRAFT













[60] G.H. Golub and V.H. Van Loan, Matrix Computations, The Johns Hopkins University Press, 1989.

[61] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, .. .. Recipes in C, Cambridge Uni-

versity Press, 1992.


August 20, 1998


DRAFT