Document Sample
Jerry Powered By Docstoc
					Fitting Solid Meshes to Animated Surfaces using
Linear Elasticity
Jaeil Choi1 and Andrzej Szymczak2

Computing correspondence between time frames of a time-dependent 3D surface is essential for
the understanding of its motion and deformation. In particular, it can be a useful tool in compres-
sion, editing, texturing or analysis of the physical or structural properties of deforming objects.
However, correspondence information is not trivial to obtain for experimentally acquired 3D ani-
mations, such as time dependent visual hulls (typically represented as either a binary occupancy
grid or as a sequence of meshes of varying connectivity).
   In this paper we present a new non-rigid fitting method that can compute such correspondence
information for objects that do not undergo large volume or topological changes, such as living
creatures. Experimental results show that it is robust enough to handle visual hull data, allowing
one to convert it into a constant connectivity mesh with vertices moving in time. Our procedure
first creates a rest state mesh from one of the input frames. That rest state mesh is then fitted to
the consecutive frames. We do that by iteratively displacing its vertices so that a combination of
surface distance and elastic potential energy is minimized. A novel rotation compensation method
enables us to obtain high quality results with linear elasticity, even in presence of significant
Categories and Subject Descriptors: I.3.5 [Computer Graphics]: Computational Geometry and
Object Modeling; I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism; I.4.8
[Image Processing and Computer Vision]: Scene Analysis
General Terms: time-dependent surfaces, tracking, elasticity
Additional Key Words and Phrases: deformation, finite element methods, fitting

In recent years, there has been rapid progress in the area of 3D scanning and
reconstruction of static geometry from scanned data. 3D scanning technologies have
found numerous applications in medicine, art, education and entertainment. At the
same time, the need to reconstruct, process and analyze time-dependent geometry
has been growing. Such tasks often require knowing how parts of the object move
and deform over time, or the ability to track individual points through time. This
problem has been attacked in many ways: by putting markers on the objects (as in

1 GVU  center, Georgia Institute of Technology, TSRB 85th 5th Street NW, Atlanta, GA 30324;
2 Department of Mathematical & Computer Sciences, Colorado School of Mines, Golden, CO


Permission to make digital/hard copy of all or part of this material without fee for personal
or classroom use provided that the copies are not made or distributed for profit or commercial
advantage, the ACM copyright/server notice, the title of the publication, and its date appear, and
notice is given that copying is by permission of the ACM, Inc. To copy otherwise, to republish,
to post on servers, or to redistribute to lists requires prior specific permission and/or a fee.
 c 20YY ACM 0000-0000/20YY/0000-0001 $5.00

                                           ACM Journal Name, Vol. V, No. N, Month 20YY, Pages 1–16.
2    ·     J. Choi and A. Szymczak

motion capture), finding features in each frame and matching them across frames
or fitting an existing model to consecutive time frames. All these approaches are
under extensive research in a variety of contexts. The approach of fitting a model
to observed data has been used in 2D and 3D applications in computer vision
and medical image registration, but previous approaches were restricted only to
small deformations, such as facial expressions or left ventricular motion in cardiac
   This paper describes a new non-rigid fitting approach. We think of the deforming
object as an elastic body, and fit a rest state mesh of fixed connectivity (constructed
from one of the input frames) to the consecutive frames, minimizing a combination
of the surface distance and elastic potential energy for each frame (Figure 1). We
use rotation compensation in order to be able to handle globally large deformations
while using linear elasticity. In contrast to [M¨ ller et al. 2002], we perform rotation
compensation on per-tetrahedron basis rather than on per-vertex basis. This makes
the procedure simlper and naturally supports computing the elastic potential energy
by summing the contributions from each tetrahedron in the mesh (Section 5.2).
   Our procedure can establish physically meaningful correspondence between all
frames of the animation, allowing to track individual points in time (Fidget 2).
Since we use a fixed connectivity mesh and model deformation using elasticity,
our approach requires two fundamental assumptions : (1) The underlying physical
object does not change topology, (2) The volume of that object does not change
too much. These assumptions are not valid for arbitrary deformations, but they
do hold for motions of living creatures, including humans. Let us stress that the
assumptions need not hold for the input 3D animation (in fact, they do not hold
for typical visual hull data), but only for the underlying physical object.
   Using a solid model (triangular mesh in 2D and tetrahedral mesh in 3D) and
elasticity simulation has many advantages. Along with the increased stability and
precision of continuum, it allows to control and measure the material properties of
the mesh elements. The representation provided by our fitting procedure can make
certain animation processing tasks such as compression, network transmission or
streaming, texturing or editing much easier than other, less structured represen-
tations like binary volumes or sequences of meshes whose connectivity changes in

Our approach is inspired by prior work on physically based deformation and energy
minimizing fitting methods in computer vision. Physically based deformation was
introduced to computer graphics in the pioneering papers of [Terzopoulos et al.
1987; Terzopoulos and Witkin 1988] and, since then, has been an active research
topic. In this paper, we use the classical linear elasticity formulation with enhance-
ments allowing to handle large rotations inspired by the stiffness warping method
described in [M¨ ller et al. 2002]. We postpone a more detailed discussion until
section 5.2.
  The approach of fitting a physically based deformable model to recover 3D shape
and nonrigid motion from volumetric data without any prior knowledge of structure
has been used in computer vision for a long time [Terzopoulos et al. 1988; Petland
ACM Journal Name, Vol. V, No. N, Month 20YY.
                                        Fitting Solid Meshes using Linear Elasticity       ·    3

Fig. 1. The result of fitting a tetrahedral mesh to a sequence of visual hull surfaces (shown as
transparent orange surfaces). The initial mesh (bottom row, left) has been constructed from one
of the frames that was chosen by user. Stills from four of the twelve videos used to compute the
visual hulls are shown in the top row.

                   (a) original shape   (b) target shape      (c) fitted result

Fig. 2. Our fitting approach can be used to find a complete correspondence between two shapes.
The texture of the original shape was transfered to the target shape using linear interpolation of
the 2D triangular mesh created by our method. In this rather extreme example, we imposed an
upper bound on the distance between a vertex and its target (section 5.1).

and Horowitz 1991; McInerney and Terzopoulos 1993]. With a lower dimensional
model, such as thin-plate surface with tension, this approach was successfully ap-
plied to small deformations, like facial expressions or left-venticular motion in car-
diac imaging. Both in dynamic simulation setting and static optimization setting,
these fitting procedures are driven by three components : topological constraints,
internal forces obtained from energy defined by deformation potentials, and feature
conformity with the data. Our fitting approach follows the same principles, but is
based on solid deformable model instead of surface model, and designed for much
more complex deformations that involve large rotations.
   Several recent research results aim to achieve goals similar to ours: to build
correspondence between frames of a 3D animation or to track features through time.
The method described in [Starck et al. 2002] allows to fit constrained surface model
of triangle meshes to time-dependent surface data. A skeleton-based approach
aiming to track points and features through time was introduced in [Brostow et al.
2004]. Interesting results of combining a skeleton computed from motion capture
and profiles from multiview video to reconstruct detailed surface geometry of a
                                                  ACM Journal Name, Vol. V, No. N, Month 20YY.
4    ·     J. Choi and A. Szymczak

moving human subject was described in [Sand et al. 2003]. However, all of the work
we are aware of are based on lower dimensional representations of the deforming
object. While lower dimensional representations have numerous advantages (like
computational efficiency and simplicity) and are sufficient in many applications, we
believe that volume based approaches will ultimately turn out to be more complete,
accurate and stable, since it can capture the interaction through interior material
under the surface. Elasticity simulated on tetrahedral meshes also provides an
elegant way to keep the local volume changes of the fitted tetrahedral mesh under
control, which is consistent with the way most real objects deform. As a result,
it allows us to deal with noise and artifacts commonly present in computer vision
data (like visual hull data [Matusik et al. 2000]) effectively.

The algorithm proposed in this paper takes a 3D animation as its input. By a 3D
animation we mean a sequence of 3D watertight surfaces (which we call frames).
Generally, these surfaces can be represented in any possible way, for example using
a binary voxel grid or a boundary representation (watertight triangle mesh). At
startup, we convert the frames into signed distance volumes which are more suitable
for our purposes. In particular, we use fast marching method [Sethian 1996] for
the calulation of the signed distance field. A typical application scenario is to use
visual hull data or other data reconstructed from multiview video as the input 3D
animation. As the output, our algorithm generates a different representation of the
input 3D animation: a fixed connectivity tetrahedral mesh whose vertices move in
  Our procedure works by first constructing the rest state mesh from one of the vol-
umes in the input sequence (ideally, one of good quality). This process is described
in Section 4. Then, this mesh is fitted to consecutive frames, both forward and
backward in time, by iteratively minimizing a combination of the surface distance
and the elastic potential energy for each frame. The details of the fitting procedure
are described in 5. Finally, we describe experimental results on synthesized and
real-world 3D motion sequences in section 6, and discuss future research directions
in section 7.

Like all finite element methods, our algorithm requires a high quality rest state
mesh in order to produce useful results in an efficient manner. The rest state mesh
is created automatically from a user-selected frame. Clearly, this frame should not
contain any large-scale topological artifacts, like two limbs meeting and forming
a topologically nontrivial loop. Expectedly, while most topologically clean frames
tend to work well for this purpose, the ones containing no large-scale geometric
artifacts and no extreme deformations (i.e. the one that is similar to the ‘average
pose’) lead to best results.
   Creating finite element meshes suitable for physical simulation has been a topic
of extensive research for a number of years. The major approaches to 3D mesh gen-
eration include the Delaunay triangulation and refinement [Shewchuk 1998; Edels-
brunner and Guoy 2002], advancing front [L¨hner 1988; Rassineux 1998; Sh¨berl o
ACM Journal Name, Vol. V, No. N, Month 20YY.
                                      Fitting Solid Meshes using Linear Elasticity   ·    5

1997] and structured grid based approaches [Yerry and Shephard 1984; Nielson and
Sung 1997]. In addition, tetrahedral meshes created from these approaches can be
further improved by the combination of small simplification and refinement tech-
niques [Cutler et al. 2004]. After comparing several approaches, we found that the
structured grid algorithm described in [Molino et al. 2003], using Body-Centered
Cubic lattice, works best for our purposes. Its simplicity and ability to retriangu-
late the input surface as well as being based on a structured grid seems to give it an
advantage over other algorithms for noisy input surfaces. Clearly, this advantage
comes at the price of losing some of surface details (because of the retriangulation),
but we were always able to recover the shape correctly on the scale we are interested

           (a) input surface          (b) initial lattice         (c) subdivision

         (d) exterior removal          (e) compression            (f) result mesh

          Fig. 3.   Generation of Finite Element Mesh. (b)-(e) are cut-away views.

   In order to build the rest state mesh, we followed the tetrahedralization proce-
dure of [Molino et al. 2003], illustrated in Figure 3. First, the domain space of
input surface (a) is divided using an uniform BCC lattice (b), and then the lattice
is adaptively subdivided (c) by “red”/“green” subdivision. After the adaptive sub-
division, exterior tetrahedra are removed (d), and the rest is ‘compressed’ to fit the
given input surface (e). In the compression step, the distribution of interior vertices
are optimized to improve the quality (aspect ratio of tetrahedra) of the mesh using
simple pattern searching approach.

This section describes the fitting procedure in more detail. We assume that the rest
state mesh R is available and it has been fitted to f -th frame. We shall denote this
fitted mesh by Mf . The goal is to fit this mesh to the (f + 1)-th frame (obtaining
the mesh Mf +1 ) by iteratively moving its vertices. We do that by treating Mf
                                                ACM Journal Name, Vol. V, No. N, Month 20YY.
6    ·     J. Choi and A. Szymczak

                                0                                             i
as an initial approximation Mf +1 of Mf +1 . Having an approximation Mf +1 , we
                           i+1                                  i
construct a better one Mf +1 by displacing the vertices of Mf +1 so that a certain
objective function is minimized. This process is repeated until convergence, yielding
the Mf +1 . The objective function (determined by R and Mf +1 ) is a combination of
the surface distance term and the elastic potential energy term, which are described
below. Both terms are quadratic functions of the coordinates of the vertices of Mf +1
(in particular, this means there are 3k variables, where k is the number of vertices
in the rest state mesh). Recall that all meshes we deal with here have the same
connectivity: they differ only in vertex coordinates.

5.1 Surface Distance Term
Since our goal is to make Mf +1 approximate the surface given by the (f + 1)-th
frame (which we denote by Sf +1 ), it is natural to introduce a term that grows
together with the distance between the unknown mesh Mf +1 and Sf +1 into the
objective function.
   Let B be the set of boundary vertices of Mf +1 . Consider a vertex v ∈ B. We
first compute its target tv i.e. the point on Sf +1 that is closest to the (known)
vertex of Mf +1 corresponding to v. With the signed distance function d from Sf +1
available, one can approximate the target of v by following the gradient of the
distance function if d(v) < 0 and following it backward for d(v) > 0.
   We would like to form the distance term D by combining the squared distances
between vertices in B and their targets, i.e. use

                               D=          wv dist2 (v, tv ).
                                                  v                              (1)

In order to make D approximate the distance to the surface Sf +1 rather than just
a combination of vertex-target distances, we use the anisotropic squared distance
function dist2 that penalizes movement in a tangent direction less than movement
along the normal direction:
                                             2      
                                              s 0 0
                dist2 (v) = (Rv (v − tv ))T  0 s2 0  Rv (v − tv ),
                                              0 0 1

where Rv is a rotation that maps the normal vector to Sf +1 at tv into a vector
pointing along the z-axis and s ∈ [0, 1] defines the amount of anisotropy. Ideally,
one could vary the anisotropy according to the differential properties of the input
surface. However, we found out that estimating them reliably for noisy data we
would like to use as the input to our algorithm is hard and therefore we simply
use s = 0.1 in all the experiments described in this paper. This creates penalties
along tangential direction that are a tenth of the penalties along the surface nor-
mal direction, which subsequently allows the fitted model to slide along the target
  Since our input animations often contain relatively large motions between con-
secutive frames, the contributions of vertices of B to the surface distance term need
to be weighted properly. This is done by varying the wv coefficients (which can be
ACM Journal Name, Vol. V, No. N, Month 20YY.
                                        Fitting Solid Meshes using Linear Elasticity       ·    7

Fig. 4. The signed distance field is shown in red (positive) and blue (negative). Blue triangles
represent the elements in the mesh that is being fitted to the isosurface of the signed distance
function corresponding to isovalue of zero. Red, blue, and cyan arrows are vertex normals, gra-
dients of distance function, and vectors pointing toward the target, respectively. Notice that the
targets of v0 and v1 are far away from the place that v0 and v1 would go to if the blue mesh was
deformed into the isosurface. Our procedure assigns zero confidence to these vertices.

thought of as confidence measures) from vertex to vertex. We define wv by
                        wv = max(0, cos θ) = max(0, nv · ∇d(v)),
where θ be the angle between vertex normal nv and gradient of distance function
d at v. This is illustrated in Figure 4.
5.2 Elastic Potential Energy Term
Assume we are given an elastic body modeled using a tetrahedral mesh and that
we apply a deformation to it by displacing the vertices. Let x0 and x be the 3k-
dimensional (where k is the number of vertices) vectors obtained by concatenating
the coordinates of all vertices in the rest state and after a deformation is applied
(respectively). Linear elasticity computes the vector of forces acting on the vertices
as a result of the deformation by
                                       F = K(x − x0 ),                                         (2)
where K is a 3k × 3k stiffness matrix defined by K =      (where u = x − x0 ), which
can be built from element-wise strain-stress relationships. The elastic potential
energy is given by
                               E=   (x − x0 )T K (x − x0 ).                            (3)
  Since linear elasticity is based on linearization of stress/strain tensors, it is useful
as a basis for efficient simulation methods. However, it is not accurate enough to
be suitable for large deformations (in particular, it is not invariant under rotations)
which makes it inadequate in most engineering applications. In graphics, where
interactivity and efficiency are of higher value, it has found wider use, but its
problems have also been well recognized [Zhuang and Canny 1999]. In an attempt
to alleviate the disadvantages of linear elasticity, [M¨ ller et al. 2002] introduced the
                                                  ACM Journal Name, Vol. V, No. N, Month 20YY.
8    ·     J. Choi and A. Szymczak

              initial state                1st iteration         after convergence

Fig. 5. Compensating rotations during fitting: Gray mesh represents rest state mesh (x0 ), blue
represents the current state (x), and the green triangles are rotated reference elements (xr ).

stiffness warping method to deal with lack of rotation invariance in linear elasticity.
Their approach can also be viewed as separating deformable and rigid components
as in [Terzopoulos and Witkin 1988], independently for each mesh element. Our
procedure is motivated by these results and proceeds as follows. We calculate pure
rotation components of the transformation from the tetrahedra in the rest state
into the corresponding tetrahedra of the currently available approximation Mf +1 .
These rotations are used to create a set of imaginary new reference elements that
reflect only element-wise rotations of current state and cannot be integrated into
a consistent mesh. From these reference elements, we can calculate the rotation
compensated elastic potential energy term for use in our procedure.
   The idea is illustrated in Figure 5. The green triangles in the figure represent
rotated new reference elements at particular iteration step. Note that the new
reference elements does not form a consistent mesh. The pure rotation component
of a transformation of a tetrahedron T into a tetrahedron T ′ is calculated in the
same way as in [Alexa et al. 2000]. Let A be the transformation that takes each
vertex of T into the corresponding vertex of T ′ . Let A = U SV T be the singular
value decomposition of A. We compute the pure rotation component as RT = U V T .
   The elastic potential energy is computed by adding up elastic potential energies
of individual tetrahedra, computed using Equation (3). However, we apply the
rotation transformation associated with a tetrahedron T to all of its vertices in the
rest state before computing the displacement vector (x− x0 ) to be used in Equation
   Suppose that the rest state mesh has m tetrahedra T1 , T2 , . . . , Tm and n vertices
v1 , v2 , . . . , vn . For a tetrahedron Tj , let x0 (xj ) be the 12-dimensional column
vector obtained by concatenating the coordinates of the vertices of Tj in the rest
state mesh (respectively, the variables corresponding to the unknown locations of
vertices of Mf +1 ). Let Rj be the rotation component of the mapping of Tj in
                                                             i         ¯
the rest state mesh into the same tetrahedron in Mf +1 . Let Rj be the 12 × 12
block diagonal matrix having four 3 × 3 blocks Rj along the diagonal. The elastic
potential energy can be re-written as
                                    1       ¯ j      ¯        ¯ j
                     E=               (xj − Rj x0 )T Kj (xj − Rj x0 ) ,                    (4)

      ¯                                                   ¯
where Kj is the stiffness matrix for the tetrahedron Tj at Rj x0 . Note that E is
ACM Journal Name, Vol. V, No. N, Month 20YY.
                                         Fitting Solid Meshes using Linear Elasticity       ·     9

                 initial state               1st iteration                   2nd

                        5th                  10th             30th          150th

Fig. 6. Iterative fitting result for large rotation. Slow convergence is due to the fact that we only
use local information, and such information is very limited in this case because of the dramatic
difference between input surface and initial mesh.

a quadratic expression in the unknowns (which appear as coordinates in the xj

5.3 Iterative Optimization for Fitting
Using equations (1) and (4), we define the objective function as a combination of
the distance and elastic potential energy terms:

                              E = αD + βE =: xT Ax − 2bT x + c,

where A, b, and c represent its quadratic, linear, and scalar terms, respectively,
and x is the vector of coordinates of the vertices of the unknown mesh Mf +1 . This
objective function E is minimized when Ax = b. Since A is sparse, symmetric and
positive definite, this linear system can be solved efficiently using the Conjugate
Gradient method. Note that the part of A corresponding to the elastic potential
energy (E) itself is singular, because strain/stress tensors are invariant under rigid
transformation. However, in all practical cases, the distance term (D) makes our
optimization problem non-degenerate, since it penalizes any translation of a vertex
with nonzero confidence away from the target.
   This simple optimization of quadratic function is iterated until convergence: after
Mf +1 is computed, we recompute A and b (recall that they depend on the rest
state mesh and the current approximation of Mf +1 ) and solve the linear system
Ax = b to obtain the locations of vertices of Mf +1 and so on until the difference
between the coordinates of the current and previous approximations is below a
certain threshold. At this point, we terminate the process and declare the current
approximation as Mf +1 . In all our experiments on 3D data discussed in Section
6, the optimization is set to terminate after 10 iterations or when the reduction of
surface distance term is less than 0.1% of initial value (whichever happens first).
For the majority of frames, only 6 − 7 iterations are performed.
   A few stages of this iterative process for a simple (although requiring an abnor-
                                                    ACM Journal Name, Vol. V, No. N, Month 20YY.
10     ·     J. Choi and A. Szymczak

                                 (a) without collision detection

                                   (b) with collision detection

Fig. 7. Collision detection. The confidence of the target association of the vertices on both body
and the hoof is decreased to zero due to the collision (indicated by red arrows).

mally large number of iterations to converge because of dramatic shape difference)
2D example are shown in Figure 6. Note that these results are similar to those
from [Alexa et al. 2000]. While their method creates continuous interpolation of
deformation when full correspondence is given, our method calculate the optimized
correspondence between two shapes when local deformations are not too big.

5.4 Avoiding collisions
If the difference between consecutive frames is small enough, the association of a
target vertex tv with a boundary vertex v described in Section 5.1 constitutes a
reasonable approximation of the correspondence between points on the surfaces
provided by the two frames. However, in typical data, the motion that occurs
between two consecutive frames may be relatively large when compared with the
size of the surface features. As a result, the association of targets to boundary
vertices v may lose continuity. This could have disastrous effect on the result of
our procedure. An example is shown in Figure 7 (a): the targets for vertices on the
horse’s leg are attracted to the body which causes the leg to fold over itself.
   We deal with this problem by preventing vertices which are far away from each
other in terms of the geodesic distance on the boundary of the rest state mesh
from being moved toward close targets. First, we estimate the geodesic distance
between all pairs of external vertices of the rest state mesh. We will call pairs of
vertices {v, v ′ } that are more than a user-defined distance d away from each other
distant. For each distant pair {v, v ′ } whose targets are close to each other (closer
than user-specified distance T ), we ‘invalidate’ the targets by setting the confidence
values wv and wv′ to zero. We found out that the procedure works well for a wide
range of values of d and T . It also does not require accurate distance computation:
in all examples discussed in this paper we use the least number of edges of a path
connecting two vertices u and w as the estimate of the geodesic distance between u
and v, thus the choice of d does not depend on the metric scale of the model or its
resolution. We consider two external vertices are ‘far away’, if the number of edges
between them is bigger than 20. As T , we use the average length of an external
ACM Journal Name, Vol. V, No. N, Month 20YY.
                                       Fitting Solid Meshes using Linear Elasticity      ·     11

                                    model size      rotated    setting eq.    solving eq.     total
                                     (vert/tet)        ¯
                                                       K        (sec(%))       (sec(%))       (sec)
                                       original       Yes      0.05 (71.4)    0.02 (28.6)     0.07
                                      (84/121)        No       0.02 (50.0)    0.01 (25.0)     0.04
                                   subdivided 1       Yes      0.27 (37.5)    0.44 (61.1)     0.72
                                    (370/968)         No       0.07 (10.0)    0.60 (85.7)     0.70
                                   subdivided 2       Yes      1.86 (15.1)   10.42 (84.5)     12.34
                                   (2035/7744)        No       0.57 (3.5)    15.58 (96.2)     16.2
                                   subdivided 3       Yes      15.15 (5.5)   260.26 (94.5)   275.41
                                  (13125/61952)       No       4.57 (1.2)    376.52 (98.8)   381.09

Fig. 8. Comparison between the performances with/without the rotation compensation, applied
to the problem of fitting a simple 3D bar to a bended shape. Left image shows (from left to right)
the result of rotation-compensated fitting (5 iterations), the result after the first iteration, and
the result without rotation-compensation (5 iterations).

edge in the rest state mesh.
   Intuitively, the above approach allows the elastic potential energy term to dom-
inate over the distance term in areas where the the target positions deduced from
the distance field are not correct. As the optimization iterates, the surface being
fitted gets closer to the target surface. Consequently, the vertex-target association
improves, the number of invalidated vertices tends to decrease fast and the fitted
surface settles into the correct result as shown in Figure 7 (b). In our experiment
data sets, most of invalid vertex-target associations were resolved in less than 3
iterations. Note that the heuristic described above comes without a guarantee: the
time-dependent surface produced by our algorithm may still have self-intersections.
In most cases, this happens if different body parts make contacts (or get close
to each other). We do not attempt to simulate dynamics of contact points and
calculate the hidden surface.

We implemented our fitting approach in C++ using a Conjugate Gradient solver
on a sparse matrix representation, and have conducted a simple experiment on 2.8
GHz, 2GB RAM Linux desktop for the problem of fitting a simple 3D bar to a
bended shape. The result is shown in Figure 8.
   We believe our result without the rotation compensation is comparable to that
of standard deformation of FEM with linear elasticity (except for the surface dis-
tance term), even though our implementation left many general optimization issues
regarding linear system solver untouched. As we can see from the result, the com-
putational cost of rotation compensation is linear to the number of the elements in
the model, and does not increase the total cost (“setting eq.” in the table of Figure
8) significantly. On the other hand, deforming a model without rotation compen-
sation not only produces incorrect fitted results, but also degrades the convergence
of the Conjugate Gradient solver.
   For more realistic situation, we also have tested our system on the following
three-dimensional datasets:
Camel :      A visual hull dataset obtained from 30fps video sequences of a moving
                                                  ACM Journal Name, Vol. V, No. N, Month 20YY.
12    ·       J. Choi and A. Szymczak

 dataset     # of frames       Signed Distance Field                  initial tet. mesh
                             (avg. size / sec per frame)   (# of vert. / # of tet. / sec to build)
  Camel          340            (184×162×140) / 0                     3908 / 16143 / 0
  Horse           48             (70×159×230) / 0                     3868 / 13409 / 0
  Dance          596            (100×126×111) / 0                      2465 / 8989 / 0

           Fig. 9.   The size of each data and processing time in the preparation steps.

           dataset      avg. # of    targetting   setting eq.   solving eq.   total time
                        iterations    (sec(%))     (sec(%))      (sec(%))        (sec)
            Camel           10        0.1(1.2)     8.1(77.3)     2.1(20.3)       10.5
            Horse           10        0.2(2.5)     6.6(58.2)     4.5(39.3)       11.4
            Dance           10        0.2(2.5)     4.5(70.4)     1.7(27.1)        6.4

Fig. 10. The number of iterations per frame and the average time (seconds per frame) spent for
each step in our fitting scheme.

Fig. 11. Results for the Horse data. The rest state mesh (upper-left) was built using the first
frame, and fitted to the subsequent time frames.

              camel puppet taken from 14 different viewpoints.
Horse:        A synthetic animation created by a designer using modeling software.
Dance:        A simulated visual hull data constructed from a 30fps motion capture
              dataset. We first produced a volumetric model by skinning the motion
              capture data and then used synthetic renderings of this surface from 6
              viewpoints as the input to the visual hull algorithm.
  For the Dance dataset, the rest state mesh was built from original synthetic
surface that was used to compute the visual hulls. For the other two datasets, the
rest state meshes were created from one of the frames. After choosing a frame for
the rest state mesh and the parameters for optimization and collision detection, our
procedure can run without intervention throughout the entire animation sequences.
The statistics on the size and the performance of each data are given in the Figure
9 and 10.
  The results of our fitting approach are shown in Figure 1, 11 and 12. In the
ACM Journal Name, Vol. V, No. N, Month 20YY.
                                      Fitting Solid Meshes using Linear Elasticity     ·     13

Fig. 12. The Dance data: the (simulated) data collection setup for the visual hull (upper-left),
the rest state mesh (upper-right), and the result of fitting (bottom row).

figures, input surfaces are illustrated as transparent orange surfaces, along with
the fitted tetrahedral meshes. To visualize the connectivity of the model, each
tetrahedra in the model was scaled down. Note how the tetrahedral model follows
the deformation of orange input surface, despite of large rotations of body parts,
noise, and false surface regions in visual hull data. Rotation compensation allows
us to use linear elasticity to address large rotations in the deformation efficiently,
and elasticity effectively make the optimization stable and robust against noise.
   Although using a physically-based solid model increases the computational cost of
the fitting significantly, we believe that our approach enables us to achieve the level
of deformation that cannot be achieved by surface model-based approaches. The
result in Figure 2 is a good example. In that example, the curvature of the surface
in the deformed area is completely different between the shapes before and after.
In fact, the curvature on either sides were inverted. If we use surface model only,
then it wouldn’t be possible to fit the initial model to the target shape without any
information on intermediate deformations. But with solid model, we can improve
the fitting iteratively, because the elasticity in the interior of the model holds all
the points together, and guide the fitting through the extreme deformation (as
illustrated in Figure 6).
   Although our procedure can track large scale motions very well, closer examina-
tion of the results reveals a few limitations and problems, as illustrated in Figure
13. Small scale deformations might be incorrectly interpreted as large scale ones:
notice that the the fetlock joint shown in Figure 13 seems to be twisted around
the leg’s axis rather than bent toward the back. We have also self-intersections in
the regions where body parts contact. Furthermore, in some cases, we noticed that
our algorithm fits a limb to a spurious geometric feature present in the visual hull
                                                 ACM Journal Name, Vol. V, No. N, Month 20YY.
14     ·     J. Choi and A. Szymczak

Fig. 13. Problems we have encountered: a twisted fetlock joint, self-intersection, a limb trapped
in the false region of visual hull, and discrepancy between surface and fitted mesh.

Fig. 14. The error for the Dance data: In this case, we have the original surface model that we
can measure the exact volume and the error of fitted model. The error measure is the percentage
of the volume of the fitted model that remained outside of the original model.

data. The last problem is the imperfect fit between input surface and the result
model. These problems originate from simple approximation of linear elasticity,
the assumption of uniform material properties, and the inevitably complex solution
space of the optimization due to complex shape and motion. Still, we can conclude
that our approach handles large scale deformation very well and is robust enough
to be applicable to realistic computer vision data. Figure 14 shows the relative
amount of error of the fitted model for each frame in the Dance dataset. The errors
are closely related to a specific choice of Young’s modulus in the stiffness matrix k.
(Small modulus means more flexibility and less robustness.)

We have presented a new fitting method to fit solid meshes to noisy and unstruc-
tured 3D animation data. We showed that complex motions that involve large
rotations of body parts can be tracked by a simple iterative optimization using
linear elasticity with rotation compensation. The results provide the correspon-
dence of any particular point on a moving 3D object throughout all the frames,
which is useful for texture mapping, compression, and further analysis of motion
and structure.
ACM Journal Name, Vol. V, No. N, Month 20YY.
                                       Fitting Solid Meshes using Linear Elasticity       ·     15

Fig. 15. Visualization of the maximum magnitude of the internal force over the entire animation
sequence at some points near the medial axis. (Colored from red to green to blue, in the decreasing
order of the magnitude of the internal deformation force.) This result demonstrates the good
possibility of recovering rigid parts and joints from a particular sequence of deformation.

   We believe our approach can be improved in many ways. In particular, we
would like to explore analyzing the structure of the deforming body based on our
algorithm’s output. Such analysis can be performed by examining the deformation
of each mesh element and could lead to detection of limbs, rigid parts and other
important features. Preliminary experiments show that such analysis can be done
by examining the internal forces computed by our algorithm (Figure 15). Apart
from being of interest of its own, taking advantage of such an analysis to improve
the rest state mesh (for example, by using different material properties for different
tetrahedra) is an intriguing possibility.
   Finally, it would also be interesting to replace the mesh-based technique described
above with a meshless approach. The advantage of meshless approaches is that
they do not require building a high quality rest state tetrahedral mesh. Initial
results based on element-free elasticity simulation are presented in [Choi et al.
2006]. Another option could be to use a graph-based technique such as [Zhou et al.

This work was supported by NSF grant DMS-CARGO-0138420. We also would
like to thank Byungmoon Kim for discussions on Finite Element Methods, Gabriel
Brostow for providing voxel sets of the Camel data, and Robert W. Summer for
sharing the Galloping Horse data.

Alexa, M., Cohen-Or, D., and Levin, D. 2000. As-rigid-as-possible shape interpolation. In
  SIGGRAPH Conference Proceedings. ACM Press, 157–164.
Brostow, G. J., Essa, I., Steedly, D., and Kwatra, V. 2004. Novel skeletal representation for
  articulated creatures. In ECCV04. Vol III: 66–78.
Choi, J., Szymczak, A., Turk, G., and Essa, I. 2006. Element-free elastic models for volume
  fitting and capture. In Proceedings of CVPR’2006. 2245–2252.
                                                  ACM Journal Name, Vol. V, No. N, Month 20YY.
16    ·     J. Choi and A. Szymczak

Cutler, B., Dorsey, J., and McMillan, L. 2004. Simplification and improvement of tetrahedral
  models for simulation. In Eurographics Symposium on Geometry Processing. 103–114.
Edelsbrunner, H. and Guoy, D. 2002. An experimental study of sliver exudation. Engineering
  with Computers, Special Issue on ‘Mesh Generation’ (10th IMR 2001) 18, 3, 229–240.
Lohner, R. 1988. Generation of three-dimensional unstructured grids by the advancing front
  algorithm. Int. Journal for Numerical Methods in Fluids 8, 1135–1149.
Matusik, W., Buehler, C., Raskar, R., Gortler, S. J., and McMillan, L. 2000. Image-based
  visual hulls. In Proceedings of ACM SIGGRAPH 2000. 369–374.
McInerney, T. and Terzopoulos, D. 1993. A finite element model for 3d shape reconstruction
  and nonrigid motion tracking. In ICCV Proceedings. Berlin, Germany.
Molino, N., Bridson, R., Teran, J., and Fedkiw, R. 2003. A crystalline, red green strategy for
  meshing highly deformable objects with tetrahedra. In 12th Int. Meshing Roundtable. 103–114.
Muller, M., Dorsey, J., McMillan, L., Jagnow, R., and Cutler, B. 2002. Stable real-
  time deformations. In Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on
  Computer animation. ACM Press, 49–54.
Nielson, G. M. and Sung, J. 1997. Interval volume tetrahedralization. In IEEE Visualization
  ’97. 221–228.
Petland, A. and Horowitz, B. 1991. Recovery of nonrigid motion and structure. IEEE Trans.
  Pattern Anal. Mach. Intell. 13, 7, 730–742.
Rassineux, A. 1998. Generation and optimization of tetrahedral meshes by advancing front
  technique. Int. Journal for Numerical Methods in Engineering 41, 651–674.
Sand, P., McMillan, L., and Popovic, J. 2003. Continuous capture of skin deformation. ACM
  Transactions on Graphics 22, 3, 578–586.
Sethian, J. A. 1996. A fast marching level set method for monotonically advancing fronts.
  Applied Mathematics 93, 4 (February), 1591–1595.
Shewchuk, J. R. 1998. Tetrahedral mesh generation by delaunay refinement. In Symposium on
  Computational Geometry. 86–95.
Shoberl, J. 1997. Netgen - an advancing front 2d/3d mesh generator based on abstact rules.
  Computing and Visualization in Science 1, 1, 41–52.
Starck, J., Hilton, A., and Illingworth, J. 2002. Reconstruction of animated models from
  images using constrained deformable surfaces. In Proc. DGCI 2002, A. Braquelaire, J.-O.
  Lachaud, and A. Vialard, Eds. 382–391.
Terzopoulos, D., Platt, J., Barr, A., and Fleischer, K. 1987. Elastically deformable models.
  In Computer Graphics. ACM Press, 205–214.
Terzopoulos, D. and Witkin, A. 1988. Physically based models with rigid and deformable
  components. IEEE Comput. Graph. Appl. 8, 6, 41–51.
Terzopoulos, D., Witkin, A., and Kass, M. 1988. Constraints on deformable models: Recov-
  ering 3d shape and nonrigid motion. Artificial Intelligence 35.
Yerry, M. A. and Shephard, M. S. 1984. Automatic three-dimensional mesh generation by the
  modified octree technique. Int. Journal for Numerical Methods in Engineering 20, 1965–1990.
Zhou, K., Huang, J., Snyder, J., Liu, X., Bao, H., Guo, B., and Shum, H.-Y. 2005. Large
  mesh deformation using the volumetric graph laplacian. ACM Trans. Graph. 24, 3, 496–503.
Zhuang, Y. and Canny, J. 1999. Real-time simulation of physically realistic global deformation.
  In SIGGRAPH Conference Proceedings. ACM Press, 270–273.

Received Month Year; Revised Month Year; accepted Month Year

ACM Journal Name, Vol. V, No. N, Month 20YY.

Shared By: