University of Southern California Dissertations and Theses

Correspondence-based analysis of 3D deformable shapes: methods and applications

(USC Thesis Other)

Correspondence-based analysis of 3D deformable shapes: methods and applications

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content CORRESPONDENCE-BASED ANALYSIS OF 3D
DEFORMABLE SHAPES: METHODS AND APPLICATIONS
by
Ilya Eckstein
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(COMPUTER SCIENCE)
December 2007
Copyright 2007 Ilya Eckstein
Acknowledgements
A few people deserve to be mentioned here—people without whom this disserta-
tion would have been unlikely to see the light of day. I would like to thank both
my advisors, Prof. Mathieu Desbrun and Prof. C.-C. Jay Kuo, who generously
supported me, each in his own way, all along this journey. I am also grateful to
my collaborators: Jean-Philippe Pons, for infecting me with his ease of manipu-
lating math, and Yiying Tong, for always putting time and effort into sharing his
knowledge with me. Going further back in time to my Technion days, I should
also mention Prof. Craig Gotsman, my master’s thesis advisor, and Prof. Gershon
Elber, who was the first to woo me with his mind-boggling lectures on geometric
modeling. Last but not the least, I am deeply grateful to my parents for bearing
with me and my embarkment on this adventure. Thank you all!
ii
Table of Contents
Acknowledgements ii
List Of Tables v
List Of Figures vi
Abstract xi
Preface xiii
Chapter 1 : Introduction 1
1.1 SignificanceoftheResearch ...... ........... ..... 1
1.2 ResearchMotivationandTechnicalChallenges ....... ..... 3
1.3 ContributionsoftheResearch ..... ........... ..... 7
1.4 OutlineoftheThesis .......... ........... ..... 9
Chapter 2 : Compression of Time-Varying Isosurfaces 10
2.1 Introduction .... ........... ........... ..... 11
2.1.1 Background ........... ........... ..... 13
2.1.2 Contributions .......... ........... ..... 14
2.2 TheStaticIsosurfaceCompressionMethod ......... ..... 15
2.3 Videovs. Dynamic3DSurfaceCoding ........... ..... 17
2.3.1 VideoCoding .......... ........... ..... 18
2.3.2 FromVideotoDynamic3DIsosurfaces ....... ..... 18
2.3.3 GeneralApproach........ ........... ..... 21
2.4 CompressionthroughBlock-basedPrediction........ ..... 22
2.4.1 IsosurfaceMacroblocks ..... ........... ..... 22
2.4.2 MotionEstimation ....... ........... ..... 23
2.4.3 MotionCompensation ..... ........... ..... 27
2.4.4 Context-basedArithmeticCoding .......... ..... 29
2.5 Results....... ........... ........... ..... 30
2.6 DiscussionandFutureWork ...... ........... ..... 32
iii
Chapter 3 : Generalized Geometric Flows and Applications to Shape
Matching 40
3.1 Introduction .... ........... ........... ..... 41
3.2 MathematicalBackgroundonGradientFlows ....... ..... 43
3.2.1 ContinuousDefinitionofGradientFlows ...... ..... 44
3.2.2 DiscreteGradientFlows .... ........... ..... 47
3.3 GeneralizedGradientFlows ...... ........... ..... 49
3.3.1 DesigningNewInnerProducts. ........... ..... 49
3.3.2 GeneralizedGradients ..... ........... ..... 50
3.3.3 TimeIntegration ........ ........... ..... 51
3.4 Application: RobustShapeMatching . ........... ..... 53
3.4.1 Background ........... ........... ..... 53
3.4.2 HausdorffGradientFlowforMatching ....... ..... 55
3.4.3 ConstructingPrior-basedInnerProducts ...... ..... 56
3.4.3.1 Quasi-RigidDeformationPrior ...... ..... 56
3.4.3.2 Quasi-ArticulatedDeformationPrior ... ..... 61
3.4.3.3 Quasi-IsometricPrior ........... ..... 65
3.4.4 ResultsofPrior-basedShapeMatching ....... ..... 68
3.5 ImplementationDetails......... ........... ..... 71
3.5.1 Pseudo-HausdorffDistanceanditsGradient .... ..... 71
3.5.2 MultiresolutionFlows...... ........... ..... 75
3.5.2.1 GeneralApproach . . ........... ..... 75
3.5.2.2 MultiresolutionShapeRepresentation .. ..... 76
3.6 GeneralizedDifferentiableDistancePotential........ ..... 81
3.6.1 BeyondVertexPositions .... ........... ..... 81
3.6.2 AlternativestotheHausdorffDistance ....... ..... 83
3.7 ApplicationsBeyondShapeMatching . ........... ..... 88
3.7.1 RegularizedGradientFlows . . ........... ..... 88
3.7.2 Volume-ControlledImplicitFairing ......... ..... 91
3.8 ConclusionandFutureWork...... ........... ..... 94
Chapter 4 : Matching Cortical Surfaces with Surface Flows 97
4.1 Introduction .... ........... ........... ..... 98
4.2 Backround ..... ........... ........... ..... 101
4.3 Method ...... ........... ........... ..... 102
4.3.1 ChoiceofaDeformationPrior . ........... ..... 103
4.3.2 MatchingCorticalSurfaces... ........... ..... 106
4.4 Results....... ........... ........... ..... 109
4.5 DiscussionandFutureWork ...... ........... ..... 111
Chapter 5 : Conclusions and Future Directions 114
Bibliography 117
iv
List Of Tables
2.1 Comparing compression performance of motion prediction vs. frame-
by-framemode. Ratesareshowninbytes. ......... ..... 32
4.1 Quality of match between deformed template and subject brains as
average L
2
distancesbetweencorrespondingsulcalcurves. ..... 110
v
List Of Figures
2.1 Static Isosurface Compression: Each leaf of the inside/outside oc-
tree (left) intersects the isosurface, and the position of the vertex
(obtained through Dual Contouring) is encoded using a sign-based
localframeasdetailedin[74]. .... ........... ..... 15
2.2 Illustrating differential prediction vs. block-based prediction in MPEG:
instead of recovering a video frame (top left) from the previous one
(top middle) with a difference image (top right), a more efficient ap-
proach is to partition the new frame in blocks (bottom left), and track
these blocks onto the previous frame to their best matching position
(exaggerated in that figure for illustrative purposes): the encoding
of the few block motion vectors and the resulting difference image
(with much less range) requires less bits due to the strong temporal
coherencepresentinusualvideosequences. ........ ..... 19
2.3 Left, for a surface slightly shifted (here, a disk moved to the right),
the difference contour between two time steps is more verbose than
each of the steps. Right, shown is a local (i.e., block-sized) surface
motion between two subsequent frames, and the result of a local
registrationprocedure. ......... ........... ..... 26
2.4 Predicting vertex position. The grey vertex is extracted from the
predicted distance function. The likelihood ellipsoid is aligned with
the local frame, describing the expected error distribution. Using the
motion vector, the position is mapped back to the previous frame and
the prediction is refined using the nearest original vertex within the
likelihood ellipsoid. .......... ........... ..... 37
2.5 Fluid simulation with falling droplets generating ripples and splashes. 37
vi
2.6 Another simulation of splashing fluid, with high topological complex-
ity. ........ ........... ........... ..... 38
2.7 MRI headscan sequence (the visible human dataset sampled at differ-
ent iso-values). Note the amount of noise and topological complexity
inthemiddleframe. .......... ........... ..... 38
2.8 Rate/distortion curves for the two fluid sequences and the MRI head-
scan sequence. We compare compression with motion prediction to
frame-by-frame encoding for different quantizations. The x-axis rep-
resents the number of bytes, while the y-axis represents the L
2
dis-
tance between the decompressed version and the original sequence
(obtained by summing eachL
2
distance for each frame) as measured
byMESH[17]. . . ........... ........... ..... 39
3.1 Articulated matching of an Armadillo template mesh to a very dif-
ferent pose, using the Hausdorff gradient flow with a quasi-articulated
prior. Our approach can embed such shape priors into conventional
surface flows to provide, e.g., robust and efficient shape matching
with no markers. Notice that although a template arm initially co-
incides with a target mesh leg (resulting in a strong local minimum),
thefinalposeissuccessfullyrecovered. ........... ..... 42
3.2 Quasi-Articulated flow: matching two hands by minimizing their
Hausdorff distance with a gradient flow incorporating an articulated
prior (using 16 regions: 3 per finger + palm). The flow automatically
morphs the red hand into the blue one, first mostly through a global
rotation, then by adjusting each finger. Color mix on the final frame
is due to the fact that the hands do not have the same sampling and
connectivity. ... ........... ........... ..... 46
3.3 Incorrect matching between two fingers obtained through a straight-
forward L
2
gradient flow of their Hausdorff distance. The absence
of rigidity priors leads to significant spurious sliding and mismatch.
Theintermediatepositionsalsoexhibitlargedeformation. ..... 53
vii
3.4 GMDS correspondences (courtesy of [23]) between two differ-
ent Armadillo poses, for a subset of 100 vertices, with the corre-
sponding Voronoi regions shown in matching colors. While most
regions are mapped adequately (e.g., the tail, shown with a black
arrow), hands and legs still have many inaccurate correspondences
(e.g., those shown with red arrows), which will result in unnatural
deformationsiffurtherusedforregistration. ........ ..... 57
3.5 Quasi-Rigid Prior using Multiresolution: Matching two full
body meshes with different sampling (∼ 28K vertices each) by min-
imizing their Hausdorff distance with a quasi-rigid gradient flow.
To speed up computations, a multi-resolution scheme was employed
(total time: 1 min). The coarsest resolution is shown in the middle.
Closeups show the final alignment (template mesh in wireframe). . . 58
3.6 Quasi-Articulated Prior: Matching two arms with the Hausdorff
gradient flow incorporating a quasi-articulated prior (using 18 sub-
parts: 16 for the hand, plus the upper and lower arm). Again, the
mesh sampling of both arms is very different. Despite the pres-
ence of fingers in different positions (adding more local minima to
the distance-based energy), our algorithm finds the right matching
through a combination of global and local rotations, with a small
amount of non-rigid deformation (for instance at the elbow). Final
alignmentofthefingersisshownincloseupontheright. . ..... 62
3.7 Quasi-Isometric Prior: Matching two arms by minimizing their
Hausdorff distance with a quasi-isometric gradient flow. Here, with-
out any prior knowledge on the leg structure, the quasi-isometric
gradient flow finds the appropriate matching. The last frame shows
both meshes, to demonstrate the quality of the match despite the
samplingdisparitybetweenthetwomodels. ........ ..... 66
3.8 Non-rigid Shape Matching: Here, an articulated template model
of a lamb (red) is automatically deformed to fit a different pose
(blue), by directly minimizing their Hausdorff distance with a gradi-
ent flow based on a quasi-articulated prior. To handle this example
(∼ 22000 vertices for each mesh) in less than two minutes, a multi-
resolution scheme was used. Notice that our prior automatically
makes the lamb move mostly rigidly first to best fit the target pose,
then adjust the limbs non-rigidly. The color mix on the final frame
shows that, even if both lambs do not have the same sampling and
connectivity,theirgeometrymatcheswell. ......... ..... 69
viii
3.9 Partial Matching of an incomplete scanned 3D mesh (in blue, cour-
tesy of U. of Washington) with a template arm in a different pose.
Despite the sparsity of the scan, a quasi-articulated prior (19 sub-
parts, including fingers) on our Hausdorff-based flow gets a good
match. ....... ........... ........... ..... 71
3.10 Illustrating asymmetric proximity: x
j
is among the nearest neighbors
of y
j
,butnotviceversa. ........ ........... ..... 73
3.11 Human body mesh, at two levels of detail: the size of a full resolution
surface(left)isreducedto10%oftheoriginal. ....... ..... 77
3.12 Local pyramid coordinates. A vertex v along with its one-ring neigh-
bors {v
i
} are projected onto the local tangent plane, where mean
value coordinates {w
i
} of v are computed. The vertical offset h en-
codedseparately. . ........... ........... ..... 78
3.13 Self-intersecting polygons (left) do not guarantee valid mean value
coordinates. Concavities in simple polygons (right) are more likely
to create self-intersections under a deformation. Thus, the concave
greyvertexisdiscardedduringconstruction. ........ ..... 80
3.14 Gaussian curvature enhanced distance functional is used to
register two extremely different rigid poses of a full body shape. A
strong local minimum is avoided with the Gaussian curvature incor-
poratedintotheHausdorffdistancefunctional........ ..... 84
3.15 Implicit fairing with local volume control. Top row: a noisy
horse model (left) becomes degenerate after aggressive implicit fair-
ing (middle), while a comparable amount of smoothing with local
volume control avoid shrinking (right); closeup shows overlaid com-
parison. Bottom: several iterations of implicit fairing of a noisy fe-
line model (far left), with local volume control (above) and without
it (below), shown for results with a comparable degree of smooth-
ness. With the original scheme, a strong smoothing (far right) causes
an extreme shrinkage and degeneracies, whereas with volume control
(λ=10
−6
), the original proportions and overall shape are retained.
The inset illustrates the non-shrinking effect for the tail region, with
thecorrespondingsurfacesoverlaid. . . ........... ..... 96
ix
4.1 White matter cortical surfaces of two different subjects. The com-
plexity and anatomical variability of folding patterns make identifi-
cation of correspondences challenging even for the human eye. . . . 99
4.2 A white matter cortex and its corresponding Partially Flattened Sur-
face (a), where the major sulci are still clearly present. A basic
Hausdorff gradient flow applied on a rotated, smoother version (b)
creates spurious deformations (c) due to lack of coherence and local
minima, while with a quasi-rigid deformation prior (d), it success-
fullyrecoversthetransformation. ... ........... ..... 103
4.3 Automatically matching a template (grey) to the subject cortex
(blue). Partially flattened representations of both surfaces are iter-
atively aligned using a Hausdorff flow with a smoothing prior. The
obtained alignment yields a correspondence between the original sur-
faces. The final color mix is due the fact that the surfaces lie on each
other. ....... ........... ........... ..... 106
4.4 a. Automatic matching of PFSs yields a close alignment for most
sulcal curves. b. Constraining only 7 out of the 23 available curves
reduces most misalignments, further improved by using the full set of
constraints (c). d. Corresponding sulcal alignment for the original
cortical surfaces. For clarity, a single cerebral hemisphere is shown. 109
4.5 A cortical surface (far left) is gradually smoothed with near volume
preservation (VP)—above, and without it—below. Far right: color-
coded Hausdorff distances (blue=smallest, red=largest) between the
original and the final results, as measured by the EPFL Mesh tool [17].
The result of VP fairing is much closer to the original for a similar
smoothingeffect. . ........... ........... ..... 113
x
Abstract
People have been studying shapes since the ancient times, using geometry to
model those that already existed in the world and to create new ones. Today,
with the power of modern computers, the task of capturing the world digitally and
recreating it virtually is no longer science fiction. This challenging, but extremely
important problem spans applications ranging from medical modeling and scientific
visualization to computer-aided design, art and entertainment. However, to achieve
that goal, we need better algorithms that can ”understand” geometry.
This research focuses on applications involving geometry analysis that relies on
the notion of shape correspondence. First, the problem of compressing time-varying
isosurfaces is studied, and a region-based motion estimation method is employed to
predict surface geometry from one frame to the next. Fundamentally, the suggested
approach reduces the motion estimation problem to local rigid body registration.
Next, to extend the context to the wider class of deformable bodies, a more general
shape matching problem is studied. The proposed solution presents a general yet
flexible deformation-based approach, with articulated shape matching used as the
test bed application. Finally, we apply the aforementioned deformation framework
xi
to the challenging problem of cortical surface matching, one of the key problems
in the field of brain mapping. The value of this approach over existing methods is
demonstrated both for automatic and user-guided registration scenarios.
xii
Preface
As we were initially looking for a new problem to sink our teeth into, it became
obvious that the lion’s share of the shape processing research was still focusing on
the analysis of static objects – leaving most dynamic scenarios as future work. That
motivated us to start our work with animated 3D data to see where it would take
us and what the real research issues were.
Since we were particularly interested in geometric data compression, the prob-
lem of compressing sequences of unstructured, time-varying geometry (a 3D equiv-
alent of the video content, if you will) seemed both important and exciting. We
decided to start with dynamic isosurfaces as a specific application, and set out to
work. Soon enough, it became clear that, data encoding aspects aside, the most
difficult problem was essentially about establishing correspondences between the
highly deformable surfaces we had to deal with. It also became clear that the
problem of matching deformable shapes is of significance well beyond the context
of compression, or even processing of time-varying geometry in general. In fact, it
arises any time a geometric association between two shapes is sought. The aim of
this PhD thesis is to explore how to go about it in a reliable and principled manner.
xiii
Chapter 1
Introduction
- What manner o’ thing is your crocodile?
- It is shaped, sir, like itself; and it is as broad as it hath
breadth: it is just so high as it is, and moves with its own
organs: it lives by that which nourisheth it; and the
elements once out of it, it transmigrates.
- What colour is it of?
-Ofitowncolourtoo.
-’Tisastrangeserpent.
W. Shakespeare, Antony and Cleopatra
1.1 Significance of the Research
Generally speaking, geometry processing algorithms deal with two fundamental
categories of problems: synthesis of new shapes and analysis of given ones. The
first category focuses mainly on applications such as shape modeling and physical
simulation of materials. On the other hand, shape reconstruction and recognition,
parametrization and discretization, segmentation and matching, compression and
approximation - all involve geometric data analysis. Research presented in this
thesis is concerned with the latter category.
1
Most analytical tasks inherently involve some notion of context. In this chap-
ter’s epigraph, for instance, the description of a crocodile is hardly informative,
since it fails to convey the animal’s shape, motion and color in relation to other
species’ characteristics. Those who fancy precise formulations over Shakespearean
metaphors may simply think of this example as a model fitting problem: given an
animal (crocodile) that is known a priori to belong to some family S of species
(reptiles), express its shape in terms of the anatomical model parameters of S.By
definition of the problem, the given animal’s geometric body features should be an-
alyzed in the context of S. Specifically, one might use a generic parametrized body
shape (a reptilian anatomical template), so that the problem reduces to fitting the
template to the given specimen.
As another example, we may consider the problem of inverse kinematics: given
a series of independent observations of a moving object over time, we seek to derive
the corresponding motion function. To find a solution, we must treat each instance
of the shape in the context of other samples. To understand how the object is
moving, we need to know how each instance corresponds to the next one along the
time axis.
In both examples given above, the analytical tasks involve solving a shape cor-
respondence (or matching) problem. The question is whether there is a general
solution to the problem that can be applied in different settings in a standard man-
ner. This appears to be a difficult and open problem since shape matching may
2
vary greatly depending on the geometric representation of the object, its physical
properties, specifics of the application, etc.
In this research, we focus on certain application contexts yet strive towards a
general and flexible framework for correspondence-based shape analysis. In the
beginning, we were concerned with the encoding (or an efficient representation) of
dynamic isosurfaces as a specific application. To remove the temporal redundancy
of isosurfaces, it turns out that we have to address a more fundamental problem:
namely, to establish correspondences between highly deformable surfaces in consec-
utive time instances. Yet the need for a technique of matching deformable shapes
is far more general, arising way in a broad spectrum of applications. As a result,
shape matching becomes the focal point of this work.
1.2 Research Motivation and Technical Challenges
The notion of shape correspondence has been utilized in a wide range of geometry
processing applications. Image-based matching (both for 2D and 3D images) is a
closely related subject in computer vision. In the following discussion, we will stay
in the context of three-dimensional objects and occasionally refer to image-based
methods. Rather than attempt to give an exhaustive survey of shape correspon-
dence techniques, this review serves as a background to the contributions of our
research (stated in the next section).
3
One classic correspondence problem in geometry processing and computer vision
is that of rigid body registration, where partial scans of an object in different coor-
dinate frames are required to be aligned into a single shape representation. Even
though only rigid transformations are considered in this case, the problem becomes
non-trivial due to incomplete and noisy input data as well as sampling differences
across scans. When a coarse alignment is given, the standard approach for high-
accuracy registration is the Iterative Closest Point algorithm [16,33], spawning a
legacy of variants and improvements. The global registration problem (i.e.,when
the whole space of rigid transformations has to be considered) is even harder, and
heuristic landmarks are typically used to drive the matching process. We refer the
reader to Planitz et al. [96] for a detailed survey on rigid body registration.
However, objects are not rigid in general, neither are they necessarily static. Ac-
quiring an accurate model of a deformable object often requires a mapping between
several different configuration (or poses). Such a mapping is important not only for
reconstruction purposes, but also for shape classification and segmentation, struc-
ture inference, inverse kinematics, texture transfer, animation compression and
many other related analytical tasks. Yet finding a meaningful mapping between
deformable surfaces without any a priori information can be extremely difficult,
or even infeasible. Hence, in practice, most solutions are model-specific and only
semi-automatic. For instance, articulated body reconstruction from multiple poses
is usually done by using a template (or skeleton) with markers added during the
4
acquisition process [6,9]. The use of markers is also typical in facial motion capture
systems [94]. In absence of markers, one may resort to user-specified correspon-
dences. However, those often amount to tedious manual labor, requiring repeated
user interaction for each new input. Statistics and machine learning techniques [8]
may be an alternative that is powerful, yet computationally expensive. For anima-
tion sequences with limited deformation from frame to frame, one may also apply
other approaches adapted from video tracking, such as optical flow [10,89].
Whether an initial set of sparse correspondences is given or not, a global map-
ping is ultimately sought. To this end, two major categories of methods can be
identified:
• those based on direct transformation (deformation) of one object or embed-
ding space onto another [34,81,99]; and
• those using the concept of cross-parametrization [72,73,100].
In its simplest form, the latter approach uses a shared parameter domain (e.g., a
sphere or a disk) for both surfaces. Thus, no plausible or minimal deformation
between the two surfaces is sought, as correspondences are defined by juxtapos-
ing their 2D maps instead, reducing the problem to the two-dimensional domain.
While this strategy is attractive for some applications, such as surface morphing
scenarios [3], it also has several drawbacks. First, both shapes are required to have
a parametrized representation, which may not even be realistic for topologically
complex objects. Second, in many cases, a meaningful deformation is the key to
5
obtaining a correct shape alignment— hence the need for deformation-based meth-
ods. The shortcoming of most techniques in that category is their limitation to a
specific deformation model (splines offer one popular choice [34,99]), which may not
fit the nature of the input data well. Furthermore, that model is often hard-coded
into the objective function, making the result less then optimal, as discussed in
Chapter 3.
Notably, for the class of near-isometric deformations (articulated motion is
one example), the method of canonical, ”bending invariant” forms has been pro-
posed [44]. This technique can be seen as a variation of the cross-parametrization
paradigm, as an intermediate pose-invariant representation is used to obtain a map-
ping between two surfaces. Yet in practice, each surface yields a somewhat different
canonical form, so the two corresponding forms still need to be aligned. To address
that shortcoming, a recent generalization of the same method presents a framework
for direct near-isometric embeddings of one surface onto another [22]. Both tech-
niques still involve high computational costs, although a multigrid approach has
been suggested for the former case [25].
Shape correspondence methods play an important role in medical applications,
say, for reconstruction and analysis of human organ models acquired by various
imaging methods such as MRI, CT, etc. For instance, matching of human cortical
surfaces is one of the central topics in the field of brain mapping [113]. This problem
will be more closely examined in Chapter 4.
6
1.3 Contributions of the Research
The contributions of this work include the following:
• Compression of Time-Varying Isosurfaces
Sequences of complex time-varying isosurfaces, as generated by medical imag-
ing procedures or physical simulations, often require vast storage and a long
transmission time. Yet the task of compressing them can be extremely chal-
lenging due to strong deformations and topology changes in the input surfaces,
making most of the previous techniques for compression of dynamic 3D con-
tent inefficient or impractical. We propose a new approach which exploits
temporal coherence of the data by adopting the paradigm of block-based mo-
tion prediction from video coding and extending it using local surface registra-
tion in Chapter 2. The resulting prediction errors across frames are treated as
a static isosurface and encoded progressively using an adaptive octree-based
scheme. To the best of our knowledge, this is the only compression method
tailored to time-varying isosurfaces.
• Generalized Geometric Flows and Applications to Shape Matching
We propose a generalized approach to surface evolution in geometry process-
ing in Chapter 3. The concept of surface flow based on functional minimiza-
tion have been used in numerous contexts. However, such gradient flows are
7
nearly always based on the canonical L
2
inner product of vector fields. Chang-
ing this inner product provides a powerful, simple, and untapped approach to
extend current flows. Specifically, we demonstrate the power of such a norm
alteration in the context of shape matching, where deformation priors, from
rigid motion to more complex cases such as articulated motion and beyond,
can be built in a gradient flow to improve results significantly. We also present
differentiable approximations of several distance metrics between meshes and
as well as the use of implicit integration to prevent spurious oscillations. Fur-
thermore, our solution is made scalable by building it into a multiresolution
framework. The result is a general yet flexible shape matching toolbox that
can be tailored to a variety of applications.
• Deformable Brain Atlases and Cortical Surface Matching
Finally, in Chapter 4 we extend the framework of generalized geometric flows
to the problem of cortical surface matching, which is extremely important
e.g. in the field of brain mapping. This task can be remarkably complex due
to extreme anatomic variability of cortical surfaces and their highly complex
geometry. Coupled with a multiscale strategy to better cope with these issues,
our flow-based approach achieves results comparable in quality to the state-
of-the-art in terms of registration accuracy, with a speedup of up to two
orders of magnitude. In addition, we show how to incorporate optional curve
matching constraints into the basic method, and demonstrate the value of the
8
outcome over existing methods in both automatic and user-guided registration
scenarios.
1.4 Outline of the Thesis
The rest of this text is organized as follows. Compression of time-varying isosurfaces
is examined in Chapter 2. The formalism generalized gradient flows and its use in
shape matching scenarios is discussed in Chapter 3, and is leveraged further in
Chapter 4 to address the problem of cortical surface matching. Finally, concluding
remarks are given in Chapter 5, along with potential extensions and future work.
9
Chapter 2
Compression of Time-Varying Isosurfaces
- How would you fit an elephant into a fridge?
- In three simple steps: open the fridge, put the elephant
there, close the fridge.
From Russian folklore
Compressing sequences of complex time-varying surfaces as generated by medi-
cal instrumentations or complex physical simulations can be extremely challenging:
repeated topology changes during the surface evolution render most of the previ-
ous compression techniques for time sequences of surfaces inefficient or impractical.
In order to provide a viable solution, we propose a new approach based upon an
existing isosurface compression technique designed for static surfaces. We exploit
temporal coherence of the data by adopting the paradigm of block-based motion
prediction developed in video coding and extending it using local surface registra-
tion. The resulting prediction errors across frames are treated as a static isosurface
and encoded progressively using an adaptive octree-based scheme. We also exploit
local spatiotemporal patterns through context-based arithmetic coding. Fine-grain
10
geometric residuals are encoded separately with user-specified precision. The other
design choices made to handle large datasets are detailed.
2.1 Introduction
The realm of dynamic surface compression, i.e., encoding and decoding of surface
motion in space and time, has so far been enjoying much less attention than its
static counterpart. In contrast, the field of signal processing has witnessed a natural
evolution from digital images to digital video over the last decade, and as a con-
sequence, digital encoding of video sequences is now common practice. Now that
dynamic 3D contents are generated at a fast growing rate in various applications
ranging from characters in movies and electronic games to scientific simulation, the
need for dynamic surface compression is becoming greater.
Not all dynamic 3D data need specific compression techniques. Computer-
generated movies, while extensively using 3D animated models and gaining an
unprecedented popularity, are treated as video once produced: the specific camera
position and lighting conditions of each scene are parts of the artistic value of
the movie, so this content is rarely encoded in 3D. Another source of 3D data is
distributed electronic games, where characters are typically animated using motion
parameters (either precomputed or generated at runtime) applied to a reference
object: here again, there is no specific need for dynamic surface compression.
11
There are, however, a number of other sources for which such compression is vi-
tal. Scientific simulation data typically requires significant computational resources
to be generated (for fluid mechanics simulation for instance), and a vast amount
of space to store the results. Similarly, acquired dynamic data such as 3D medical
images and sequences of range scans demand a great deal of storage space as well.
Compressing this type of animated 3D content is, in these cases, key to their proper
handling.
Although some work has been done on compression of fixed-topology mesh an-
imation and time-varying volume data, we will focus on a specific type of dynamic
3D data, namely time-varying isosurfaces. Isosurfaces are a useful representation
for visualizing volumetric content, be it in scientific simulation or medical applica-
tions; this is particularly true for time-varying volumetric data, where one is often
interested in observing the evolution of one or more specific iso-contours (corre-
sponding to particular values of intensity, density, etc.), rather than the volume in
general. Unfortunately, the size of the data often becomes the bottleneck, even once
the volume dataset has been converted to a polygonal mesh. In particular, for appli-
cations involving remote data monitoring and analysis, high-resolution isosurfaces
become impractical without efficient compression. The challenge in compressing
isosurfaces is to cope with the level of geometric and topological complexity they
often have, as deformation may be arbitrary and topology can change over time.
12
Consequently, none of the existing dynamic mesh coding methods are applicable
to this case as we explain next.
2.1.1 Background
Existing methods for compression of dynamic polygonal meshes [4,21,55,62,66,75]
rely on known inter-frame vertex correspondences—typically, the mesh connectiv-
ity is fixed in time. This specificity renders those methods unsuitable in our con-
text, where not only the surface sampling but also the surface topology may vary
considerably between frames. Anuar and Guskov [10] proposed to solve the cor-
respondence problem by converting the input meshes to adaptive signed distance
volumes, followed by motion estimation using optical flow. Consequently, motion
is estimated as a dense vector field, impractical for compression purposes.
Other related work has been done in volumetric data compression. Mascarenhas et
al. [80] suggest a streaming scheme for volumetric grids, where voxels are ordered to
facilitate extraction of isosurfaces for subsequent isovalues, a very valuable feature
for data analysis. However, their approach does not support encoding of specific
(subset of) isosurface(s) and their motion in time, rendering it inefficient for our
purpose. Guthe and Straßer [56], and Sohn et al. [104] compress time-varying
volumetric data using wavelet transforms and motion estimation, but do not report
compression rates of single isosurfaces. Another recent work by Gregorsky et al.
[51] is concerned with the extraction of time-varying isosurfaces from compressed
13
volume data. There again, compression is applied to the whole raw data, rather
than just a single isosurface. Finally, Ibarria et al. [61] suggest a method for out-
of-core compression of n-dimensional scalar fields, mostly applied to 3D and 4D
datasets. Their main contribution is the so-called Lorenzo prediction, which is a
generalization of the parallelogram predictor to an arbitrary dimension. While this
predictor could in theory be used to predict corresponding inside/outside grid val-
ues for two subsequent frames, they provide no evaluation for values with a binary
range.
2.1.2 Contributions
In this work we extend the concept of block-based motion estimation, traditionally
used in video coding, for compression of dynamic isosurfaces extracted from a series
of uniformly sampled volumetric grids. Leveraging the progressive encoding method
of Lee et al. [74] for static isosurface compression, we provide a codec allowing rapid,
progressive display of the sequence during decoding. The unstructured nature of
time-varying isosurfaces allows our algorithm to deal naturally with surfaces of high
geometric complexity and dynamic topology, unlike any of the existing codecs.
14
Figure 2.1: Static Isosurface Compression: Each leaf of the inside/outside octree
(left) intersects the isosurface, and the position of the vertex (obtained through
Dual Contouring) is encoded using a sign-based local frame as detailed in [74].
2.2 The Static Isosurface Compression Method
Before delving into the details of the dynamic case, we first review the static iso-
surface compression scheme that we will build upon, and establish terminology. A
more detailed coverage can be found in the original text [74].
Input We are given a 3D grid of uniformly sampled volume data of size (2
n
+1)
3
,
along with an isovalue which determines the desired isosurface to encode.
Octree Generation and Surface Extraction Based on the input grid values,
a sign bitmap of the given volume (i.e., a binary grid of inside/outside values) is
generated. An octree representation of this sign bitmap is then constructed bottom-
up: whenever all eight child cells are homogeneous (i.e., contains no sign change),
they are deleted. This results in a tree has two kinds of leaf nodes: (a) coarse level
leaves in the homogeneous regions; and (b) heterogeneous leaves, corresponding to
the original grid cells intersected by the isosurface (the so-called narrow band). For
surface extraction, the Dual Contouring method [65] is applied, assigning a vertex to
each leaf cell, contrasting with the Marching Cubes algorithm [78] which positions
15
vertices at edge intersections instead. Dual Contouring produces watertight meshes
for any sparse or full octree description of the original bitmap and can even handle
sharp features if necessary. In the remainder of this chapter, we will call connectivity
the adjacency structure of the octree, and geometry the positions of the vertices
stored at the heterogeneous leaves of this octree as they roughly correspond to the
connectivity and geometry of a usual polygonal mesh (see Figure 2.1).
Connectivity Encoding At this stage, the octree representation of the sign
bitmap is encoded in a breadth-first fashion. Eight children of each cell correspond
to 27 grid corners, each of them encoded using a context based arithmetic coder.
For each homogeneous child, a leaf bit is encoded in a separate stream indicating
whether the corresponding cell is indeed a leaf node—to let the decoder know that
no subsequent data is needed for that particular subtree of the sign bitmap.
Geometry Encoding The sign bitmap alone is often sufficient to approximate
the original surface, by simply applying Dual Contouring after assigning +/-1 val-
ues to outside/inside grid corners. A useful feature of this approximation is that
the distortion (i.e., error introduced by the encoding process) is bounded by the
grid cell size. For some applications, this level of precision may prove sufficient
(surface smoothing techniques are often applied to improve the visual quality and
the quantization error in many cases). Otherwise, local geometric residuals are en-
coded for each leaf cell to provide the final vertex positions. Each residual vector
is encoded in a local coordinate frame, as the normal component typically bears
16
more information than the tangential components, and therefore needs a better
quantization accuracy. At each leaf cell, the local frame is given by the binary
sign configuration of that cell (see Figure 2.1, right), leading to 256 possible cases,
precomputed offline.
Discussion This algorithm provides, to the authors’ knowledge, the most ef-
ficient compression of static isosurfaces. As such, it can be used on time-varying
sequences on each successive isosurface. Obviously, this straightforward approach
is rather limited as no temporal coherence in the sequence can be either detected
or taken advantage of to reduce the bitrate. A modification of this approach to
render the compression method more efficient is presented next.
2.3 Video vs. Dynamic 3D Surface Coding
Recall that the input to our algorithm is a series of 3D volumetric frames{F
i
},all
sampled uniformly at the same resolution 2
n
+ 1, paired with isovalues {α
i
} de-
termining the desired isosurfaces. The na¨ ıve choice of frame-by-frame compression
mentioned above does not take advantage of the temporal coherence. We thus pro-
pose to rely on the paradigm of block-based motion prediction used in video coding
to exploit both geometrical and temporal coherence.
17
2.3.1 Video Coding
In video coding (e.g., in the MPEG format), a common approach to encode a single
frame at time t in a given video sequence is as follows (Figure 2.2):
• predict frame t basedonframe t-1 through motion prediction of small pixel
clusters called macroblocks;
• encode the estimated motion of each block;
• encode the prediction error image (i.e., the error between the prediction and
the actual frame).
The prediction error image is encoded using JPEG (i.e., block-based DCT) or
JPEG2000 (i.e., wavelets), two of the most efficient algorithms in image compres-
sion. Note that if two successive frames are too different as it happens in movies
when a scene ends and another one starts, the motion prediction and correction
can cost more than a direct encoding of the new frame; to deal with this case, each
frame is labeled as an I-frame or a P-frame: I-frames are Individually encoded,
while P-frames are Predicted from previous frames.
2.3.2 From Video to Dynamic 3D Isosurfaces
We now wish to mimic the approach used in video compression by first splitting
the surface in small block-like parts, estimating local motion of the isosurface parts
from the previous isosurface, then encoding the difference between the resulting
18
Figure 2.2: Illustrating differential prediction vs. block-basedpredictioninMPEG:
instead of recovering a video frame (top left) from the previous one (top middle)
with a difference image (top right), a more efficient approach is to partition the new
frame in blocks (bottom left), and track these blocks onto the previous frame to
their best matching position (exaggerated in that figure for illustrative purposes):
the encoding of the few block motion vectors and the resulting difference image
(with much less range) requires less bits due to the strong temporal coherence
present in usual video sequences.
prediction and the exact shape in a compact bitstream. The most delicate step
in this encoding process is the motion estimation between two consecutive frames.
Indeed, compression rates will be directly affected the quality of the estimation.
Additionally we have to deal with 3D shape matching instead of pixel comparisons,
rendering the issue significantly more complex. One simple strategy is to restrict
our attention to changes in the sign bitmap over time, thus reducing the problem
to motion search in three-dimensional binary images. Then, for each macroblock
in the predicted frame, we can search for the best matching block in the previous
19
frame, within a limited search range. Matching error can be measured as the sum
of absolute differences (SAD) over all the values in a macroblock, just as it is done
in video compression. Alas, in general this approach is not very effective. To see
why a more complex solution is needed, we have to understand the challenges that
are specific to our problem:
• Lack of Correspondence: Since each isosurface is extracted independently,
there are no correspondences between vertex positions of two consecutive
shapes. As a result, even for a simple global translation between two sub-
sequent frames, new vertex positions never reproduce the underlying motion
exactly. Moreover, connectivity often changes from frame to frame, which
makes tracking vertices impossible in practice.
• Frame Difference: In the case of 2D video, the difference between two almost
similar frames F
t−1
and F
t
will typically produce an error image with a lower
dynamic range compared to the original F
t
, that is, an image much cheaper
to encode. Unfortunately, this is no longer the case in our setting. Bear in
mind that when encoding a isosurface, we are basically working with a binary
inside/outside function. Consequently reduction of dynamic range is simply
impossible: the difference between two frames F
t
− F
t−1
is equivalent to the
exclusive “or” between the two. In fact, the na¨ ıve approach of subtracting
two binary images and encoding the resulting contour typically increases the
amount of information to be encoded, as shown in Figure 2.3 (left).
20
• Temporal Resolution: Unlike most video codecs which assume a high frame
rate (typical standards are 24-30 frames/second), here we have to deal with
arbitrary time step, thus arbitrary deformation.
• Improving the geometric rate/distortion curve: Although the sign bitmap
data provides a bounded-error approximation of the shape, overall compres-
sion rates are dominated by the geometric residuals (about 90%). Therefore,
to achieve a significant compression gain, we need to focus on the latter com-
ponent.
2.3.3 General Approach
With the above considerations in mind, we propose the following principles for
compression of time-varying isosurfaces:
• Motion estimation needs to be as accurate as possible, and not limited by
the grid resolution. We therefore cast this task as a local surface registration
problem. As we are only interested in high quality prediction, low-confidence
estimates will be discarded, and static encoding will be used in those regions
instead. Details are discussed in section 2.4.2.
• Instead of attempting to track the isosurface directly, we adopt an indirect
approach in which we track the surface distance field. Then for each voxel
in the narrow band, we use the estimated distance values on its corners to
predict a vertex position. See section 2.4.3 for details.
21
• As motion prediction errors tend to follow a different distribution from the
static residuals, different categories of errors are quantized and encoded sep-
arately, as shown in section 2.4.4.
2.4 Compression through Block-based Prediction
We begin describing our contribution by providing a simplified overview of the
encoding/decoding scheme:
ALGORITHM:
Given two consecutive isosurfacesS
i
andS
i−1
Estimate motion vectorsM fromS
i−1
toS
i
(encoder only)
Encode/decodeM
PredictS
i
based onS
i−1
andM
Encode/decode prediction residuals.
We now describe the various steps of our encoding/decoding algorithm, pointing
out the choices we made to render the technique as scalable as possible.
2.4.1 Isosurface Macroblocks
We chose to treat each frame as an array of macroblocks of user-defined size (e.g.,
8× 8× 8 voxels in most of our tests), and encode one set of motion parameters
per macroblock. Note that while using simple blocks for motion prediction may
22
be very suboptimal in some 3D animation scenarios, it becomes a practical choice
when dealing with time-varying topology.
2.4.2 Motion Estimation
While video coding algorithms typically use pixel-based translational motion of
blocks (i.e., displacements with granularity defined by a fixed fraction of the pixel
size), this fixed granularity is too limiting in our context as our type of data carries
also sub-voxel information. To achieve better motion estimation through accurate
geometric fitting, we change the approach as follows.
Splitting the Grid in Blocks Given an isosurface S
i
to be predicted at time step i,
let B
i
be the set of heterogeneous macroblocks intersected by S
i
. We first partition
the surface S
i
into block-sized regions P
j
such as

j
P
j
= S
i
and P
j
∩ P
k
= ∅ for
any j= k. Thatis,foreachmacroblock b
j
∈ B
i
, we define P
j
the subset of vertices
of S
i
that are in b
j
.
Block Backtracking The goal of the next step is to find a region of the previous
surface S
i−1
(that both the coder and the decoder know at this point) that best
approximate the shape present in a given block. As the deformation field between
two successive frames can be arbitrary, we cannot expect to find a “perfect fit”
between the shape of the surface within a block and a part of the previous surface.
We will therefore only look for a good-enough fit; in other words, we wish to
find a translation T
j
for each block b
j
that minimizes the squared distance metric
23
between T
j
P
j
and a region of S
i−1
(see Figure 2.3). Note that we restrict our
search to translational transformation for the same reasons that video compression
techniques use translation only: this simple temporal domain prediction technique
(known as motion predictive coding in the context of video) is adopted because of
its computational simplicity as well as its coding gain (very few parameters need to
be sent). Since prediction is applied to local blocks in temporally adjacent frames,
the union of translations provides a reasonable approximation.
Block-based Shape Matching To find the best local fit for each P
j
,weemploy
a local registration approach where each macroblock patch is registered against
the previous isosurface through the geometric optimization procedure proposed by
Mitra et al. [85]. This method is a generalization of the Iterative Closest Point algo-
rithm [33], offering improved stability and convergence rate. Although the original
technique is designed for registration of point clouds, we adapted the algorithm to
handle our isosurface shape matching as follows:
1. Given the previous isosurface S
i−1
(as known to the decoder), construct an
approximation

d
2
of the square of the distance to S
i−1
. Toachieve ahigh-
quality approximation, we sample the squared distance to S
i−1
at twice the
original grid resolution.
2. Given the set of vertices P
j
in the macroblock b
j
, find a translation T that
minimizes

d
2
(T(P
j
),S
i−1
). This is done by a Gauss-Newton optimization
procedure, employing a gradient descent step at each iteration [85]. Note
24
that

d
2
gives us a closed form quadratic error potential for every point in
space: as a result, a gradient descent direction can be easily found for every
point q∈ P
j
, and a general minimization step can be computed by solving a
least-squares linear system (see Appendix A for more details).
Algorithmic Details To compute and store the

d
2
field, we use the d2Tree octree-
based structure [85], where each cell stores the polynomial coefficients of the locally
fitted quadratic approximant. The fitting process involves solving a local linear
least-squares system for each relevant cell, making the evaluation of

d
2
the most
computationally demanding step in the encoding process. However, we render the
process efficient by restraining the computations to within a fixed distance from the
block: a perfect match with the previous surface is most unlikely to be at a long
range, but rather nearby the block’s spatial position. Moreover, we implemented a
lazy evaluation of

d
2
, i.e., approximants are evaluated on demand only, thus avoiding
unnecessary computations for cells that are not used during the registration process.
Finally, the least-squares matrix is the same for all non-boundary cells (one other
advantage of only using translations!), and therefore we precompute its inverse
once and for all. These optimizations save a significant amount of computation and
memory, making the approach practical yet still robust.
Checking Quality of Match The best fit is only valid if it is a good-enough esti-
mate; otherwise, sending the motion estimate followed by a large error correction
25
may be require more bits than sending the right shape in the first place. Conse-
quently, in cases when the final registration error of a block exceeds an acceptable
threshold, we discard the prediction and issue a “no-motion” flag instead. The cor-
responding macroblock is then treated by the decoder as an I-block, i.e., encoded
by the static method.
Handling Fast Motion So far we have assumed only moderate inter-frame defor-
mations, i.e., low-speed motion. Naturally, the above algorithm can handle faster
motion by simply increasing the approximation coverage around the surface. In
practice, however, this change turns out to be too computationally extensive. An
alternative approach is to find a coarse motion vector by a simple voxel-based search
(MPEG-style) over the distance field which needs to be evaluated anyway, then re-
fine the motion estimate step using the registration-based local search. Our fast
voxel-based motion search extends the method of Chalidabhongse and Kuo [28] to
the 3D case.
Figure 2.3: Left, for a surface slightly shifted (here, a disk moved to the right), the
difference contour between two time steps is more verbose than each of the steps.
Right, shown is a local (i.e., block-sized) surface motion between two subsequent
frames, and the result of a local registration procedure.
26
2.4.3 Motion Compensation
From the previous frame and the motion vectors, the decoder must now start de-
ducing the new frame: this process is typically referred to as motion compensation
in video coding. We now review the main steps involved in our 3D extension.
Connectivity Encoding/Decoding We first proceed by sending the connectivity
(i.e., the inside/outside function of all grid nodes) of the new frame using the static
encoding described in [74]. Our tests has shown that there is very little to gain
to try to predict connectivity from the motion vector, as the encoding of the full-
blown information only requires a few bits (connectivity typically represents 10%
of the total bitstream). This design choice also allows rapid progressive display
of the decoded isosurface, which is not possible when connectivity and geometry
bitstreams are interleaved.
Vertex Position Prediction Given the correct connectivity, the decoder knows
that each heterogeneous cell needs to find an inside vertex position. As stated
earlier, due to connectivity changes and sampling artifacts we cannot in general
trace vertex positions directly. We therefore proceed on a cell-by-cell basis, trying
to make a good prediction of each vertex position given the previous frame, the
new motion vectors, and the new, correct connectivity. Two cases can happen:
• First, a prediction on the vertex position of a heterogeneous cell is used only
if the predicted distance values in the narrow band agree with the actual
connectivity. Otherwise, we know that any temporal prediction of the new
27
vertex position is pointless: too much deformation has happened locally. Note
that this case happens quite often for fast motion. For those unpredictable
positions, the decoder simply expects the coder to send the correct position
using the spatial encoding used in [74], i.e., using a normal and tangential
components from a prediction purely deduced from the connectivity.
• Otherwise, the decoder first extracts a vertex position x
1
from the predicted
distance function. The final predicted position is obtained as follows. Using
the motion vector, we map x
1
back to the previous frame, and find the vertex
x
2
closest to it within the likelihood ellipsoid. The ellipsoid is defined in
the local frame (as defined by the voxel connectivity) with its radii being
parameters specified a-priori. These correspond to typical standard deviations
of the error in the normal/tangential directions. If there is no such vertex,
prediction is discarded. Otherwise, x
2
is the predicted position. This heuristic
has been, in all our tests, the most reliable way to filter out predictions that
do not help reduce the total entropy. See Figure 2.4 for an illustration.
In fact, whether or not a prediction is available, our tests has shown that the
tangential components of the prediction are mostly noise (therefore, not helpful
for prediction): the decoder thus only uses the normal component of the predic-
tion. Finally, we encode, send, and decode the geometric residuals, i.e., geometry
prediction error, using the static approach from [74] (as we have already exploited
28
temporal coherence), which completes the process. The next section focuses on
arithmetic coding of different data components.
2.4.4 Context-based Arithmetic Coding
Further reduction in bit rates are obtained using context-based coding as explained
next.
Motion Vector Encoding
As motion vectors represent the least dominant component in the overall bit rate,
we use a straightforward encoding scheme. For each macroblock, a single bit flag is
first generated, indicating if the block has a motion vector. Then for each motion
vector, X, Y and Z components are uniformly quantized and encoded separately,
using a context-based arithmetic coder. The context for each component is based
on the values of the already encoded motion vectors of the immediate neighbors.
Connectivity Encoding
Our connectivity encoding scheme uses exactly the same context as the one for
static isosurfaces [74]: based on our tests, this context appears to be superior to
others.
Geometry Encoding
Geometric information, i.e., the positions of the vertices inside heterogeneous cells
can also benefit from context encoding. The results reported in here have been
generated by separating the encoding of the residuals for predicted vs. unpredicted
29
vertices (see Section 2.4.3). Indeed, they show very different distributions, and are
therefore better encoded using two streams. The stream for predicted vertices is
best encoded using a Lloyd quantization on the residual, as the error distribution
shows a significant peak with few outliers. Finally, the context used in both cases
is the connectivity-based, 8-bit context used in [74].
2.5 Results
We are faced with the same issues as the authors from [74] reported: we are not
aware of any existing compression techniques for time-varying isosurfaces, and com-
paring our compression performance with methods for different data modalities
(e.g., volumes) is bound to be biased. This makes the comparisons of our results
quite challenging. For a fair assessment of the performance of our motion-prediction
technique (MP), we compare it to the frame-by-frame approach (FF) where each
frame uses the progressive encoding of [74] with the best published rates for an
isosurface encoder. We focused our tests on highly dynamic models (fast motion
and/or evolution of the surfaces) as it is the most challenging type of motion com-
pression: obviously, a simple translation motion of an object can be achieved with
extremely low rates, but such encoding is not as interesting.
To demonstrate our results, we used three challenging datasets:
• Droplets Dimensions: 257
3
× 20. A fluid simulation with topology changes
(see Figure 2.5).
30
• FluidSplash Dimensions: 257
3
× 20. Another fluid simulation, this time of
a splash with a high topological complexity (see Figure 2.6).
• HeadScan Dimensions: 257
3
× 36. An MRI head scan, sampled at different
isovalues to visualize different tissues (see Figure 2.7, left).
Figure 2.8 illustrates the rate/distortion (R-D) performance of both methods.
Table 2.1 provides some concrete rates for the median quantization parameters.
The results of our dynamic surface encoder on such fast-pace animation exhibit
a typical gain of 12% to 18% over the best frame-by-frame encoding method. It
can be seen from the R-D curves that the highest gain is achieved for lower rates,
reflecting the well known fact that fine detail is hard to encode. On the other
hand, it is precisely the high frequency data that dominates the bit rate, hence the
moderate overall gain.
This leads us to the important question of the most appropriate distortion
measure for isosurfaces. Minimizing the L
2
distance to the originally extracted
isosurface (a usual quality criterion) turns out to be quite costly. The reason is that
we are forced to encode the exact vertex coordinates, as chosen by the extraction
procedure. The actual ground truth, however, is the original volume itself, and
various contouring algorithms exist to extract the actual surface. Therefore, using
the volume data for error assessment may be the right choice for many applications.
Designing a codec with such error in mind would allow more freedom in vertex
placements—thus greatly reducing the bitstream size.
31
Finally, we note that our compression technique is also progressive “in space”:
as more bits come in, the decoder can visualize a progressive refinement of the next
frame. It was our initial goal to send the animation frame by frame to get a nicely
progressive codec; however, we believe that using fully space-time compression,
where the user cannot visualize the whole result until the whole sequence is sent
but for which more freedom in the way we compress the space-time data is available,
could give dramatically better compression. Unfortunately, this approach does not
lend itself to streaming, while our presented method does.
Dataset # frames FF MP Gain
Droplets 20 1528932 1343131 12%
FluidSplash 20 2196948 1841577 16%
Headscan 36 7185869 5935626 18%
Table 2.1: Comparing compression performance of motion prediction vs. frame-by-
frame mode. Rates are shown in bytes.
2.6 Discussion and Future Work
We have presented an approach to encoding time varying isosurfaces inspired by
the simple framework used in video encoding. Mixed with the progressive encoding
method of Lee et al. [74] for static isosurface compression, we provided a codec
allowing rapid, progressive display of the sequence during decoding. The unstruc-
tured nature of time-varying isosurfaces allows our algorithm to deal naturally with
32
surfaces of high geometric complexity and dynamic topology, unlike any of the ex-
isting dynamic 3D geometry codecs. Although beyond the scope of this work, the
method may also prove useful to encoding of additional surface properties, such as
normals, color, texture, or in the context of fluid simulation, density and velocity
fields: the latter is of particular interest here, as it eliminates the need for motion
estimation.
Further potential improvements are numerous: nearly all optimization ideas
that have been proposed for MPEG and H.264, such as the use of bi-directional
prediction (B-frames), variable macroblock size, etc. can all be added to our frame-
work. Our algorithm could also gain from improving the efficiency of the local
registration method. Finally, extending our method to other non-conformal repre-
sentations of time-varying geometry, such as moving point sets and free-viewpoint
video could also be of interest. In the case of point-based geometry, for instance,
although a different, format-specific method is required to encode I-frames, the
motion estimation module of the proposed encoder could still be used as is. Also,
compared to isosurface data, the related content often implies locally more rigid
motion (e.g., articulated deformations) and thus could exhibit favorable prediction
rates. Moreover, in the case where a specific deformation model is known (or can
be inferred), the general block-based motion estimation approach could be replaced
33
with a more efficient model-specific registration scheme. This brings us to the gen-
eral problem of deformable shape matching which is addressed in a wider context
in the next chapter.
Acknowledgements
We are grateful to Haeyoung Lee, Sanjit Patel and Yiying Tong for their insights and
to Santiago Lombeyda, Andreas S¨ oderstr¨ om, Adam Bargteil and Richard Keiser
for their data sets. The research has been funded in part by the Integrated Me-
dia Systems Center (a National Science Foundation Engineering Research Center,
Cooperative Agreement No. EEC-9529152), and by DOE (DE-FG02-04ER25657)
and NSF (CCR-0133983, DMS-0453145) grants.
Appendix A
For completeness, we provide a derivation of the geometric optimization procedure
used for macroblock registration. Due to limited space, we only explain the process
in 2D; details can be found in the original paper [85].
For an isosurface S
i
, the local squared distance approximant

d
2
at point p=(x, y)
is given by

d
2
(p)= Ax
2
+ Bxy + Cy
2
+ Dx + Ey + F,
34
where A, B, C, D, E, F are the approximant coefficients, valid locally around p.
Note that for any point q on S
i
,

d
2
(q) = 0 by definition. Now, given a set of points
{p
j
} on isosurface S
i+1
, we find a translation T that aligns them best with S
i
by
iteratively minimizing

j

d
2
(p
j
+ T). In each iteration (x
j
,y
j
)→ (x
j
+ t
x
,y
j
+ t
y
):
if for a small translation (t
x
,t
y
),

d
2
(p
j
) stays approximately fixed, the residual error
for (p
j
)isgiven by
ε(t
x
,t
y
)=

x
j
+ t
x
y
j
+ t
y
1

W
p
j

x
j
+ t
x
y
j
+ t
y
1

t
,
where W
p
j
is a matrix representing the approximant

d
2
around p
j
. Setting ∇ε to
zero, we obtain the linear system:
⎡ ⎢ ⎢ ⎣

j
⎛ ⎜ ⎜ ⎝ 2AB
B 2A
⎞ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ t
x
t
y
⎤ ⎥ ⎥ ⎦ =
−
⎡ ⎢ ⎢ ⎣

j
⎛ ⎜ ⎜ ⎝ 2Ax
j
+ By
j
+ D
Bx
j
+2Cy
j
+ E
⎞ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎦ This system is solved at each iteration of the gradient descent process, until a
convergence criterion is satisfied or the maximum number of iterations is reached.
The approximant

d
2
(.) is fitted locally by solving another least-squares system. The
35
corresponding matrix does not depend on a global coordinate system and therefore
can be precomputed for most cases.
36
predicted
position
final
position
motion
vector
t-1
t
Figure 2.4: Predicting vertex position. The grey vertex is extracted from the
predicted distance function. The likelihood ellipsoid is aligned with the local frame,
describing the expected error distribution. Using the motion vector, the position is
mapped back to the previous frame and the prediction is refined using the nearest
original vertex within the likelihood ellipsoid.
Figure 2.5: Fluid simulation with falling droplets generating ripples and splashes.
37
Figure 2.6: Another simulation of splashing fluid, with high topological complexity.
Figure 2.7: MRI headscan sequence (the visible human dataset sampled at different
iso-values). Note the amount of noise and topological complexity in the middle
frame.
38
Figure 2.8: Rate/distortion curves for the two fluid sequences and the MRI head-
scan sequence. We compare compression with motion prediction to frame-by-frame
encoding for different quantizations. The x-axis represents the number of bytes,
while the y-axis represents the L
2
distance between the decompressed version and
the original sequence (obtained by summing each L
2
distance for each frame) as
measured by MESH [17].
39
Chapter 3
Generalized Geometric Flows and Applications to
Shape Matching
Errors like straws upon the surface flow:
Who would search for pearls must dive below.
John Dryden
Geometric flows are ubiquitous in geometry processing. Curve and surface evo-
lutions based on functional minimization have been used in the context of surface
diffusion, denoising, shape optimization, minimal surfaces, and geodesic shortest
paths to mention a few. However, such gradient flows are nearly always (and of-
ten implicitly) based on the canonical L
2
inner product of vector fields. In this
chapter, we point out that changing this inner product provides a powerful, sim-
ple, and untapped approach to extend current flows. We demonstrate the value of
such a norm alteration in the context of shape matching, where deformation priors
(ranging from rigid motion to more complex cases such as articulated motion) can
be built in a gradient flow to significantly improve results. Implementation details,
40
including a differentiable approximation of the Hausdorff distance between meshes
and the use of implicit integration to prevent spurious oscillations, are presented
later in this chapter. Finally, to illustrate the versatility of the method in discrete
surface processing in general, we also show additional applications outside of shape
matching.
3.1 Introduction
Geometric flows have been extensively used in mesh processing. In particular,
surface flows based on functional minimization (i.e., evolving a surface so as to
progressively decrease an energy functional) is a common methodology in geom-
etry processing with applications spanning surface diffusion [109], denoising of
scanned meshes [38], shape optimization and surface design [19,71,120], minimal
surfaces [95], (geodesic) shortest paths [18], and animation [53]. Such gradient
flows even extend to image processing where they appear, e.g., in the foundations
of active contours [27] and shape spaces [70,122]. While a large part of the image
processing community implements gradient flows using an Eulerian methodology
(typically, with level sets [90,101]—see also [41,87] in graphics), Lagrangian rep-
resentations of surfaces based on triangle meshes are most common in graphics.
In this Lagrangian setting, discretization of continuous flows is usually achieved
through the use of discrete differential operators and/or using finite element tech-
niques [43, 82, 95]. In this chapter, we point out that most (if not all) surface
41
gradient flows used in mesh processing are based on the canonical L
2
inner product
of vector fields—although this choice is rarely mentioned or acknowledged. In fact,
changing this inner product is shown to be a powerful, versatile, and untapped
approach to design novel generalized Lagrangian gradient flows of discrete surfaces.
Figure 3.1: Articulated matching of an Armadillo template mesh to a very differ-
ent pose, using the Hausdorff gradient flow with a quasi-articulated prior.Our
approach can embed such shape priors into conventional surface flows to provide,
e.g., robust and efficient shape matching with no markers. Notice that although a
template arm initially coincides with a target mesh leg (resulting in a strong local
minimum), the final pose is successfully recovered.
Outline We first give a brief exposition of the mathematical background involved
in the design of conventional gradient flows in Section 3.2, stressing the (often im-
plicit) choice of the L
2
inner product over the space of vector fields on meshes. In
42
Section 3.3, we then propose to alter this L
2
inner product through the introduction
of an additional operator to offer a simple, versatile, and powerful way to design
generalized gradient flows of discrete surfaces—an idea recently introduced in [30]
and [107] in the Eulerian setting. We come back to the applications in Section 3.4,
where we show how deformation priors render generalized flows particularly relevant
to shape matching: we introduce a series of prior-based inner products to develop
a multi-resolution shape matching procedure between two meshes through a gener-
alized pseudo-Hausdorff distance gradient flow. Details related to the mesh-based
implementation are elaborated in Section 3.5.1.
Finally, we show additional applications of the generalized flows approach, be-
yond shape matching. In Section 3.7.1, we demonstrate that the use of Sobolev
inner products is an effective way to regularize any explicit gradient flows, similar
to implicit fairing for the special case of the mean curvature flow, that we enhance
to offer local volume control (Section 3.7.2).
3.2 Mathematical Background on Gradient Flows
We start by giving a brief exposition of the mathematical background involved in
the design of Lagrangian gradient flows. Note that we will refer to the Euclidean
(R
3
) metric as |.|, and to its associated inner product as the conventional dot
product “·” notation. Norms and inner products of vector fields are denoted as
.

and., . instead, to avoid confusion.
43
3.2.1 Continuous Definition of Gradient Flows
Let S be a surface inR
3
(assumed with no boundary for simplicity), and suppose
that we are given an energy functional E as a measure of surface “quality”, i.e.,
the lower the energy, the “fairer” the surface is. A gradient flow is a motion of the
surface that follows the steepest descent of the functional, i.e., the motion follows
the ”gradient” of the functional. However, in the continuous setting, the definition
of the gradient of the energy functional is not straightforward: one cannot use
the traditional concept of partial derivatives along each dimension since we are
dealing with an infinite-dimensional space of deformation (thereafter denoted by
F). Using the calculus of variations, we can however define a theory of curve and
surface evolution [47] using the weaker notion of Gˆ ateaux derivative. The Gˆ ateaux
derivative of E at S (denoted by DE(S)[v] for its value in the direction v ∈ F)
is defined to be the rate of variation of the energy when the surface undergoes a
deformation along a vector field v:
DE(S)[v]
def
=
dE(S+v)
d

=0
= lim
→0
E(S+v)−E(S)

. (3.1)
Inner Product of Vector Fields Recall now that we wish to define a surface
evolution that decreases the energy the fastest. However, the very notion of steep-
ness requires a metric, i.e, a way to measure distances and angles. To that effect,
we must equip our space F of all possible deformation fields with an inner product
44
,
F
. With this additional structure on the space of deformations of a surface, we
can now formally define the gradient ofE relative to the inner product,
F
, denoted
∇
F
E(S), as the vector field such as:
∀v∈ F, DE(S)[v]=∇
F
E(S),v
F
, (3.2)
Here, we assume existence and uniqueness of this gradient for simplicity—these
topics are outside the scope of this article; but be aware that these properties
intimately depend on the choice of energy functional.
Gradient Flow Now that the gradient of the energy functional is properly de-
fined, shape motion via functional minimization can be performed by evolving an
initial surfaceS
0
in the direction of the negative gradient to best reduce the energy,
yielding:
⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ S(0) =S
0
dS
dt
=−∇
F
E(S)
(3.3)
Once again, we assume that this flow is well defined without further analysis or
justification. We note that if this evolution exists, it does decrease the energyE as
expected:
dE(S)
dt
= lim
→0
E(S(t)− ∇
F
E(S))−E(S(t))

= ∇
F
E(S),−∇
F
E(S)
F
=−
∇
F
E(S)

2
F
≤ 0 .
45
Figure 3.2: Quasi-Articulated flow: matching two hands by minimizing their
Hausdorff distance with a gradient flow incorporating an articulated prior (using
16 regions: 3 per finger + palm). The flow automatically morphs the red hand into
the blue one, first mostly through a global rotation, then by adjusting each finger.
Color mix on the final frame is due to the fact that the hands do not have the same
sampling and connectivity.
Choice of Inner Product In the construction of a gradient flow as we just
reviewed, we remained agnostic vis-a-vis the inner product on F. In fact, every
inner product yields its own gradient, and thus, a different gradient flow. However
it is widely assumed, overtly or covertly, that the deformation space F is the Hilbert
space of square integrable velocity fields L
2
(S,R
3
) equipped with the canonical L
2
inner product onS:
u,v
L
2
=

S
u(x)· v(x) dx , (3.4)
where dx is the area element of the surface. Even with this choice of space for F
(the only deformation space that we will consider in this paper), picking a different
inner product is a degree of freedom for gradient flows that is rarely exploited in
geometry processing—to the point where most classical gradient flows reported in
the literature (mean curvature flow, geodesic active contours, etc) improperly refer
to the energy gradient, implicitly assuming the L
2
inner product.
46
3.2.2 Discrete Gradient Flows
We now assume that the surface S is represented by a triangle (piecewise-linear)
mesh X, and the position of its vertices is denoted by x
i
. A piecewise linear vector
field on S is a collection V ={v
i
} of vectors at the vertices, interpolated linearly
over each triangle: v =

i
v
i
φ
i
, where the linear, interpolating basis functions
φ
i
satisfy: ∀x∈S,

i
φ
i
(x) = 1. Note that for this particular setup, a discrete
surface flow exactly coincides to a motion of each vertex x
i
along v
i
.
L
2
Inner Product and L
2
-Gradient With this choice of discrete setup, the
canonical L
2
inner product of two discrete vector fields U and V becomes particu-
larly simple to express:
U,V
L
2
= U
T
MV , (3.5)
where M ={M
ij
} is a symmetric definite positive matrix (called the mass matrix
in the Finite Elements literature) defined through 3×3blocks:
M
ij
=Id
3

S
φ
i
(x)φ
j
(x) dx .
Note that M is sparse, but not diagonal. A classical approximation, simplifying
computations considerably with limited accuracy loss, is to use mass lumping which
turns M into a diagonal matrix where M
ii
is the area of the (Voronoi or barycentric)
dual cell of x
i
times the identity matrix Id
3
[37,82].
47
Discrete Gradient Flow If we call X the vector of all the positions x
i
of the
vertices of the mesh, the discrete counterpart of the L
2
-gradient of a surface energy
functionalE is:
∇
L
2E(X)= M
−1
∂E
∂X
(X) ,
where ∂E/∂X is the derivative of the energy with respect to the position of the
vertices. The discrete L
2
gradient descent then follows straightforwardly:
dX
dt
=−M
−1
∂E
∂X
(X) . (3.6)
Example The most canonical example of surface flow is certainly the mean cur-
vature flow [47]. Its discrete treatment on meshes was first studied in [95], and
more recently in [38], where the flow was formulated as
dx
i
dt
=−
1
2A
i
dA
i
dx
i
, (3.7)
with x
i
denoting the i-th vertex of the mesh, and A
i
denoting the 1-ring surface
area of this vertex. This last term being proportional to the area of the barycentric
dual of x
i
(and thus, proportional to the lumped mass matrix diagonal term M
ii
),
one can see that this discrete flow corresponds to the L
2
-minimization of the area of
the surface—the exact analog of its continuous counterpart, called Laplace-Beltrami
flow.
48
3.3 Generalized Gradient Flows
We now present how to design other inner product than L
2
on the space F of vector
fields before showing how it results in a family of generalized gradient flows.
3.3.1 Designing New Inner Products
Recently, Charpiat et al. [30] have shown for the Eulerian setting that shape
optimization through gradient descent can greatly benefit from the design of an
application-dependent inner product on the deformation space. In particular, they
noted that for any self-adjoint positive-definite linear operator L : F → F,a new
inner product can be defined by
u,v
L
=u,Lv
L
2
=Lu,v
L
2
. (3.8)
This is evidently a very special type of inner product, as it is defined with respect
to the L
2
norm. The advantage is however that if∇
L
2E(S) exists, then∇
L
E(S) also
exists and can be expressed as:
∇
L
E(S)=L
−1
∇
L
2E(S) . (3.9)
49
Indeed, one can verify that:
∀v∈ F, DE(S)[v]=∇
L
2E(S),v
L
2
=

LL
−1
∇
L
2E(S),v

L
2
=

L
−1
∇
L
2E(S),v

L
.
In other words, this procedure is of significant practical interest because it allows
to modify any existing L
2
gradient flow. We now show how to adapt this idea in
the mesh setting.
3.3.2 Generalized Gradients
The discrete counterpart of the linear operatorL is a matrix that we will denote as
L ={L
ij
}. Then, using Eq.(3.8), the L inner product of two discrete vector fields
U and V is expressed through
U,V
L
=LU,V
L
2
= U
T
MLV. (3.10)
Thus theL-gradient of the energy can be easily expressed as
∇
L
E(X)=(ML)
−1
∂E
∂X
(X)= L
−1
∇
L
2E(X) .
Note that we made use of the self-adjointness of the operator L,which in our
discrete setting means that the ML matrix is symmetric (so that U,LV
L
2
=
50
U
T
(ML)V = U
T
(L
T
M)V = LU,V
L
2
). In fact, the discrete L matrix is not
necessarily symmetric; however, left-multiplying by the mass matrix M will make
the final product symmetric. The reader can verify that all the inner products
proposed in this text necessarily satisfy this property although their corresponding
matrices L may not be symmetric by themselves.
3.3.3 Time Integration
Once our generalized gradient is defined, the associated generalized flow is obtained
via time integration. A direct implementation of an explicit Euler time stepping
scheme for a generalized gradient flow yields:
X
t+dt
= X
t
− dt (ML)
−1
∂E
∂X
(X
t
) . (3.11)
We will see that for most of our proposed generalized flows, the inverse of ML is
given analytically (see, e.g., Eq.(3.15)), rendering the computational cost of time
integration negligible. The only exception is the prior in Section 3.4.3.3, where
(ML)
−1
has no closed form expression. Even in this case, integration can still be
done efficiently: the matrix ML is sparse, meaning that multiplying it by a vec-
tor is not a quadratic, but a linear operation. Thus, a preconditioned conjugate
gradient method (PCG) [13] will be, in practice, very efficient, as is the case in
Section 3.4.3.3. We can also use implicit integration for our generalized gradient
flows through a linearization of
∂E
∂X
with respect to the position of the vertices.
51
This modification, tantamount to keeping the metric of the surface constant over
the integration step, has been advocated for the mean curvature flow in [38] as an
easy way to remove the usual time step constraints inherent to explicit integration
techniques and add robustness to mesh flows. We will suppose that such a lineariza-
tion exists (as will be the case in all of our flows) such that:
∂E(Xt)
∂X
= A
t
X
t
+ B
t
,
where A is diagonal, positive definite matrix (the same argument can be made for
sparse SPD matrices). Then, Eq.(3.11) can be rewritten by using the position of
the vertices at the next time step in the evaluation of the gradient:
X
t+dt
= X
t
− dt (ML)
−1
(AX
t+dt
+ B) .
That is, one can find the next mesh X
t+dt
by solving the following linear system:
(ML + dtA)X
t+dt
= MLX
t
− dtB . (3.12)
Interestingly, as we will show in Section 3.7.1, this implicit integration of flows can
be seen as a special case of the explicit integration of a regularized gradient flow.
52
Figure 3.3: Incorrect matching between two fingers obtained through a straight-
forward L
2
gradient flow of their Hausdorff distance. The absence of rigidity
priors leads to significant spurious sliding and mismatch. The intermediate posi-
tions also exhibit large deformation.
3.4 Application: Robust Shape Matching
3.4.1 Background
A particularly interesting application of surface flows is shape matching (or registra-
tion). In fact, due to the notorious multitude of local minima in object registration
scenarios, this task can be seen as a ”stress test” for any geometric optimization
method. One way to formulate is problem is as follows. Suppose that we have two
3D models (a template and an instance) of a same subject (be it an object, a whole
human body, or just an organ) in different poses with no known correspondence
(this is typically the case for most models acquired from 3D scanners and medical
imaging). The shape matching problem amounts to recovering a relevant mapping
(i.e., induced by a plausible deformation) between these two different poses—see
53
Fig. 3.2 for an example of properly established correspondence. For shapes given
as meshes or point clouds, matching can be done using only global registration if
the two shapes are quite similar [16, 32, 50, 86] —see Planitz et al. [96] for a de-
tailed survey—or using more complex approaches allowing for the matching of less
similar shapes [34, 39, 122]. For widely different shapes, an alternative is to rely
on Markov networks, approximate inferences, and geodesic nearness constraints to
automatically compute Correlated Correspondences [8], and use these established
correspondences as markers for more robust registration [9]. When no training set
is available, a common procedure is to evolve the template shape, often initially
positioned by the user, in the direction of steepest descent of an energy measur-
ing the distance to the instance shape. However, due to the highly non-convex and
non-linear nature of most energy functionals, an unsupervised gradient descent flow
gets trapped in irrelevant states corresponding to local minima, seriously limiting
its applicability and efficiency (see Figure 3.3). The use of markers as in [5,7,106]
can help tremendously, but requires (potentially heavy) user supervision. Alter-
natively, for the class of near-isometric deformations (such as articulated motion),
generalized multidimensional scaling (GMDS) has been proposed [23,24] based on
the MDS method [44]. A mapping is formulated as a parametrization of one surface
onto the other such that the distortion of pairwise geodesic distances is minimized.
Yet without any prior knowledge about the shape structure, an optimal mapping
can be extremely hard to recover, as shown in Fig. 3.4.
54
Our Prior-based Approach We now explore an alternative where our general-
ized gradient flows are beneficially applied to embed priors (such as quasi-rigidity,
quasi-isometry, or quasi-articulation) in the modeling of the gradient flow. The use
of priors, often used in object tracking or shape fairing/completion [40], consists
in introducing a measure of the plausibility of the different configurations, and bi-
asing the optimization process towards the most probable ones —thus improving
robustness of the estimation.
3.4.2 Hausdorff Gradient Flow for Matching
A typical shape matching algorithm tries to deform the template to match the
instance by reducing a geometric distance. Because a L
2
-type distance is known
to be too forgiving in comparing two shapes (as well as prone to numerous local
minima), we resort to an approximation of the Hausdorff distance instead, based
on a method recently introduced in [29] in the context of level sets. The exact
symmetric Hausdorff distance d(S,T ) between surface S andT is
d(S,T)=max

max
p∈S
min
q∈T
|p− q| , max
q∈T
min
p∈S
|p− q|

Since this expression is not differentiable, we use instead a pseudo-Hausdorff dis-
tance d
H
(X,Y) between two distinct meshes X = {x
i
}
i=1..P
and Y = {y
j
}
j=1..Q
.
We refer the reader to Section 3.5.1 for the actual derivation of this distance, as
it significantly differs from the original Eulerian description in [29] to allow for a
55
robust treatment of irregular sampling. The L
2
gradient flow based on our approx-
imate Hausdorff distance is then written as
dX
dt
=−M
−1
∂d
H
∂X
(X,X
∞
) , (3.13)
where X
∞
is the instance mesh, and the mesh X is initialized as the template
mesh. As discussed previously, such a naive minimization is unable to yield relevant
correspondences between two dissimilar shapes, as the energy landscape is too
complex and non-linear to avoid getting stuck in one of the numerous local minima
(see Fig. 3.3). While other similarity measures could be added to the Hausdorff
distance for better matching (see the use of the normal field in [35,98,123]), we will
refrain from discussing them here for simplicity.
3.4.3 Constructing Prior-based Inner Products
To increase efficiency and robustness of shape matching as well as further exemplify
the versatility of generalized flows, we present in this section three different possible
alterations of a gradient flow.
3.4.3.1 Quasi-Rigid Deformation Prior
The case of quasi-rigid deformation prior, where a L
2
-gradient flow is altered by
penalizing non-rigid motion, was introduced in Charpiat et al. [30], and we adapt
it here to our Lagrangian setting due to its relevance to our goals. To that effect,
56
Figure 3.4: GMDS correspondences (courtesy of [23]) between two different Ar-
madillo poses, for a subset of 100 vertices, with the corresponding Voronoi regions
shown in matching colors. While most regions are mapped adequately (e.g., the
tail, shown with a black arrow), hands and legs still have many inaccurate cor-
respondences (e.g., those shown with red arrows), which will result in unnatural
deformations if further used for registration.
we denote byA the surface area of a given shapeS,andby c =
1
A

S
x dx its center
of mass. With these two definitions, we decompose the space of rigid motions into
the subspaces of translations, of instantaneous rotations around the center of mass,
and of global scaling,defined as
T =

x→ t| t∈R
3

,
R =

x→ ω ω ω× (x− c)| ω ω ω∈R
3

,
S =

x→ s

x− c

| s∈R

.
These subspaces are mutually orthogonal for the L
2
inner product. We finally
57
Figure 3.5: Quasi-Rigid Prior using Multiresolution: Matching two full body
meshes with different sampling (∼ 28K vertices each) by minimizing their Haus-
dorff distance with a quasi-rigid gradient flow. To speed up computations, a multi-
resolution scheme was employed (total time: 1 min). The coarsest resolution is
shown in the middle. Closeups show the final alignment (template mesh in wire-
frame).
denote by N the orthogonal complement of these subspaces, i.e. the space of
strictly non-rigid motions. Thus, the space F of all deformations satisfies:
F = T ⊥ R⊥ S ⊥ N.
A new inner product, altering the respective ”strength” of the different compo-
nents of a deformation, can now be expressed by following the general procedure
of Eq.(3.8) and introducing:
L = λ
T
Π
T
+ λ
R
Π
R
+ λ
S
Π
S
+Π
N
, (3.14)
58
where Π
T
,Π
R
,Π
S
,Π
N
denote orthogonal projections on different subspaces, and
λ
T
, λ
R
, λ
S
> 0 are penalty factors of translations, rotations, scalings, respectively.
The lower a penalty factor is, the stronger its corresponding type of motion will be
in the new gradient flow. Naturally, when all the penalty factors are equal to 1,
L is the identity and the new inner product coincides with the L
2
inner product.
Note also thatL fulfills our requirements of being self-adjoint positive definite and
invertible, guaranteeing a proper generalized gradient flow.
Mechanical Analogy Interestingly, this inner product has a straightforward
mechanical interpretation. We start by defining the following quantities for a given
shape S and a deformation field v (we adopt the naming conventions used in me-
chanics for these quantities):
• Discrete Moment of Inertia:
I =

S
!

x− c

2
I
3
− (x− c)(x− c)
T
"
dx
• Discrete Angular Momentum:
Ω =

S
(x− c)× v(x) dx ,
where I
3
is the 3x3 identity matrix. Note that these quantities depend on the current
shape and deformation field, although we have dropped the time dependency in our
59
notations for conciseness. Note also that the integrals become discrete weighted
sums in our case, since v is piecewise-linear.
Now, the squareL-norm of a velocity field v equals to:

v

2
L
= λ
T
A|t|
2
+ λ
R
ω ω ω
T
Iω ω ω + λ
S
tr(I)
3
s
2
+
Π
N
v

2
L
2
,
where t, ω ω ω and s represent the translational, rotational, and scaling part of the
motion, respectively. The above expression is best interpreted as a weighted sum
of kinetic energies. In particular, the first term is the translational kinetic energy
associated to the motion of the centroid, the second term is the rotational kinetic
energy, and the last term is the kinetic energy of the purely non-rigid residual
motion.
Expression of New Gradient The new gradient ∇
L
E(S) can be obtained
from the L
2
-gradient using Equation (3.9), where:
L
−1
=
1
λ
T
Π
T
+
1
λ
R
Π
R
+
1
λ
S
Π
S
+Π
N
=Id+
1− λ
T
λ
T
Π
T
+
1− λ
R
λ
R
Π
R
+
1− λ
S
λ
S
Π
S
. (3.15)
Thus, the newL-gradient is now simply a linear combination of the L
2
-orthogonal
projections of the L
2
-gradient on T, R and S. Each of these projections are found
60
through the following equations, with the integrals becoming weighted sums inth
discrete case:
∀v∈ F, (Π
T
v)(x)=
1
A

S
v(y) dy ,
(Π
R
v)(x)= I
−1
Ω × (x− c) ,
(Π
S
v)(x)=

S
v(y)· (y− c) dy

S
|y− c|
2
dy
(x− c) .
A similar procedure to isolate the instantaneous affine deformation in order to favor
affine motions is also presented in [30].
3.4.3.2 Quasi-Articulated Deformation Prior
We now extend the idea introduced above to motion of articulated bodies with
quasi-rigid links. This prior will be particularly useful when used on mapping
the hand in Figure 3.2 onto another deformed hand. The bone structure within
each finger severely restricts the possible spectrum of deformation that a hand can
usually undergo, and building this knowledge into the gradient flow will significantly
improve the robustness of the matching procedure.
We will assume that a segmentation of S =

k=1..K
S
k
into K parts is known,
and that the articulation connectivity graph is a tree: G ={(S
i
,S
j
)|S
i
is a predecessor of S
j
}.
Assume also that for each S
k
a corresponding center of rotation (joint) c
k
is given
(in practice, we estimate joint locations trivially by computing centroids of bound-
aries between connected parts). We will denote by R, A and N respectively the
61
Figure 3.6: Quasi-Articulated Prior: Matching two arms with the Hausdorff
gradient flow incorporating a quasi-articulated prior (using 18 subparts: 16 for the
hand, plus the upper and lower arm). Again, the mesh sampling of both arms is
very different. Despite the presence of fingers in different positions (adding more
local minima to the distance-based energy), our algorithm finds the right matching
through a combination of global and local rotations, with a small amount of non-
rigid deformation (for instance at the elbow). Final alignment of the fingers is
shown in closeup on the right.
orthogonal subspaces ofRigid,Articulated andNon-rigid motions of the full object.
Similarly, R
k
, A
k
and N
k
denote the subspaces of rigid, articulated and non-rigid
motions of the k
th
articulated part. Notice that since we are dealing with connected
links, the translation and scaling parts of R
k
(described in the previous section)
are irrelevant, and only the rotation part R
k
is needed. Note also that each A
k
depends solely on rotations—not only about its own center c
k
, but also about its
predecessors’ joints. Thus, the projection Π
A
k
on the articulated motion subspace
A
k
is best computed via a recursive formula through the tree:
⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Π
A
1
=Π
R
1
Π
A
k
=Π
A
pred(k)
+Π
R
k
(Id− Π
A
pred(k)
)
, (3.16)
62
where pred(k) denotes the predecessor of link k in the tree structure and the root
index is defined as 1. We can now obtain a weighted inner product by first isolat-
ing the rigid component of a deformation, then applying to the non-rigid residual
an {articulated ⊕ non-rigid} decomposition of each part S
k
of the object taken
separately. Thus, the square L-norm of a deformation field v is defined as:

v

2
L
= λ
Π
R
v

2
L
2
+

k

μ
k

Π
R
k
Π
N
v

2
L
2
+
Π
N
k
Π
N
v

2
L
2

,
where λ (resp. μ
k
) is the penalty factor of the rigid motion of the full object (resp.
of part S
k
). These penalty factors typically fulfill ∀k, λ < μ
k
< 1, so that rigid
motion is favored over articulated motion, while the latter is favored over non-rigid
motion. Finally, the linear operatorL corresponding to the norm in Equation (3.17)
becomes:
L = λ Π
R
+

k

μ
k
Π
A
k
Π
N
+Π
N
k

(3.17)
and the corresponding closed form inverse is:
L
−1
=
1
λ
Π
R
+

k

1
μ
k
Π
A
k
Π
N
+Π
N
k

. (3.18)
63
Extension to Hierarchical Linked Structures We remark that the above ideas
can be easily generalized to a hierarchical structure, i.e., a multi-level decomposition
of the object, with different penalty factors both for the different parts and for the
different levels. This approach allows to prioritize global motions over local ones,
to obtain an ”as-rigid-as-possible” articulated deformation. For instance, a hand
can be a subtree in a body hierarchy, represented by single component on level k,
but split into a palm and fingers on level k + 1. This can be achieved with another
recursive construction:

v

2
L,S
= λ(S)
#
#
Π
A(S)
v
#
#
2
L
2
+

S

⊂S
#
#
Π
N(S)
v
#
#
2
L,S
,
whereA(S) is the space of articulated motion on a partS containing subpartsS

k
,
and λ(S) is the corresponding penalty factor. To end the recursion correctly, we
define:

v

2
L,S
= λ(S)
#
#
Π
A(S)
v
#
#
2
L
2
+
#
#
Π
N(S)
v
#
#
2
L
2
for the leaf parts of the structure (i.e., the quasi-rigid deformation prior of Section
3.4.3.1).
Handling Joints A shortcoming of the above articulated priors is that they do
not directly penalize disagreement between the motion of different parts at the joints
(i.e., the nodes of the tree). As a result, a gradient flow derived from these priors
may be discontinuous across parts (i.e., at the joints). To overcome this problem,
64
we use a simple partition of unity of the shapeS rather than a strict partitioning:
by modeling the points close to the joints as simultaneously belonging to several
parts with different probabilities (”blending”), we recover the smoothness of the
deformation. Note that while in our case the probabilities are obtained simply by
a local blending of the parts around the joints, this approach generalizes to any
probabilistic segmentation if available.
3.4.3.3 Quasi-Isometric Prior
If no a-priori knowledge is available on the range of deformation that a particular
shape can undergo, a reasonable assumption valid for a a wide class of models is
quasi-isometry, i.e., that a surface deforms such that Euclidean distances between
pairs of surface points, sampled from the whole surface or from a neighborhood
of a certain size, are nearly preserved (this assumption implies a near-preservation
of both angles and areas at the scale of the chosen neighborhood). For instance,
articulated body motion is approximately isometric, but not vice versa.An isometric
deformation is simple to define mathematically. Let T
x
S denote the tangent plane
to S at x. The condition that a deformation field v must fulfill to preserve length
on the shape (i.e., to be an isometry) is:
∀x∈S, ∀τ ∈ T
x
S,τ·
∂v
∂τ
=0 .
65
Figure 3.7: Quasi-Isometric Prior: Matching two arms by minimizing their
Hausdorff distance with a quasi-isometric gradient flow. Here, without any prior
knowledge on the leg structure, the quasi-isometric gradient flow finds the appro-
priate matching. The last frame shows both meshes, to demonstrate the quality of
the match despite the sampling disparity between the two models.
A continuous inner product could be derived from this condition. For simplicity,
we present a discrete version directly. This prior bears similarity in its goals to the
“semi-local rigidification” from Charpiat et al. in [31] in the context of level sets,
but it differs from this previous work in its foundations and realization. A closely
related approach was recently and independently suggested in [69].
Gradient for Quasi-Isometric Prior First, we note that the length preserva-
tion of an edge e
ij
=(x
i
,x
j
) is expressed as
(v
i
− v
j
)· (x
i
− x
j
)=0 .
66
From this condition, we can derive the following inner product of two discrete vector
fields U and V directly:
U,V
L
=U,V
L
2
+ λ
iso

e
ij
w
ij
[(u
i
− u
j
)· (x
i
− x
j
)] [(v
i
− v
j
)· (x
i
− x
j
)] ,
where λ
iso
is a penalty factor, and{w
ij
} are the edge-associated cotangent weights,
commonly used to define a discrete surface metric [95]. In simple terms, this weight-
ing is required to ensure that the above discrete sum converges to a continuous
surface integral as the sampling is refined. The associated operator is:
(MLV)
i
= M
ii
v
i
+ λ
iso

j∈N
1
(i)
w
ij
[(v
i
− v
j
)· (x
i
− x
j
)] (x
i
− x
j
) ,
where N
1
(i) denotes the one-ring neighborhood of vertex i,and M
ii
is, again, the
corresponding lumped mass matrix diagonal term. This prior can be optimized by
adapting the size of the neighborhood and the penalty factor λ
iso
to the desired
degree and scale of local rigidity, such that the ML matrix is still sparse and thus
can be inverted rapidly using PCG. We then obtain a computationally efficient
and versatile prior, which encompasses articulated motion without the need for
a segmentation of the shape. Note that constraining local Euclidean distance is
generally not sufficient to enforce a quasi-rigidity of the shape. Thus, one can add
67
fake edges (joining a few opposite vertices in each limb) with our same isometric
prior when added rigidity is desired.
3.4.4 Results of Prior-based Shape Matching
We show in a series of numerical experiments that upgrading the gradient flow of
the approximate Hausdorff distance with our deformation priors of Section 3.4.3
yields an efficient method for shape matching. Note that in all the following ex-
amples, the sampling of the template and of the query meshes are (purposely)
completely different. Unless specified otherwise, the meshes shown in this paper
were obtained using E-Frontier’s Poser. To accelerate convergence by several or-
ders of magnitude, we implemented a multiresolution approach for the articulated
examples. The algorithm precomputes coarsified versions [48] of both the tem-
plate and the instance mesh and run the Hausdorff gradient flow on those simpler
meshes, before refining them back to the original resolution - using pyramid co-
ordinates [103] - for final alignment, as illustrated by figure 3.5 (see Section 3.5.2
for further details). This multiresolution implementation results in reasonable tim-
ings: most flows shown in this paper are done in a minute or less, the fastest being
the hand in Figure 3.2 (10s.), while aligning the whole human body in Figure 3.5
takes 1 minute for 28K vertices. Note also that we used implicit integration for
stability. Figs.3.8-3.6 demonstrate the ability of our quasi-articulated prior of Sec-
tion 3.4.3.2 to generate a relevant matching between body parts in different poses
68
Figure 3.8: Non-rigid Shape Matching: Here, an articulated template model of
a lamb (red) is automatically deformed to fit a different pose (blue), by directly min-
imizing their Hausdorff distance with a gradient flow based on a quasi-articulated
prior. To handle this example (∼ 22000 vertices for each mesh) in less than two
minutes, a multi-resolution scheme was used. Notice that our prior automatically
makes the lamb move mostly rigidly first to best fit the target pose, then adjust the
limbs non-rigidly. The color mix on the final frame shows that, even if both lambs
do not have the same sampling and connectivity, their geometry matches well.
(using λ
R
= λ
T
varying from 5e-3 to 5e-4, and the μ’s decreasing by a factor of
10 per level in the hierarchy). This prior can be used as soon as a segmentation
(even a coarse one) of the template shape is available. In Figure 3.7 (λ
iso
= 50),
we show that our quasi-isometric can nonetheless accommodate articulated bodies
with different poses without known segmentation, by enforcing the rigidity of the
deformation only at a local scale. Finally, in Figure 3.9, we show the applicability
of our method to noisy and incomplete meshes, with the example of a full arm and
shoulder. Note that, although we obtain satisfying results in this experiment, the
symmetric Hausdorff distance is not naturally suited to partial matches. Yet, our
approach was able to resolve this matching procedure.
69
Comparison to Related Work The results in this chapter are intended to show
a ”stress test” for the generalized flows methodology rather than the panacea of
shape matching. Nevertheless, a closer look at shape matching literature points out
the significant benefits that generalized flows bring. Several techniques [5,26,92]
essentially employ various non-rigid variants of ICP; one known shortcoming of
ICP is its reliance on closest point-pairs of estimated correspondences, which are
often erroneous unless the two meshes are pre-aligned. A common way to improve
robustness is to add a relaxation term, affecting the quality of the final match.
Our flow approach bypasses these issues due to the global nature of the Hausdorff
distance functional, without sacrificing the quality of the final match. Allen et al.
focused on human body matching [5] with fairly similar poses, using a combination
of locally-affine transformations (similarly to [45]) while enforcing a smoothly vary-
ing deformation field. To our knowledge, this setup does not have a simple way
to build prior knowledge into the deformation, and rigid motion cannot be clearly
separated (the same holds for [92]), necessitating the use of markers for each new
object to avoid local minima. By comparison, our quasi-articulated prior only re-
quires a one-time effort of segmenting the template mesh. Another popular way to
constrain the deformation is via thin-plate splines [26,34]. Interestingly, adding a
thin-plate energy term is similar to using the high-order Sobolev-like gradient (Sec-
tion 3.7.1) of the original energy, making it a special case of our proposed setup.
We will come back to a related deformation prior in the next chapter, where we
70
Figure 3.9: Partial Matching of an incomplete scanned 3D mesh (in blue, courtesy
of U. of Washington) with a template arm in a different pose. Despite the sparsity of
the scan, a quasi-articulated prior (19 subparts, including fingers) on our Hausdorff-
based flow gets a good match.
focus specifically on the problem of cortical surface matching and show how the
proposed framework can be leveraged effectively to tackle this challenging task.
3.5 Implementation Details
3.5.1 Pseudo-Hausdorff Distance and its Gradient
We define our approximation of the Hausdorff distance between two distinct meshes
X ={x
i
}
i=1..P
and Y ={y
j
}
j=1..Q
as follows:
d
H
(X,Y)=
$
1
P

i
M
x
ii
f
−1
i
+
1
Q

j
M
y
jj
g
−1
j
% 1
2α
− (3.19)
with f
i
=
1
Q

j
M
y
jj
d
−α
ij
,g
j
=
1
P

i
M
x
ii
d
−α
ij
and d
ij
=|x
i
− y
j
|
2
+
2
,
where M
x
ii
and M
y
jj
are elements of the diagonal mass matrices of X and Y, respec-
tively (required to compensate for irregular sampling), >0and α ≥ 0. Proving
that the above expression converges to the Hausdorff distance between the two
meshes when α → +∞ and the sampling of the two meshes increases follows the
71
continuous proof of [29]. An additional argument necessary to the proof is the
elementary fact that for any set of positive values{u
i
}
i=1..N
,
lim
α→+∞
&
1
N
N

i=1
u
α
i
'
1/α
=max
1≤i≤N
u
i
.
Note however that our Hausdorff approximation does not satisfy the triangular
inequality. Nevertheless, it is symmetric with respect to the two meshes (while only
sided distances are often used in practice for other methods, see for instance [93]),
and only takes positive values. Finally, it has a decisive advantage over the real
Hausdorff distance expression: it is differentiable with respect to the position of
the vertices of the two meshes. Specifically, we have:
∂d
H
(X,Y)
∂x
i
=
(d
H
(X,Y)+ )
1−2α
P · Q
M
x
ii

j
x
i
− y
j
d
α+1
ij
M
y
jj
(f
−2
i
+ g
−2
j
) . (3.20)
Computational Issues The complexity of computing d
H
or its gradient is O(P·
Q), which can be prohibitive when using large datasets. In practice, we restrict
the sums in f
i
and g
i
(Eq.(3.19)) to only the -nearest neighbor pairs, reducing the
complexity to a more tractable O(P)+ O(Q) for coarsely aligned shapes, without
a noticeable loss of accuracy. The key observation from Eq. (3.19) is that for each
x
i
∈ X, the term f
i
is mostly determined by the set of points y
j
∈ Y closest to
72
x
i
(a similar claim holds for g
j
as well). Formally, we can define this set as the
-nearest neighbors of x
i
in Y, i.e.,
Y
i
={j|y
j
∈ Y and
x
i
− y
j

≤ (1 + )min
y
k
∈Y

x
i
− y
k

},
where is a ”padding” factor. The larger α is, the smaller can be (in our imple-
mentation, we used α∈ [2..4] and ∈ (0..1]). Informally, we will say that Y
i
is the
set of all points in Y that affect x
i
. Yet, to design a proper approximation of
(
d
H
,
we need to distinguish between two different kinds of neighborhoods:
− →
Y
i
= {j|y
j
∈ Y affects x
i
},
← −
Y
i
= {j|y
j
∈ Y is affected by x
i
}.
Note that in general,
− →
Y
i
=
← −
Y
i
, as seen on Figure 3.10: for a small , x
i
affects by
y
j
but is not affected by y
j
.
Thus, in Eq. (3.19), the term f
i
−1
approximating min
y
j
∈Y

x
i
− y
j

can now be
x
i y
j
d
Figure 3.10: Illustrating asymmetric proximity: x
j
is among the nearest neighbors
of y
j
, but not vice versa.
73
written as
f
i
≈
1
A(|
− →
Y
i
|)

j∈
− →
Y
i
M
y
jj
d
ij
−α
,
whereA(.) denotes the surface area associated with the give vertex set. Doing the
same for g
j
, we can rewrite Eq. (3.19) as
∂
(
d
H
(X,Y)
∂x
i
=
(
(
d
H
(X,Y)+ )
1−2α
$
M
x
ii
A(P)
f
−2
i
∂f
i
∂x
i
+
1
A(Q)
Q

j=1
M
y
jj
g
−2
j
∂g
j
∂x
i
%
=
(
(
d
H
(X,Y)+ )
1−2α
M
x
ii
⎡ ⎣ f
−2
i
A(P)

j∈
− →
Y
i
M
y
jj
A(|
− →
Y
i
|)
x
i
− y
j
d
α+1
ij
+
1
A(Q)
Q

j=1
M
y
jj
A(|
− →
X
j
|)
g
−2
j
x
i
− y
j
d
α+1
ij
⎤ ⎦ .
At first glance, the computational complexity for the second term is stillO(Q), i.e.,
O(P · Q) overall complexity to compute the gradient. Yet actually, that summation
now becomes sparse, as
∂g
j
∂x
i
is non-zero only for j∈
← −
Y
i
. And hence,
∂
(
d
H
(X,Y)
∂x
i
=
(
(
d
H
(X,Y)+ )
1−2α
M
x
ii

f
−2
i
A(P)

j∈
− →
Y
i
M
y
jj
A(|
− →
Y
i
|)
x
i
−y
j
d
α+1
ij
+
1
A(Q)

j∈
← −
Y
i
M
y
jj
A(|
− →
X
j
|)
g
−2
j
x
i
−y
j
d
α+1
ij

.
For bounded inter-surface distances (such as when the two shapes are in a coarse
alignment), the size of the -neighborhood for each point is also bounded, and
therefore the above expression has a constant complexity, amounting toO(P + Q)
time to compute the full gradient. On the other hand, the worst case complexity
(when all points are far away from the opposite surface) is still O(P · Q). Further
74
elimination of the computational bottleneck is made possible through the use of
multi-resolution, as described next.
3.5.2 Multiresolution Flows
3.5.2.1 General Approach
The surface matching approach presented so far suffers from one practical drawback:
evaluating the approximate Hausdorff energy gradient may still be ofO(P · Q)com-
plexity in the worst case (see Section 3.5.1), thus becoming a bottleneck in complex
scenarios. A more general observation concerning any flow is that excessive geomet-
ric detail often introduces a multitude of local minima in the energy minimization
process. A natural solution, therefore, is to apply the method in a multiresolution
manner, both for robustness and computational efficiency. Intuitively speaking,
while the main bulk of vertices in any given mesh merely convey high-frequency
detail, this detail is typically not needed to obtain an approximate solution. Specif-
ically in the shape matching setting, our modified algorithm can be summarized as
follows:
ALGORITHM:
Given two shapes to be registered,
1. Faithfully simplify both shapes.
2. Find a coarse alignment.
75
3. Refine shapes.
4. Improve alignment.
Steps (2) and for (4) are performed using the basic method. Steps (1) and (3)
require a multiresolution shape representation that allows one to reduce the given
surface complexity to a desired level of detail, then reconstruct the shape back to the
original scale once some deformation has been applied to the coarse configuration.
These details are discussed in the following section.
3.5.2.2 Multiresolution Shape Representation
Simplification Strategy Since we use polygonal meshes to represent shapes in
this work, our choice for mesh simplification method is the well-known Garland and
Heckbert scheme [48]. Basically, a given mesh is simplified to a desired degree by a
series of edge collapse operations, where a quadric-based error metric is employed
to minimize distortion. Since our surface matching procedure strives to minimize
the point-based Hausdorff distance between two given surface, it is preferable that
the coarse mesh vertices be a subset of the original vertex set, i.e., vertices that are
retained in the simplification process do not change their positions (a.k.a. half-edge
collapse simplification). Conveniently, the Garland and Heckbert scheme supports
half-edge collapse operations just as well. After each (half-)edge collapse, the local
connectivity change is stored in the simplification history.
Depending on the specific deformation prior, the simplification process may need
76
to be further adapted: e.g., when dealing with an articulated template mesh, it
is important indeed to keep enough vertices at the ”joints”, so that the coarsened
shape can still bend like an articulated object. In that case, we boost the priority of
those vertices by an importance factor: vertices that belong to more than a single
shape component become less likely to be decimated.
Figure 3.11: Human body mesh, at two levels of detail: the size of a full resolution
surface (left) is reduced to 10% of the original.
Shape Reconstruction Once an initial alignment has been found, a full-scale
shape has to be reconstructed from the deformed coarse mesh and the simplifica-
tion history. In this case, a valid shape reconstruction method should satisfy the
following requirements:
• Local Support: The position of each decimated vertex should be encoded
locally, with respect to nearby vertices. Locality is an important property for
any hierarchical representation.
77
• Rigid Transformation Invariance: whenever the coarse mesh undergoes merely
a rigid transformation, the reconstructed shape should reflect the same trans-
formation with respect to the original fine-scale configuration.
• Interpolation Property: The deformation of each reconstructed vertex should
be an interpolation of the coarse deformation field of the control mesh.
With these properties in mind, we chose a variant of pyramid coordinates [103]
for hierarchical shape representation. Essentially, pyramid coordinates locally cap-
ture angles and lengths by representing each vertex with respect to its neighbors
in the mesh. This representation consists of mean value coordinates [59] in the
tangential dimension and a scalar offset in the normal direction. As our construc-
tion scheme slightly differs from the original one, the complete algorithm is stated
below. See Figure 3.12 for an illustration.
V
i
V
h
V'
i
V'
w
i
Figure 3.12: Local pyramid coordinates. A vertex v along with its one-ring neigh-
bors{v
i
} are projected onto the local tangent plane, where mean value coordinates
{w
i
} of v are computed. The vertical offset h encoded separately.
78
Algorithm: Given a mesh vertex v be inR
3
,
1. Let N
1
(v)={v
0
, ..., v
k−1
} be the set of its neighboring vertices, ordered with
a consistent orientation, and let p
b
be the barycenter of N
1
(v).
2. The local normal n =(n
x
,n
y
,n
z
) is defined as the area-weighted normal of
the triangle set {(v
i
,v
i+1
,p
b
)}
k−1
i=0
.
3. We then define the local projection plane as P = n
x
x + n
y
y + n
z
z + d ,
d =−n· p
b
.
4. Let v

and{v

0
, ..., v

k−1
} denote the projections of v and{v
0
, ..., v
k−1
}, respec-
tively, on P,andlet h=(v− v

)· n be the offset in the normal direction.
5. Finally, let l =
1
k

k−1
i=0
|v

i
− v

| and β = h/l.
6. Then the pyramid coordinates ψ(v) are defined as
ψ(v)={(w
0
, ..., w
k−1
),β},
where (w
0
, ..., w
k−1
) is the vector of normalized mean value coordinates of v

w.r.t. {v

0
, ..., v

k−1
} (computed using the algorithm from [59]).
During the reconstruction step, the position of v is restored as follows:
v =
&
k−1

i=0
w
i
v

i
'
+(βl)n,
79
where {v

i
}, l and n are recomputed according to their above definitions.
The above scheme still has a caveat. Although mean value coordinates are one of
the most general barycentric forms, they are only valid if the neighborhood polygon
of v does not have self-intersections once projected on a plane. In our setting, self-
intersections in the projection plane can occur (a) during the simplification step
or (b) during the reconstruction step, after the coarse mesh has been deformed.
The former case is easy: as ”problematic” vertices are detected, we can simply
disallow their decimation. The latter case is more complicated, as the coarse mesh
may be deformed arbitrarily after the hierarchy has been constructed. We treat
this problem heuristically as follows. Although mean value coordinates can be
validly defined for simple non-convex polygons, local concavities are more likely
to degenerate into self-intersection under a deformation (see Figure 3.13). We
thus enforce convex configurations in the construction process: concave vertices are
discarded and assigned zero weights.
Figure 3.13: Self-intersecting polygons (left) do not guarantee valid mean value
coordinates. Concavities in simple polygons (right) are more likely to create self-
intersections under a deformation. Thus, the concave grey vertex is discarded
during construction.
80
3.6 Generalized Differentiable Distance Potential
So far we have merely considered the point-based Hausdorff distance as the potential
function. Now we will show how the distance potential can be generalized to adapt
to various contexts.
3.6.1 Beyond Vertex Positions
Discrete shapes are more than just point sets, they also possess higher order proper-
ties. Some of those have been advocated for use in shape similarity measures: from
differential properties such as surface normals [35,98] and curvature [45,49] to more
complex local features such as shape contexts [15], spin images [63] and geodesic
fans [123]. However, incorporating such information in differential-geometric op-
timization problems may be a non-trivial task. Consider, for instance, the vertex
normal: since it is a deformation-dependent property, its derivative with respect to
perturbations of the related mesh vertices must be computed. While a closed form
discrete approximation of that derivative can be accommodated (in fact, we have
derived one in our experiments), the resulting numerical behavior depends on the
mesh quality and resolution, and may lead to instabilities for arbitrary meshes.
Since the definition of a geometric similarity often depends on the application,
we will confine our further discussion to articulated shape matching. In this con-
text, the Gaussian curvature of a surface becomes an attractive candidate to be
a part of a shape similarity metric. Unlike the mean curvature which reflects the
81
embedding of a surface in the ambient space, the Gaussian curvature is intrin-
sic, i.e., invariant under isometric deformations. Since articulated body motion is
approximately isometric, this assumption allows us to avoid differentiation of the
Gaussian curvatures altogether! In the following, we will show how to incorporate
that Gaussian curvatures into the distance potential function.
Recall that our approximation d
H
(X,Y) of the Hausdorff distance between
two meshes X = {x
i
}
i=1..P
and Y = {y
j
}
j=1..Q
is given by Eq. (3.19), and the
corresponding partial derivatives can be written (Eq. (3.20), with some abuse of
notation) as
∇
x
i
d
H
(X,Y)=
(d
H
(X,Y)+ )
1−2α
P · Q
M
x
ii

j
∇
x
i
d
ij
d
α+1
ij
M
y
jj
(f
−2
i
+ g
−2
j
) .
Now we would like to generalize the definition of d
ij
to
d
ij
=|x
i
− y
j
|
2
[1 + γ(K
G
(x
i
)− K
G
(y
j
))
2
]+
2
,
where K
G
denotes the discrete Gaussian curvature operator [82], and γ is the rel-
ative weight given to the curvature component. Thus, under the assumption of
isometric motion,
∇
x
i
d
ij
=2(x
i
− y
j
)[1 + γ(K
G
(x
i
)− K
G
(y
j
))
2
].
82
In practice, estimated discrete Gaussian curvature can be noisy for an arbitrary
mesh. Therefore, we use Laplacian-smoothed curvature values instead:

K
G
(X) = ((1− β)Id + βΔ)K
G
(X)
where Δ is the discrete Laplace-Beltrami operator [82] and β is a smoothing factor.
This generalization helps disambiguate the basic point-based distance and thus
reduce the amount of local minima: Figure 3.14 illustrates the added power of the
generalized distance, with a rigid registration case is chosen for simplicity. Likewise,
most complex examples of non-rigid matching in this chapter were also generated
using this generalized metric.
3.6.2 Alternatives to the Hausdorff Distance
So far, we have shown that Hausdorff distance minimization can be a powerful tool
for enforcing a complete match. Yet, such a match may not always exist. For
instance, with range scans, it is often the case that only a partial match between
two surfaces is possible. Under such scenarios, the symmetric Hausdorff distance
may not be the natural choice of a potential function. Indeed, this distance can be
large as long as there is at least one unaligned point on either surface, even if that
point is, in fact, an outlier.
One alternative is the (squared) L
2
distance, which has commonly been used in
object registration. Whenever point correspondences are available, the L
2
distance
83
Figure 3.14: Gaussian curvature enhanced distance functional is used to
register two extremely different rigid poses of a full body shape. A strong local
minimum is avoided with the Gaussian curvature incorporated into the Hausdorff
distance functional.
84
is typically the most natural choice with a simple closed form solution. That, how-
ever, is no longer true in absence of correspondences, in which case the symmetric
L
2
distance d
L
2(S,T ) between two surfaces S andT is given by
d
L
2(S,T)=

p∈S
min
q∈T
|p− q| +

q∈T
min
p∈S
|p− q| .
In practice, whether a complete or a partial match is sought, the choice between the
Hausdorff and the L
2
distance depends on the application. Notably, when matching
an incomplete scan X to a complete model Y, an asymmetric analog of either
distance can (and should) be used. In our optimization toolbox, we have included
both functionals for a side-by-side evaluation (the two happen to be quite similar in
terms of implementation—see details below). Overall, most our experiments have
indicated that, when coupled with an adequate deformation prior, the Hausdorff
distance has favorable convergence properties, even the case of partial matching!
This can be explained by the fact that the deformation prior largely eliminates the
notorious problem of small outliers, while the resulting Hausdorff flow is still very
“unforgiving” with respect to large unaligned regions. For the sake of completeness,
a derivation of the L
2
distance minimizing flow is shown below.
85
By analogy to our definition of the approximate Hausdorff distance (Eq. (3.19)),
we can define a differentiable approximation of the symmetric L
2
distance between
two meshes X ={x
i
}
i=1..P
and Y ={y
j
}
j=1..Q
as follows:
(
d
L
2(X,Y)=

i
f
1
2α
i
+

j
g
1
2α
j
− (P + Q)
where, just like before f
i
=
1
Q

j
d
−α
ij
,g
j
=
1
P

i
d
−α
ij
,
and again, > 0and α ≥ 0. Note that we have have omitted the mass matrices
here, assuming a regular sampling for the sake of clarity. Differentiating the above
expression with respect to a vertex position yields:
∂
(
d
L
2(X,Y)
∂x
i
=−

j
x
i
− y
j
d
α+1
ij
(
1
Q
f
i
−
1+2α
2α
+
1
P
g
j
−
1+2α
2α
) .
In terms of the computational complexity, the above expression is very similar to
the Hausdorff distance gradient (Eq. (3.20)), with a slightly higher cost due to the
fractional exponents.
In the particular case of minimizing an (approximated) asymmetric L
2
distance
of an incomplete surface X to a complete model Y, the corresponding partial
derivative with respect to a vertex position becomes
∂
− →
d
L
2(X,Y)
∂x
i
=−
1
Q
f
i
−
1+2α
2α

j
x
i
− y
j
d
α+1
ij
,
86
and the asymmetric Hausdorff distance counterpart can be derived in a similar
fashion. Note that the above expression is now an implicit vector function of Y
in R
3
, which is the target surface that remains static during the minimization
process. In other words, the value of this function for each point x∈R
3
is defined
solely by Y, as long as it remains fixed. Given sufficient memory resources, this
implicit function can be precomputed approximately over a discretized domain (e.g.,
a regular 3D grid), resulting in a major speedup of the minimization process (see
e.g. [86] for a similar idea). This property may become particularly attractive when
the execution speed is a major bottleneck: even in the case of complete matching,
accuracy may be traded for efficiency by approximating the symmetric distance
with the asymmetric one. Keep in mind, however, that for the Hausdorff distance
such an approximation is in general not accurate. Therefore, the asymmetric L
2
distance may often become a practical preference in this particular scenario. Finally,
this strategy may no longer valid in traditional template matching applications,
whenever a complete template is matched to a partial surface and not vice versa.
87
3.7 Applications Beyond Shape Matching
3.7.1 Regularized Gradient Flows
As extensively discussed previously, most existing gradient flows used in geometry
processing are based on L
2
, they all share the consequences induced by this partic-
ular choice of inner product. First, since the L
2
norm of a vector field completely
disregards its spatial coherence, a conventional gradient flow may (and will, see
Figure 3.3) produce highly non-rigid deformation, thus degenerate meshes. Sec-
ond, one can also show that such gradient flows are potentially very susceptible to
noise, and often get stuck in local minima of the energy. To remediate these flaws,
Sundaramoorthi et al. [107] proposed recently to use a regularizing inner product,
namely, a Sobolev norm. They demonstrated the multiple advantages of this choice
in the context of Eulerian (Level Sets) active contours.
Sobolev Inner product Our Lagrangian set-up can accommodate the same
change of norm. For meshes, the Sobolev norm H
1
derives from the following inner
product:
u,v
H
1
=

S
u(x)· v(x)dx + λ

S
∇u(x)·∇v(x)dx .
It is a simple matter of integration by part to show that, for closed surfaces, this
inner product corresponds to the linear operator L(u)=u− λΔu, where Δ is the
88
Laplace-Beltrami operator (discretized for instance in [82,95]), and λ is an arbitrary
positive constant.
Sobolev Gradient Flow With this inner product, we can define the H
1
-gradient
of an energyE, and perform explicit integration of the corresponding gradient flow:
X
t+dt
= X
t
− dt (Id− λΔ)
−1
∇
L
2E(X
t
).
Since ∇
L
2E is linearized in each iteration, a step of Sobolev gradient flow can be
performed by solving the following linear system:
(Id− λΔ)X
t+dt
=(Id− λΔ− dt∇
L
2E)X
t
. (3.21)
The solution of this linear system now couples the motion of each vertex to the mo-
tion of the other vertices. This exemplifies the regularization effect: for an explicit
integration of a L
2
-gradient, vertices move basically independently of each other,
while for H
1
, vertices move in concord. Finally, we can extend the regularization
scheme to higher Sobolev-like norms. For instance, a higher-order inner product
can be defined through L(u)= u− λΔu + μΔ
2
u, leading to an additional term
on each side of the linear system to solve. This regularization is very reminiscent
of a well-known smoothing method: it turns out that the explicit integration of
the H
1
-gradient flow of the surface area of a mesh is known as implicit fairing [38].
89
Indeed, we mentioned in Section 3.2.2 that the L
2
-gradient of the area corresponds
to the pointwise Laplace-Beltrami operator. Substituting this term in Eq. (3.21)
and choosing λ = dt, we get the usual expression:
(Id− dtΔ) X
t+dt
= X
t
.
Seeing implicit fairing as simply being a regularized explicit mean curvature flow
is another way of understanding its notable stability. The reader will have no-
ticed that, the generalized mean curvature flow defined by the higher-order inner
product corresponds once again to an implicit fairing, this time with higher-order
operators—a known way to improve the filtering quality of this smoothing flow [38].
A similar regularizing effect can be used for any flow using Eq. (3.21) with the ad-
ditional Δ
2
=Δ◦Δterms.
Discretization Remark Here again, care must be taken to ensure that the dis-
crete case is consistent with its continuous counterpart. This implies the use of
the discrete pointwise analog of the Laplace-Beltrami operator (denoted here by
K, following the notation of Meyer et al. [82]), as opposed to non-normalized dis-
crete surface Laplacian L
d
[95] that acts on area-weighted mesh vertices. In this
case, K(X) directly corresponds to the discrete mean curvature flow (Eq. (3.7)),
and the relation to implicit fairing carries over exactly in the discrete case. Finally,
to dissipate potential confusion, we will note that the two different discretization
90
of the Laplace-Beltrami operator are closely related by K = M
−1
L
d
,where M is
the lumped mass matrix. For a detailed discussion on pointwise vs. integrated
Laplacians, the reader is referred to [116].
3.7.2 Volume-Controlled Implicit Fairing
Now we take the above regularized scheme one step further by modifying it to
allow a direct control over the local volume change during the flow. Once again, a
simple modification of the inner product is all it takes. To better understand and
appreciate the problem, let us first briefly review the related work.
Background of Volume Preserving Surface Fairing Surface fairing is a cen-
tral geometry processing tool, routinely used for denoising and smoothing applications—
see [20], Chap. 7 for an overview of the extensive literature on the subject. A noto-
rious shortcoming of basic fairing schemes is non-uniform volume shrinkage. This
calls for non-uniform local volume compensation treatment which, unlike global vol-
ume control [38], becomes a non-trivial problem. Existing remedies often suffer from
sensitivity to choice of parameters [109], are subject to unnatural global effects [57],
or resort to higher-order, computationally intensive surface evolutions [19]. Alter-
natively, constrained optimization methods involving either local volume [77] or
distance [58] constraints cannot always copewithoutliers. Instead, hereweshow
a simple, efficient and unconditionally stable smoothing scheme that implements a
locally volume-controlled flow, as described below.
91
Quasi-Volume Preserving Flows Let VP denote the subspace of all volume-
preserving deformations, with NVP being its L
2
-orthogonal complement, and let
Π
VP
and Π
NVP
be the projection operators on the corresponding subspaces. We
define a volume-controlled norm of a field v simply as:

v

2
L
VP
=
Π

2
VP
v +
1
λ

Π
NVP
v

2
,
where λ is a penalty factor on volume-preserving motion: a smaller λ will make
such motion ”cheaper”, thus favored. Following the procedure from Eq.(3.8), we
get:
L
VP
=Π
VP
+
1
λ
Π
NVP
, thus:L
VP
−1
= λId + (1− λ)Π
VP
,
which can be used to obtain a new gradient flow, as in Eq.(3.9). For instance, a
volume-controlled implicit fairing is achieved through solving:
(L
VP
− λΔ)X
t+dt
=L
VP
X
t
. (3.22)
Finally, the projection operator Π
VP
on the subspace of locally volume-preserving
motions can be implemented as follows. A deformation field v of a discrete mesh
X = {x
i
} with a normal field N = {n
i
} and a lumped mass matrix M,causesa
local volume change δV
i
= v
i
· M
ii
n
i
around each vertex x
i
,where M is assumed
to be diagonal, i.e., the total area is A =

i
M
ii
.(M
ii
n
i
is akin to the usual
92
area-weighted normal at a vertex.) Then averaging δV
i
over the one-ring area of x
i
yields:
δV
i
=
δV
i

j∈N
+
(i)
M
jj
,
where N
+
(i) denotes the vertex i and its immediate neighbors. Thus, the local
volume-changing velocity component of each x
i
becomes:
¯ v
i
= n
i

j∈N
+
(i)
δV
j
,
yielding: Π
VP
(v
i
)= v
i
− ¯ v
i
. In the matrix form, theL
VP
−1
is given by
L
VP
ij
−1
=
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1+(λ− 1)m
i
n
i
n
i
T
j = i
(λ− 1)m
j
n
i
n
j
T
j ∈ N
+
(i)
0otherwise
,
where m
i
=
M
ii

j∈N
+
(i)
M
jj
. With this procedure, λ is a“knob”that makes the flow
more or less volume-preserving—it will become fully volume-preserving as λ goes
to zero. Note that controlling the volume locally also ensures global volume control.
Also, note that the extra linear operator is zero for all pairs of vertices that are
not immediate neighbors in the mesh, and therefore does not affect the sparsity of
the original linear system in implicit fairing. Fig. 3.15 illustrates the obvious non-
shrinking effect of this generalized flow compared to the original implicit fairing for
a same amount of smoothing. Lastly, let us remark that the above deformation
93
prior is not limited to fairing: it could also be applied in other contexts, such
as registration of deformable anatomical organs, which often emphasizes volume
preservation.
3.8 Conclusion and Future Work
We have extended the notion of gradient flows of triangle meshes by pointing out
that one can tailor the inner product on vector fields instead of using the canonical
L
2
. We gave details on implementation of such generalized geometric flows using
explicit and implicit integration techniques. We applied this formalism to improve
conventional approaches to shape matching, by first defining an approximate Haus-
dorff gradient flow, then by applying deformation priors (such as articulated bodies
or quasi-isometry) to render the flow significantly more robust (i.e., less prone to
local minima and erroneous matching). Finally, we use a multi-resolution approach
to speed up convergence and add scalability. Our future work includes extending
the toolbox of energies and priors for the particular case of shape matching. For
instance, landmarks are often used in practice to alleviate the registration task
(e.g. in [50]). Incorporating landmarks into our deformation energy is simple
and could combine the strengths of both approaches. We also plan on further
investigating the power of altering the norms in mesh flows in the more general
context of shape representation and analysis. Promising results have been shown
94
in this area [69,83,121], where properly defining continuous morphs (homotopies)
and geodesic distances in the space of curves requires altering the L
2
metric.
Acknowledgements We are grateful to Xiaohan Shi for sharing the deformed Armadillo
model and to Michael and Alex Bronstein for testing their GMDS algorithm on it. This work
waspartiallyfunded byNSF (CAREERCCR-0133983andITR DMS-0453145),DOE(DE-FG02-
04ER25657), and Caltech’s Center for the Mathematics of Information.
95
Figure 3.15: Implicit fairing with local volume control. Top row: a noisy horse
model (left) becomes degenerate after aggressive implicit fairing (middle), while a
comparable amount of smoothing with local volume control avoid shrinking (right);
closeup shows overlaid comparison. Bottom: several iterations of implicit fairing
of a noisy feline model (far left), with local volume control (above) and without
it (below), shown for results with a comparable degree of smoothness. With the
original scheme, a strong smoothing (far right) causes an extreme shrinkage and
degeneracies, whereas with volume control (λ=10
−6
), the original proportions and
overall shape are retained. The inset illustrates the non-shrinking effect for the tail
region, with the corresponding surfaces overlaid.
96
Chapter 4
Matching Cortical Surfaces with Surface Flows
The human brain, then, is the most complicated
organization of matter that we know.
Isaac Asimov
Despite being routinely required in medical applications, deformable surface
registration is notoriously difficult due to large intersubject variability and complex
geometry of most medical datasets. We present a general and flexible deformable
matching framework based on generalized surface flows that efficiently tackles these
issues through tailored deformation priors and multiresolution computations. The
value of our approach over existing methods is demonstrated for automatic and
user-guided registration of cortical surfaces.
97
4.1 Introduction
Matching (or registration) of deformable surfaces is a fundamental problem in med-
ical image analysis and computational anatomy. One particularly challenging in-
stance of the problem arises in the field of human brain mapping, where deformable
registration of two cortical surfaces (e.g., segmented from MRI volumetric data) is
required for intersubject comparisons and intrasubject analysis of neuroanatomi-
cal surface data. In these applications, cortical matching becomes an important
tool for understanding brain anatomy and function. For instance, anatomical la-
beling of brain structures from 3D imaging of the patient’s brain is one of the
key challenges in brain mapping [113]. Typically, this task is tackled by matching
a previously labeled canonical atlas of the brain to the subject anatomy. Other
related studies include progression of disorders such as Alzheimer’s disease, brain
growth patterns, genetic influences [110] and the effects of drug abuse on the struc-
ture and function of the brain [88]. The challenge in registering two cortices lies
in the wide inter-subject variability and the convoluted geometry of the cortical
surface, thus becoming a real ”stress test” for any general deformable registration
technique, as shown in Figure 4.1. To overcome this challenge, various landmark-
based and landmark-free methods have been developed, e.g. [2,12,46,60,112,114],
as elaborated in the background section below.
In this chapter we show how our framework of generalized surface flows intro-
duced in Chapter 3 can be leveraged for the case of cortical matching. As in other
98
Figure 4.1: White matter cortical surfaces of two different subjects. The complex-
ity and anatomical variability of folding patterns make identification of correspon-
dences challenging even for the human eye.
shape matching cases discussed earlier, the proposed method involves an iterative
deformation of the template surface to match the target, until convergence criteria
are met. As a result, the deformation trajectory is available as a by-product for
evaluation and determination of the best fit. This property is particularly impor-
tant in computational anatomy applications [14, 52, 64, 84, 122], where finding an
optimal deformation between two related shapes becomes the central issue.
The use of geometric flows in medical image analysis is not new. For instance,
active contours (or snakes) and deformable models have been widely applied to
reconstruct surfaces from volumetric images (e.g., [119]). The main drawback of
these methods is their sensitivity to local minima, which can become particularly
severe when matching of geometrically complex objects is sought. Here again, we
take advantage of the generalized flows formalism to systematically deal with this
99
issue using problem-specific prior knowledge. Specifically, the contributions of this
chapter are as follows:
• We reduce the cortical surface registration problem to a deformable alignment
of the corresponding partially flattened surfaces, generated rapidly via an
implicit fairing procedure.
• We show how to achieve an automatic alignment of such partially flattened
surfaces using a generalized pseudo-Hausdorff flow equipped with a Sobolev-
like deformation prior.
• We demonstrate how to incorporate curve correspondence constraints into the
surface matching framework, to allow optional user control over the registra-
tion process and achieve a finer alignment of the cortical folds.
• Our approach applies multiscale strategies to reduce both geometric and
computational complexities: geometrically—using implicit fairing of the in-
put surfaces to simplify the matching task, and computationally—employing
coarse-to-fine mesh representations to optimize performance.
A qualitative comparison of the method to a state-of-the-art in brain matching
technique shows a comparable quality of sulcal alignment, but at much faster (near-
interactive) rates.
100
4.2 Backround
Brain matching and mapping has been an active research area for the last twenty
years. Here we will only sketch a short outline of existing approaches, with focus on
the recent work. For a comprehensive, although somewhat dated, survey of the sub-
ject, the reader is referred to Thompson and Toga [91]. Among existing methods,
two major categories can be identified. In the first category, parameterization-based
techniques start by finding mappings between each cortical surface and shared pa-
rameter domain, such as a plane or a sphere. These maps are then treated as
2D images in the parameter domain and aligned using geometric features such as
normals and curvature [42,46,54,88,111,112] or sulcal landmarks [2,91]. The of-
ten large change in metric due to the mapping needs to be accounted for while
performing the alignment process in the parameter domain [76, 91], e.g., by em-
ploying covariant derivatives [2], adding to the computational costs and numerical
instabilities. Wang et al. [115] avoid the use of an intermediate parameter domain
by computing a direct cross-parametrization of the two surfaces, requiring a large
number of manually placed correspondences whose identification is often laborious
and error-prone. Another class of techniques operates directly in the ambient space
by finding a 3D warping field that aligns the cortical features. Most of these meth-
ods are volume-based, aiming to align image features rather than surfaces, with
similarity measures based on voxel intensities [11, 117, 118], invariant geometric
101
moments [36] or mutual information [97]. While these methods are useful for cross-
modality (MRI, CT, PET, etc.) registration as well as for coarse-level matching,
their inherent inability to preserve structural information presents a disadvantage
for cortical matching. To overcome this limitation, Liu et al.proposed a hybrid
approach [108], enhancing the initial volumetric warp [36] with a deformable cortex
alignment step. Still, important details of the latter remain not entirely clear, and
the overall computation time amounts to several hours. In the following, we present
a purely surface-based method that can align cortices up to an arbitrary precision
in a few minutes.
4.3 Method
Casting cortical registration as a deformable template matching problem allows us
to take advantage of the computational framework introduced in Chapter 3 by it-
eratively minimizing the Hausdorff distance between the template and the target
cortical surfaces. Recall that a na¨ ıve minimization of the Hausdorff distance for
two dissimilar shapes is unlikely to yield relevant correspondences, as the energy
landscape is too complex and non-linear to avoid getting stuck in one of the nu-
merous local minima (this holds even for low-complexity surfaces, as illustrated by
Figure 4.2, (b-d)). Thus, to apply the method effectively, we still need to supply
an adequate deformation prior.
102
Figure 4.2: A white matter cortex and its corresponding Partially Flattened Surface
(a), where the major sulci are still clearly present. A basic Hausdorff gradient flow
applied on a rotated, smoother version (b) creates spurious deformations (c) due
to lack of coherence and local minima, while with a quasi-rigid deformation prior
(d), it successfully recovers the transformation.
4.3.1 Choice of a Deformation Prior
One fundamental difficulty arising in intersubject anatomical matching stems from
the fact that a ground-truth deformation model is rarely available, as this process
does not mimic any natural phenomenon. Therefore in practice, smoothness of
the template-to-subject mapping typically becomes the default requirement. To
enforce such a smoothness constraint, we can take advantage of the regularizing
H
1
prior introduced in Section 3.7.1. Recall that this prior is defined as L(u)=
u− λΔu, where Δ is the Laplace-Beltrami operator, λ is the desired smoothness
factor, and the corresponding surface flow is given in its general form by Eq. (3.21).
Thus, the Hausdorff H
1
flow can be obtained simply by substituting the Hausdorff
distance d
H
(X,Y) between the template mesh X and the target Y as the objective
functional:
(Id− λΔ)X
t+dt
=(Id− λΔ− dt
∂d
H
∂X
)X
t
.
103
The only remaining step is to replace the continuous operator L with its discrete
analog, given by ML
H
1 = M− λL
d
, where L
d
is the discrete non-normalized
surface Laplacian, as discussed in section 3.7.1. Thus, each iteration of the discrete
Hausdorff H
1
flow can be computed by solving the following linear system:
(M− λL
d
)X
t+dt
=(M− λL
d
− dt
∂d
H
∂X
)X
t
, (4.1)
or, with implicit time integration:
(M− λL
d
+ dt
∂d
H
∂X
)X
t+dt
=(M− λL
d
)X
t
. (4.2)
For a stronger regularization effect, one can extend to above scheme to higher
order Sobolev-type norms. For instance, we can define a higher-order prior L(u)=
u + μΔ
2
u,where Δ
2
=Δ◦ Δ. With a somewhat higher computational cost, the
resulting scheme is tantamount to regularizing the instantaneous deformation with
a thin-plate spline energy term (see e.g. [34]). In practice, we stick to the sparser
H
1
prior in this work.
As-Rigid-As-Possible Smooth Deformations Since the two input shapes are
generally given in separate coordinate frames, it is typically desired to first bring
them into a rigid alignment. For that purpose we can reuse the quasi-rigid deforma-
tion prior L
R
from Section 3.4.3.1. Furthermore, it is easy to combine both priors
104
by first recovering and boosting the rigid component of a given motion field factor,
followed by a regularization of the non-rigid residual. The proper combination of
the two can be derived through the following square norm of a deformation field u:

u

2
R,H
1
= λ
R

Π
R
u

2
L
2
+ λ

S
∇˜ u(x)·∇˜ u(x)dx ,
where Π
R
projects the given deformation field to the rigid motion subspace, λ
R
is
the corresponding cost factor, and ˜ u(x)=(Id−Π
R
)u(x). Thus, the corresponding
combined deformation prior is given in its discrete form by:
ML
R,H
1 = λ
R
MΠ
R
− λ(Id− Π
R
)
T
L
d
(Id− Π
R
) .
Recall that, even though Π
R
is a full matrix, the actual procedure has a linear time
complexity since several mechanical quantities of the shape and the associated
deformation field (e.g., surface area, mean translation, total angular momentum,
etc.) can be precomputed. However, it is less clear at first glance clear how to apply
Π
T
R
in a similar fashion. And yet, since M and M(Id− Π
R
) are both symmetric,
we get: (Id− Π
R
)
T
= M(Id− Π
R
)M
−1
, which is again a linear-time operator. To
sum up, by applying this combined prior we obtain an as-rigid-as-possible smooth
surface flow that can adequately and seamlessly resolve both the rigid and the
non-rigid phases of the alignment process in a single procedure.
105
Figure 4.3: Automatically matching a template (grey) to the subject cortex (blue).
Partially flattened representations of both surfaces are iteratively aligned using a
Hausdorff flow with a smoothing prior. The obtained alignment yields a correspon-
dence between the original surfaces. The final color mix is due the fact that the
surfaces lie on each other.
4.3.2 Matching Cortical Surfaces
Basic Algorithm We are now ready to apply the Hausdorff flow approach to
match a template cortical surface (e.g., a digital atlas of the cortex) to an instance
surface, e.g., segmented from a MRI scan. One na¨ ıve solution would be to perform
the minimization directly on the input surfaces, combined with the H
1
deformation
prior for regularization. This process is still likely to get stuck in a local minimum
due to the highly convoluted geometry of the cortex. However, even if we managed
to get the two surfaces into a complete alignment, the result would hardly be ade-
quate, as intercortex correspondence is in general not well-defined due to extreme
anatomical variability of the cortical structures. In practice, quality of match is
measured by the alignment of the major sulcal patterns that can be consistently
106
identified in all normal brains. Thus, minimum intersurface distance alone is not a
sufficient condition for an acceptable solution. In view if this problem, the use of
Partially Inflated Surfaces (PFS) has been advocated for cortical matching [1,112].
The idea is to filter out excessive surface detail through, e.g., Mean Curvature
Smoothing [38]), applying a limited amount of smoothing in order to facilitate
matching while preserving the principal sulcal patterns—see Figure 4.2(a) for an
illustration. We adopt this approach, with one important difference. While corre-
spondence between two PFSs is typically computed by matching their maps in a
common parameter domain, we eliminate the need for these intermediate mappings
by aligning the PFSs directly. Our strategy is summarized below, and illustrated
in Figure 4.3:
ALGORITHM:
1. Partially flattenS andT , obtaining S

andT

, respectively.
2. Apply Generalized Hausdorff Flow to achieve an arbitrarily close alignment
ofS

withT

, yielding a correspondence map ϕ between the two.
3. ReturnS→S
ϕ
− →T

→T as a bijective map between S andT .
The first step can be done rapidly using Mean Curvature Smoothing (MSC), with
implicit time integration allowing an arbitrarily large time step. Note that MSC
is a classical example of gradient flow, so our whole approach fits nicely into the
flow-based methodology. The crux of the algorithm lies in the second step, where
107
thetemplatePFS S

(which can be precomputed for repeated use) is iteratively
deformed to match T

. To regularize the flow, we use the L
R,H
1 operator from
Section 4.3.1. In practice, once the rigid component of the motion vanishes, L
R,H
1
can be replaced with a simpler H
1
prior for efficiency. As the surfaces get closer,
we switch to implicit time integration (Eq. 4.2) to avoid oscillations and accelerate
convergence. Finally, the process is further optimized for high-resolution models
through a multiresolution strategy, as described in Section 3.5.2, yielding a speedup
of several orders of magnitude.
Adding Constraints As shown in Figure 4.4(a), the above procedure manages to
automatically align most sulci, but cannot guarantee a fully correct match for all of
them when a strong sulcal variability is present. A common remedy is to incorporate
constraints, i.e., expert-specified sulcal curves, to control the mapping. In our case,
adding constraints to the pseudo-Hausdorff energy is quite straightforward. Indeed,
matching of two curves on opposite surfaces is just another distance minimization
problem—this time, between sets of surface points that lie on the two curves. Thus,
we can reuse the same Hausdorff distance approach, assigning a separate but similar
weighted energy term (one for each pair of corresponding curves) to those mesh
vertices that are incident on the curves—instead of the global Hausdorff potential.
Adding point constraints, if needed, is even simpler. Note also that the constrained
deformation is still kept smooth due to the use of the H
1
prior.
108
Figure 4.4: a. Automatic matching of PFSs yields a close alignment for most sul-
cal curves. b. Constraining only 7 out of the 23 available curves reduces most
misalignments, further improved by using the full set of constraints (c). d. Cor-
responding sulcal alignment for the original cortical surfaces. For clarity, a single
cerebral hemisphere is shown.
4.4 Results
The proposed cortical matching algorithm was tested with a dataset of six subject
brains, segmented from MRI scans using the BrainSuite tool [102], each supple-
mented with a set of sulcal curves marked by an expert according to the LONI
Sulcal Tracing Protocol [110]. As illustrated by Figure 4.3, the algorithm automat-
ically computes a near zero-distance alignment for two partially flattened cortical
109
surfaces, effectively yielding an intercortex correspondence. It results in a reason-
ably close alignment for most sulcal curves, further improved through the addition
of constraints. Figure 4.4 shows that most sulci could be matched automatically,
and constraining only a subset of the sulcal curves is sufficient, thus significantly
reducing the amount of manual effort required.
Table 4.1 summarizes a limited evaluation of our algorithm (GHF) based on six
pairs of subject brain images, compared to HAMMER [36] which is a state-of-the-
art, publicly available brain matching tool. Although the two methods operate on
different modalities, distances between corresponding subject and deformed tem-
plate sulcal curves can be measured in both cases. Even without resorting to con-
straints (to make a fair comparison to the landmark-free HAMMER), our method
demonstrates a comparable quality of match, with clearly superior computation
times: under 5 min on a standard PC, as opposed to several hours. Note also that
for PFS sulci, registration error is even lower, which illustrates the quality of the
core deformable matching procedure.
Method / Data MeanL
2
Distance Per Case (mm) Total Average
HAMMER / Original Sulci 4.67 4.62 4.79 5.05 5.13 4.90 4.87
GHF / Original Sulci 5.49 5.02 4.56 5.16 4.97 4.53 4.96
GHF / PFS Sulci 4.13 3.87 3.54 4.12 4.07 3.32 3.84
Table 4.1: Quality of match between deformed template and subject brains as
average L
2
distances between corresponding sulcal curves.
110
4.5 Discussion and Future Work
We have presented a practical and flexible multiresolution framework for deformable
registration of anatomical surface data, based on generalized geometric flows. In
the case of cortical matching, evaluations indicate quality comparable to state of
the art methods, with near-interactive computation times.
Still, the presented solution is not without limitations. For instance, self-
intersections may occur during the deformation. The likelihood of such self-intersections
and “pinching” effects may increase in presence of curve constraints, if their cor-
responding weighting factor is not chosen properly (in fact, one can easily design
constraints that do not allow any intersection-free solution). Even though the choice
of an adequate weighting factor did not present much of a challenge in our experi-
ments, this issue is still a limitation in the general case. In the future, we wish to
investigate seamless ways to combine surface and curve distance terms in a single,
parameter-free energy functional. Another potential strategy to address this short-
coming is through a dedicated deformation prior assigned to the energy constraint
term. For instance, a prior that boosts tangential motion of constrained vertices
will cause the sulcal curves to “slide” along the surface rather then penetrating it.
We are also exploring ways to generalize the definition of a geometric distance and
design new priors to improve automatic matching of sulcal features. The central ob-
servation is that, in its current realization, the automatic matching procedure may
unfold an existing sulcus (or gyrus), even when a proper sulcal alignment can be
111
obtained merely by a local surface stretching. Therefore, it is a promising strategy
to model the instantaneous non-rigid deformation as a combination of stretching
and bending terms (as it is often done e.g. in thin shell simulations in computer
graphics [53,116]), downplaying the bending component in the deformation prior.
So far, some progress has been done in this direction, and further experiments are
currently underway.
Yet another potential improvement is related to the generation of partially flat-
tened representation of the cortex, which is currently done via the standard Mean
Curvature Smoothing [38]. Table 4.1 shows that the quality of sulcal alignment is
slightly degraded at the final level, compared to the PFS level. This observation
suggests that with a surface smoothing procedure that offers better local shape
preservation, the final alignment could be improved. Furthermore, we already have
one promising candidate for such a procedure, in the volume-controlled surface fair-
ing method, described in Section 3.7.2. Indeed, Figure 4.5 shows that smoothing
the cortex with near volume preservation results in a surface much closer to to the
original in terms of the Hausdorff distance. We are currently investigating the effect
of such a modification on our cortical matching scheme.
Finally, using generalized flows to compute continuous morphs that follow geodesics
in shape spaces [14,69,84] is another exciting avenue of future work, with important
applications in computational anatomy.
112
Figure 4.5: A cortical surface (far left) is gradually smoothed with near volume
preservation (VP)—above, and without it—below. Far right: color-coded Hausdorff
distances (blue=smallest, red=largest) between the original and the final results,
as measured by the EPFL Mesh tool [17]. The result of VP fairing is much closer
to the original for a similar smoothing effect.
113
Chapter 5
Conclusions and Future Directions
In this research we focused on several analytical problems in discrete geometry
processing, with all the proposed solutions relying on the notion of shape cor-
respondence. Starting with dynamic geometry compression, where a geometric
optimization approach was leveraged for local motion estimation, through gener-
alized surface flows—first explored in a wider context, then harnessed in a flex-
ible and powerful shape matching framework, to applications in computational
neuroanatomy—the unifying philosophy in all these studies amounts to approach-
ing the related surface correspondence problems from a geometric optimization
perspective. For many of these problems, research presented in this text merely
marks the first step in that direction, paving the way to numerous directions of
future work. Let us briefly mentioned some of them.
• Template-free matching and automatic structure inference of articulated shapes
114
One interesting problem is extending the matching method for near-articulated
shapes to a completely template-free algorithm. Here, it is a promising di-
rection to reuse the quasi-isometric prior from Section 3.4.3.3. As such, this
prior strives to preserve isometry for local regions on the surface, which is not
enough to ensure articulation: local isometry does not prevent bending and
strong variations in distance between points that are located far apart on the
surface, but nearby in the Euclidean domain. However, one could extend the
notion of quasi-isometry to a truly three-dimensional prior by tetrahedralizing
the shape interior and enforcing length preservation for the tetrahedra edges.
Alternatively, combining local volume preservation (Section 3.7.2) with the
quasi-isometric prior could provide a tradeoff between generality and sim-
plicity. Also, one particularly interesting application context is 3D shape
tracking, which brings us back to the case of time-varying surfaces, discussed
in Chapter 2. In this context, a denser temporal sampling leads to a limited
amount of deformation between consecutive shape instances, thereby allow-
ing for a simpler deformation model. Lastly, automatic segmentation of the
tracked object into articulated parts is another motivating problem.
• Shape warping that follows geodesics in shape spaces
Gradient descent methods are greedy: their generated minimization paths to
the global minimum (if found) are not guaranteed to be the shortest. Conse-
quently, shape warping done via generalized gradient flows may not follow a
115
globally optimal trajectory (with the definition of optimality usually depend-
ing on the chosen metric of vector fields), which is an important requirement
in several applications. In a recent work, Kilian et al. [69] proposed an al-
gorithm to generate such optimal surface evolutions in the context of shape
morphing. However, the technique is limited to pairs of input meshes with
given point correspondences. Computational anatomy [52] is another area
where global optimality of the warp is of major significance. There, it is a
common practice to express the global warp φ as a minimizer of the following
”energy” functional:

1
0

v
t

2
L
dt + U(φ(X
0
),X
1
) , (5.1)
where the second term is the potential energy of the warped template shape
X
0
(e.g., its Hausdorff distance to the target) and the first term can be
thought of as path integral of the kinetic energy, where the metric is de-
fined with respect to L. Finally, φ is obtained through a time integration
of the instantaneous velocity fields {v
t
},from t =0 to t =1. Note that
this formulation bears a similarity to generalized gradient flow approach. Yet
the important distinction is that, unlike with gradient flows, here the global
minimum corresponds to a trajectory that is globally optimal, and therefore
needs to be solved via a global minimization whose complexity often becomes
prohibitive for real-life cases. The natural question is, therefore, whether we
116
could still use a greedy algorithm to find the (near-)optimal solution. Inter-
estingly, this problem bears a similarity to mechanics, where body dynamics
are often described through the least action principle. This principle, stating
that a dynamical system always finds the optimal course from one position
to another, can be enforced in discrete Lagrangian mechanics through the
use of variational time integration [67,79,105] which can be realized through
a greedy algorithm! This suggests a clear payoff in looking at the optimal
warping problem from the Lagrangian dynamics point of view. While some
of these ideas require further exploration, they lead to an exciting avenue of
future work.
• Extension to deformable surface-based brain atlases
Finally, research presented in Chapter 4 focuses specifically on cortical match-
ing, and could be used to fit a deformable cortical atlas to the subject cortex.
It is therefore a natural generalization to extend the approach to the interior
of the brain. Specifically, given a generic brain atlas represented as a collection
of manifold surfaces corresponding to the anatomical brain structures [68], a
generalized method for fitting such an atlas to patient-specific data will have
important applications in brain mapping and neurosurgical planning.
117
Bibliography
[1] Cachia A., Mangin J., Riviere D., Kherif F., Boddaert N., Andrade A.,
Papadopoulos-Orfanos D., Poline J., Bloch I., Zilbovicius M., Sonigo P.,
Brunelle F., and Regis J. A primal sketch of the cortex mean curvature:
a morphogenesis based approach to study the variability of the folding pat-
terns. IEEE Transactions on Medical Imaging, 22(6):754–765, 2003.
[2] Joshi A.A, Shattuck D.W., Thompson P.M., and Leahy R.M. A framework
for registration, statistical characterization and classification of cortically con-
strained functional imaging data. IPMI, 2005.
[3] M. Alexa. Recent Advances in Mesh Morphing. Computer Graphics Forum,
21(2):173–196, 2002.
[4] Marc Alexa and Wolfgang M¨ uller. Representing animations by principal com-
ponents. Computer Graphics Forum, 19(3), 2000.
[5] B. Allen, B. Curless, and Z. Popovi´ c. The space of human body shapes:
reconstruction and parameterization from range scans. ACM SIGGRAPH,
pages 587–594, 2003.
[6] Brett Allen, Brian Curless, and Zoran Popovi´ c. Articulated body deformation
from range scan data. ACM Trans. on Graphics (SIGGRAPH), 21(3):612–
619, 2002.
[7] Brett Allen, Brian Curless, and Zoran Popovi´ c. Articulated body deforma-
tion from range scan data. ACM Transactions on Graphics (SIGGRAPH),
21(3):612–619, 2002.
[8] D. Anguelov, P. Srinivasan, H.-C. Pang, D. Koller, S. Thrun, and J. Davis.
The correlated correspondence algorithm for unsupervised registration of non-
rigid surfaces. In Advances in Neural Information Processing Systems, 2004.
[9] Dragomir Anguelov, Praveen Srinivasan, Daphne Koller, Sebastian Thrun,
Jim Rodgers, and James Davis. Scape: shape completion and animation
of people. ACM Trans. on Graphics (SIGGRAPH), 24(3):408–416, August
2005.
118
[10] Nizam Anuar and Igor Guskov. Extracting animated meshes with adaptive
motion estimation. In Vision, Modeling, and Visualization, pages 63–71,
2004.
[11] J. Ashburner, P. Neelin, DL Collins, A. Evans, and KJ Friston. Incorporating
prior knowledge into image registration. Neuroimage, 6(4):344–352, 1997.
[12] M. Bakircioglu, U. Grenander, N. Khaneja, and M. I. Miller. Curve matching
on brain surfaces using frenet distances. Human Brain Mapping, 6:329–333,
1998.
[13] R. Barret, M. Berry, T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Ei-
jkhout, R. Pozo, C. Romine, and H. van der Vonst. Templates for the Solution
of Linear Systems: Building Blocks for Iterative Methods. SIAM, 1994.
[14] M.F. Beg, M.I. Miller, A. Trouv´ e, and L. Younes. Computing Large De-
formation Metric Mappings via Geodesic Flows of Diffeomorphisms. IJCV,
61(2):139–157, 2005.
[15] S. Belongie, J. Malik, and J. Puzicha. Shape matching and object recogni-
tion using shape contexts. Pattern Analysis and Machine Intelligence, IEEE
Transactions on, 24(4):509–522, 2002.
[16] P.J. Besl and N.D. McKay. A method for registration of 3D shapes. IEEE
Transactions on PAMI, 14(2):239–256, 1992.
[17] MESH : Measuring Error between Surfaces using the Hausdorff distance.
http://mesh.epfl.ch.
[18] S. Bischoff, T. Weyand, and L. Kobbelt. Snakes on triangle meshes. Bildver-
arbeitung fr die Medizin, pages 208–212, 2005.
[19] Alexander I. Bobenko and Peter Schr¨ oder. Discrete willmore flow. In Sym-
posium on Geometry Processing, pages 101–110, July 2005.
[20] M. Botsch, M. Pauly, C. Rossl, S. Bischoff, and L. Kobbelt. Geometric model-
ing based on triangle meshes. International Conference on Computer Graph-
ics and Interactive Techniques, 2006.
[21] Hector M. Briceo, Pedro V. Sander, Leonard McMillan, Steven Gortler, and
Hugues Hoppe. Geometry videos: a new representation for 3d animations.
In SCA ’03: Proceedings of the 2003 ACM SIGGRAPH/Eurographics Sym-
posium on Computer animation, pages 136–146, 2003.
[22] A.M. Bronstein, M.M. Bronstein, and R. Kimmel. Generalized multidimen-
sional scaling: A framework for isometry-invariant partial surface matching.
Proceedings of the National Academy of Sciences, 103(5):1168–1172, 2006.
119
[23] A.M. Bronstein, M.M. Bronstein, and R. Kimmel. Calculus of non-rigid
surfaces for geometry and texture manipulation. IEEE Transactions on Vi-
sualization and Computer Graphics, to appear, 2007.
[24] A.M. Bronstein, M.M. Bronstein, and R. Kimmel. Joint intrinsic and extrin-
sic similarity for recognition of non-rigid shapes. Technical Report CIS-2007-
01, Dept. of CS, Technion, 2007.
[25] M.M. Bronstein, A.M. Bronstein, R. Kimmel, and I. Yavneh. A multigrid ap-
proach for multi-dimensional scaling. Proceedings of Copper Mountain Con-
ference on Multigrid Methods, 2005.
[26] BJ Brown and S. Rusinkiewicz. Non-rigid range-scan alignment using thin-
plate splines. Int. Symp. on 3D Data Processing, Visualization and Trans-
mission, pages 759–765, 2004.
[27] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours. The Inter-
national Journal of Computer Vision, 22(1):61–79, 1997.
[28] Junavit Chalidabhongse and C.-C. Jay Kuo. Fast motion vector estimation
using multiresolution-spatio-temporal correlations. IEEE Trans. on Circuits
and Systems for Video Technology, 7(3):477–488, 1997.
[29] G. Charpiat, O. Faugeras, and R. Keriven. Approximations of shape metrics
and application to shape warping and empirical shape statistics. Foundations
of Computational Mathematics, 5(1):1–58, 2005.
[30] G. Charpiat, R. Keriven, J.-P. Pons, and O. Faugeras. Designing spatially
coherent minimizing flows for variational problems based on active contours.
In International Conference on Computer Vision, volume 2, pages 1403–1408,
2005.
[31] G. Charpiat, P. Maurel, J.-P. Pons, R. Keriven, and O. Faugeras. Generalized
gradients: Priors on minimization flows. Technical Report 06-23, ENPC-
CERTIS, march 2006.
[32] Y. Chen and G. Medioni. Object modeling by registration of multiple range
images. In IEEE Conference on Robotics and Automation, volume 3, pages
2724–2729, 1991.
[33] Y. Chen and G. Medioni. Object modelling by registration of multiple range
images. Image Vision Comput., 10(3):145–155, 1992.
[34] H. Chui and A. Rangarajan. A new point matching algorithm for non-rigid
registration. Computer Vision and Image Understanding, 89(2-3):114–141,
2003.
120
[35] D. Cohen-Steiner, P. Alliez, and M. Desbrun. Variational shape approxima-
tion. ACM Transactions on Graphics (TOG), 23(3):905–914, 2004.
[36] Shen D. and Davatzikos C. Hammer: Hierarchical attribute matching mech-
anism for elastic registration. IEEE Transactions on Medical Imaging, 21(8),
2002.
[37] Gilles Debunne, Mathieu Desbrun, Marie-Paule Cani, and Alan H. Barr.
Adaptive simulation of soft bodies in real-time. In Computer Animation,
pages 133–144, May 2000.
[38] M. Desbrun, M. Meyer, P. Schr¨ oder, and A.H. Barr. Implicit fairing of ir-
regular meshes using diffusion and curvature flow. ACM SIGGRAPH, pages
317–324, 1999.
[39] T.K. Dey, J. Giesen, and S. Goswami. Shape segmentation and matching
with flow discretization. In Workshop on Algorithms and Data Structures,
pages 25–36, 2003.
[40] James R. Diebel, Sebastian Thrun, and Michael Br¨ unig. A Bayesian method
for probable surface reconstruction & decimation. ACM Trans. Graphics,
25(1):39–59, 2006.
[41] M. Droske and M. Rumpf. A Levelset Formulation for Willmore Flow. In-
terfaces & Free Boundaries, 6:361–378, 2004.
[42] Meier D.S. and Fisher E. Atlas-based anatomic labeling in neurodegener-
ative disease via structure-driven atlas warping. Journal of Neuroimaging,
15(1):16–26, 2005.
[43] Gerhard Dziuk, Ernst Kuwert, and Reiner Schatzle. Evolution of elastic
curves inR
n
: Existence and computation. SIAM J. Math. Anal., 33(5):1228–
1245, 2002.
[44] A. Elad and R. Kimmel. On bending invariant signatures for surfaces. IEEE
Transactions on Pattern Analysis and Machine Intelligence, 25(10):1285–
1295, 2003.
[45] J. Feldmar and N. Ayache. Rigid, affine and locally affine registration of
free-form surfaces. International Journal of Computer Vision, 18(2):99–119,
1996.
[46] B. Fischl, M. I. Sereno, R. B. H. Tootell, and A. M. Dale. High-resolution
inter-subject averaging and a coordinate system for the cortical surface. Hu-
man Brain Mapping, 8:272–284, 1998.
121
[47] M. Gage and R.S. Hamilton. The heat equation shrinking convex plane
curves. Journal of Differential Geometry, 23:69–96, 1986.
[48] M. Garland and P.S. Heckbert. Surface simplification using quadric error
metrics. ACM SIGGRAPH, pages 209–216, 1997.
[49] T. Gatzke, C. Grimm, M. Garland, and S. Zelinka. Curvature maps for local
shape comparison. Proc. Shape Modeling International, 2005.
[50] N. Gelfand, N.J. Mitra, L.J. Guibas, and H. Pottmann. Robust global regis-
tration. In Eurographics Symposium on Geometry Processing, pages 197–206,
2005.
[51] Benjamin Gregorski, Joshua Senecal, Member-Mark A. Duchaineau, and
Member-Kenneth I. Joy. Adaptive extraction of time-varying isosurfaces.
IEEE Transactions on Visualization and Computer Graphics, 10(6):683–694,
2004.
[52] U. Grenander and M.I. Miller. Computational anatomy: An emerging disci-
pline. Quarterly of Applied Mathematics, 56(4):617–694, 1998.
[53] Eitan Grinspun, Anil N. Hirani, Mathieu Desbrun, and Peter Schr¨ oder. Dis-
crete shells. In ACM SIGGRAPH Symposium on Computer Animation, pages
62–67, 2003.
[54] X. Gu, Y. Wang, TF Chan, PM Thompson, and S.T. Yau. Genus zero surface
conformal mapping and its application to brain surface mapping. Medical
Imaging, IEEE Transactions on, 23(8):949–958, 2004.
[55] Sumit Gupta, Kuntal Sengupta, and Ashraf A. Kassim. Compression of dy-
namic 3d geometry data using iterative closest point algorithm. Comput. Vis.
Image Underst., 87(1-3):116–130, 2002.
[56] Stefan Guthe and Wolfgang Straßer. Real-time decompression and visualiza-
tion of animated volume data. In VIS ’01: Proceedings of the conference on
Visualization ’01, pages 349–356, Washington, DC, USA, 2001. IEEE Com-
puter Society.
[57] G. Hermosillo, O. Faugeras, and J. Gomes. Unfolding the Cerebral Cortex
Using Level Set Methods. Int. Conf. on Scale-Space Theories in Comp. Vi-
sion, 8:12–28, 1999.
[58] K. Hildebrandt and K. Polthier. Constraint-based fairing of surface meshes.
Proceedings of the 2007 Eurographics/ACM SIGGRAPH Symposium on Ge-
ometry processing, 2007.
122
[59] K. Hormann and M. S. Floater. Mean value coordinates for arbitrary planar
polygons. ACM Transactions on Graphics, mar 2006.
[60] M. K. Hurdal, K. Stephenson, P. L. Bowers, D. W. L. Sumners, and D. A.
Rottenberg. Coordinate system for conformal cerebellar flat maps. NeuroIm-
age, 11:S467, 2000.
[61] Lawrence Ibarria, Peter Lindstrom, Jarec Rossignac, and Andrzej Szymczak.
Out-of-core compression and decompression of large n-dimensional scalar
fields, 2003.
[62] Lawrence Ibarria and Jarek Rossignac. Dynapack: space-time compression
of the 3d animations of triangle meshes with fixed connectivity. In SCA
’03: Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on
Computer animation, pages 126–135, 2003.
[63] AE Johnson and M. Hebert. Using spin images for efficient object recogni-
tion in cluttered 3Dscenes. Pattern Analysis and Machine Intelligence, IEEE
Transactions on, 21(5):433–449, 1999.
[64] SC Joshi and MI Miller. Landmark matching via large deformation diffeo-
morphisms. Image Processing, IEEE Transactions on, 9(8):1357–1370, 2000.
[65] Tao Ju, Frank Losasso, Scott Schaefer, and Joe Warren. Dual contouring of
hermite data. ACM Trans. Graph., 21(3):339–346, 2002.
[66] Zachi Karni and Craig Gotsman. Compression of soft body animation se-
quences. Computers & Graphics, 28:25–34, 2004.
[67] L. Kharevych, W. Yang, Y. Tong, E. Kanso, JE Marsden, P. Schr¨ oder, and
M. Desbrun. Geometric, variational integrators for computer animation. Pro-
ceedings of the 2006 ACM SIGGRAPH/Eurographics symposium on Com-
puter animation, pages 43–51, 2006.
[68] R. Kikinis, M.E. Shenton, D.V. Iosifescu, R.W. McCarley, P. Saiviroon-
porn, H.H. Hokama, A. Robatino, D. Metcalf, C.G. Wible, C.M. Portas,
et al. A digital brain atlas for surgical planning, model-driven segmentation,
and teaching. IEEE Transactions on Visualization and Computer Graphics,
2(3):232–241, 1996.
[69] Martin Kilian, N. J. Mitra, and H. Pottmann. Geometric modeling in shape
space. In ACM Transactions on Graphics (SIGGRAPH), 2007.
[70] Eric Klassen, Anuj Srivastava, Washington Mio, and Shantanu H. Joshi.
Analysis of planar shapes using geodesic paths on shape spaces. IEEE Trans.
on Pattern Analysis and Machine Intelligence, 26(3):372–383, 2004.
123
[71] Leif Kobbelt. Discrete fairing and variational subdivision for freeform surface
design. The Visual Computer, 16(3-4):142–158, 2000.
[72] V. Kraevoy and A. Sheffer. Cross-parameterization and compatible remeshing
of 3D models. ACM Transactions on Graphics (TOG), 23(3):861–869, 2004.
[73] V. Kraevoy and A. Sheffer. Template-based mesh completion. Eurographics
Symposium on Geometry Processing, pages 13–22, 2005.
[74] Haeyoung Lee, Mathieu Desbrun, and Peter Schr¨ oder. Progressive encoding
of complex isosurfaces. ACM Trans. Graph., 22(3):471–476, 2003.
[75] Jerome Edward Lengyel. Compression of time-dependent geometry. In SI3D
’99: Proceedings of the 1999 symposium on Interactive 3D graphics, pages
89–95, 1999.
[76] Nathan Litke, Marc Droske, Martin Rumpf, and Peter Schr¨ oder. An image
processing approach to surface matching. In Proceedings of the Symposium
on Geometry Processing, 2005.
[77] X. Liu, H. Bao, H.Y. Shum, and Q. Peng. A novel volume constrained
smoothing method for meshes. Graphical Models, 64(3/4):169–182, 2002.
[78] William E. Lorensen and Harvey E. Cline. Marching cubes: A high reso-
lution 3d surface construction algorithm. In Computer Graphics (Proc. of
SIGGRAPH), volume 21, pages 163–169, 1987.
[79] JE Marsden and M. West. Discrete mechanics and variational integrators.
Acta Numerica, 10(00):357–514, 2003.
[80] Ajith Mascarenhas, Martin Isenburg, Valerio Pascucci, and Jack Snoeyink.
Encoding volumetric grids for streaming isosurface extraction. In 3DPVT
’04: Proceeding 2nd International Symposium on 3D Data Processing, Visu-
alization and Transmission, 2004, pages 665–672, 2004.
[81] T. McInerney and D. Terzopoulos. Deformable models in medical image
analysis. Mathematical methods in biomedical image analysis, pages 171–180,
1996.
[82] M. Meyer, M. Desbrun, P. Schr¨ oder, and A.H. Barr. Discrete differential-
geometry operators for triangulated 2-manifolds. In Proceedings of the Int.
Workshop on Visualization and Mathematics, May 2002.
[83] P.W. Michor and D. Mumford. Riemannian geometries of space of plane
curves, 2005. Preprint.
124
[84] M.I. Miller, A. Trouv´ e, and L. Younes. Geodesic Shooting for Computational
Anatomy. Journal of Mathematical Imaging and Vision, 24(2):209–228, 2006.
[85] Niloy J. Mitra, Natasha Gelfand, Helmut Pottmann, and Leonidas Guibas.
Registration of point cloud data from a geometric optimization perspective.
In SGP ’04: Proceedings of the 2004 Eurographics/ACM SIGGRAPH sympo-
sium on Geometry processing, pages 22–31, New York, NY, USA, 2004. ACM
Press.
[86] N.J. Mitra, N. Gelfand, H. Pottmann, and L. Guibas. Registration of point
cloud data from a geometric optimization perspective. In Eurographics Sym-
posium on Geometry Processing, pages 23–32, 2004.
[87] Ken Museth, David E. Breen, Ross T. Whitaker, and Alan H. Barr. Level
set surface editing operators. ACM Transactions on Graphics (SIGGRAPH),
21(3):330–338, 2002.
[88] G. G. Nahas and T. F. Burks, editors. Drug Abuse in the Decade of the Brain.
IOS Press, 1997.
[89] Mitra N.J, Flory S., Ovsjanikov M., Gelfand N., Guibas L., and Pottmann H.
Registration of point cloud data from a geometric optimization perspective.
In Proceedings of the 2007 Eurographics/ACM SIGGRAPH Symposium on
Geometry processing, 2007.
[90] Stanley Osher and Ronald Fedkiw. Level Set Methods and Dynamic Implicit
Surfaces, volume 153 of Applied Mathematical Sciences. Springer-Verlag, New
York, 2003.
[91] Thompson P. and Toga A. A framework for computational anatomy. Com-
puting and Visualization in Science, 5:1–12, 2002.
[92] M. Pauly, N.J. Mitra, J. Giesen, M. Gross, and L.J. Guibas. Example-based
3D scan completion. Symposium on Geometry Processing, pages 23–32, 2005.
[93] Jianbo Peng, Daniel Kristjansson, and Denis Zorin. Interactive modeling of
topologically complex geometric detail. ACM Trans. Graph., 23(3):635–643,
2004.
[94] F. Pighin, R. Szeliski, and D. Salesin. Resynthesizing facial animation
through 3d model-based tracking. Proceedings, International Conference on
Computer Vision, 1999.
[95] U. Pinkall and K. Polthier. Computing discrete minimal surfaces and their
conjugates. Experimental Mathematics, 2(1):15–36, 1993.
125
[96] BM Planitz, AJ Maeder, and JA Williams. The correspondence framework
for 3D surface matching algorithms. Computer Vision and Image Under-
standing, 97(3):347–383, 2005.
[97] JPW Pluim, JBA Maintz, and MA Viergever. Mutual-information-based reg-
istration of medical images: a survey. Medical Imaging, IEEE Transactions
on, 22(8):986–1004, 2003.
[98] S. Rusinkiewicz and M. Levoy. Efficient variants of the ICP algorithm. Proc.
3DIM, pages 145–152, 2001.
[99] S. Sandor and R. Leahy. Surface-based labeling of cortical anatomy using
a deformable atlas. Medical Imaging, IEEE Transactions on, 16(1):41–54,
1997.
[100] J. Schreiner, A. Asirvatham, E. Praun, and H. Hoppe. Inter-surface mapping.
ACM Transactions on Graphics (TOG), 23(3):870–877, 2004.
[101] J. A. Sethian. Level Set Methods and Fast Marching Methods, volume 3 of
Monographs on Appl. Comput. Math. Cambridge University Press, Cam-
bridge, 2nd edition, 1999.
[102] D.W. Shattuck and R.M. Leahy. BrainSuite: an automated cortical surface
identification tool. Med. Image Anal, 6(2):129–142, 2002.
[103] A. Sheffer and V. Kraevoy. Pyramid coordinates for morphing and defor-
mation. 3D Data Processing, Visualization and Transmission, 2004. 3DPVT
2004. Proceedings. 2nd International Symposium on, pages 68–75, 2004.
[104] Bong-Soo Sohn, Chandrajit Bajaj, and Vinay Siddavanahalli. Feature based
volumetric video compression for interactive playback. In VVS ’02: Pro-
ceedings of the 2002 IEEE symposium on Volume visualization and graphics,
pages 89–96, Piscataway, NJ, USA, 2002. IEEE Press.
[105] A. Stern and M. Desbrun. Discrete geometric mechanics for variational time
integrators. International Conference on Computer Graphics and Interactive
Techniques, pages 75–80, 2006.
[106] Carsten Stoll, Zachi Karni, Christian R¨ ossl, Hitoshi Yamauchi, and Hans-
Peter Seidel. Template deformation for point cloud fitting. In Symposium on
Point-Based Graphics, pages 27–35, 2006.
[107] G. Sundaramoorthi, A. Yezzi, and A. C. G. Mennucci. Sobolev active con-
tours. IJCV, 73(3):345–366, 2007.
126
[108] Liu T., Shen D., and Davatzikos C. Deformable registration of cortical struc-
tures via hybrid volumetric and surface warping. NeuroImage, 22:1790–1801,
2004.
[109] Gabriel Taubin. A signal processing approach to fair surface design. In
Proceedings of ACM SIGGRAPH, pages 351–358, August 1995.
[110] Paul M. Thompson, Michael S. Mega, Christine Vidal, Judith Rapoport, and
Arthur W. Toga. Detecting disease-specific patterns of brain structure using
cortical pattern matching and a population-based probabilistic brain atlas.
In Proc. 17th IPMI2001, Davis, CA, USA, pages 488–501, 2001.
[111] D. Tosun and J.L. Prince. Cortical Surface Alignment Using Geometry
Driven Multispectral Optical Flow. Inf Process Med Imaging, 19:480–92,
2005.
[112] Duygu Tosun, Maryam E. Rettmann, and Jerry L. Prince. Mapping tech-
niques for aligning sulci across multiple brains. Medical Image Analysis,
8(3):295–309, 2005.
[113] Toga A. W. and Mazziotta J. C.˙ Brain Mapping: The Methods.Academic
Press, second edition, 2002.
[114] Y. Wang, X. Gu, K.M. Hayashi, T.F. Chan, P.M. Thompson, and S.T. Yau.
Brain surface parameterization using riemann surface structure. In MICCAI
2005, pages 657–665, 2005.
[115] Y. Wang, B.S. Peterson, and L.H. Staib. 3D Brain surface matching based on
geodesics and local geometry. Computer Vision and Image Understanding,
89(2-3):252–271, 2003.
[116] M. Wardetzky, M. Bergou, D. Harmon, D. Zorin, and E. Grinspun. Discrete
quadratic curvature energies. Computer Aided Geometric Design, to appear,
2007.
[117] R. P. Woods, S. T. Grafton, C. J. Holmes, S. R. Cherry, and J. C. Mazziotta.
Automated image registration: I. J. Comp. Assist. Tomogr., 22:139–152,
1998.
[118] R.P. Woods, S.T. Grafton, J.D.G. Watson, N.L. Sicotte, and J.C. Mazziotta.
Automated image registration: II. Intersubject validation of linear and non-
linear models. J Comput Assist Tomogr, 22(1):153–165, 1998.
[119] C. Xu and JL Prince. Snakes, shapes, and gradient vector flow. IEEE TIP,
7(3):359–369, 1998.
127
[120] Guoliang Xu, Qing Pan, and Chandrajit L. Bajaj. Discrete surface model-
ing using partial differential equations. Computer Aided Geometric Design,
23(2):125–145, 2006.
[121] A.J. Yezzi and A. Mennucci. Conformal metrics and true ”gradient flows”
for curves. In International Conference on Computer Vision, pages 913–919,
2005.
[122] Laurent Younes. Optimal matching between shapes via elastic deformations.
Image and Vision Computing, 17(5-6):381–389, 1999.
[123] S. Zelinka and M. Garland. Similarity-based surface modelling using geodesic
fans. Proceedings of the 2004 Eurographics/ACM SIGGRAPH symposium on
Geometry processing, pages 204–213, 2004.
128

Abstract (if available)

Abstract People have been studying shapes since the ancient times, using geometry to model those that already existed in the world and to create new ones. Today, with the power of modern computers, the task of capturing the world digitally and recreating it virtually is no longer science fiction. This challenging, but extremely important problem spans applications ranging from medical modeling and scientific visualization to computer-aided design, art and entertainment. However, to achieve that goal, we need better algorithms that can "understand" geometry.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Feature-preserving simplification and sketch-based creation of 3D models

PDF

Rapid creation of photorealistic large-scale urban city models

PDF

Articulated human body deformation from in-vivo 3d image scans

PDF

Hybrid mesh/image-based rendering techniques for computer graphics applications

PDF

Automatic image matching for mobile multimedia applications

PDF

Interactive rapid part-based 3d modeling from a single image and its applications

PDF

Machine learning methods for 2D/3D shape retrieval and classification

PDF

A transparent framework of trust-based collaborative decision-making

PDF

Line segment matching and its applications in 3D urban modeling

PDF

Accurate image registration through 3D reconstruction

PDF

Hybrid methods for music analysis and synthesis: audio key finding and automatic style-specific accompaniment

PDF

Context-based information and trust analysis

PDF

Combining object recognition and tracking for augmented reality

PDF

Geometric methods of image registration and analysis

PDF

3D deep learning for perception and modeling

PDF

Facial gesture analysis in an interactive environment

PDF

Deformable geometry design with controlled mechanical property based on 3D printing

PDF

Part based object detection, segmentation, and tracking by boosting simple shape feature based weak classifiers

PDF

Software quality analysis: a value-based approach

PDF

Anatomically based human hand modeling and simulation

Asset Metadata

Creator Eckstein, Ilya (author)

Core Title Correspondence-based analysis of 3D deformable shapes: methods and applications

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Computer Science

Publication Date 10/18/2007

Defense Date 09/04/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag computer graphics,geometry processing,medical image analysis,OAI-PMH Harvest

Language English

Advisor Kuo, C.C. Jay (committee chair), Desbrun, Mathieu (committee member)

Creator Email ilyaeck@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m875

Unique identifier UC1313830

Identifier etd-Eckstein-20071018 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-560694 (legacy record id),usctheses-m875 (legacy record id)

Legacy Identifier etd-Eckstein-20071018.pdf

Dmrecord 560694

Document Type Dissertation

Rights Eckstein, Ilya

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu