Mutual Information-based 3D Sur

Document Sample
 Mutual Information-based 3D Sur Powered By Docstoc
					    Mutual Information-based 3D Surface Matching with Applications to Face
                       Recognition and Brain Mapping

         Yalin Wang                         Ming-Chang Chiang                          Paul M. Thompson
    Mathematics Department               Laboratory of Neuro Imaging              Laboratory of Neuro Imaging
           UCLA                           UCLA School of Medicine                  UCLA School of Medicine                       

                       Abstract                               1. Introduction
   Face recognition and many medical imaging applica-             Many face recognition algorithms have been proposed in
tions require the computation of dense correspondence vec-    the last few decades [24]. Several approaches (e.g., “eigen-
tor fields that match one surface with another. In brain       face” methods) encode patterns of geometric and intensity
imaging, surface-based registration is useful for tracking    variation between faces, and compute metrics to determine
brain change, and for creating statistical shape models of    the degree of differences between individual faces. Related
anatomy. Based on surface correspondences, metrics can        work has focused on image matching for tracking facial
also be designed to measure differences in facial geome-      features in video images. However, all 2D (image-based)
try and expressions. To avoid the need for a large set of     face recognition systems are somewhat sensitive to facial
manually-defined landmarks to constrain these surface cor-     expressions and illumination conditions. 3D geometric sur-
respondences, we developed an algorithm to automate the       face matching can solve these problems and may offer better
matching of surface features. It extends the mutual infor-    recognition performance.
mation method to automatically match general 3D surfaces          Surface models are also widely used in medical imaging
(including surfaces with a branching topology). We use dif-   to assist in data visualization, nonlinear image registration,
feomorphic flows to optimally align the Riemann surface        and surface-based signal processing or statistics. Surface
structures of two surfaces. First, we use holomorphic 1-      models are often generated in computational anatomy stud-
forms to induce consistent conformal grids on both sur-       ies to support computations, e.g., when statistically com-
faces. High genus surfaces are mapped to a set of rectan-     bining or comparing 3D anatomical models across subjects,
gles in the Euclidean plane, and closed genus-zero surfaces   or mapping functional imaging parameters onto anatomical
are mapped to the sphere. Next, we compute stable geo-        surfaces. Often the comparison of data on two anatomi-
metric features (mean curvature and conformal factor) and     cal surfaces is required, and a correspondence field must be
pull them back as scalar fields onto the 2D parameter do-      computed to register one surface nonlinearly onto the other.
mains. Mutual information is used as a cost functional to     Multiple surfaces can be registered nonlinearly to construct
drive a fluid flow in the parameter domain that optimally       a mean shape for a group of subjects, and deformation map-
aligns these surface features. A diffeomorphic surface-to-    pings can encode shape variations around the mean. This
surface mapping is then recovered that matches surfaces in    type of deformable surface registration has been used to de-
3D. Lastly, we present a spectral method that ensures that    tect developmental and disease effects on brain structures
the grids induced on the target surface remain conformal      such as the corpus callosum and basal ganglia [26], the hip-
when pulled through the correspondence field. Using the        pocampus [7], and the cortex [27]. Nonlinear matching of
chain rule, we express the gradient of the mutual informa-    brain surfaces can also be used to track the progression
tion between surfaces in the conformal basis of the source    of neurodegenerative disorders such as Alzheimer’s dis-
surface. This finite-dimensional linear space generates all    ease [7], to measure brain growth in development [26], and
conformal reparameterizations of the surface. Illustrative    to reveal directional biases in gyral pattern variability [19].
experiments apply the method to face recognition and to the       Surface registration has numerous applications, but a
registration of brain structures, such as the hippocampus     direct mapping between two 3D surfaces is challenging
in 3D MRI scans, a key step in understanding brain shape      to compute. Often, higher order correspondences must
alterations in Alzheimer’s disease and schizophrenia.         be enforced between specific anatomical points, curved
                                                              landmarks, or subregions lying within the two surfaces.
One common way to achieve this is to first map each of          1.1. Previous Work
the 3D surfaces to canonical parameter spaces such as a
sphere [11, 3] or a planar domain [20]. The surface cor-           Some researchers [17, 21] incorporate a 3D model in
respondence problem can then be addressed by computing         face recognition research. Bronstein et al. [5] propose a
a flow in the parameter space of the two surfaces [26, 9],      3D face recognition approach based on geometric invari-
which induces a correspondence field in 3D. Furthermore,        ants to compute isometric deformations. Several variational
correspondences may be determined by using a minimum           or PDE-based methods have been proposed for match-
description length (MDL) principle, based on the compact-      ing surfaces. Surfaces may be represented by parametric
ness of the covariance of the resulting shape model [10].      meshes [11], level sets, or both representations [20]. An-
Anatomically homologous points can then be forced to           genent et al. [2, 1] represent the Laplace-Beltrami operator
match across a dataset. Recently, Twining et al. [28] pro-     as a linear system and implement a finite element approxi-
posed a theoretical framework to unify groupwise image         mation for parameterizing brain/colon surfaces via confor-
registration and average model construction. In their ap-      mal mapping. Gu et al. [13] found a unique conformal map-
proach, an information-based model of the correspondences      ping between any two genus zero manifolds by minimizing
among a group of images becomes a part of the registration     the harmonic energy of the map. Gu and Vemuri [12] also
process.                                                       matched 3D shapes by first conformally mapping them to
    By the Riemann uniformization theorem, all surfaces        a canonical domain and aligning their 2D representations
can be conformally embedded in a sphere, a plane or a          over the class of diffeomorphisms. They demonstrated their
hyperbolic space. The resulting embeddings form special        algorithm on genus zero closed surfaces.
groups. Using holomorphic 1-forms and critical graphs,             The mutual information (MI) method [30, 29] measures
global conformal parameterization [15] can also be used to     the statistical dependence of the voxel intensities between
conformally map any high genus surface (i.e., a surface with   two images. This measure of agreement can be used to
branching topology) to a set of rectangular domains in the     tune the parameters of a registration transform such that
Euclidean plane. In this paper, we show how to use con-        MI is maximal when the two images are optimally aligned.
formal parameterizations to assist in the matching of arbi-    The MI method has been successful for rigid [31] and non-
trary 3D face and anatomical surfaces. Mutual information      rigid [22, 25] image registration. Here, we generalize it to
is used to drive a diffeomorphic fluid flow that is adjusted     match 3D surfaces. For MI to work, a monotonic map-
to find appropriate surface correspondences in the param-       ping in grayscales between images is not required, so im-
eter domain. In this study, we chose the mean curvature        ages from different modalities can be registered [18]. Her-
and the conformal factor of the surfaces as the differetial    mosillo et al. [16] adopted linear elasticity theory to reg-
geometric features to be aligned, as these are intrinsic and   ularize the variational maximization of MI. D’Agostino et
stable. These choices are purely illustrative. In fact, any    al. [8] extended this approach to a viscous fluid scheme al-
scalar fields defined on the surfaces could be matched, e.g.,    lowing large local deformations, while maintaining smooth,
cortical thickness maps, or even functional imaging signals    one-to-one topology [6].
or metabolic data. Since conformal mapping and fluid reg-
istration techniques generate diffeomorphic mappings, the      1.2. Basic Idea
3D shape correspondence established by composing these
mappings is also diffeomorphic (i.e., provides smooth one-        Suppose
                                                                               is an oriented surface. The map from       to
to-one correspondences).                                       a local coordinate
                                                                                            plane is a conformal map when
                                                                                                ¦¤        ©
   We also present a spectral approach for ensuring that the   the first fundamental form satisfies:
                                                                                                                       %  !
                                                                                                                        &© $#¤ ©  "§ ¦¤ © 
                                                                       !               !
grid induced on the target surface by the correspondence            © 
                                                                    . Here      ¦¤ is called the conformal factor, a func-
field, remains conformal. Grid orthogonality is advanta-        tion that scales the metric at each point
                                                                    ¨¡                                            . We say           #¤ ¡
geous for accurate numerical discretization of PDEs or for     ©¥ ¦¤  is a conformal coordinate of
                                                                         ¥                                  . Locally, each                 (
signal processing on the resulting surface meshes. For high    surface patch is covered by a conformal coordinate chart.
genus surfaces, the global conformal parameterization is not   For high genus surfaces, the local conformal parameteriza-
unique and all the conformal parameterizations form a lin-     tion can be extended to cover the whole surface. By the
ear space. The degrees of freedom in this space of con-        Riemann-Roch theorem and the circle-valued Morse theo-
formal grids are tuned to maximize the mutual information      rem, a high genus surface (         ) can be completely cov-
                                                                                                                2 0
energy of features between the two surfaces. Because the       ered by a set of non-overlapping segments. Each segment
conformal structure is intrinsic and the conformal parame-     can be conformally mapped to a rectangle. With the Gauss
terization continuously depends on the Riemannian metric       and Codazzi equations, one can prove that a closed surface
                                                                      "§                                                         !
on the surface, our method is also stable and computation-     ¦¤ 4    in    with conformal parameter
                                                                               6                                 is uniquely
                                                                                                                      ¦¤ '§
ally efficient.                                                 determined by its conformal factor                        #¤ 
                                                                                                           and its mean cur-
 !          ¤
   "§  ¦98 "§  , up to a rigid motion. We call a tuple of
       ¤ !
 ¦@8 "§  ¦¤  and   a conformal representation of the sur-
face         . We can solve the surface registration problem
            ¦¤ 4
by computing intrinsic geometric features from the confor-
mal mesh, and aligning them in the parameterization do-
main. To align these scalar fields, we use a fluid registration
technique in the parameter domain that is driven by mutual
information. With conformal mapping, we essentially con-
vert the surface registration problem into an image registra-
                                                                                                                                                         Figure 1. Illustrates surface conformal struc-
tion problem, for which MI methods are especially advan-
                                                                                                                                                         ture. (a) shows the cut we introduce on a
tageous. Finally, by invoking the surface partitioning tech-
                                                                                                                                                         face surface. The face becomes an open
nique induced by holomorphic 1-forms, our surface-based
                                                                                                                                                         boundary genus one surface. (b) shows a
mutual information method works on general surfaces with
                                                                                                                                                         global conformal parameterization of the sur-
arbitrary topologies.
                                                                                                                                                         face. (c) shows the horizontal trajectory. (d)
                                                                                                                                                         shows the rectangle to which the face sur-
2. Theoretical Background                                                                                                                                face is conformally mapped.

2.1. Global Conformal Parameterization

   Suppose      ,     are two surfaces. Locally they can be ¡                                         ¡©            ¡                         map it to a square as in (d) and get its conformal parameter-
                                               © ¥ § ¦© 4  © ¥ § BA ¡ 4 A                                          §¡                       ization.
represented as             ,           , where
                                                        are  ¥¨¤ ¡                                       ¥
                                                                                                          #¤     © ¥ ¦¤¥
                    6H5 FG© 5ED© 4 4
their local coordinates, and                     are vector-
                                                               §                                                                                  An atlas is a collection of consistent coordinate charts on
valued functions. The first fundamental form of
    r a§ A  q§ih"§ 'cW aYWX USQ
               p `Y                                       is                                       S Q USQ                                    a manifold, where transition functions between overlapping
          2            , where
                       cW fdccW
                       gb e b
                                  `                        .     )                                     V $ ) TSRP ¡© 
                                                                                                      ¥ ¥               Q I                  coordinate charts are smooth.
Similarly, the first fundamental form of
                          ©@A S Q TSQ          is defined in                                                                                       We treat           © 5
                                                                                                                                                                as isomorphic to the complex plane, where
the same way:                        $ V &USQ I  ©©  ¡
                                     ¥ ¥ s)
                                          . Define a map-
                                                                                                                                                                       "§                 v     !
                                                                                                                                                                                                 qh %           ! w
                                                                                                                                             the point         is equivalent to
                                                                                                                                                                  ¦¤       !h
                                                                                                                                                                            z1w               , and         is        #¤
ping                     between two surfaces. Using lo-                                       © BA F                  A C
                                                                                                                        vut                   equivalent to     yv
                                                                                                                                                                     §    . Let be a surface in
                                                                                                                                                                                       §             with an            6
   © 5 F@© fPt
cal coordinates, can be represented as                     ,                                                                                                                                                  F |~„CG¨v
                  5 C                                                               § ¡ © § t© § ¡ ¡                                          atlas              , where
                                                                                                                                                                   ~¤
                                                                                                                                                                ’€¨v |z}          is a chart, and
                                                                                                                                                                                   ƒ¨v ‚q¤                             …
                               . Then any tangent vector               " © ¥ #¤ t ¡  ¥ #¤ q¤ § ¡  t
                                                                                           ¥                          ¥ t                         …
                                                                                                                                                 maps an open set             to the complex plane .
                                                                                                                                                                            G ~
  tw        on     will be mapped to a tangent vector                                                            A       © $ $#¤
                                                                                                                             ¥         ¥         An atlas is called conformal if (1). each chart
                                                                                                                                                                                                           €¨v ‡z¤
on     ,                       y 9x                                                                               €y
                                                                                                                  9x          © BA
                     ƒ G`‚ b y                     `‚‚ ’ `‚ ’3ˆ†„`cy
                                                       ‘“                     ‘i ‡‰ … ƒ ‚ €                                                  is a conformal chart. Namely, on each chart, the first funda-
                             b                 R‚ ’ ` ’
                                             ”• ‘i                             “
                                                                              ‘i                                                       (1)   mental form can be formulated as        ‰  v   v   
                                                                                                                                                                                      Šv ƒ¨ © ¨#¤ uˆ© s¡
                                                                                                                                                                                        y                    ;
                                                                                                                                              (2). the transition maps
                                                          i ` “
                                                          “                                                                                                                            F £~ 1 q¤  C ƒ1cv
                                                                                                                                                                                           ‹  ~ v          Ž v Œ ‹
                                      ™ – &–
 The length of    is                    . We use the length
                                           t w § w¡t  s ) ™ ˜ˆI w
                                                                    ™                — –                        t                                  ‹  ~¤ ‹
                                                                                                                                                  ƒE~ ’ z‘v are holomorphic.
of to define the length of               ¥© V $#¤
                                       . Namely, we define          ¥                                                             t w            A chart is compatible with a given conformal atlas if
a new metric for
t                      which is induced by the mapping
                                                                                                              A                               adding it to the atlas again yields a conformal atlas. A con-
and the metric on       . We call this metric the pull-back
                        –                                                                                                                     formal structure (Riemann surface structure) is obtained by
metric, and denote it by   w
                           t  . Replacing       in the above             ©© fet 9A
                                                                                                                                              adding all compatible charts to a conformal atlas. A Rie-
equation by (1), we get the analytic formula for the pull-                                                                                    mann surface is a surface with a conformal structure.
back metric,r $k j k € j €                                                                                                                        One coordinate chart in the conformal structure in-
                 y yq                                  qq
                                    W W b — ` b m m p ljv… ‚
                                                                                  €                      i m k i ‚ ’a¨€hyg
                                                                                                                                              troduces a conformal parameterization between a surface
                b b cW cW                                                                                                               (2)
                               gb db                                                   Tgd
                                                                                             n                                                patch and the image plane. The conformal parameteriza-
                                                                                                                                              tion is angle-preserving and intrinsic to the geometry, and
 We call a conformal mapping, if there exists a positive                                                                                      is independent of the resolution and triangulation.
                         t               §¡                                             ¡© su © ¥ § ¡ ¦¤ t © sfet
scalar function      § ¡             © ¥ ¦¤ 
                             , such that    ¥                .                                       ¥  © d                                  Locally, a surface patch is covered by a conformal coor-
where            © ¥ #¤ 
                    is called the conformal factor.
                           ¥                                                        ¡                                                         dinate chart. For high genus surfaces, the local conformal
    Intuitively, all the angles on       are preserved on    .                A                                                    © A        parameterization can be extended to cover the whole sur-
Figure 1 shows a conformal mapping example. Figure 1                                                                                          face except at several points. These exceptional points are
shows a face surface. We introduce a cut on the nose (a) and                                                                                  called zero points. By the Riemann-Roch theorem, there are
                                                                                                                                                      r w r
change it to a genus one surface. We illustrate the conformal                                                                                        zero points on a global conformal structure of a genus
parameterization via the texture mapping of a checkerboard                                                                                    )  closed surface. By the circle-valued Morse theorem, the
in (b). By tracing a horizontal trajectory (c), where the ini-                                                                                iso-parametric curves through the zero points segment the
tial tracing point was manually selected, we conformally                                                                                      whole segment the whole surface to patches, where each
patch is either a topological disk, or a cylinder. The seg-                                    tion for a high genus surface is fundamentally determined
mentation is determined by the conformal structure of the                                      by the choice of the holomorphic 1-form.
surface and the choice of the global conformal parameteri-                                        Figure 3 shows two different conformal parameteriza-
zation.                                                                                        tions of a lateral ventricle surface of a HIV/AIDS patient
   Figure 2 shows an example of the conformal parameter-                                       subject. (a) shows a uniform parameterization result and
ization of a dog surface model. We introduce 3 cuts on the                                     (b) shows a nonuniform parameterization result. Note that
surface and change it to a genus 2 surface. The computed                                       although both of these are conformal, one has greater area
conformal structure is shown in (a). (b) shows a partition of                                  distortion than the other.
the dog surface, where each segment is labeled by a unique
color. (d) shows the parameterization domain. Each rectan-
gle is the image, in the parameterization domain, of a sur-
face component in (c).

                                                                                                  Figure 3. Shows a uniform (left panel) and
                                                                                                  a non-uniform (right panel) global conformal
                                                                                                  parameterization for the same surface, the
    Figure 2. Illustrates conformal parameteriza-                                                 lateral ventricles of the human brain.
    tion for a high genus surface. (a) shows the
    conformal structure for a model of a dog. Af-
    ter introducing cuts between both ears and
    on the bottom, we turn the dog model into an                                               2.3. Conformal Representation of a General Sur-
    open boundary genus 2 surface. (b) shows a                                                       face
    partition of this surface, with a unique color                                                                                                                         
                                                                                                  For a general surface , we can compute conformal co-
    labeling each part. (c) shows the parameter-                                               ordinates      ¦¤
                                                                                                                to parameterize . Based on these coordi-
    ization domain. Each surface component in                                                  nates, one can derive scalar fields including the conformal
    (b) is conformally mapped to a rectangle in
                                                                                                                   '§                                                                                   !
                                                                                               factor,#¤    , and mean curvature,
                                                                                                                                          , of the surface     !
                                                                                                                                                                '§           ¦98¤                            
    (c). The color scheme shows the association                                                position vector         :                              #¤                                                     
    between elements in (b) and (c).
                                                                                                                                                                       q m iq Vv…
                                                                                                                                                                           — ™ — m ¨ ŠW ŠW           ¤     ¤
                                                                                                                                                                         § ¥ §‚¥                ‚ cW cW                      (3)
                                                                                                                                                                                                    § ¦ ¥
                                                                                                                 « q
                                                                                                             ‚ ¨ ¬… — m ª
                                                                                                                                                                      « q — m aq ‚ W ®‚ W m q ¡ m
                                                                                                                     § ¥                                                  § ¥ Y cW                       —                    (4)
                                                                                                                                                                                        §         ¥
                                                                                                                                                                                                  W § ¥
2.2. Optimal Global Conformal Parameterization                                                  We can regard the tuple           as the conformal represen-                                  D8  ¤
                                                                                               tation of         . We have the following theorem [14].
                                                                                                                   ¦¤                                                              ! '§  
    Given a Riemann surface        §  with a conformal atlas                                       Theorem: A closed surface              in     with confor-        65  #¤
            , a holomorphic 1-form is defined by a fam-
                     § i §  v  q}  ~¤             “                                        mal parameter          is uniquely determined by its confor-
                                                                                                                                           !'¤§                                                                   !
ily                 , such that (1).
              i  “  v  q}   ~¤                    , where    u  #¤  t   “
                                                                    v  v                      mal factor               #¤ 
                                                                                                                    and its mean curvature             up to                                                     '§ #98
     is holomorphic on         , and (2). if
                                                                ” •  v                                                                                                                                              ! ¤
   ~                     aiy
                          œ›              t            v
                                                    ‹#¤e‹                                  rigid motions. A simply connected surface               with a                                                      "#§ ¤ 4
is the coordinate transformation on
                          ›                                          – ¤ ‹ 
                                                              —˜’eE~  ~
                                                          , then                               boundary in        and conformal parameter
                                                                                                                             6                       is deter-                                                          ¦¤
       ‹ š¤ ‹ t       , i.e., the local representation of the
          v                       Ššl™t
                                   v¤                                                       mined by its conformal factor
                                                                                                             "§                          and its mean curva-                            #¤         
differential form satisfies the chain rule.
               “                                                                               ture¦98
                                                                                                    ¤        and the boundary position.
    For a Riemann surface with genus
                                                 , all holomor-
                                                             Ÿ ž
                                                               )                                  Clearly, various fields of scalars or tuples could be used
phic 1-forms on form a complex -dimensional vector
                               r                              ¡   )                            to represent surfaces in the parameter domain. Because the
space ( real dimensions), denoted by
                             )                      . The con-
                                                           ¤ ¡                                conformal structure is intrinsic and independent of the data
formal structure of a higher genus surface can always be                                       resolution and triangulation, we use the conformal repre-
                                                                                                                                  "§                                           !
represented in terms of a holomorphic 1-form basis, which
                  r                                                                            sentation,         ¦¤ 
                                                                                                                   and          , represent the 3D surfaces.
is a set of     functions
                      )                  F ¡ ‡C Q “
                                             ¢             '§ © 5
                                                                .      2
                                                                           i§         a§
                                                                                           )   This representation is stable and computationally efficient.
Any holomorphic 1-form is a linear combination of these
                                                 “                                             Figure 4 illustrates computed conformal factor and mean
functions. The quality of a global conformal parameteriza-                                     curvature indexed by color on a hippocampal surface.
                                                                                                          Q qm ´ ³ Í                    ¡ "h ”               h§¡
                                                                                                                                                             "™"h            q ‚ ‚ m —Q                                              h
                                                                                                    £`´ ³ ` Q m ³ ) #Ì      c© »Ï ”  Ï ¡
                                                                                                                             ¤             ¤           £©      ¤ Ï
                                                                                                                                                                  ƒ7Î             Q
                                                                                                               Ò Ñ
                                                                                                               e¨r           Vr Ô © xw
                                                                                                                             Ò h               Ò Ñ
                                                                                                                                                l¨r           Ò Ô h
                                                                                                                                                              $r ¡© xw
                                                                                                        Ž           Ð
                                                                                                                    q¤   ©            ¤ ¹ Ó Ž
                                                                                                                                   © £V¥ 'e ’      z¤ e  ©
                                                                                                                                                    Ð                    ¢$¥ ae
                                                                                                                                                                        ¤ ¹ Ó
                                                                                                    as the Parzen kernel, and “ ” denotes convolution.    Õ

                                                                                                    3. The Surface Mutual Information Method
                                                                                                       for an Arbitrary Genus Surface

   Figure 4. Illustrates the computed conformal                                                         Next, suppose we want to match two high genus surfaces
   factor and mean curvature on a hippocampal                                                       (i.e., surfaces with the same branching topology). To ap-
   surface. On the left are two views of the hip-                                                   ply our surface mutual information method piecewise, we
   pocampal surface, colored according to the                                                       first compute the conformal representations of the two sur-
   conformal factor. The right two are colored                                                      faces based on a global conformal parameterization. Mu-
   by mean curvature.                                                                               tual information driven flows are then applied to align the
                                                                                                    computed conformal representations, while enforcing con-
                                                                                                    straints to guarantee continuity of the vector-valued flow at
                                                                                                    the patch boundaries. When the chain rule is used, we can
2.4. Mutual Information for Surface Registration
                                                                                                    further optimize the mutual information matching results by
    We now describe the mutual information functional used
                                                !                 '§
                                                                  !                                 optimizing the underlying global conformal parameteriza-
to drive the scalar fields       and      ¦¤ into correspon-
                                                            #98                                    tion.           ¡
dence, effectively using the equivalent of a 2D image reg-                                              Let     and    be two surfaces we want to match and the
                                                                                                                                          ©                                             ¡                              ¡
istration in the surface parameter space (i.e., in conformal                                        conformal parameterization of        is , conformal param-
                                                                                                                                                            ¡   ¡                                                          Ö
                              ¡                                                                     eterization for     is ,           and         are rectangles
coordinates). Let and be the target and the deforming
                          ¯                ©¯                                                                                                  ©                  ©
                                                                                                                                                            ¤ Ö £Ö                                              © ¤$Ö©               ¡              
template images respectively, and                      . Let
                                                            F 5 C © §
                                                            ˆ© ‡9£¯ Š¡ ¯         5                  in   © 5. Instead of finding the mapping from            to                                                      ”                                          ©
© 5 † ¡   be the common parameter domain of both surfaces                                           directly, we can use mutual information method to find a
                                                                                                                                                         ¡ ¡
(if both are rectangular domains, the target parameter do-                                          diffeomorphism                                ‰× F eRÖ
                                                                                                                                     , such that the diagram be-
                                                                                                                                                ©          × C
main can first be matched to the source parameter domain                                             low commutes:                             ”Ù ¡ x†ØŒ Ž © Ö
                                                                                                                                        . Then the map can be
                                                                                                                                                     Ö Œ Ö                                                                                ”
using a 2D diagonal matrix). Also, let be a deformation                                             obtained from the following commutative diagram,
vector field on . The MI of the scalar fields (treated as 2D
                     ¡                                                                                                                                                                Ü              ‚
images) between the two surfaces is defined by
                           q                                                                                                                                         `¤           ¤ Û
               ‚ Q y Q y q ‚ ‚ m Q — ` Q q m £m ³ ¸¶q ‚ Q — Q m fƒ‚ £… q m °
                                             ´ ·µ                   ²±                                                                                                                                       ‚
                     ` Q £³ ` Q ³ p ´                       ` ´³          ¥                   (5)                                                          `ޗ                                           Þ                                                 (8)
                                                                                                                                                                     ‚Ý                         ‚Ý
       h      w               h      ¡
                                      "h         ¡ º    ¡ "h                                                                                                     Ú                Ú Û ßÞ
where             ¥¤ © ¯ º 
           ‘ h #$fš¤ w R c© ¤ ¥ ¹uh   »#¤ #¤ t '¨"¢¹
                             ,   ¡   ¥¡ ¯ R h § ¡ ¤ h                                    
and          ©    ‘ ¦f¯
                         ¥¤ ©  and       ¥ ¯
                                         '»¦¤ š¤  .
                                                   º £©
                                                            ¤ ¥ ¹                                               ¡                                ¡                                  ¡
   We adopted the framework of D’Agostino et al. [8] to
                                                                                                     Ž © &àxŒ Ö ˆ”
                                                                                                         Ö Œ Ö                    . Because , and                                   Ö Ö                            £Ö
                                                                                                                                                                                                                   ©         are all diffeomor-
                                                                                                    phisms,    ”         is also a diffeomorphism.
maximize MI with viscous fluid regularization. Briefly,
the deforming template image was treated as embedded
in a compressible viscous fluid governed by Navier-Stokes                                            3.1. Mutual Information Contained in Maps be-
equation for conservation of momentum [6], simplified to a                                                tween High Genus Surfaces
linear PDE:
               §       ! ¿ ¿ ½ %       ! ¿ ½ ! ¼                                                       A global conformal parameterization for a high genus
         Ÿ  ‘ ¦¤ @%  e ƒ¤ ‘ 9† ¤ % © e¾
              ¥ Á       À À                                                                  (6)   surface can be obtained by integrating a holomorphic 1-
                ½                           !                                                      form . Suppose                               is a holomor-
                                                                                                                                                                   "§ Q £}
                                                                                                                                                                      h                  § r
                                                                                                                                                                                         „a§                        r
   Here is the deformation velocity, and and are the                                                            “                                                        “          2                £e
                                                                                                                                                                                                     ee                      
viscosity constants. Following the derivations in [8], we                                           phic 1-form basis, where an arbitrary holomorphic 1-form
                                                                                                                                                       Q Q  ¡ p … Q© I 
take the first variation of      with respect to , and use the
                                            ¤                                                      has the formula                      . Assuming the target
                                                                                                                                                         “                  “
Parzen window method [23] to estimate the joint probability
                                           h§¡                                            §         surface’s parameterization is fixed, the mutual information
density function (pdf)           . The driving force
                                       c©      ¤ ¥ ¹                            ‘ ¦¤ Á
                                                                                  ¥                energy between it and the source surface’s parameteriza-
that registers features in the 2D surface parameter space is                                        tion is denoted by        , which is a function of the lin-
given by                                                                                                                                                                                  Q 
                                        q‚ m ¡ … q m                                                ear combination of coefficients . The necessary condi-
       m ‚ ##q
           Êq    m ‚ — q m É m È Ç‚       Q— Q ´      —
                                                                                                    tion for the optimal holomorphic 1-form is straightforward,                                                                                  ‚
  ¥ Ž b ° ¥ Ž b ° b z° ®eQ cWÆ W d           ` Å HÃ
                                                 Ä   ¥ b                                     (7)      h '§   ⨠                             § r
                                                                                                                                              „i§              r
                                                                                                    2  Ÿ  d cW W             . If the Hessian matrix   £e
                                                                                                                                                         ee                   )                                                           g ¨ â W d W¨ W ¤
       where                "¨"h ¼ §
                            h§¡   is      the      Ë     area     of       the       param-         is positive definite, then will reach the minimum. If the              á
           domain        ©
                                   ¤ ¥ ¡                                                      %    Hessian matrix is negative definite, will be maximized.                                       á
    Our surface mutual information method depends on the                                                                                                                      9. If                  , return t. Otherwise, assign
                                                                                                                                                                                           ¨ñ ó° á ° å á
                                                                                                                                                                                          á û ÷å                                                          á
selection of holomorphic 1-form . To get an optimal sur-                                             “                                                                  to       and repeat steps through .
                                                                                                                                                                             ÷ á              ü                           ý
face mutual information matching result, we need find the
optimal holomorphic 1-form for mutual information metric.                                                                                                                  Currently, we use the following numerical scheme in
                                                                                            Q Q  ¡ p … ©Q Ù
Suppose a holomorphic function                   , our goal                                   “ aä丧 I "§ “ Q
                                                                                                 r §ããã         h                                                      step :     ü
is to find a set of coefficients                  that maxi-                                     )            2                                                             1. Compute               and
                                                                                                                                                                                                                ! Ô°
                                                                                                                                                                                                                ,  å    á
                                                                                                                                                                                                                                              "§ Q s Ô $
                                                                                                                                                                                                                                        $ Ô ° å  á
mize the mutual information energy,        . We can solve                                               ° Há
                                                                                                          å                                                             2
                                                                                                                                                                              § ã § r
                                                                                                                                                                                  ;       )
this optimization problem numerically as follows:                                                                                                                                                                         r §ããã § r        h
                                                                                                                                                                                                                          aä丄a§  "§ Q Ô !
                                                                                                                   UrUrr                                                   2. Compute                       ;         )"§             2          s 
                       y         y                                                                                                             ¨y                          3. Compute           )
                                                                                                                                                                                                    § ã r
                                                                                                                                                                                                                     h Q Ô°
                                                                                                                                                                                                                              s ˜$  ó×
                                                                                                                                                                                                                                     å á           ò
      ´ ëlq
     uêìÉ` ê ê é ‡æ y ⠗ ‡æ y â G…m      y                                                  ‚ ´ uêì ê            rUUrr                `‚ UrUr¨ r y ‡ï‰    q‚ ´ ê"ê
                                                                                                                                                           í ì
                                                                                                                                                                        Equation 9.
                 è        è           ‡çâ
                                     è æ                                                    ‚ É uêì ê                           ðq”• z‚ ¨ y             Rq‚É ê"ê
                                                                                                                                                        î í ì
    ` uêì          §        ¥                                                                                                        í
                           "§                                           y
 where        is the conformal coordinate.
   Once we compute                            , we can use
                                                                             "§ ‚æ ¨ yâ
                                                                              h è
                                                                                                                 r §ããã § r
steepest descent to optimize the resulting mutual informa-
tion. A complete description of the surface mutual informa-
tion method follows.

Algorithm 1 Surface Mutual Information Method (for sur-
faces of arbitrary genus)                                  ¡
   Input (mesh          and     ,step length , energy differ-
                                                     A                       © A                                                ò
ence threshold      ),                           ¡ ¨ñ
   Output(                    ) where minimizes the surface
                             A C
                             |óò                                   F    © A                              ò
mutual information energy.
   1. Compute global conformal parameterization of two
               aä丄a§  '„i§  q§ Q Q ¡ p … ©Q ô S
               r §ããã § r   hõ r            p
surfaces,    ) ¡          2        2            “       I
                                                       , where
                                                             “                                                                                                                Figure 5. Matching surface features in 2D pa-
  is the surface genus number of two surfaces        and      ,                                                                                            © A                rameter domains using Mutual Information.
)                    A               i¸ä¸„i§  "„i§  z§ QS
                                     r §ããã § r     h§ r       p                                                                                                              Geometric features on 3D hippocampal sur-
and                              )    are the coefficients of a
                                                 2         2     
linear combination of holomorphic function basis elements.                                                                                                                    faces (the conformal factor and mean curva-
The steps include computing the homology basis, cohomol-                                                                                                                      ture) were computed and compound scalar
ogy basis, harmonic 1-form basis and holomorphic 1-form                                                                                                                       fields (e.g., 8xconformal factor + mean cur-
basis [15].                                                                                                                                                                   vature) were mapped to a 2D square by
   2. Compute holomorphic flow segmentation of the target                                                                                                                      conformal flattening. In the 2D parameter
                          © A
              , from the global conformal parameterization,                                                                                                                   domain, data from a healthy normal sub-
© “, which conformally maps the 3D surface to a set of rect-                                                                                                                  ject (the template, leftmost column) was reg-
angles in the Euclidean plane.                                                                                                                                                istered to data from several patients with
   3. Compute 2D conformal representation for the target                                                                                                                      Alzheimer’s disease (target images, second
                                            '§                                     !
                                                                                   '§                                           !
                                                                                                                                "§                                            column). Each mapping can be used to ob-
surface,           and       ©
                          #¤     , where      is the confor-
                                                                      #¤ © 8                                             ¦¤        
mal coordinate;                                                                                                                                                               tain a reparameterization of the 3D surface of
   4. Compute holomorphic flow segmentation of the source
                                ¡                                                                                                                   ¡                         the normal subject, by convecting the origi-
                                                                                                                                                             "§               nal 3D coordinates along with the flow. The
surface,    !  , and 2D conformal representation
              '§ A ¡                                                                                                                              ¦¤              
and             ; #¤ 8                                                                                                                                                       deformed template images are shown in the
   5.      Apply the mutual information method to op-                                                                                                                         third and fourth (gridded) columns. The grids
timize the correspondence between two surfaces,                                                                                                                               show how the fluid transform expands some
    '§ ¡ 
    !                       ! ¡
                      F " "§ ¦¤ 8 §                                             q§ "§
                                                                                 p !      § "§ ©
                                                                                            !                                                         C öò i§
                                                                                                                                                           r                  highly curved features to match similar fea-
         "§ #¤ S ¤
         !               i§
                          r       p
                                  q§                      and                2  " ¦¤ © 8  ¦¤  ¤
            ¦¤ 8             2  
                      ; and compute mutual information en-                                                                                                                    tures. Importantly, there are some consis-
ergy  ’á
      ÷å    ;                        °                                                                                                                                        tent 3D geometric features that can be re-
   6. Compute derivative .                                                    ó×
                                                                              ò                                                                                               identified in the 2D parameter domain; e.g.,
   7. Update the global conformal parameterization of          ¡                                                                                         ¡                    bright areas (arrows) correspond to high cur-
source surface,        , by changing the coefficients   A                                                                                            !¤                     vature features in the head of the hippocam-
        !                                                                                                                                                                     pus.
    òñ ¤ò .
    Šu fó×
   8. Compute mutual information energy , with steps
       a"§ ø
       ú§ ù                                                                                                                                  á
4. Experimental Results                                         by the joint distribution of features lying in both surfaces.
                                                                This is a natural idea, in that it uses conformal parameteri-
    To make the results easier to illustrate, we chose to en-   zation to transform a surface matching problem into an im-
code the profile of surface features using a compound scalar     age registration problem. Whether or not this approach pro-
                 ! "§  ÿ ! "§           !
                                         '§                     vides a more relevant correspondences than those afforded
function  #98 %  ¦¤ $P  ¦‰þ
           ¤                   ¤           . We linearly nor-
malized its dynamic range to the pixel intensity range 0 to     by other criteria (minimum description length, neural nets,
255. Several examples are shown, mapping one face to the        or hand landmarking) requires careful validation for each
other face and one hippocampal surface to another. Face         application. Optimal correspondence depends more on util-
mapping is useful for face recognition. The deformable sur-     ity for a particular application than on anatomical homol-
face registration of hippocampus is important for tracking      ogy. Because different correspondence principles produce
developmental and degenerative changes in the hippocam-         different shape models, we plan to compare them in future
pus, as well as computing average shape models with ap-         for detecting group differences in brain structure, and indi-
propriate boundary correspondences. Figure 5 shows the          vidual differences in face recognition applications. . If sta-
matching fields for several pairs of hippocampal surfaces,       tistical power is increased in group comparisons, this would
establishing correspondences between distinctive features.      support the use of correspondence fields established by mu-
                   !                                            tual information on surfaces.
The velocity field in Eqn 6 was computed iteratively by
convolution of the force field with a filter kernel derived by
Bro-Nielsen and Gramkow [4]. The viscosity coefficients
   and were set to 0.9 and 6.0 respectively. The deforma-  !
tion field in the parameter domain ( ) was obtained from
by Euler integration over time, and the deformed template
image was regridded when the Jacobian determinant of the
deformation mapping at any point in  ¥       was smaller than
0.5 [6]. At each step, the joint pdf was updated and the
MI re-computed. Iterations were stopped when MI was no
longer monotonically increasing or when the number of it-
erations reached 350. The Parzen parameter was set to 10
for smoothing the joint pdf. In Figure 6, we show face sur-
faces and hippocampal surfaces to be matched, where the
face surface was built with a high resolution, real-time 3D
face acquisition system [32] and hippocampal surfaces were
built from 3D MRI scans of the brain. Specifically, (a), (b)        Figure 6. Surface matching results from our
and (e), (f) show two surfaces to be matched; (c), (d) and         method. Panels (a), (b), (c) and (d) show the
(g) and (h) show 3D vector displacement map, connecting            surfaces being matched. (a) and (b) are two
corresponding points on the two surfaces. (c) and (g) are be-      face surfaces. (e) is a normal subject’s hip-
fore and (d) and (h) after reparameterization of the source        pocampal surface and (d) is an Alzheimer’s
surface using a fluid flow in the parameter domain. These            disease patient’s hippocampal surface. We
more complex 3D vector fields store information on geo-             flow the surface from (a) to (b) and from (e)
metrical feature correspondences between the surfaces.             to (f), respectively. Panels (c), (d), (g) and (h)
                                                                   show the 3D vector displacement map, con-
                                                                   necting corresponding points on the two sur-
5. Conclusions and Future Work                                     faces, (c) and (g) before and (d) and (h) af-
                                                                   ter reparameterization of the source surface
   We extended the mutual information method to match              using a fluid flow in the parameter domain.
general surfaces. This is useful for face recognition and has      After reparameterization, a leftward shift in
numerous applications in medical imaging. Our examples             the vertical isocurves adds a larger tangen-
of matching various hippocampal surfaces are relevant for          tial component to the vector field. Even so,
mapping how degenerative diseases affect the brain, as well        the deformed grid structure remains close to
as building statistical shape models to detect the anatomical      conformal. These more complex 3D vector
effects of disease, aging, or development. The face and hip-       fields store information on geometrical fea-
pocampus are used as specific examples, but the method is           ture correspondences between the surfaces.
general and is applicable in principle to other surfaces.
   Surface-based mutual information automates the match-
ing of surfaces by computing a correspondence field guided
References                                                           [18] B. Kim, J. Boes, K. Frey, and C. Meyer. Mutual informa-
                                                                          tion for automated unwarping of rat brain autoradiographs.
 [1] S. Angenent, S. Haker, R. Kikinis, and A. Tannenbaum.                NeuroImage, 5(1):31–40, 1997.
     Nondistorting flattening maps and the 3D visualization of        [19] G. Kindlmann, D. Weinstein, A. Lee, A. Toga, and
     colon CT images. IEEE TMI, 19:665–671, 2000.                         P. Thompson. Visualization of anatomic covariance tensor
 [2] S. Angenent, S. Haker, A. Tannenbaum, and R. Kikinis.                fields. In Proc. IEEE EMBS, San Francisco, CA, 2004.
     Conformal geometry and brain flattening. MICCAI, Lecture         [20] A. Leow, C. Yu, S. Lee, S. Huang, R. Nicolson, K. Hayashi,
     Notes in Computer Science, 1679:271–78, 1999.                        H. Protas, A. Toga, and P. Thompson. Brain structural map-
 [3] M. Bakircioglu, S. Joshi, and M. Miller. Landmark matching           ping using a novel hybrid implicit/explicit framework based
     on brain surfaces via large deformation diffeomorphisms on           on the level-set method. NeuroImage, 24(3):910–27, 2005.
                                                                     [21] N. Mavridis, F. Tsalakanidou, D. Pantazis, S. Malassiotis,
     the sphere. In SPIE Medical Imaging, volume 3661, pages
                                                                          and M. Strintzis. The HISCORE face recognition applica-
     710–15, 1999.
 [4] M. Bro-Nielsen and C. Gramkow. Fast fluid registration of             tion: Affordable desktop face recognition based on a novel
     medical images. In Visualization in Biomedical Computing             3D camera. In Proc. ICAV3D, pages 157–60, 2001.
                                                                     [22] C. Meyer, J. Boes, B. Kim, P. Bland, K. Zasadny, P. Kison,
     (VBC’96)., volume 1131, pages 267–76. Springer, 1996.
                                                                          K. Koral, K. Frey, and R. Wahl. Demonstration of accu-
 [5] A. Bronstein, M. Bronstein, and R. Kimmel. Expression-
                                                                          racy and clinical versality of mutual information for auto-
     invariant 3D face recognition. In Proc. Audio & Video-based
                                                                          matic multimodality image fusion using affine and thin plate
     Biometric Person Authentication (AVBPA), Lecture Notes in
                                                                          spline warped geometric deformation. Medical Image Anal-
     Comp. Science, volume 2688, pages 62–69. Springer, 2003.
 [6] G. Christensen, R. Rabbitt, and M. Miller. Deformable                ysis, 1(3):195–206, 1997.
                                                                     [23] E. Parzen. On the estimation of probability density function.
     templates using large deformation kinematics. IEEE TIP,
                                                                          Annals of mathematical statistics, 33:1065–76, 1962.
     5(10):1435–47, 1996.                                            [24] A. Pentland. Looking at people: sensing for ubiquitous and
 [7] J. Csernansky, S. Joshi, L. Wang, J. Haller, M. Gado,
                                                                          wearable computing. IEEE TPAMI, 22(1):107–119, 2000.
     J. Miller, U. Grenander, and M. Miller. Hippocampal mor-        [25] D. Rueckert, L. Sonoda, C. Hayes, D. Hill, M. Leach, and
     phometry in schizophrenia via high dimensional brain map-            D. Hawkes. Nonrigid registration using free-form defor-
     ping. Proc. Natl. Acad. Sci., 95:11406–11, 1998.                     mations: application to breast MR images. IEEE TMI,
 [8] E. D’Agostino, F. Maes, D. Vandermeulen, and P. Suetens.             18(8):712–21, 1999.
     A viscous fluid model for multimodal non-rigid image reg-        [26] P. Thompson, J. Giedd, R. Woods, D. MacDonald, A. Evans,
     istration using mutual information. Medical Image Analysis,          and A. Toga. Growth patterns in the developing brain de-
     7(4):565–75, 2003.                                                   tected by using continuum-mechanical tensor maps. Nature,
 [9] C. Davatzikos. Spatial normalization of 3D brain images              404(6774):190–3, Mar. 2000.
     using deformable models. J. Comp. Assisted Tomography,          [27] P. Thompson, R. Woods, M. Mega, and A. Toga. Mathemat-
     20(4):656–65, 1996.                                                  ical/computational challenges in creating population-based
[10] R. Davies, C. Twining, T. Cootes, J. Waterton, and C. Taylor.        brain atlases. In Human Brain Mapping, volume 9, pages
     A minimum description length approach to statistical shape           81–92, Feb. 2000.
     modeling. In IEEE TMI, volume 21, pages 525–37, 2002.           [28] C. Twining, T. Cootes, S. Marsland, V. Petrovic,
[11] B. Fischl, M. Sereno, R. Tootell, and A. Dale. High-                 R. Schestowitz, and C. Taylor. A unified inforamtion-
     resolution inter-subject averaging and a coordinate system           theoretic approach to groupwise non-rigid registration and
     for the cortical surface. In Human Brain Mapping, volume 8,          model building. In IPMI, 2005.
     pages 272–84, 1999.                                             [29] P. Viola and W. M. Wells. Alignment by maximization of
[12] X. Gu and B. Vemuri. Matching 3D shapes using 2D con-                mutual information. In Proc. ICCV, pages 16–23, Washing-
     formal representations. In MICCAI, Lectures Notes in Com-            ton, DC, USA, 1995.
     puter Science, volume 3217, pages 771–80. Springer, 2004.       [30] W. Wells, P. Viola, H. Atsumi, S. Nakajima, and R. Kikinis.
[13] X. Gu, Y. Wang, T. Chan, P. Thompson, and S.-T. Yau.                 Multi-modal volume registration by maximization of mutual
     Genus zero surface conformal mapping and its application             information. Medical Image Analysis, 1(1):35–51, 1996.
     to brain surface mapping. IEEE TMI, 23(8):949–58, Aug.          [31] J. West, J. Fitzpatrick, M. Wang, B. Dawant, C. Maurer,
     2004.                                                                R. Kessler, R. Maciunas, C. Barillot, D. Lemoine, A. Col-
[14] X. Gu, Y. Wang, and S.-T. Yau. Geometric compression                 lignon, F. Maes, P. Suetens, D. Vandermeulen, P. van den
     using Riemann surface. Communication of Information and              Elsen, S. Napel, T. Sumanaweera, B. Harkness, P. Hemler,
     Systems, 3(3):171–82, 2005.                                          D. Hill, D. Hawkes, C. Studholme, J. Maintz, M. Viergever,
[15] X. Gu and S.-T. Yau. Computing conformal structures                  G. Malandain, X. Pennec, M. Noz, G. Maguire, M. Pollack,
     of surfaces. Communication of Information and Systems,               C. Pelizzari, R. Robb, D. Hanson, and R. Woods. Compar-
     2(2):121–46, 2002.                                                   ison and evaluation of retrospective intermodality brain im-
[16] G. Hermosillo. PhD thesis, Universit de Nice (INRIA-


                                                                          age registration techniques. J. Comp. Assisted Tomography,
     ROBOTVIS), Sophia Antipolis, France, 2002.                           21(4):554–68, 1997.
[17] J. Huang, V. Blanz, and B. Heisele. Face recognition using      [32] S. Zhang and P. Huang. High-resolution, real-time 3D shape
     component-based SVM classification and morphable mod-                 acquisition. In IEEE CVPR Workshop, volume 3, page 28,
     els. In S.-W. Lee and A. Verri, editors, SVM, LNCS 2388,             Washington DC, USA, 2004.
     pages 334–41. Springer-Verlag, 2002.