An Intrinsic Geometric
Framework for Simultaneous
NonRigid Registration and
Segmentation of Surfaces *
N. Lord, J. Ho, and B.C. Vemuri
email: [nal, jho, vemuri]@cise.ufl.edu
UFCISE Technical Report TR06001
Computer and Information Sciences Department
University of Florida
Bldg. CSE, Room 301
Gainesville, Florida 32611
April 18, 2006
* This research was supported in part by the grant NIHR01NS04812 to BCV and the
UF Alumni Fellowship to NAL. Data sets were provided by Dr. Christiana Leonard
of the UF McKnight Brain Institute.
Abstract. In clinical applications where structural asymmetries between
homologous shapes have been correlated with pathology, the questions
of definition and quantification of 'asymmetry' arise naturally. When not
only the degree but the position of deformity is thought relevant, asym
metry localization must also be addressed. Asymmetries between paired
shapes can and have been formulated in the literature in terms of (non
rigid) diffeomorphisms between the shapes. For the infinity of such maps
possible for a given pair, we define optimality as the minimization of
total distortion, where 'distortion' is in turn defend as deviation from
isometry. We thus propose a novel variational formulation for segmenting
asymmetric regions from surface pairs based on the minimization of a
functional of both the deformation map and the segmentation boundary,
which controls gradient discontinuity of the map. This minimization is
achieved via a quasisimultaneous evolution of the map and curve. Our
formulation is inherently intrinsic and parameterizationindependent. We
present examples using both synthetic data and pairs of left and right
hippocampal structures, hippocampus malformation being linked with
such neurological disorders as !il and schizophrenia.
1 Introduction
Many neurological disorders have been clinically correlated with shape abnormal
ities in specific brain structures. For example, hippocampal shrinkage is associ
ated with epilepsy and schizophrenia, among other conditions. A precise analysis
of such abnormalities, which are usually considered to be failures of symmetry
of the left and right halves of the structure in question, could unlock the abil
ity not only to track the progression of existent disease, but to identify atrisk
individuals for preventative treatment. However, devising a framework in which
to conduct the necessary shape comparison is a nontrivial matter. The short
comings of raw volumetric comparison are obvious. Landmarkbased methods
are inapplicable to problems involving many anatomical structures (including
the hippocampus), because many biological structures do not exhibit the sorts
of readily and consistently identifiable local features on which such methods de
pend. The medial surface representation of Gerig et al. [1][2] represents a more
sophisticated take on volumetric analysis that includes some localization of vol
ume disparity, but the economy of the method leads both to noise sensitivity and
an acknowledged inability to distinguish between different sorts of deformation
(e.g. certain forms of elongation vs. bending). This is due to the expression of
complicated behaviors in terms of simple distances which discard some amount
of directionality.
More recently, Wang et al.[3][4][5] presented a technique involving the align
ment of surfaces based on the formation of a 2D parametricc) diffeomorphic
map from Riemannian surface structure information, with this map representing
the nonrigid registration between the surfaces. This approach arguably makes
use of what surface I. ,ilre" information is there to be had from relatively
nondescript shapes, without imposing an illdefined demand for landmark iden
tification. However, in forming the correspondence by maximizing the mutual
information of the chosen surface characteristics over the entirety of both sur
faces, these methods ignore key considerations. For one, while the conceptual
division of the body into anatomical structures is hardly arbitrary in a general
sense, for the purpose of deformation analysis, it may very well be. That is to
say, we may wish to identify a selected portion of the hippocampus as struc
turally abnormal relative to its correspondent, rather than regarding the entire
hippocampus as deformed. Furthermore, the MI of the domain and target surface
features does not in and of itself restrict the differential properties of the map,
only what properties should be found at given locations on the target surface.
Therefore two regions of the target surface that exhibit the same Riemannian
characteristics cannot be differentiated between, despite the fact that one choice
may be physically unreasonable as a map target. Thus this framework demands
a separate regularizing constraint on the map evolution.
To address these concerns, we present a novel approach to the nonrigid reg
istration of anatomically correspondent surfaces which rests on the fundamental
assumptions that (a) the relative deformity of anatomical structures for which
symmetry is expected can be intuitively and precisely quantified as the devia
tion from isometry of the deformation map between their surfaces, and (b) since
deformities may well be locally confined to certain regions of the given surfaces,
the evolution of the global correspondence must allow for a partial disconnect
between normal and abnormal regions, as these are by definition not expected
to exhibit the same deformation patterns. As such, surface segmentation rests
at the heart of our approach, which can be viewed not only as a registration
tool but also as an asymmetry localizer and quantifier. Our definition of system
energy in terms of isometry inherently balances the consideration that the map
should accurately match Riemannian surface characteristics with the idea that
the map should incur as little deformation as possible in doing so. This gives
way to a simple, elegant formulation relatively free of heuristics. Since we here
restrict our consideration to surfaces that are topologically cylindrical (appropri
ate for hippocampi among many other structures), we are able to parameterize
with a single patch (periodic at a pair of ends) for each surface in the pair, and
as the formulation is intrinsic, we can confine evolution to these 2D parametric
domains, retaining elegance. The evolution is driven quasisimultaneously by the
energies of both the segmenting curve and the deformation field, in the spirit of
the works of Yezzi et al. and Unal et al.[6] [7], where the former group restricted
its work to rigid image registration and the latter generalized the idea to the
nonrigid deformation case. In [8], Vemuri et al. presented a registrationassisted
image segmentation technique. Of the many proposals in the literature which
followed, however, none addressed the issue of joint manifold segmentation and
registration, which is our focus. Our work differs notably in that we evolve the
segmenting curve and the registration map on a general 2D manifold (rather than
a flat image plane), and rather than registering and segmenting image data, we
register and segment the shape characteristics of the very surface on which the
processes are operating.
The rest of the paper is organized as follows: Section 2 contains the mathe
matical formulation and numerical algorithm used to solve it. Section 3 contains
experimental results on both synthetic and real data sets. In section 4, we draw
conclusions.
2 Simultaneous Segmentation and Registration
Algorithm
Let S1 and S2 be two surfaces in IR3. The Euclidean metric in IR3 induces Rie
mannian metrics gl and g2 on S1 and S2, respectively. The goal of the algorithm
is to segment regions in S1 and their corresponding regions in S2 such that the
metric structures of the corresponding regions match optimally. More specifi
cally, the algorithm computes a homeomorphism f between Si and S2 and a
set of closed curves 71 on Si. Let 72 denote the closed curve in S2 that is the
image of i7 under f. The restriction of f to the complement Sl\7y is a diffeo
morphism between S1\7y and S2\y2. If U1, , U, and V1,  V, denote the
collections of open components of S1\7i and S2\72, respectively, then f maps
Ui diffeomorphically onto Vi for each i and f matches the Riemannian structures
between that of Ui and of Vi = f(U). See Figure 1. We solve the simultaneous
St . s2
Is ii
Fig. 1. Left: The two given surfaces, S1 and S2. Right: The desired result of simul
taneous segmentation and registration of these two surfaces. The map f is a homeo
morphism between S1 and S2. f maps the segmenting curve 71 on Si to a curve 72
on S2, and f establishes correspondences between connected components of Si\7i and
S2\72. U1 and U2 are two disjoint open sets in Si\7i, and V1 and V2 are two disjoint
open sets in S2\72. f maps Ui diffeomorphically onto V for each i. The restriction of
f to Ui optimally matches the Riemannian structures of Ui and V.
segmentation and registration problem outlined above using a variational frame
work. The energy functional E is defined as a functional of a pair (f, y7), where
71 is a set of closed curves on S1 and f : S1 S2 is a homeomorphism which is
C on Sl\yi. Let f*g2 denote the pullback metric (the first fundamental form)
on gi, and S(f, i7) is defined as
'(f, ~) I f*g2 g2dA+ df(lQ(t))/ /' ./ (1)
JS1\7i .
The second integral computes the length of the segmenting curve 72 on S2, while
the first integral computes the matching (registration) cost between the two Rie
mannian structures on the complement, Si\7i. By computing the length of the
segmenting curve on S2 instead of Si, the energy functional tightly couples the
two somewhat disparate processes, segmentation and registration. The integrand
of the first term above, If*g2 gfl2, is the usual L2norm between the two ten
sors g1 and f*g.9 In local coordinates, it is given as the usual L2norm between
the two matrices Gi and JtG2J:
IJtG2J Gi 9, (2)
where J is the Jacobian of f expressed in the local coordinates system, and Gi
and G2 are the usual 2 x 2 positive definite matrices expressing the metric tensors
gi and g2 in the given local coordinates on S1 and the induced one (using f) on
S2.
At this point, a comparison between our approach and the usual Mumford
Shah framework [9] is called for. Given an image I, the energy functional for
MumfordShah is
MS(I' ) j 1I' 12dA+ VIIdA+ [ (t)Idt.
In the above, D denotes the image and 7 the segmenting curve on D. The third
term, which measures the length of 7, corresponds directly to the second term
of C in 1. The first two terms of MumfordShah, representing data fidelity and
smoothness, respectively, may seem to lack analogous entries in E. However,
we argue below that the first integral in the definition of C serves both as the
fidelity and smoothing term for I, and we can thus consider C as a generalization
of the usual MumfordShah functional to the problem of simultaneous shape
registration and segmentation.
For our problem, we want to compare two surfaces B1 and B2 via their
respective Riemannian structures gi and g2, and to extract a meaningful seg
mentation from the comparison. Instead of an image I' as in the usual Mumford
Shah framework, our variable is the map f between Si and S2. In our case, the
data fidelity requirement appears as the penalization of f for failure to match
the Riemannian structures of the surfaces, while the smoothing term appears
as the penalization of f for large firstorder derivative magnitudes. A close ex
amination of Equation 1 reveals that the first integral indeed contains both
requirements. First, the norm between the two tensors f*g9 and gi as defined
in Equation 2 clearly measures the quality of match between the two metric
tensors. Second, the Jacobian J contains all of the firstorder derivatives of f,
and the squarednorm in Equation 2 indirectly measures the magnitudes of the
firstorder derivatives of f using both metrics gl and g2. In sum, we can interpret
the first term in Equation 1 as the integration of the first two terms in the usual
MumfordShah functional MS. In MS, these two terms are necessarily separate
because there is no natural connection between matching pixel intensities and
restricting the gradient VI. In our problem, the connection between matching
and smoothing is natural: since we are matching Riemannian structures, which
are tensors defined on tangent spaces, any matching between these two struc
tures has to involve some kind of (implicit) matching between tangent spaces,
and the appearance of derivatives in defining the matching is then unavoidable.
2.1 Minimizing the Energy Functional
In implementing the energy minimization as defined above, we rely on alter
nating iterations between a ChanVese level set curve segmenter [10], extended
to the general 2D manifold domain, which produces a step towards the best
estimate of the segmentation of the energy field given the current map, and a
sparse Hessianexploiting quasiNewton optimization scheme acting on the non
rigid map between the parametric domains, which provides a step towards the
minimumenergy map between the surfaces as constrained by the current po
sition of the segmenting curve. The manner in which the curve separates the
evolution occurring on its exterior from that occurring on its interior relates to
a Neumann condition imposed across the boundary, akin to what is suggested
in the MumfordShah curve evolution implementation in [11]. The curve length
minimizing term implementation is straightforward, though curvature computa
tion must take nonuniform surface length elements into account. Hippocampus
surface data used in this investigation was brought into rigid alignment by ICP
in a preprocessing step, then 2D sliced and manually segmented by a trained
neuroanatomist into a 40x21 grid periodic in one direction. This grid in effect
forms a parametric domain for the corresponding surface, and the Riemannian
surface characteristics are thus calculated by numerical methods and stored in
said domain.
3 Results
We here present examples of the algorithm in action, validating its behavior on
synthetic cases and demonstrating results obtained on real hippocampus pairs.
As a starting point, we consider the trivial case of a pair of cylinders, where
one member of the pair has had its surface distorted according to an outward
normal vector field of magnitude dictated by a Gaussian distribution. That is
to say, its surface has been bumped. In Fig. 2, we see a comparative case in
volving a cylinder with two such distortions versus its unaltered counterpart.
The intuition behind our method asserts that we should "blame" the applied
bumps for observed differences between the Riemannian surface characteristics
of the two shapes. Obviously, then, the desired result is the segmentation of the
two clearly visible bumps on the left cylinder, which is what we observe. (Red
'hot spots,' the curve interiors, represent regions of highest deformation energy
under the current map and segmentation. The curve itself is found at the border
between the red and blue regions.) Note the textured quality of the surfaces,
which occurs as a result of the addition of a SD=0.02 zeromean Gaussian noise
field to the (X,Y,Z) coordinates of the (unit radius) cylinders prior to output,
for demonstration of robustness in the face of noise in surface scanning.
As indicated previously, the case of the cylinder is a trivial one in this context
in that there exists a "natural" parameterization of the surface which results in
a constant metric tensor field over the parametric domain. Thus, we further val
idate using synthetic alterations to real data exhibiting arbitrary shape, within
the confines of cylindrical topology. Since the hippocampus meets these criteria
and forms part of our clinical motivation, it seems a natural choice. Selecting
a single hippocampal horn, we again apply a small outward normal Gaussian
deformation to its surface, whose location we hope to see the segmenter recover.
However, rather than illustrating invariance to surface point cloud noise, we now
seek to additionally validate invariance of the solution to reparameterizations of
the given surfaces. To pause for conceptual refreshment, recall that in compar
ing surfaces we are in fact comparing two parametric domains based on the two
tensorvalued functions corresponding to each of those domains, those functions
being the Riemannian surface structure information of each parameterized sur
face. But just as two surfaces (or their parametric domains, as we do in our
implementation) can be mapped to one another nonrigidly in infinitely many
ways, so too can each parametric domain map to its topologicallyequivalent pa
rameterized surface in infinitely many ways. The Riemannian structure tensor
fields themselves depend upon which of these infinitely many parameterizations
are chosen. When synthetically producing one surface from another, we pos
sess a trivial alignment between the domains under which the two tensor fields
are identical, excepting the patch corresponding to the synthetic bump. In real
cases, of course, data will not have this sort of property. The validation is thus
incomplete until we can show that under a suitably deformed reparameteriza
tion of one of the surfaces being compared, a result comparable to the above is
achieved. We can reparameterize simply by randomly generating a distortion of
the parametric domain's regular grid, then interpolating a new point cloud from
the given data based on the new locations of the gridpoints. It is completely
appropriate to confine said distortion to the order of a pixel in each direction:
after all, we have assumed close rigid alignment and sensibly spaced data points
from our preprocessing, and if these constraints are violated, they can always
be reinstated by repeating the algorithms that enforced them in the first place.
The purpose of the nonrigid registration evolution is to further refine the surface
correspondence following bulk rigid alignment of the shapes based on the sur
face characteristics observed locally, not to register raw point clouds. We cannot
show the parameterization warp used due to space constraints: it is of 0.1 SD in
each coordinate.
It is visually evident that the segmentation succeeds despite the change of
parameterization of one of the surfaces being compared. The deformation map
II ,.
Fig. 2. Left column: Segmentation of prominent distortions (bumps) from cylindrical
surface, based on comparison to homologous cylinder lacking bumps. Surfaces are both
distorted with distinct random noise fields (seen as textured effect). Top: naive initial
ization; Middle: 3 iterations; Bottom: 10 iterations. Right column: Segmentation of
synthetic bump from surface of real hippocampus, based on comparison to same hip
pocampus without bump. The two hippocampi are under different parameterizations.
Top: naive initialization; Middle: 5 iterations; Bottom: 20 iterations.
between the parametric domains (not shown) undergoes considerable evolution
through this process. This evolution can be best understood as the method's
attempt to invert the reparameterization, to achieve the optimal alignment of the
identical surface patches external to the segmented bump, while at the same time
trying to find the minimum energy description of the acknowledged deformation
of the patch on the curve's interior.
In presenting the following real example, we issue the caveat that the method's
success lies in the eye of the beholder, because there is of course no ground truth
associated with this sort of problem. For this reason we have selected a case
with an ".i ,. i which can easily be agreed upon. Fig. 3 shows the left and
right hippocampi clinically segmented from an MRI scan of an epilepsy patient.
We take it as selfevident that there exists a localizable region of high distor
tion between the left and right shapes, and assume that one will agree that our
method has in fact converged to a segmentation of that region. Illustrations of
the segmented energy field and parametric warp evolutions have been withheld
due to space limitations.
Fig. 3. Segmentation of distortion between real hippocampi of epileptic. Top left: naive
initialization; Top right: 1 iteration; Bottom left: 5 iterations; Bottom right: 20 itera
tions.
4 Conclusion
We have presented a novel scheme for simultaneous nonrigid registration and seg
mentation of homologous shapes, wherein the interdependent registration and
segmentation processes are driven by intrinsic geometric characteristics of the
shapes themselves. As such, we are able not only to identify surface pairs rep
resenting large deformations, but also to specify which subregions of the shapes
appear most likely to be involved in the deformation, to produce an estimate
of the implicit deformation field, and to quantify the deformation energy of the
segmented subregions separately from that of the remaining surface patches.
Given the number of considerations interplaying, the variational principle driv
ing the evolution is notably compact. We have seen that the method's results
are reasonably invariant to local parameterization changes and added surface
noise, and results on real data (while inherently lacking in ground truth for
comparison) appear promising. Future efforts could include the application of
the method to extensive real datasets that have been clinically preclassified ac
cording to the presence or absence of pathology, as well as an attempt to correlate
the locations and intensities of the segmented deformed regions with the given
pathology classifications. It should also be noted that the tworegion ChanVese
segmentation method (which is often used for its strong balance of effective
ness and relative simplicity of implementation) is merely a good "first draft" of
a segmentation scheme for this sort of problem. A more robust segmentation
method, such as the general MumfordShah functional extended to arbitrary
surfaces, would doubtless prove beneficial, as our energy fields will not always
neatly obey the twomean framework. Finally, it goes without saying that all
evolution/optimizationbased methods (apart from the most trivial cases) suffer
from some amount of initializationdependence. Intelligent initialization schemes
could be expected to further boost performance.
References
1. Gerig, G., Muller, K.E., Kistner, E.O., Chi, Y.Y., Chakos, M., Styner, M., Lieber
man, J.A.: Age and treatment related local hippocampal changes in schizophrenia
explained by a novel shape analysis method. In: MICCAI. (2003) 653660
2. Styner, M., Lieberman, J.A., Gerig, G.: Boundary and medial shape analysis of
the hippocampus in schizophrenia. In: MICCAI. (2003) 464471
3. Wang, Y., Gu, X., Hayashi, K.M., Chan, T.F., Thompson, P.M., Yau, S.T.: Surface
parameterization using riemann surface structure. In: Proc. Int. Conf. on Computer
Vision. (2005) 10611066
4. Wang, Y., ('C I.... M.C., Thompson, P.M.: Mutual informationbased 3d surface
matching with applications to face recognition and brain mapping. In: Proc. Int.
Conf. on Computer Vision. (2005) 527534
5. Wang, Y., C('I .... M.C., Thompson, P.M.: Automated surface matching using
mutual information applied to riemann surface structures. In: MICCAI. (2005)
777674
6. A.Yezzi, Zollei, L., Kapur, T.: A variational framework for joint segmentation and
registration. In: CVPR MMBIA. (2001)
7. Unal, G.B., Slabaugh, G.G.: Coupled pdes for nonrigid registration and segmen
tation. In: CVPR. (2005) 168175
8. Wang, F., Vemuri, B.C.: Simultaneous registration and segmentation of anatomical
structures from brain mri. In: MICCAI. (2005) 1725
9. Mumford, D., Shah, J.: Optimal approximations by piecewise smooth functions
and associated variational problems. Commun. Pure Appl. Math 42(4) (1989)
10. Chan, T.F., Vese, L.A.: Active contours without edges. IEEE Transactions on
Image Processing 10(2) (2001) 266277
11. Tsai, A., Yezzi, A.J., Willsky, A.S.: Curve evolution implementation of the
mumfordshah functional for image segmentation, denoising, interpolation, and
magnification. IEEE Transactions on Image Processing 10(8) (2001) 11691186
