intensities in the kth bin, $n_k$ is the number of values in $Y_k$ with $N = \sum_k n_k$. This cost is optimized to find a 12-parameter affine transformation.

In addition, we compute the inverse transformation to take the subject image back to its native space at the end of the segmentation.

2. We find the location or boundary of the skin with the skin organism. Its initial shape is of the skin of our template image found using the marching cubes method. We dilated this substantially to ensure we encompass the head of any subject registered to the template. Our template-fitted mesh needs to be further refined to fit our subject. To do this, we analyze the subject's intensities and apply a threshold to create a binary image masking out the head. We also use behavior 1 to sense the threshold image and evolve our skin organism's mesh to find this perimeter. We iterate the deformations dictated by behavior 1 (applying smoothing at every step) until there is no significant movement of the surface or we reach a maximum iteration bound. We specify this bound based on images being reasonably aligned to the template, an approach we used for all our deformations. The adjacent eyes are handled in a similar manner.

3. Our eye organisms find the eye boundary by sensing the 3-means classified image. We initialize the eye organisms' meshes by fitting them to our template and eroding them to make sure they lie within a subject image's eyeballs. The eyeball locations are found by the organisms sensing the 3-means image with behavior 2, which chooses a label found in the eyeball region.

4. Our next step locates the cerebrospinal fluid (CSF) that surrounds the brain. We achieve this using our skin organism by contracting its mesh into the head through the skin and skull. The skin and skull locations are roughly classified in our 2-means image and we apply behavior 3 to sense it and find the CSF boundary. To further refine this boundary the behavior also makes the skin organism deform through the eye organisms because there is more information about the CSF location by the eyes.

5. This step finds the tissues surrounding the eyeballs that need to be excluded from the brain delineation. We attain this goal by expanding the eye organisms further by sensing the 2-means image along with behavior 2 again, this time behavior 2 looks for a different label in the classified image, one that gives an understanding of tissues around the eyes. The eyes now furnish a better understanding to locate the brain.

6. We complete our plan by finding the brain using our brain organism. Every step in the plan supports of this step and all of our knowledge up to this point will come into play. Our brain organism begins by sensing the 3-means image and gradient image with behavior 4. The behavior is cognizant of the other organisms' locations and uses them to constrain its evolution. With the completion of behavior 4 we further refine the boundary by sensing the 3-means image again with behavior 5, which results in the brain being encapsulated by our brain mesh. The mesh is then converted by the 3D rasterization scheme to a binary volume to which we apply our inverse
transform from FLIRT to bring the subject's delineation back into its native space to complete the segmentation.

<table><caption>Table 1: Plan for each organism to skull-strip an MRI image.</caption><thead><tr><th></th><th>Skin Organism</th><th>Eye Organism</th><th>Brain Organism</th></tr></thead><tbody><tr><th>Initialization</th><td></td><td></td><td></td></tr><tr><th>Deformation 1</th><td></td><td></td><td></td></tr><tr><th>Deformation 2</th><td></td><td></td><td></td></tr></tbody></table>

We make a binary volume of the brain organism and apply the inverse transformation back to the subject image space, completing the segmentation.

## 3 Error Learning

We are able to learn the types of errors our method makes, by comparing the masks it generated with expert manual delineations. [33] introduced an algorithm using Adaboost [3] to learn a weighting of a set of features used to classify if a voxel has been correctly classified by a prior algorithm. This ‘Adaboost wrapper’ algorithm uses a set of corresponding automated and manual segmentations, as well as intensity images to find features that lie in regions that the first-pass method incorrectly classifies. We use this algorithm to learn situations in which our method makes errors and thereby improve the segmentation.

## 4 Validation

The masks from the deformable organism method are compared with the gold standard manual segmentations using standard distance, overlap, and error met-
rics. We use the Hausdorff distance measure [7] to find the distance from the furthest point in the deformable organism method mask to the closest point of the mask in the manual delineation. We also compare the expert and automated masks using the Jaccard coefficient (union overlap), Dice coefficient (mean overlap), false negative rate, and false positive rate described in [11].

We compared the differences in metric values across methods using paired-sample *t*-tests to understand if their results were statistically different.

# 5 Experiments

We tested our deformable organism skull stripping method on our set of 630 manually-labeled subject images. In addition, we also ran BET (BET2, FSL 4.1.5, default parameters), BSE (BSE 10a, default parameters), and the Watershed algorithm (Freesurfer 5.0.0, default parameters), and assessed their errors using standard distance-based, overlap, and error metrics.

Typically patterns of error in our method were learned by selecting a subset of our segmentation results and using the error learning algorithm (the 'Adaboost wrapper' approach). We then segmented a new subset of images with the error classified and corrected. We repeated this experiment 10 times, using 10 random images from our results to train and 10 random (but non-overlapping) images to test. Masks were then compared to expert ground truth before and after error correction; note that the test set of images was independent of those used for training the error correction step.

# 6 Results

## 6.1 Subject Data

Our subject data consisted of 630 T1-weighted magnetic resonance (MR) images from healthy young adults, between 20 and 30 years of age. These images are from Australian twins, and have been used in numerous prior analyses [10].

Each of the images had been manually skull-stripped by a neuroanatomically trained expert. These manual labels were used as the gold standard to compare with automatic segmentation results of our method and the other 3 widely-used methods. The subjects were scanned with a 4-Tesla Bruker Medspec whole-body scanner. 3D T1-weighted images were acquired using a magnetization-prepared rapid gradient echo sequence, resolving the anatomy at high resolution. Acquisition parameters were: inversion time ($T_I$)/repetition time ($T_R$)/echo time ($T_E$)=700/1500/3.35 ms, flip angle=8°, slice thickness=0.9 mm with a 256 × 256 × 256 acquisition matrix.

In addition, we used one of the 40 images from the LONI Probabilistic Brain Atlas (LPBA40) [26]. Each image had 56 different structures manually labeled, including a mask of the brain.

Table 2 shows the distance, overlap, and error metrics for the automated skull-stripping algorithms compared to manual segmentation. We compare BET,
BSE, and the Watershed method. Average metrics over the 630 subject images are shown. Our paired-sample *t*-tests showed the metric values across all the methods were significantly different. These test results showed that the deformable organisms approach was statistically better than the others in the Hausdorff distance and false negative error for our dataset.

The deformable organisms method took a few minutes to run on the subject images we used on a machine with dual 64-bit 2.4 gigahertz AMD Opteron 250 CPU with 8 gigabytes of memory.

We list average results of deformable organisms with and without error correction versus manual training in Table 3. Random samples of 20 images from the 630 were selected, using 10 to train and 10 to test the error correction. We repeat this 10 times and average the results. These average results were used for two-sample *t*-test that found statistically significant improvement in Jaccard, Dice, and false positive error metrics.

<table><caption>Table 2: Distance, overlap, and error metrics comparing automated results with manual skull-stripping.</caption><thead><tr><th></th><th>Hausdorff Distance</th><th>Jaccard Overlap</th><th>Dice Overlap</th><th>False Negative Error</th><th>False Positive Error</th></tr></thead><tbody><tr><td>Deformable Organisms</td><td>36.3475±24.3842</td><td>0.8478±0.0242</td><td>0.9175±0.0143</td><td>0.0253±0.0115</td><td>0.1328±0.0280</td></tr><tr><td>BET</td><td>41.7997±24.7972</td><td>0.8860±0.0183</td><td>0.9395±0.0104</td><td>0.0711±0.0218</td><td>0.0491±0.0189</td></tr><tr><td>BSE</td><td>64.7451±25.7714</td><td>0.8348±0.1295</td><td>0.9040±0.0842</td><td>0.1045±0.1456</td><td>0.0720±0.0323</td></tr><tr><td>Watershed</td><td>67.4672±8.7859</td><td>0.3650±0.0599</td><td>0.5321±0.0617</td><td>0.4511±0.0613</td><td>0.4831±0.0648</td></tr></tbody></table>

<table><caption>Table 3: Distance, overlap, and error metrics comparing the deformable organisms segmentation with and without error correction versus manual skull-stripping.</caption><thead><tr><th></th><th>Hausdorff Distance</th><th>Jaccard Overlap</th><th>Dice Overlap</th><th>False Negative Error</th><th>False Positive Error</th></tr></thead><tbody><tr><td>Basic DO</td><td>35.4780±2.7890</td><td>0.8485±0.0009</td><td>0.9178±0.0005</td><td>0.0256±0.0006</td><td>0.1318±0.0010</td></tr><tr><td>DO with correction</td><td>35.5194±3.1299</td><td>0.8858±0.0013</td><td>0.9393±0.0007</td><td>0.0253±0.0004</td><td>0.0929±0.0012</td></tr></tbody></table>

# 7 Discussion

The metrics in Table 2 suggest that the performance of our deformable organisms approach is comparable to that of other widely used methods. It has the lowest Hausdorff distance average between its automatic results and the gold standard delineations. It did have a higher false positive error, meaning it may include slightly more voxels inside the boundary of the brain.

Table 3 shows that learning errors improves the segmentation results; all metrics examined were improved, especially the Jaccard coefficient. This additional training step may be useful if a large data set needs to be segmented, making it reasonable to segment some images manually for error correction. The method could be trained on a small subset of the manual and automatically segmented data, in a first pass, and the error corrections learned could be useful to segment the rest of the dataset.
Deformable organisms provide an adaptable framework to perform segmen-
tation. They can encode a high-level plan into deformable models, to help them
work together to accomplish segmentation tasks. The different control layers
may be adapted to fit any type of segmentation tasks.

References

1. Boesen, K., Rehm, K., Schaper, K., Stoltzner, S., Woods, R., Luders, E., Rottenberg, D.: Quantitative comparison of four brain extraction algorithms. NeuroImage 22(3), 1255–1261 (2004)

2. Desbrun, M., Meyer, M., Schroder, P., Barr, A.: Implicit fairing of irregular meshes using diffusion and curvature flow. In: Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques. pp. 317–324. ACM Press/Addison-Wesley Publishing Co. (1999)

3. Freund, Y., Schapire, R.: A desicion-theoretic generalization of on-line learning and an application to boosting. In: Computational Learning Theory. pp. 23–37. Springer (1995)

4. Gering, D., Nabavi, A., Kikinis, R., Hata, N., ODonnell, L., Grimson, W., Jolesz, F., Black, P., Wells III, W.: An integrated visualization system for surgical planning and guidance using image fusion and an open MR. Journal of Magnetic Resonance Imaging 13, 967–975 (2001)

5. Gharachorloo, N., Gupta, S., Sproull, R., Sutherland, I.: A characterization of ten rasterization techniques. ACM SIGGRAPH Computer Graphics 23(3), 355–368 (1989)

6. Hamarneh, G., McIntosh, C.: Physics-based deformable organisms for medical image analysis. In: Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. vol. 5747, pp. 326–335 (2005)

7. Huttenlocher, D., Klanderman, G., Rucklidge, W.: Comparing images using the Hausdorff distance. IEEE Transactions on Pattern Analysis and Machine Intelligence pp. 850–863 (1993)

8. Ibanez, L., Schroeder, W., Ng, L., Cates, J., et al.: The ITK software guide. Kitware Inc (2005)

9. Jenkinson, M., Smith, S.: A global optimisation method for robust affine registration of brain images. Medical Image Analysis 5(2), 143–156 (2001)

10. Joshi, A., Lepore, N., Joshi, S., Lee, A., Barysheva, M., Stein, J., McMahon, K., Johnson, K., de Zubicaray, G., Martin, N., et al.: The contribution of genes to cortical thickness and volume. NeuroReport 22(3), 101 (2011)

11. Klein, A., Andersson, J., Ardekani, B., Ashburner, J., Avants, B., Chiang, M., Christensen, G., Collins, D., Gee, J., Hellier, P., et al.: Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration. NeuroImage 46(3), 786–802 (2009)

12. Lorensen, W., Cline, H.: Marching cubes: A high resolution 3D surface construction algorithm. ACM Siggraph Computer Graphics 21(4), 163–169 (1987)

13. McInerney, T., Hamarneh, G., Shenton, M., Terzopoulos, D.: Deformable organisms for automatic medical image analysis. Medical Image Analysis 6(3), 251–266 (2002)

14. McInerney, T., Terzopoulos, D.: Deformable models in medical image analysis: a survey. Medical Image Analysis 1(2), 91–108 (1996)
15. McIntosh, C., Hamarneh, G.: I-DO: A Deformable Organisms framework for ITK. Medical Image Analysis (2006)

16. McIntosh, C., Hamarneh, G.: Spinal crawlers: Deformable organisms for spinal cord segmentation and analysis. Medical Image Computing and Computer-Assisted Intervention–MICCAI 2006 pp. 808–815 (2006)

17. McIntosh, C., Hamarneh, G.: Vessel crawlers: 3D physically-based deformable organisms for vasculature segmentation and analysis. In: Computer Vision and Pattern Recognition, 2006 IEEE Computer Society Conference on. vol. 1, pp. 1084–1091 (2006)

18. McIntosh, C., Hamarneh, G., Mori, G.: Human limb delineation and joint position recovery using localized boundary models. In: Motion and Video Computing, 2007. WMVC’07. IEEE Workshop on. pp. 31–31. IEEE (2007)

19. Ohtake, Y., Belyaev, A., Bogaevski, I.: Polyhedral surface smoothing with simultaneous mesh regularization. In: Geometric Modeling and Processing 2000. Theory and Applications. Proceedings. pp. 229–237. IEEE (2000)

20. Pineda, J.: A parallel algorithm for polygon rasterization. ACM SIGGRAPH Computer Graphics 22(4), 17–20 (1988)

21. Prasad, G., Joshi, A., Thompson, P., Toga, A., Terzopoulos, D., Shattuck, D.: Skull-stripping with deformable organisms. In: International Symposium on Biomedical Imaging, ISBI 2011. pp. 1662–1665. IEEE (2011)

22. Resnick, S., Pham, D., Kraut, M., Zonderman, A., Davatzikos, C.: Longitudinal magnetic resonance imaging studies of older adults: a shrinking brain. Journal of Neuroscience 23(8), 3295–3301 (2003)

23. Roche, A., Malandain, G., Pennec, X., Ayache, N.: The correlation ratio as a new similarity measure for multimodal image registration. Medical Image Computing and Computer-Assisted Intervention–MICCAI 1998 p. 1115 (1998)

24. Segonne, F., Dale, A., Busa, E., Glessner, M., Salat, D., Hahn, H., Fischl, B.: A hybrid approach to the skull stripping problem in MRI. NeuroImage 22(3), 1060–1075 (2004)

25. Shattuck, D., Leahy, R.: BrainSuite: an automated cortical surface identification tool. Medical Image Analysis 6(2), 129–142 (2002)

26. Shattuck, D., Mirza, M., Adisetiyo, V., Hojatkashani, C., Salamon, G., Narr, K., Poldrack, R., Bilder, R., Toga, A.: Construction of a 3D probabilistic atlas of human cortical structures. NeuroImage 39(3), 1064–1080 (2008)

27. Shattuck, D., Prasad, G., Mirza, M., Narr, K., Toga, A.: Online resource for validation of brain segmentation methods. NeuroImage 45(2), 431–439 (2009)

28. Smith, S.: Fast robust automated brain extraction. Human Brain Mapping 17(3), 143–155 (2002)

29. Steels, L.: The artificial life roots of artificial intelligence. Artificial Life 1(1.2), 75–110 (1993)

30. Terzopoulos, D., Platt, J., Barr, A., Fleischer, K.: Elastically deformable models. ACM Siggraph Computer Graphics 21(4), 205–214 (1987)

31. Thompson, P., Hayashi, K., De Zubicaray, G., Janke, A., Rose, S., Semple, J., Herman, D., Hong, M., Dittmer, S., Doddrell, D., et al.: Dynamics of gray matter loss in Alzheimer’s disease. Journal of Neuroscience 23(3), 994–1005 (2003)

32. Thompson, P., Hayashi, K., Sowell, E., Gogtay, N., Giedd, J., Rapoport, J., De Zubicaray, G., Janke, A., Rose, S., Semple, J., et al.: Mapping cortical change in Alzheimer’s disease, brain development, and schizophrenia. NeuroImage 23, S2–S18 (2004)
33. Wang, H., Das, S., Suh, J., Altinay, M., Pluta, J., Craige, C., Avants, B., Yushkevich, P.: A Learning-Based Wrapper Method to Correct Systematic Errors in Automatic Image Segmentation: Consistently Improved Performance in Hippocampus, Cortex and Brain Segmentation. NeuroImage (2011)

34. Woods, R., Dapretto, M., Sicotte, N., Toga, A., Mazziotta, J.: Creation and use of a Talairach-compatible atlas for accurate, automated, nonlinear intersubject registration, and analysis of functional imaging data. Human Brain Mapping 8(2-3), 73-79 (1999)
# Selecting Scales by Multiple Kernel Learning for Shape Diffusion Analysis

Umberto Castellani<fnref n="1" />*, Aydın Ulaş<fnref n="1" />, and Vittorio Murino<fnref n="1" />,<fnref n="2" />,
Marcella Bellani<fnref n="3" />, Gianluca Rambaldelli<fnref n="3" />, Michele Tansella<fnref n="3" />, Paolo
Brambilla<fnref n="3" />,<fnref n="4" />

**Abstract.** Brain morphological abnormalities can typically be detected by advanced geometrical shape analysis techniques. Recently, shape diffusion methods have proved to be very effective in providing useful descriptions for brain classification purposes. In particular, they allow the analysis of such shapes at multiple scales, but the selection of the correct range of scales remains an open issue heavily affecting the performance of methods, and it needs to be estimated adaptively for different classes of shapes. In this paper, we focus on the diffusion scale selection in order to define a robust shape descriptor for brain classification. To this end, geometric features are extracted for each scale and the best feature combination is selected by employing *multiple kernel learning* (MKL). In the presented experiments, we compare the shape of Thalamic regions in order to discriminate between normal subjects and schizophrenic patients. We demonstrate that MKL allows to obtain classifiers which are more accurate with respect to other competing algorithms for schizophrenia detection. Moreover, using the weights computed by the MKL algorithm, we can select at which scale the features are more effective for schizophrenia classification.

**Keywords:** multiple kernel learning, schizophrenia, heat kernel, spectral shape analysis, support vector machines

## 1 Introduction

Recent advances in geometric shape analysis have led to a larger diffusion of computational anatomy methods, aimed at characterizing or modeling the morphological variations of biological shapes. One of the typical applications is analyzing the anatomy of organs, known as being possibly affected by abnormalities due to a certain disease, of several persons in order to discriminate between normal and pathological subjects [11,1]. To this aim, effective shape analysis

<footnote n="1">¹ Department of Computer Science, University of Verona, Italy</footnote>

<footnote n="2">² Istituto Italiano di Tecnologia, Genova, Italy</footnote>

<footnote n="3">³ Department of Medicine and Public Health, University of Verona, Italy</footnote>

<footnote n="4">⁴ ICBN Center, University of Udine and Verona, Italy</footnote>

<footnote>* Corresponding author.</footnote>
techniques are crucial to extract geometric features with high discriminant prop-
erties. A wide class of methods are based on the encoding of the deformation
which aligns a pair of subjects, but such an approach requires the solution of
a complex problem due to the non-linear registration between different shapes.
More recently, new methods have been proposed to encode the shape geometric
properties into a *descriptor* which compactly represents the whole shape, and
performing the comparison by computing similarities in the descriptor space
without any registration procedure. Among the shape analysis methods, dif-
fusion geometry approaches are very promising since they are able to capture
*intrinsic* characteristic of the shape. More specifically local geometric proper-
ties are encoded by the so-called *Heat Kernel* [16] which exploits heat diffusion
characteristics at a given scale. The general idea consists of gaining information
about the neighborhood of a point on the shape by recording the dissipation of
heat over time from that point onto the rest of the shape. The fixed time is very
important since it allows to capture different kinds of information: *local* shape
characteristics are highlighted through the behavior of heat diffusion over short
time periods, and, conversely, *global* shape properties are observed while con-
sidering longer periods [16,10]. So doing, simply varying a single parameter (the
time), it is possible to characterize the properties of a shape at different scales.
In particular, the so called *Heat Kernel Signature (HKS)* [16] has been proposed
to encode simultaneously the contribution of local features for a fixed set of
scales into a single shape descriptor. This general approach has been successfully
applied for object retrieval [4] and brain classification [6]. However, the choice
of the range of the time periods to be evaluated (i.e., the *scales*) is critical and
depends on the considered shape. In fact, for a particular shape, some scales may
be highly discriminative, while some other scales should encode useless informa-
tion. In this paper, we propose a new approach for integrating and selecting the
contribution of geometric features collected at different scales by utilizing a Mul-
tiple Kernel Learning (MKL) approach. In general, MKL algorithms can learn
a weighted combination of different kernel functions able to exploit information
coming from multiple sources. In our case, the different sources are the features
extracted at different scales. Therefore, several kernels are computed (i.e., one
kernel per scale) and a set of weights are estimated for the kernel combination.
In this fashion, we can choose the most discriminative scales by selecting those
associated to the highest weights, and viceversa. Moreover, kernel combination
leads to a new similarity measure which increases the classification accuracy. It
is important to note that in our approach we aim at selecting the best shape
characteristics for classification purposes, hence, our selection is driven by the
performance of a Support Vector Machine (SVM) classifier. We have applied
our method for brain classification in schizophrenic subjects: we have adopted
a Region of Interest (ROI)-based method by analysing the shape of the Thala-
mic region, employing a *volumetric*-heat kernel computed for each voxel of the
MRI scan at different scales, as described in our previous work [6]. This paper
improves [6] for both methodological aspects, by proposing the automatic scale
selection procedure and promising results. The rest of the paper is organized as
follows. In Section 2, the basics on shape diffusion procedures are reported. Section 3 describes the Multiple Kernel Lerning strategy, and the proposed method is reported in Section 4. Results are shown in Section 5 and conclusions are finally drawn in Section 6.

## 2 Shape analysis by heat diffusion

Considering a shape $M$ as a compact Riemannian manifold [5], the heat diffusion on shape<fnref n="5" /> is defined by the heat equation:

$$ (\Delta_M + \frac{\partial}{\partial t})u(t, \mathbf{m}) = 0; \qquad (1) $$

where $u$ is the distribution of heat on the surface, $\mathbf{m} \in M$, $\Delta_M$ is the Laplace-Beltrami operator which, for compact spaces, has discrete eigendecomposition of the form $\Delta_M = \lambda_i \phi_i$. In this way, the heat kernel has the following eigendecomposition:

$$ h_t(\mathbf{m}, \mathbf{m}') = \sum_{i=0}^{\infty} e^{-\lambda_i t} \phi_i(\mathbf{m}) \phi_i(\mathbf{m}'), \qquad (2) $$

where $\lambda_i$ and $\phi_i$ are the $i^{th}$ eigenvalue and the $i^{th}$ eigenfunction of the Laplace-Beltrami operator, respectively. The heat kernel $h_t(\mathbf{m}, \mathbf{m}')$ is the solution of the heat equation with initial point heat source in $\mathbf{m}$ at time $t = 0$, and heat value in ending point $\mathbf{m}' \in M$ after time $t$. The heat kernel is isometric invariant, it is informative, and stable [16].

In the case of volumetric representations, the volume is sampled by a regular Cartesian grid composed by voxels, which allows the use of standard Laplacian in $R^3$ as the Laplace-Beltrami operator. We use finite differences to evaluate the second derivative in each direction of the volume. The heat kernel on volumes is invariant to volume isometries, in which shortest paths between points inside the shape do not change. Note that in real applications exact volume isometries are limited to the set of rigid transformations [15], however, also non-rigid deformations can faithfully be modelled as approximated volume isometries in practice. It is also worth noting that, as observed in [16,15], for small $t$ the autodiffusion heat kernel $h_t(\mathbf{m}, \mathbf{m})$ of a point $\mathbf{m}$ with itself is directly related to the scalar curvature $s(\mathbf{m})$ [15]. More formally:

$$ h_t(\mathbf{m}, \mathbf{m}) = (4\pi t)^{-3/2}(1 + \frac{1}{6}s(\mathbf{m})). \qquad (3) $$

In practice, Equation 3 states that the heat tends to diffuse slower at points with positive curvature, and viceversa. This gives an intuitive explanation about the geometric properties of $h_t(\mathbf{m}, \mathbf{m})$, and suggests the idea of using it to build a shape descriptor [16].

<footnote n="5">⁵ In this section, we borrow the notation from [16,5]</footnote>
## 3 Multiple Kernel Learning

The main idea behind kernel methods [17] is to transform the input feature space to another space (eventually with a larger dimension) where the classes are linearly separable. In particular, by employing the SVM classifier, the discriminant function after the training phase becomes $f(x) = \langle w, \Phi(x) \rangle + b$, where $w$ and $b$ are the parameters of the hyperplane which separates two classes, and $\Phi(\cdot)$ is the mapping function. Using the dual formulation and the kernel trick, one does not have to define this mapping function explicitly and the discriminant function can be written as

$$f(\mathbf{x}) = \sum_{i=1}^{N} \alpha_{i} y_{i} k(\mathbf{x}_{i}, \mathbf{x}) + b \quad (4)$$

where $k(\mathbf{x}_i, \mathbf{x}_j) = \langle \Phi(\mathbf{x}_i), \Phi(\mathbf{x}_j) \rangle$ is the kernel function that calculates a similarity metric between data instances.

More recently, MKL methods have been proposed [3,13] for learning a combination $k_\eta$ of several kernels:

$$k_{\eta}(\mathbf{x}_i, \mathbf{x}_j; \eta) = f_{\eta}(\{k_m(\mathbf{x}_i^m, \mathbf{x}_j^m)\}_{m=1}^P; \eta) \quad (5)$$

where the combination function $f_\eta$ forms a single kernel from $P$ base kernels using the parameters vector (i.e., weights) $\eta$. Such new kernel must be a valid kernel$^6$ [9] and can be plugged in Equation 4 for classification purposes. Different kernel functions correspond to different notions of similarity and instead of searching which works best, the MKL method does the picking for us, or may use a combination of kernels. MKL also allows us to combine different representations, possibly coming from different sources or modalities.

There is significant work on the theory and application of MKL, and most of the proposed algorithms differ among them by the optimization method employed to estimate the weights and by the used combination rule [3,13,14]. In this paper we focus on linear-MKL methods [3,13], whose general formulation is defined as:

$$k_{\eta}(\mathbf{x}_i, \mathbf{x}_j; \eta) = \sum_{m=1}^{P} \eta_m k_m(\mathbf{x}_i^m, \mathbf{x}_j^m) \quad (6)$$

with $\eta_m \in \mathbb{R}$. As a simplest combination approach, the so called fixed rules [9] use the combination function in Eq. (6) with all weights equally set to $\eta_m = 1$. Similarly, the mean-rule takes the mean of the kernels by setting all $\eta_m = 1/P$. Indeed, in the most general case the weights $\eta_m$ are automatically estimated by a learning by example approach. More specifically, MKL methods search for a combination of kernels that maximizes a generalized performance measure (i.e, maximum margin classification errors [9]). To this aim, in the training phase, both MKL weights and SVM parameters are simultaneously estimated within the same optimization problem.

<footnote>$^6$ The validity of the kernel depends by the combination function.</footnote>
# 4 The Proposed Method

The proposed method can be summarized in the following main steps:

1. MRI data collection.

2. Feature extraction at multiple scales.

3. Learning weights and classifier by MKL.

4. Scale selection and performance evaluation.

<figure><img src="image_56.png" /><figcaption>Fig. 1. General scheme of the proposed method.</figcaption></figure>

**MRI data collection.** In order to employ a learning-by-example approach, we need a collection of samples for both healthy subjects and patients. Source data are MRI scans where shape information can be provided in terms of volumetric data.

**Feature extraction at multiple scales.** According to the shape diffusion analysis described in Section 2, for each subject geometric features are extracted at multiple scales: a set of time values $(t_1, t_2, \dots, t_n)$ are defined, and the autodiffusion value is computed for each voxel $m$, leading to:

$$H_{t_i}(M) = \{h_{t_i}(m, m), \forall m \in M\}.$$

Then, such values are accumulated into a histogram $r_i = hist(H_{t_i}(M))$. In this manner, we obtain a set of $n$ sources of shape representation $\{r_1, \dots, r_n\}$, each one encoding the global shape at a certain scale. The number of bins for each histogram is chosen as 100.
**Learning weights and classifier by MKL.** The contribution of geometric features extracted at each scale are combined by employing the MKL strategy as described in 3. Each shape representation $r_i$ is associated to a kernel $k_m$ by leading to $n = P$ kernels. Indeed, both the weights ($\eta_1 \cdots \eta_P$) and the SVM parameters are estimated. In order to obtain the best classification accuracy according to the *max-margin* paradigm an *alternating* approach is used between the optimization of kernel weights and the optimization of the SVM classifier. In each step, given the current solution of kernel weights, MKL solves a standard SVM optimization problem with the combined kernel. Then, a specific procedure is applied to update the kernel weights.

**Scale selection and performance evaluation.** Once the MKL procedure is completed, we obtain a two-fold advantage: i) we can select the best scale contributions by keeping only the scales associated to the highest weighs, and ii) we can compose a new kernel from the weighted contributions of the best scales, which can be evaluated for classification purposes.

# 5 Experiments

This section is organized in the following parts: i) data gathering, ii) experimental methodology, iii) results, and iv) discussion.

## 5.1 Data Gathering

Quantitative data collection and processing in MRI-based research implies to face several methodological issues to minimize biases and distortions. The standard approach to deal with these issues is following well-established guidelines dictated by international organizations, such as the World Health Organization (WHO), or codified by respected institutions, such as leading universities. All patients received a diagnosis of schizophrenia according to the criteria of the Diagnostic and Statistical Manual of Mental Disorders [2]. In this work, we employ a ROI-based approach [11], so only a well defined brain subpart has been considered in this study. More specifically, we focus our analysis on the left-Thalamus whose abnormal activity has been already investigated in schizophrenia[8]. ROIs have been manually traced by experts, according to well defined medical protocols. The data set used in this work is composed by MRI brain scans of 30 patients affected by schizophrenia and 30 healthy control subjects.

## 5.2 Experimental protocol

In our experiments, we apply leave-one-out (LOO) cross-validation to assess the performance of the technique. Since LOO is used as the cross validation technique, we do not report standard deviations or variances. We compare our results using *k*-fold paired *t*-test at *p* = 0.05. We collect geometric features at 11 scales generating different shape representations $r_{01}, \cdots, r_{11}$. In practice, each
representation $r_i$ is a feature vector $x_i$ which is plugged in the MKL framework. We employ the dot product as basic kernel function (i.e., linear kernel) since it avoids the estimation of free kernel parameters. Different strategies to combine the different shape representations have also been evaluated:

- Single Best Kernel (Single-best): an SVM is trained separately per each representation. Therefore, the performance of the classification are evaluated separately at each scale. So doing, we can evaluate the independent contributions coming from the different sources of information and select the best one.

- Feature concatenation (SVM-con): the contributions coming from the different sources are concatenated into a single feature vector. Then, a single SVM is employed for classification $^7$.

- Rule-based MKL (RBMKL): as baseline MKL approach, the so called rule-based method is evaluated: the kernels computed at each scale are combined by simply taking their average (i.e., $\forall m, \eta_m = 1/P$).

- Simple MKL (SimpleMKL): a simple but effective MKL algorithm is employed [14] by addressing the MKL problem through a weighted 2-norm regularization formulation with additional constraint on the weights that encourages sparse kernel combination. It is a popular approach and its code is publicly available$^8$.

- Group Lasso MKL (GLMKL): it denotes the group Lasso-based MKL algorithms proposed by [12,18]. A closed form solution for optimizing the kernel weights based on the equivalence between group-lasso and MKL is proposed. In our implementation, we used $l_1$-norm on the kernel weights and learned a convex combination of the kernels.

## 5.3 Results

The first evaluation scores are shown in Table 1, which reports the single-best kernel accuracies for all feature representations. We can observe that the best performance is obtained at 78.33 % using r02 which is shown as bold face in the table. The entries marked with "*" show the accuracies which are statistically significantly less accurate than the best algorithm using k-fold paired t-test at $p = 0.05$.

<table><caption>**Table 1.** Single-kernel SVM accuracies.</caption><tr><td>r01</td><td>r02</td><td>r03</td><td>r04</td><td>r05</td><td>r06</td><td>r07</td><td>r08</td><td>r09</td><td>r10</td><td>r11</td></tr><tr><td>75.00</td><td><b>78.33</b></td><td>76.67</td><td>76.67</td><td>73.33</td><td>*66.67</td><td>68.33</td><td>70.00</td><td>76.67</td><td>71.67</td><td>70.00</td></tr></table>

Second, concatenating the features in a single vector leads to 83.33 % accuracy.

<footnote>$^7$ We use LIBSVM software [7] to train the SVM.</footnote>

<footnote>$^8$ http://asi.insa-rouen.fr</footnote>
Third, using the proposed three different MKL algorithms, we combined the eleven kernels by introducing the weights $\eta_m$. Table 2 reports the results of the best single-kernel SVM, the accuracy of the concatenated feature set, and the three MKL-based algorithms trained. The values in parentheses show the percentage of controls classified as schizophrenia and the percentage of patients classified as healthy respectively. We achieve an accuracy of 86.67%, reached by combining eleven kernels with the SimpleMKL approach. This result is better than all other MKL settings and single-kernel SVMs. Further, GLMKL achieves 85% accuracy which is still higher than that reached by the feature concatenation method. We can also note that we cannot overcome SVM-con when we use RBMKL, as the latter gives equal weight to each kernel. In fact, if there are inaccurate representations in the given set, the overall mean combination accuracy may be less of that reached using the single best. Conversely, when the weights are automatically estimated, such as in SimpleMKL and GLMKL the selection of the most reliable information is carried out by the MKL procedure and the overall performance improves.

<table><caption>**Table 2.** MKL accuracies (false positives and negatives are reported in brackets).</caption><thead><tr><th>Single-best</th><th>SVM-con</th><th>RBMKL</th><th>SimpleMKL</th><th>GLMKL</th></tr></thead><tbody><tr><td>*78.33 (10, 11.6)</td><td>83.33 (8.3, 8.3)</td><td>*81.67 (10, 8.3)</td><td><strong>86.67</strong> (6.6, 6.6)</td><td>85.00 (8.3, 6.6)</td></tr></tbody></table>

In Figure 2, we plotted the weights of MKL for both SimpleMKL and GLMKL algorithms. Note that the estimated weights are coherent in the two algorithms. As expected, the best representation is r02, which has the highest weights. Although the other representations with high weights (r08, r11 and r05) do not provide much accurate single-kernel SVMs results, their contributions to the overall accuracy in the combination is higher than those given by the other kernels. This demonstrates that when considering combinations, even a representation which does not lead to very precise results may contribute to raise the overall combination accuracy. Moreover, we can also deduce that these four representations are the most useful in discriminating between healthy and schizophrenic subjects, and we may focus the attention on these properties only.

Using this information, we also performed the above pipeline using only these four representations, and we can observe the results in Table 3. Using this subset, we get the highest accuracy with SimpleMKL<fnref n="9" />, reaching 88.33% of accuracy. We can also observe an increase in RBMKL.

<footnote n="9">⁹ Note that in principle the same result should have been obtained automatically from MKL algorithms on all representations. In practice, this is not the case in our experiment due to the fact that the estimated solution is trapped into a local minimum.</footnote>
<figure><img src="image_57.png" /><figcaption>**Fig. 2.** Combination weights in MKL using the linear kernel. Top: using SimpleMKL, Bottom: using GLMKL.</figcaption></figure>
<table><caption>**Table 3.** MKL accuracies on the selected subset of representations (false positives and negatives are reported in brackets).</caption><thead><tr><th>SVM</th><th>SVM-con</th><th>RBMKL</th><th>SimpleMKL</th><th>GLMKL</th></tr></thead><tbody><tr><td>*78.33 (10, 11.6)</td><td>*83.33 (6.6, 10)</td><td>*83.33 (6.6, 10)</td><td><b>88.33</b> (6.6, 5)</td><td>85.00 (6.6, 8.3)</td></tr></tbody></table>

## 5.4 Discussion

In this work, we have shown in general that MKL algorithms perform better than both single-best kernel SVMs and feature concatenation strategies. We have also observed that RBMKL (which does not compute weights while combining kernels) does not outperform the feature concatenation approach. Conversely, when the kernel combination is carried out by estimating proper weights, a drastic improvement is instead obtained. The kernel weights also allow us to extract useful information: it is interesting to observe that, for both MKL algorithms with the highest accuracy, four representations have the maximum effect (i.e., the highest weights), i.e., r02, r08, r11, and r05, with r02 being the best single-kernel. We use this information to select a smaller number of representations to reduce the costs of the feature extraction phase. Finally, we can also observe that by using such subset we can reach the best accuracy overall.

## 6 Conclusions

In this paper, we focus on scale selection for anatomical shape characterization. By employing a shape diffusion approach, we extract several shape descriptors at different scales in order to discriminate between healthy subjects and patients affected by schizophrenia. We have shown that machine learning techniques can be useful to improve the shape analysis in these (biomedical) contexts. We propose a Multiple Kernel Learning algorithm for the automatic estimation of the best feature representation for classification purposes. In this way, being driven by the training data, we are able to choose the scales of the heat kernel which are more suitable to describe our kind of shapes. In particular, in our experiments addressing the Thalamic region classification, we have shown that both small and high scales are crucial. Actually, the best accuracy is observed at r02 for which very local information are collected from the shape. Nevertheless, when also higher scales are considered the performance is further improved, meaning that also global shape information is relevant.

## Acknowledgements

We acknowledge financial support from the FET programme within the EU FP7, under the SIMBAD project (contract 213250). We thank Dr. Mehmet Gönen for the implementation of the MKL algorithms.
References

1. Agarwal, N., Port, J.D., Bazzocchi, M., Renshaw, P.F.: Update on the use of MR for assessment and diagnosis of psychiatric diseases. Radiology 255(1), 23–41 (2010)

2. American Psychiatric Association: Diagnostic and statistical manual of mental disorders, DSM-IV. Washington DC, 4th edn. (1994)

3. Bach, F.R., Lanckriet, G.R.G., Jordan, M.I.: Multiple kernel learning, conic duality, and the smo algorithm. In: Proceedings of the twenty-first international conference on Machine learning, ICML ’04. pp. 41–48 (2004)

4. Bronstein, A.M., Bronstein, M.M., Ovsjanikov, M., Guibas, L.J.: Shape Google: geometric words and expressions for invariant shape retrieval. ACM Transaction on Graphics 30(1), 1–20 (2011)

5. Bronstein, A.M., Bronstein, M.M., Ovsjanikov, M., Guibas, L.J.: Shape recognition with spectral distances. IEEE Trans. Pattern Analysis and Machine Intelligence 33(5), 1065–1071 (2011)

6. Castellani, U., Mirtuono, P., Murino, V., Bellani, M., Rambaldelli, G., Tansella, M., Brambilla, P.: A new shape diffusion descriptor for brain classification. In: Medical Image Comp. Computer-Assisted Intervention (MICCAI) (2011)

7. Chang, C.C., Lin, C.J.: LIBSVM: a library for support vector machines (2001), http://www.csie.ntu.edu.tw/$\sim$sjlin/libsvm

8. Corradi-DellAcqua, C., Tomelleri, L., Bellani, M., Rambaldelli, G., Cerini, R., Pozzi-Mucelli, R., Balestrieri, M., Tansella, M., Brambilla, P.: Thalamic-insular dysconnectivity in schizophrenia: Evidence from structural equation modeling. Human Brain Mapping p. in press (2011)

9. Cristianini, N., Shawe-Taylor, J.: An Introduction to Support Vector Machines and other Kernel-based Learning Methods. Cambridge University Press (2000)

10. Gebal, K., Baerentzen, J.A., Aanaes, H., Larsen, R.: Shape analysis using the auto diffusion function. In: In SGP (2009)

11. Giuliani, N.R., Calhouna, V.D., Pearlson, G.D., Francis, A., Buchanan, R.W.: Voxel-based morphometry versus region of interest: a comparison of two methods for analyzing gray matter differences in schizophrenia. Schizophrenia Research 74(2–3), 135–147 (2005)

12. Kloft, M., Brefeld, U., Sonnenburg, S., Zien, A.: $l_p$-norm multiple kernel learning. Journal of Machine Learning Research 12, 953–997 (2011)

13. Lanckriet, G.R.G., Cristianini, N., Bartlett, P., Ghaoui, L.E., Jordan, M.I.: Learning the kernel matrix with semidefinite programming. Journal of Machine Learning Research 5, 27–72 (December 2004)

14. Rakotomamonjy, A., Bach, F., Canu, S., Grandvalet, Y.: SimpleMKL. Journal of Machine Learning Research 9, 2491–2521 (2008)

15. Raviv, D., Bronstein, A.M., Bronstein, M.M., Kimmel, R.: Volumetric heat kernel signatures. In: Workshop on 3D Object Retrieval (2010)

16. Sun, J., Ovsjanikov, M., Guibas, L.: A concise and provably informative multi-scale signature based on heat diffusion. In: Proceedings of the Symposium on Geometry Processing. pp. 1383–1392 (2009)

17. Vapnik, V.N.: Statistical learning theory. John Wiley and Sons (1998)

18. Xu, Z., Jin, R., Yang, H., King, I., Lyu, M.R.: Simple and efficient multiple kernel learning by group Lasso. In: Proceedings of the 27th International Conference on Machine Learning, ICML ’10. pp. 1175–1182 (2010)
# Probabilistic Modeling of Landmark Distances and Structure for Anomaly-proof Landmark Detection

Shouhei Hanaoka<fnref n="1" />, Yoshitaka Masutani<fnref n="1" />,², Mitsutaka Nemoto<fnref n="1" />, Yukihiro Nomura<fnref n="1" />, Takeharu Yoshikawa³, Naoto Hayashi³, Naoki Yoshioka⁴, Kuni Ohtomo<fnref n="1" />,²

**Abstract.** A combinatorial optimization algorithm for detecting multiple anatomical landmarks is presented. It can determine the positions of over 100 landmarks concurrently, taking spatial correlations of all landmark pairs into account. Provided that a set of landmark candidate lists is given by sensitivity-optimized single-landmark detectors, the proposed algorithm can find the most probable combination of them through solving a MAP estimation-based combinatorial optimization problem. Additionally, it is designed to handle subjects with "segmentation anomaly of the spinal column," a common anatomical anomaly of the spine. The proposed system was evaluated with 156 landmarks in 50 datasets, using virtually created detector output sets. In the result, the algorithm achieved 97.6% of spinal anomaly estimation accuracy even with 50 points of candidates given per landmark, as well as 96.2% of accuracy in landmark candidate selection. From these results, usefulness of the proposed algorithm for subjects with spinal anomaly was suggested.

**Keywords:** Landmark, Combinatorial optimization, MAP estimation, Anatomical anomaly, Computed tomography, Spine

## 1 Introduction

Landmark point detection algorithms are extensively researched and widely used in various medical image processing applications. However, it is a difficult task to detect a large number of landmarks correctly, because the human body includes a lot of similar regions sharing their appearances. Even most of frequently-used landmarks do not have truly unique local shape or intensity. Consequently, detection results often include a certain number of false positive candidates.

Furthermore, some sort of important anatomical landmarks can intrinsically be non-existent, not only in patients due to pathological condition, but even in healthy subjects due to anatomical anomalies (fig. 1a). Confirming such a situation is very difficult by detecting each landmark independently and sequentially.

<footnote n="1">¹ Department of Radiology, ² Division of Radiology and Biomedical Engineering, Graduate School of Medicine, ³ Department of Computational Diagnostic Radiology and Preventive Medicine and ⁴ Department of Integrated Imaging Informatics, The University of Tokyo hospital, 7-3-1 Hongo, Bunkyo-ku, Tokyo, Japan
hanaoka-tky@umin.ac.jp</footnote>
One solution is to determine the entire landmark positions simultaneously, taking their spatial relation into account. It can be done by dividing the whole problem into 2 sequential phases: the individual landmark detection phase and the combinatorial optimization phase (figs. 1b, 1c). In the former phase, each landmark is detected by a single detector. Each detector is optimized to maximize the sensitivity, not the specificity, so the detection result forms a candidate list which includes a lot of false positives. In the latter phase, the best combination of choice from all candidate lists is selected with the use of a priori knowledge on the inter-landmark relationship. Defect of any landmark can also be detected through this optimization phase.

<figure><img src="image_58.png" /><figcaption>**Fig. 1.** (a) An example of anatomical anomaly and landmark deficit in the 5th lumbar vertebra. (b) A schema of the framework of the proposed method. Note that, for the non-existent landmark (□), all detected candidates are rejected by the following combinatorial optimization.</figcaption></figure>

Among related studies, Seifert et al. reported a framework to automatically detect 19 landmarks. [1] In their approach, a belief propagation algorithm is used with prior knowledge about landmarks’ geometric relationships such as “to the right of”, “close to”, etc. Though their result was excellent, the prior knowledge used and the way to build it were not described precisely.

In this paper, we introduce a framework to determine over 100 landmark positions reliably. Especially, we focus upon our method to solve the landmark-set combinatorial optimization problem with use of maximum *a posteriori* (MAP) estimation. The method was evaluated with 156 landmarks in 50 human body CT image datasets. Virtually-created detector outputs were used in this study, rather than those of real detectors [2], in order to evaluate the ability of our framework to handle a large number of false candidates. Furthermore, handling of common anomaly (alteration in the number of thoracic/lumbar vertebrae) was also evaluated.
## 2 Methods

### 2.1. Definition

**2.1.1. Detector output.** Firstly, we defined a landmark detector generally. In our framework, each detector has to output not only a series of candidate positions, but also estimated probabilities of them. Based on this, outputs from one detector (for the m-th landmark) are defined as:

<table><tr><td>$N_m$</td><td>number of the candidates</td><td>(0 ≤ $N_m$)</td></tr><tr><td>$\mathbf{c}_m^i$</td><td>coordinates of the i-th candidate</td><td>(1 ≤ $i$ ≤ $N_m$)</td></tr><tr><td>$p_m^i$</td><td>detector-estimated probability of the i-th candidate</td><td>(0 &lt; $p_m^i$ ≤ 1)</td></tr><tr><td>$p_m^3$</td><td>probability that “none of candidates $\mathbf{c}_m^i$ is correct”</td><td>(0 ≤ $p_m^3$ ≤ 1)</td></tr></table>

The whole set of output is also defined as $I_m = \{N_m, \mathbf{c}_m^1, \dots, \mathbf{c}_m^{N_m}, p_m^1, \dots, p_m^{N_m}, p_m^3\}$.

Any output $I_m$ must satisfy $p_m^3 + \sum p_m^i = 1$.

Letting the true landmark position be $\mathbf{x}_m$, the detector-estimated probability $p_m^i$ can be interpreted as a conditional probability with a certain detector output $I_m$.

$$p_m^i = p(\mathbf{x}_m = \mathbf{c}_m^i | I_m) \qquad (1)$$

If $p_m^3 > 0$, it indicates that there are some possibilities of “no true landmark position is included in the candidate set $\{\mathbf{c}_m^1, \dots, \mathbf{c}_m^{N_m}\}$.” The proposed algorithm can consider such a situation, and handles it as one extra state which will be represented as “$\mathbf{x}_m = \mathfrak{I}$” in this paper. The state $\mathfrak{I}$ can be a “true answer” when the target landmark is out of the imaging range, or corresponding landmark does not exist anywhere (e.g., due to some anatomical anomaly or pathological condition). In this describing method, the probability $p_m^3$ can also be interpreted as a conditional probability as

$$p_m^3 = p(\mathbf{x}_m = \mathfrak{I} | I_m). \qquad (2)$$

**2.1.2. The prior probability distribution of LM positions.** Secondly, the prior probability function for all possible landmark position sets has to be defined in advance.

Let $\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_M$ be the positional vectors of total M landmarks, and $\mathbf{X} = (\mathbf{x}_1^\top \mathbf{x}_2^\top \dots \mathbf{x}_M^\top)^\top$ be the concatenated form of them. Each element $\mathbf{x}_m$ can be regarded as a stochastic variable whose domain is $\mathbf{x}_m \in \{\mathfrak{I}, \mathbf{c}_m^1, \dots, \mathbf{c}_m^{N_m}\}$. The aim is to approximate the prior probability distribution of $\mathbf{X}$ as a single function $p(\mathbf{X})$, in order to use in the following MAP estimation.

In this study $p(\mathbf{X})$ is defined as a function of the squared distances between all landmark pairs. When any $\mathbf{x}_m$ ($1 \le m \le M$) satisfies $\mathbf{x}_m = \mathfrak{I}$, however, the corresponding term of $p(\mathbf{X})$ is replaced by a conventionally-defined constant term. The calculation method of $p(\mathbf{X})$ without considering $\mathfrak{I}$ is discussed in the next chapter, which will be followed by the general definition of $p(\mathbf{X})$.
(i) When $\mathbf{x}_m \neq \mathfrak{T}, \forall m$. Let

$$d_{i,j} = |\mathbf{x}_i - \mathbf{x}_j|^2 \quad (1 \le i < j \le M) \quad (3)$$

be the squared distance between the i-th and j-th landmarks. Note that $d_{i,j}$ is defined only if both $\mathbf{x}_i$ and $\mathbf{x}_j$ are not $\mathfrak{T}$. Then, the distance is normalized by its average $E(d_{i,j})$ and variance $V(d_{i,j})$ in the training datasets. The normalized distance $g_{i,j}$ follows the equation

$$g_{i,j} = \frac{d_{i,j} - E(d_{i,j})}{\sqrt{V(d_{i,j})}} \quad \text{if } \mathbf{x}_i \neq \mathfrak{T}, \mathbf{x}_j \neq \mathfrak{T}. \quad (1 \le i < j \le M) \quad (4)$$

Let a vector $\mathbf{G} = (g_{1,2} g_{1,3} \dots g_{i,j} \dots g_{M-1,M})$ be the concatenated normalized squared distances between all landmark pairs. Note that the vector $\mathbf{G}$ has $M C_2 = \frac{M(M-1)}{2}$ of elements. Then, the prior probability distribution $p(\mathbf{X})$ was approximated by a multivariate normal distribution of the vector $\mathbf{G}$. That is,

$$p(\mathbf{X}) = \frac{1}{(\sqrt{2\pi})^{M C_2} \cdot \sqrt{|\mathbf{V}|}} \exp\left(-\frac{1}{2} \mathbf{G}' \mathbf{V}^{-1} \mathbf{G}\right) \quad (5)$$

where $\mathbf{V}$ represents the covariance matrix of $\mathbf{G}$, which is calculated from training sets, and $|\mathbf{V}|$ is the determinant of $\mathbf{V}$. The size of the matrix $\mathbf{V}$ is $M C_2 \times M C_2$.

(ii) General definition. Eq. (5) can be written in an extended form

$$p(\mathbf{X}) = \frac{1}{(\sqrt{2\pi})^{M C_2} \cdot \sqrt{|\mathbf{V}|}} \exp\left(-\frac{1}{2} \sum_{i=1}^{M-1} \sum_{j=i+1}^{M} \sum_{k=1}^{M-1} \sum_{l=k+1}^{M} g_{i,j} \cdot \{\mathbf{V}^{-1}\}_{(i,j),(k,l)} \cdot g_{k,l}\right) \quad (6)$$

where $\{\mathbf{V}^{-1}\}_{(i,j),(k,l)}$ is the corresponding element of the matrix $\mathbf{V}^{-1}$. Here, the summed term $g_{i,j} \cdot \{\mathbf{V}^{-1}\}_{(i,j),(k,l)} \cdot g_{k,l}$ is not available when any of $\mathbf{x}_i, \mathbf{x}_j, \mathbf{x}_k$ or $\mathbf{x}_l$ is $\mathfrak{T}$. In order to define $p(\mathbf{X})$ generally, we replaced them as follows:

$$p(\mathbf{X}) = \frac{1}{(\sqrt{2\pi})^{M C_2} \cdot \sqrt{|\mathbf{V}|}} \exp\left(-\frac{1}{2} \sum_{i=1}^{M-1} \sum_{j=i+1}^{M} \sum_{k=1}^{M-1} \sum_{l=k+1}^{M} t_{(i,j),(k,l)}\right) \quad (7)$$
$$t_{(i,j),(k,l)} = \begin{cases} g_{i,j} \cdot \{\mathbf{V}^{-1}\}_{(i,j),(k,l)} \cdot g_{k,l} & \text{if all of } \mathbf{x}_i, \mathbf{x}_j, \mathbf{x}_k \text{ and } \mathbf{x}_l \text{ are not } \mathfrak{T} \\ 0 & \text{otherwise, and } (i, j) \neq (k, l) \\ \gamma & \text{otherwise, and } (i, j) = (k, l) \end{cases}$$

When $(i,j) \neq (k,l)$, the summed term $t_{(i,j),(k,l)}$ evaluates how $g_{i,j}$ and $g_{k,l}$ are correlated as like that of training sets. When no prior information on $g_{i,j}$ (or $g_{k,l}$) is available, the expected value of this term should be zero. That is why it should be replaced by zero in case of $\mathfrak{T}$.

On the other hand, when $(i,j)=(k,l)$, the term $t_{(i,j),(i,j)} = \{\mathbf{V}^{-1}\}_{(i,j),(i,j)} \cdot g_{i,j}^2$ always has a positive value which evaluates how the distance $d_{i,j}$ varies from that of the training sets. Therefore, we replaced it by a positive constant $\gamma$. Because increasing this
parameter will reduce the probability $p(\mathbf{X})$ for any $\mathbf{x}_m=\mathfrak{T}$, $\gamma$ can be regarded as a parameter which controls how the algorithm avoids the state $\mathfrak{T}$, or how it prefers any detected candidate $\mathbf{c}_m^i$, for each landmark. We empirically selected $\gamma=2$ based on our preliminary experiments (the data is not shown).

**2.1.3. a posteriori probability distribution.** Once the detector outputs are given, the probability a posteriori can be calculated with Bayes' theorem. From the series of detector outputs $I_1, I_2, ..., I_M$ and the prior probability distribution $p(\mathbf{X})$, the posterior probability can be calculated by the Bayes' as:

$$p(\mathbf{X} | I_1, I_2, I_3, ..., I_M) = \frac{p(I_1, I_2, I_3, ..., I_M | \mathbf{X}) \cdot p(\mathbf{X})}{p(I_1, I_2, I_3, ..., I_M)} \quad (8)$$

The denominator is constant and independent of $\mathbf{X}$. Therefore, the maximum a posteriori estimation of landmark position set $\mathbf{X}$ is as follows:

$$\hat{\mathbf{X}} = \arg \max_{\mathbf{x}} p(\mathbf{X} | I_1, I_2, I_3, ..., I_M) = \arg \max_{\mathbf{x}} p(I_1, I_2, I_3, ..., I_M | \mathbf{X}) \cdot p(\mathbf{X}) \quad (9)$$

We assumed that all detector outputs are independent of each other, as well as independent on the positions of the other landmarks. It means that each detector output is only dependent on the corresponding landmark position. Then the term $p(I_1, I_2, I_3, ..., I_M | \mathbf{X})$ in (9) can be divided into the product of single-landmark conditional probabilities $p(I_m | \mathbf{x}_m)$. That is,

$$p(I_1, I_2, I_3, ..., I_M | \mathbf{X}) = \prod_{m=1}^{M} p(I_m | \mathbf{x}_m) = \prod_{m=1}^{M} \frac{p(\mathbf{x}_m | I_m) \cdot p(I_m)}{p(\mathbf{x}_m)} \quad (10)$$

The term $p(I_m)$ in this formula is independent of $\mathbf{X}$, so it can be ignored through the MAP estimation. The denominator $p(\mathbf{x}_m)$ is also can be ignored, because the term is the probability distribution of one single landmark position $\mathbf{x}_m$ without any prior information. Ignoring the term is equivalent to regarding it to be constant and homogeneous anywhere. The remaining term $p(\mathbf{x}_m | I_m)$ is the detector-estimated probability as described in Eqs. (1) and (2). That is,

$$p(\mathbf{x}_m | I_m) = \begin{cases} p_m^i & \text{if } \mathbf{x}_m = \mathbf{c}_m^i, \exists i \\ p_m^\mathfrak{T} & \text{if } \mathbf{x}_m = \mathfrak{T} \end{cases} \quad (11)$$

In conclusion, the MAP estimation will be performed by the following formula:

$$\hat{\mathbf{X}} = \arg \max_{\mathbf{x}} p(\mathbf{X} | I_1, I_2, I_3, ..., I_M) = \arg \max_{\mathbf{x}} \left[ \left\{ \prod_{m=1}^{M} p(\mathbf{x}_m | I_m) \right\} \cdot p(\mathbf{X}) \right] \quad (12)$$
## 2.2. Implementation

**2.2.1. Tikhonov's regularization.** In practice, the dimension of $\mathbf{G}$ in Eq. (5) (or Eq. (7)), which equals $M(C_2 = \frac{M(M-1)}{2})$, can be much greater than the number of training cases $N$. In such a case, the covariance matrix $\mathbf{V}$ will be rank-deficit and has no inverse matrix $\mathbf{V}^{-1}$. To avoid them, Tikhonov's regularization [3] was performed by replacing $\mathbf{V}$ with its regularized matrix $\mathbf{V}_{\text{reg}}$ as follows:

$$ \mathbf{V}_{\text{reg}} = \mathbf{V} + \lambda \mathbf{I} \qquad (13) $$

We empirically selected $\lambda = 1.0$ for this study.

Note that the calculation of Eq. (5) can be speeded up by using the following formula (derived from the Woodbury matrix identity),

$$ \mathbf{V}_{\text{reg}}^{-1} = (\mathbf{V} + \lambda \mathbf{I})^{-1} = (\mathbf{U}^t \mathbf{D} \mathbf{U} + \lambda \mathbf{I})^{-1} = \lambda^{-1} \mathbf{I} - \lambda^{-2} \mathbf{U}^t (\mathbf{D}^{-1} + \lambda^{-1} \mathbf{I})^{-1} \mathbf{U}. \qquad (14) $$

Here, $\mathbf{D}$ is a $(N-1) \times (N-1)$ diagonal matrix whose diagonal elements are nonzero eigenvalues of $\mathbf{V}$, and $\mathbf{U}$ is a $(N-1) \times M(C_2)$ matrix whose row vectors are corresponding eigenvectors. $\mathbf{I}$ and $\mathbf{I}'$ are identity matrices whose sizes are $M(C_2)$ and $N-1$, respectively.

**2.2.2. Combinatorial optimization.** The maximization of Eq. (12) for all combinations of $\mathbf{X} = \{\mathbf{x}_m\}, \mathbf{x}_m \in \{\mathfrak{S}, \mathbf{c}_m^1, \dots, \mathbf{c}_m^{N_m}\}$ is performed by a Gibbs' sampler-based simulated annealing algorithm reported by Geman & Geman [4]. In the algorithm, a virtual temperature $T$ is introduced to modify the probability distribution. The distribution in (12) was modified as follows:

$$ p_{\text{modified}}(\mathbf{X}; T) = \frac{1}{Z} \cdot \left[ \left\{ \prod_{m=1}^{M} p(\mathbf{x}_m | I_m) \right\} \cdot p(\mathbf{X}) \right]^{-\frac{1}{T}} \qquad (15) $$

$Z$ is a normalization factor in order to make the sum of probability 1.

When $T=1$, the modified distribution equals to the original in (12). With higher $T$, the modified distribution is almost homogeneous for any $\mathbf{X}$ in its domain. However, it becomes sharper and more local with lower $T$. In the limiting case of $T \to +0$, the modified distribution has nonzero probability only at its global maximum point. (Fig. 2)
<figure><img src="image_59.png" /><figcaption>**Fig. 2.** A simple example how the probability distribution is modified by the temperature. The green curve illustrates an original distribution (a unary mixture Gaussian distribution for example). The modified distributions with $T = 10$ and $0.1$ are illustrated with red and blue curves, respectively. Note that the distribution becomes flat with higher temperature, while it becomes sharper with lower temperature.</figcaption></figure>

In the simulated annealing, each $x_m (m = 1, 2, ..., M)$ is sequentially and repeatedly sampled from the modified distribution by a Gibbs' sampler. The sampling begins with very high $T$, which decreases gradually, and finally it gets so low that the system converges to the maximum point. In this study the simulated annealing were performed for 100 cycles with $T=1000$ (in order to cancel the effect of the initial condition; so-called *burn-in*), then it was gradually cooled down to $T=0.01$ through 1000 cycles.

**2.2.3. Handling of vertebral anomalies.** A majority of human beings have 12 thoracic and 5 lumbar vertebrae. Segmentation anomaly of the spinal column is a common anatomical anomaly in which the subject has 11 or 13 thoracic, and/or 4 or 6 lumbar, vertebrae. The prevalence is, in a report, about 9% [5]. This anomaly is very problematic in both defining and detecting vertebral landmarks (fig. 1a). The “state $\mathfrak{T}$” approach described in chapter 2.1. is not enough for them, because it is not only a local banishment of a single anatomical entity but causing a global morphological change in the spine.

To overcome this, a series of “anomaly landmark position set converters” are introduced. One converter can convert any landmark position set in a subject with a certain type of anomaly (e.g., 6 lumbar vertebrae, or “6L”) into a virtually normalized landmark position set (i.e., a landmark position set as if she or he has only 5 lumbar vertebrae). It is simply performed by replacing each landmark coordinates by an appropriate internally dividing point between two of them (fig. 3). Because 7 types of anomalies (11T, 13T, 4L, 6L, 11T+6L, 13T+4L) were considered in this study, 7 different converters were designed. Through one of these converters, abnormal spines can be converted into a “normalized” one, with which the prior probability $p(X)$ can be calculated.
<figure><img src="image_60.png" /><figcaption>**Fig. 3.** The anomaly conversion. In this example, (a) a landmark position set in a case of 4-lumbar vertebra (4L) anomaly is virtually converted into (b) a “normal” landmark set with 5 vertebrae. With this conversion, the posterior probability $p(X)$ of the given landmark position set can be calculated in the same manner as in the cases with (c) normal 5-lumbar spines.</figcaption></figure>

Because the algorithm do not know which anomaly is correct, it must estimate it. It can be done by (i) hypothesizing one anomaly (or normal), (ii) calculating the posterior probability under the hypothesis and (iii) comparing the probabilities between all hypotheses. In detail, two different strategies were evaluated:

(i) *Comparison after all optimization.* Firstly, a series of combinatorial optimization with all converters are performed. Then, the converter with the largest probability is chosen as the estimated anomaly (or normal).

(ii) *Comparison in situ.* The comparison is performed at the end of every cycle in the simulated annealing. In other words, the state “which converter is currently selected” is also dealt with as one extra variable to be optimized.

## 2.3. Evaluation

**2.3.1. Virtual detector output construction.** The detector outputs were virtually created for each landmark and for each CT dataset. If the target landmark existed in the subject's real body, it was included as a candidate, as well as 25, 50, 75 or 100 of false positive candidates. Each false positive candidate was determined randomly following a 3-D Gaussian probability distribution (the center corresponds to the true point, and the standard deviations were $\sigma_x = \sigma_y = \sigma_z = \sqrt{3} \cdot 100$ millimeters). Any false candidates within 20 millimeters from the true point were removed. If the target landmark did not exist in reality, only false positive points were added, using another adjacent landmark as the distribution center. The detector-estimated probability $p_m^3$ was fixed to be 0.05, and all of $p_m^i$ were set to be uniform.

**2.3.2. Experimental settings.** The experiments were performed on a workstation with an Intel® Core™ i7-2600 processor and one NVIDIA® Tesla™ C2050 GPU
computing processor. The latter was utilized to speed up the calculation. The processing times were 9.7 minutes per case for *comparison after all optimization* strategy and 3.6 minutes for *comparison in situ* strategy with 100 candidates per landmark.

5 times of experiments were performed for each experimental setting, using different sets of virtual detector outputs.

The statistical model of squared distances (the mean $E(d_{ij})$, variance $V(d_{ij})$ and the covariance matrix $V$) was calculated by leave-one-case out method.

## 3 Results and discussion

The proposed algorithm was evaluated with 156 bony landmark points which were manually inputted in 50 thin-slice clinical human body CT datasets. The list of landmark used is available in [2]. Among the 50 datasets, total 9 had a segmentation anomaly (three 13T, one 11T, three 6L, two 4L and none of 13T+4L or 11T+6L).

The summary of results is shown in Table 1. An example result is also shown in Figure 4. The accuracies of anomaly estimation were from 87.6 to 97.6%, varying among strategies and number of candidates. The accuracy of anomaly estimation was better in *comparison after optimization* strategy than in *comparison in situ* strategy. Though the former strategy takes approximately 3 times longer time, it seems to be useful especially for a larger number of candidates.

<table><caption>**Table 1.** The entire result of anomaly estimation and combinatorial optimization.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">No. of<br>cands</th><th colspan="3">anomaly estimation accuracy (%)</th><th colspan="3">landmark determination:<br>existent-in-reality LMs (%)</th><th colspan="2">non-existent<br>LMs (%)</th></tr><tr><th>anomaly<br>cases</th><th>normal cases</th><th>overall</th><th>TP</th><th>FN</th><th>FP<sub>candidate</sub></th><th>TN</th><th>FP<sub>existence</sub></th></tr></thead><tbody><tr><td rowspan="4">comparison<br>after opt.</td><td>25</td><td>88.9 ±13.6</td><td>99.5 ±1.1</td><td><strong>97.6 ±3.3</strong></td><td>96.60 ±2.90</td><td>0.09 ±0.13</td><td>3.31 ±2.78</td><td rowspan="4">100<br>±0</td><td rowspan="4">0<br>±0</td></tr><tr><td>50</td><td>88.9 ±7.9</td><td>99.5 ±1.1</td><td><strong>97.6 ±1.7</strong></td><td>96.25 ±0.34</td><td>0.08 ±0.05</td><td>3.67 ±0.32</td></tr><tr><td>75</td><td>86.7 ±14.5</td><td>99.0 ±1.3</td><td><strong>96.8 ±2.7</strong></td><td>94.40 ±0.39</td><td>0.09 ±0.08</td><td>5.51 ±0.36</td></tr><tr><td>100</td><td>77.8 ±11.1</td><td>94.6 ±4.0</td><td><strong>91.6 ±3.0</strong></td><td>91.33 ±1.03</td><td>0.26 ±0.09</td><td>8.42 ±1.00</td></tr><tr><td rowspan="4">comparison<br>in situ</td><td>25</td><td>84.4 ±12.7</td><td>100 ±0</td><td><strong>97.2 ±2.3</strong></td><td>97.65 ±0.36</td><td>0.09 ±0.08</td><td>2.26 ±0.32</td><td rowspan="4"></td><td rowspan="4"></td></tr><tr><td>50</td><td>73.3 ±18.6</td><td>100 ±0</td><td><strong>95.2 ±3.4</strong></td><td>95.67 ±0.70</td><td>0.17 ±0.12</td><td>4.15 ±0.59</td></tr><tr><td>75</td><td>64.4 ±19.9</td><td>99.5 ±1.1</td><td><strong>93.2 ±3.9</strong></td><td>93.37 ±1.02</td><td>0.24 ±0.15</td><td>6.40 ±0.91</td></tr><tr><td>100</td><td>48.9 ±14.9</td><td>96.1 ±2.8</td><td><strong>87.6 ±2.6</strong></td><td>89.64 ±0.61</td><td>0.40 ±0.05</td><td>9.96 ±0.57</td></tr></tbody></table>
<figure><img src="image_61.png" /><figcaption>**Fig. 4.** An example result of a case with 6 lumbar vertebra (6L) anomaly. (a) The false (green) and true (orange) candidates outputted by the virtual detector for the tip of right transverse process of the 4th lumbar vertebra. (b) The optimization result for 156 landmarks.</figcaption></figure>

The optimization result of each landmark in each case was also evaluated. Any landmark which existed in reality was classified as one of TP, $FP_{candidate}$ or FN. Within them, TP means that the algorithm correctly selected the true position candidate. $FP_{candidate}$ means that the algorithm selected any false position candidate. If the algorithm concluded that the target landmark did not exist, but in reality it existed, it was categorized as FN. In this study, the accuracy ratios (ratios of TP) for existent landmarks were varied from 89.64 to 97.65%.

Additionally, all of non-existing landmarks were classified as either TN or $FP_{existence}$. TN means that the target landmark did not exist in reality and the algorithm correctly found it (either by adopting some anomaly hypothesis which does not include the target landmark, or by selecting the state $\mathfrak{T}$ for the landmark.) If the algorithm chose any candidate other than $\mathfrak{T}$, it was classified as $FP_{existence}$. In this study all of non-existent landmarks were determined as TN, so the accuracy was 100%.

It is noticeable that, throughout all experimental condition, none of non-existing landmarks were mistaken as $FP_{existence}$. In fact, without any exception, all mistakes in anomaly estimation were underestimation of thoracic / lumbar vertebral number. It may imply that our algorithm was not sufficiently optimized for the problem, favoring anomalies having less vertebral number, and yet to be investigated in the future work. Also, we are now planning to overcome this limitation by adding some additional term, which evaluates the regularity of vertebral bones' alignment, to the MAP estimation.

# 4 Conclusion

A novel combinatorial optimization algorithm for landmark detection and anomaly estimation was presented. The proposed method showed fair results even with a large number of landmark position candidates. Additionally, its feasibility to estimate segmentation anomaly of the vertebrae, which is one of the most common and problematic anomalies in detecting bone landmarks, has been shown. Therefore, we believe that our algorithm is useful in medical image analysis such as a pre-
process for computer-assisted detection/diagnosis (CAD) applications. The future work will include improvement of accuracy in anomaly estimation and evaluation with real detector outputs with a large number of datasets.

## Acknowledgement

This study is a part of the research project "Computational Anatomy for Computer-aided Diagnosis and therapy: Frontiers of Medical Image Sciences", which is financially supported by the grant-in-aid for scientific research on innovative areas MEXT, Japan.

## References

1. Seifert, S., Barbu, A., Zhou, S.K., et. al.: Hierarchical parsing and semantic navigation of full body CT data. In: Pluim, J.P.W., Dawant, B.M. (eds.) SPIE (2009)

2. Nemoto, M., Masutani, Y., Hanaoka, S. et. al.: A unified framework for concurrent detection of anatomical landmarks for medical image understanding. In: Proc. SPIE 2011, 7962--121 (2011)

3. Crimi, A., Sporring, J., Bruijne, M., et. al.: Prior knowledge regularization in statistical medical image tasks. Probabilistic Models for Medical Image Analysis 2009, MICCAI, Imperial College London, UK (2009)

4. Geman, S., Geman, D.: Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol.6, no.6, 721--741 (1984)

5. Carrino, J. A., Campbell, P. D., Lin, D.C., et al.: Effect of spinal segment variants on numbering vertebral levels at lumbar MR imaging. Radiology, 259:203--212 (2011)
Hippocampal Morphometry Study by
Automated Surface Fluid Registration and its
Application to Alzheimer's Disease

Jie Shi¹, Yuting Wang², Paul M. Thompson³, and Yalin Wang¹

¹ School of Computing, Informatics and Decision Systems Engineering, ASU, Tempe,
AZ 85281, USA

² Mathematics Department, UCLA, Los Angeles, CA 90095, USA

³ Lab. of Neuro Imaging, UCLA School of Medicine, Los Angeles, CA 90095, USA
jie.shi@asu.edu

**Abstract.** Alzheimer's disease (AD) is a severe and growing public health crisis. Efforts are underway to look for AD early detection in an efficient manner. Among all the AD biomarkers, hippocampal atrophy assessed on high-resolution T1-weighted MRI is the best established and validated. Hippocampal morphometry is increasingly used in the AD research, with modeling the hippocampus as a 3D parametric surface mesh. However, a major question in the analysis is how to align corresponding surface regions across subjects. Here we develop a system for detecting AD symptoms on hippocampal surfaces with an automated surface fluid registration method, which is based on conformal surface representation and mutual information regularized image fluid registration. Since conformal mappings are diffeomorphic and the mutual information method is able to drive a diffeomorphic flow that is adjusted to enforce appropriate surface correspondences in the surface parameter domain, combining conformal and fluid mappings will generate 3D shape correspondences that are diffeomorphic. We also incorporate in the system a novel method to compute curvatures using surface conformal parameterization. Experimental results in three hippocampal datasets show that the new system outperformed an early similar method and the popular SPHARM tool.

# 1 Introduction

Alzheimer's disease (AD) doubles in frequency of onset every 5 years after age 60, afflicting 1% of those aged 60 to 64, and 30-40% of those over 85. Many MRI-based measures of atrophy in several structural measures, including whole-brain, hippocampus [14, 16, 19], and ventricular enlargement [19], etc., correlate closely with changes in cognitive performance, supporting their validity as markers of AD progression [8]. Among all the biomarkers, hippocampal atrophy assessed on high-resolution T1-weighted MRI is the best established and validated. As a result, detection of valid and efficient morphometry changes in hippocampus and their correlation with other cognitive functions and biomarkers becomes a
key research topic for clinical diagnosis and monitoring of patients with sus-
pected Alzheimer's disease. Although most hippocampus analysis used volume
as the atrophy measurement [10], recent researches [1, 16, 19, 26] demonstrated
that surface-based subcortical structure analysis may offer more advantages
because these methods studied patterns of hippocampal subfield atrophy and
produced detailed point-wise correlation between atrophy and cognitive func-
tions/biological markers [19].

Brain surface deformation studies typically require the computation of dense
correspondence vector fields that match one surface with another. Many brain
surface registration methods have been proposed [6, 20, 27]. Often, higher order
correspondences must be enforced between specific anatomical points and curved
landmarks lying within the two surfaces. This can be achieved by first mapping
each of the 3D surfaces to canonical parameter spaces such as a sphere [7] or
a planar domain [21]. A flow, computed in the parameter space of the two sur-
faces [4], induces a correspondence field in 3D. Artificial neural networks can
rule out or favor certain types of feature matches [15]. Correspondences may be
determined by using a minimum description length (MDL) principle, based on
the compactness of the covariance of the resulting shape model [5]. A key direc-
tion in surface registration research has been the computation of a diffeomorphic
surface map that also matches automatically identified surface features.

Using holomorphic 1-forms, a global conformal parameterization can be de-
veloped to conformally map a surface with complex topology (e.g. a surface
with branching topology) to a set of rectangular domains in the Euclidean plane.
The resulting parameterization helps in discretizing partial differential equations
(PDEs) for smoothing, denoising, or matching signals defined on the surface. The
mutual information (MI) method has been widely used to drive a diffeomorphic
flow in image registration. By adjusting the mutual information method to en-
force appropriate surface correspondences in the parameter domain, any scalar-
valued signals defined on the surfaces can also be aligned using the same flow
field. Conformal maps and fluid registration techniques can be combined to avoid
having to define a large set of manually-defined landmarks to constrain brain
surface correspondences. Since they generate diffeomorphic mappings, confor-
mal and fluid mappings together could generate 3D shape correspondences that
are diffeomorphic (i.e., smooth one-to-one correspondences). In [21, 22], Wang
et al. proposed an automated surface fluid registration method based on confor-
mal mapping and mutual information regularized image fluid registration and
applied it to register human faces and hippocampus. Here we develop a system
based on this technique for studying hippocampus in AD and incorporate a novel
method to compute surface curvatures as proposed in [12]. Our major contribu-
tions can be summarized as: (1). Introduction of a new stable method to compute
surface curvatures. (2). An automated hippocampal surface registration system
validated in three AD datasets with better performance than a previous similar
method [23] and SPHARM [18]. (3). The system will be publicly available [24].
Last, although the current system finds applications in AD detection, it is a
general method which may be applied to many other applications.
## 2 Method

### 2.1 Surface Conformal Parameterization

Let $S$ be a surface in $\mathbb{R}^3$ with an atlas $\{(U_\alpha, z_\alpha)\}$, where $(U_\alpha, z_\alpha)$ is a coordinate chart defined on $S$. The atlas thus is a set of consistent charts with smooth transition functions between overlapping charts. Here $z_\alpha: U_\alpha \to \mathbb{C}$ maps an open set $U_\alpha \subset S$ to a complex plane $\mathbb{C}$. If on any chart $(U_\alpha, z_\alpha)$ in the atlas, the Riemannian metric or the first fundamental form can be formulated as $ds^2 = \lambda(z_\alpha)^2 dz_\alpha d\bar{z}_\alpha$, and the transition maps $z_\beta \circ z_\alpha^{-1}: z_\alpha(U_\alpha \cap U_\beta) \to z_\beta(U_\alpha \cap U_\beta)$ are holomorphic, the atlas could be called conformal. Given a conformal atlas, a chart is compatible with the atlas if adding this chart still generates a conformal atlas. A conformal structure is obtained by adding all possible compatible charts to a conformal atlas. A Riemann surface is a surface with a conformal structure. One coordinate chart in the conformal structure introduces a conformal parameterization between a surface patch and the image plane. The conformal parameterization is angle-preserving and intrinsic to the surface geometry.

Let $(x_\alpha, y_\alpha)$ be the local parameter on the chart $(U_\alpha, z_\alpha)$ on $S$, a differential 1-form in $(x_\alpha, y_\alpha)$ can be defined as

$$ \omega = f(x_{\alpha}, y_{\alpha})dx_{\alpha} + g(x_{\alpha}, y_{\alpha})dy_{\alpha} \quad (1) $$

where $f, g$ are smooth functions. $\omega$ is a closed 1-form, if in each local parameter $(x_\alpha, y_\alpha)$, $\frac{\partial f}{\partial y_\alpha} - \frac{\partial g}{\partial x_\alpha} = 0$. If $\omega$ is the gradient of another function defined on $S$, it can be called an exact 1-form. An exact 1-form is also a closed 1-form. If a closed 1-form $\omega$ satisfies $\frac{\partial f}{\partial x_\alpha} + \frac{\partial g}{\partial y_\alpha} = 0$, then it is a harmonic 1-form. The gradient of a harmonic 1-form is an exact harmonic 1-form. The Hodge star operator acting on a differential 1-form gives the conjugate differential 1-form

$$ {}^*\omega = -g(x_{\alpha}, y_{\alpha})dx_{\alpha} + f(x_{\alpha}, y_{\alpha})dy_{\alpha} \quad (2) $$

Intuitively, the conjugate 1-form ${}^*\omega$ is obtained by rotating $\omega$ by a right angle everywhere. If $\omega$ is harmonic, so is its conjugate ${}^*\omega$. The holomorphic 1-form consists of a pair of conjugate harmonic 1-forms:

$$ \tau = \omega + \sqrt{-1}{}^*\omega \quad (3) $$

For a Riemann surface $S$ with genus $g > 0$, its conformal structure can always be represented in terms of a holomorphic 1-form basis, which is a set of $2g$ functions $\tau_i: K_i \to \mathbb{R}^2, i = 1, 2, \dots, 2g$ [22]. Any holomorphic 1-form $\tau$ is a linear combination of these functions. This finite-dimensional linear space generates all possible conformal parameterizations of surface $S$ and the quality of a global conformal parameterization is fundamentally determined by the choice of the holomorphic 1-form [22]. Then the conformal parameterization $\phi$ at point $p$ can be computed by integrating the holomorphic 1-form: $\phi(p) = \int_\gamma \tau$, where $\gamma$ is any path joining $p$ to a fixed point $c$ on the surface. Figure 1 (a) illustrates a brain cortical surface and its conformal parameterization to a square.
## 2.2 Surface Conformal Representation

It has been known that surface registration requires defining a lot of landmarks in order to align corresponding functional regions. Labeling features could be accurate but time-consuming. Here we show that surface conformal parame- terization could represent surface geometric features, thus avoiding the manual definition of landmarks.

For a general surface and its conformal parameterization $\phi : S \to \mathbb{R}^2$, the conformal factor at a point $p$ can be determined by the formulation:

$$ \lambda(p) = \frac{\text{Area}(B_{\epsilon}(p))}{\text{Area}(\phi(B_{\epsilon}(p)))} \qquad (4) $$

where $B_{\epsilon}(p)$ is an open ball around $p$ with a radius $\epsilon$. The conformal factor $\lambda$ encodes a lot of geometric information about the surface and can be used to compute curvatures and geodesic. In our system, we compute the surface mean curvatures only from the derivatives of the conformal parameterization as proposed in [12], instead of the three coordinate functions and the normal, which are generally more sensitive to digitization errors. Mathematically, the mean curvature is defined as:

$$ H = \frac{1}{2\lambda}\text{sign}(\phi)|\Delta\phi|, \quad \text{where } \text{sign}(\phi) = \frac{\langle \Delta\phi, \vec{N} \rangle}{|\Delta\phi|}. \qquad (5) $$

Using this formulation of $H$, we need to use the surface normal $\vec{N}$ only when computing $\text{sign}(\phi)$, which takes the value 1 or -1. Thus, the surface normal does not need to be accurately estimated and still we can get more accurate mean curvatures. Using the Gauss and Codazzi equations, one can prove that the conformal factor and mean curvature uniquely determine a closed surface in $\mathbb{R}^3$, up to a rigid motion. We call them the conformal representation of the surface. Figure 1 (b) shows the computed conformal factor (left) and mean curvature (right) on a hippocampal surface with color indices according to the values. Since conformal factor and mean curvature could represent important surface features and they are intrinsic to the surface, they may be used for surface registration.

## 2.3 Surface Fluid Registration Regularized by Mutual Information

After computing intrinsic geometric features, we align surfaces in the parameter domain with a fluid registration technique. Using conformal mapping, we essen- tially convert the surface registration problem to an image registration problem. The mutual information (MI) method has been successfully used to drive a dif- feomorphic flow in rigid [25] and non-rigid [13, 17] image registration. Image registration will be optimized when MI between two images is maximized. For MI to work, a monotonic mapping in grayscales between images is not required, so images from different modalities can be registered [11]. Hermosillo et al. [9] adopted linear elasticity theory to regularize the variational maximization of MI.
<figure><img src="image_62.png" /><figcaption>Fig. 1. Illustration of conformal parameterization (a) and geometric features (b). (c) shows the matching of geometric features in the 2D parameter domains using fluid registration with synthetic surfaces. Geometric features on 3D surfaces were computed and mapped to 2D conformal parameter domains.</figcaption></figure>

D'Agostino et al. [3] extended this approach to a viscous fluid scheme allowing large local deformations, while maintaining smooth, one-to-one topology [2]. We call this approach MI regularized fluid registration.

In [21, 22], Wang et al. proposed an automated surface fluid registration method combining conformal mapping and image fluid registration [3]. Let $I_1, I_2$ be the conformal representations of the target and the deforming template surfaces, respectively, the MI between two surfaces was defined as [21, 22]

$$I(u) = \int_{\mathbb{R}^2} p_u(i_1, i_2) \log \frac{p_u(i_1, i_2)}{p(i_1)p_u(i_2)} di_1 di_2 \quad (6)$$

where $p(i_1) = P(I_1(x) = i_1)$, $p_u(i_2) = P(I_2(x - u) = i_2)$, and $p_u(i_1, i_2) = P(I_1(x) = i_1 \& I_2(x - u) = i_2)$. Since conformal mapping and MI regularized fluid registration generate diffeomorphic mappings, a diffeomorphic surface-to-surface mapping is then recovered that matches surfaces in 3D. In our system, we adopt their methods of conformal mapping and fluid registration. However, our system differs from theirs in the computation of surface features as introduced in Sec. 2.2. The new way to compute mean curvature is more stable and less sensitive to normal computation, thus gives better representation of the surface features for registration.
# 3 Results

## 3.1 Synthetic Surface Registration Result

First, we illustrate the algorithm on synthetic surfaces. Figure 1 (c) illustrates a synthetic surface example. A pair of simple S-shape surfaces were generated. Corresponding 2D images were generated based on the sum of the local conformal factor and the mean curvature, expressed in the conformal parameterization domain. Some black horizontal lines were drawn on the surfaces to show equal distances on the surfaces and represent the differences in their shapes. The locations of the highest and lowest intensities are different (as shown by the positions of the horizontal stripes in the 2D images below). Using surface-based fluid registration, in the last image, the obtained horizontal line positions demonstrated an improved matching between features lying in the two surfaces.

In this experiment, we aim to visually explain the work flow of the algorithm. First, we conformally mapped each surface onto a planar domain and computed their conformal representation by combining conformal factor and mean curvature. Second, we scaled the conformal representation to form the 2D feature images as shown by the first two 2D images in Fig. 1 (c). As we pointed out in Sec. 2.2, the 2D images, i.e. the conformal representations of the surfaces clearly show the characteristics of the surfaces. Third, by registering the two feature images using fluid registration method, we induced a change in both the template feature image and the template surface. The last image in Fig. 1 (c) demonstrates that without changing the shape of the surface, the features on the template surface are well aligned with the target surface.

## 3.2 Hippocampal Surface Morphometry in Alzheimer's Disease

In this experiment, we test the robustness of our system by applying it to two clinical studies of hippocampal changes in Alzheimer's disease (AD). In our system, we leave two holes at the front and back of the hippocampal surface, representing its anterior junction with the amygdala, and its posterior limit as it turns into the white matter of the fornix. The resulting structure can then be logically represented as an open boundary genus-one surface, i.e., a cylinder. To better visualize the matching of surface features, we chose to encode surface features using a compound scalar function based on the local conformal factor and the mean curvature. A similar technique for AD study in hippocampus was proposed in [23], but with surface registration using constrained harmonic map. We take their method as a comparison. After the cross-subject registration is computed with one target surface selected, we examine shape differences using the multivariate tensor-based morphometry (MTBM) [23]. MTBM computes statistics from the Riemannian metric tensors that retain the full information in the deformation tensor fields, thus is more powerful in detecting surface differences than many other statistics [23].

Figure 2 illustrates the experimental results on a group of hippocampal surface models extracted from 3D brain MRI scans of 12 individuals with AD and
<figure><img src="image_63.png" /><figcaption>**Fig. 2.** The comparison of our method and a previous similar method on map of local shape differences (*p*-values), based on the multivariate TBM method with hippocampal surfaces in 12 AD versus 14 control subjects [19].</figcaption></figure>

14 matched healthy control subjects, the same dataset as in [19], where the data was manually traced by experienced neurologists. With MTBM, we can see that the significant areas detected with two methods are consistent. But with surface fluid registration, we detected larger areas with more significant group differences in the surface parameterization tensor, which is related to the relative area of regions in disease versus normality. The overall statistical significance of these group difference maps, based on permutation testing of the suprathreshold area of statistics (and therefore corrected for multiple comparisons) were *p*=0.0063 for the left hippocampal surface and 0.0298 for the right. This outperformed the previous method [23] (0.0205 for the left hippocampal surface and 0.1026 for the right). The cumulative distribution function (CDF) plot shows the comparison of resulting *p*-values for each method. While the line *y* = *x* represents null hypothesis, which is according to no difference, steeper curve shows more differences detected.

Figure 3 illustrates the experimental results on another hippocampal dataset from the Alzheimer's Disease Neuroimaging Initiative (ADNI) dataset (http://www.loni.ucla.edu/ADNI). Mild cognitive impairment (MCI) is an intermediate stage between the expected cognitive decline of normal aging and the more pronounced decline of dementia. If MCI could be found and treated, the risk of AD will be significantly reduced. However, at MCI stage, changes in brain surface are not significant thus impose more difficulty on the detection. Here we randomly chose 40 AD, 40 MCI, and 40 control subjects to test the effectiveness of our system. The hippocampus data was automatically segmented by a prior work [14]. With MTBM, we can see that, in the three experiments, our system
demonstrated better results than the previous method. Particularly, our system gave better MCI detection results when comparing with both AD and control subjects. In the experiment, all group difference *p*-maps were corrected using false discovery rate (FDR). The FDR method decides whether a threshold can be assigned to the statistical map that keeps the expected FDR below 5% (i.e., no more than 5% of the voxels are false positive findings). The CDF plots show the uncorrected *p*-values (as in a conventional FDR analysis). The *x* value at which the CDF plot intersects the *y* = 20*x* line represents the FDR corrected *p*-value or *q*-value. It is the highest statistical threshold that can be applied to the data, for which at most 5% false positives are expected in the map. In general, a larger *q*-value indicates a more significant difference in the sense that there is a broader range of statistic threshold that can be used to limit the rate of false positives to at most 5%. The use of the *y* = 20*x* line is related to the fact that significance is declared when the volume of suprathreshold statistics is more than 20 times that expected under the null hypothesis. Table 1 gives the FDR corrected *p*-values comparison.

## 3.3 Comparison with SPHARM

In this study, we tested our system on the whole ADNI baseline dataset (http://www.loni.ucla.edu/ADNI). The dataset consists of 233 healthy controls, 410 subjects with MCI, and 200 patients with AD. We excluded 1 subject from the control group and 2 subjects from the MCI group due to name duplication. For subjects with duplicated names, we retained the one which is the repeated scan. The hippocampal surfaces were automatically segmented using FIRST (http://www.fmrib.ox.ac.uk/fsl/first/index.html). FIRST is an integrated surface analysis tool developed as part of the FSL library, which is written mainly by members of the Analysis Group, FMRIB, Oxford, UK. We compared our system's results with the popular SPHARM tool [18] in surface registration on the same dataset. We adopted the suggested parameters for hippocampus in the manual of SPHARM [18]. In this experiment, 1 subject from each group (AD, MCI, control) failed the FIRST segmentation step probably due to the original images' resolution or contrast. After FIRST segmentation, we extracted the left and right hippocampi and saved each of them into a binary image. We fed the binary images as the input of both systems. As a result, 231 control, 199 AD, and 407 MCI subjects were successfully registered by our system, which shows the robustness of the system compared with SPHARM, in which 5 control, 17 AD, 14 MCI subjects failed either due to segmentation failure or parameterization failure. Fig. 4 shows the experimental results. From the *p*-map and the CDF plots, we can see that, with MTBM, our system outperformed the SPHARM

<table><caption>**Table 1.** FDR corrected *p*-values on hippocampal surfaces</caption><thead><tr><td></td><td>Our Method</td><td>Prior Method [23]</td></tr></thead><tbody><tr><td>AD-CTL</td><td>0.0476</td><td>0.0467</td></tr><tr><td>MCI-AD</td><td>0.0217</td><td>0.0199</td></tr><tr><td>MCI-CTL</td><td>0.0221</td><td>0.0185</td></tr></tbody></table>
in all three group difference studies. Table 2 gives the FDR corrected *p*-values comparison. Considering fairness, we also made comparisons by excluding those subjects that failed in SPHARM from our system in the statistical study. Table 3 gives the FDR corrected *p*-values comparison with the failed subjects excluded, i.e., 226 control, 182 AD, and 393 MCI subjects were studied in this experiment. From the table, we can see that our system still outperformed SPHARM on this dataset. Furthermore, because AD and MCI failed more than control subjects, the AD-MCI group has been more affected as expected.

# 4 Conclusion and Future Work

We develop an automated surface fluid registration system for hippocampal surface registration. Experiments on different types of AD hippocampal datasets all demonstrate our system's stronger statistical power. Ongoing work is to apply this system to automatically map lateral ventricle enlargements in Alzheimer's disease and those at risk.

# References

1. Apostolova, L.G., et al.: Subregional hippocampal atrophy predicts Alzheimer's dementia in the cognitively normal. Neurobiol Aging 31(7), 1077-1088 (2008)

2. Christensen, G.E., et al.: Deformable templates using large deformation kinematics. IEEE Trans Image Process 5(10), 1435-1447 (1996)

3. D'Agostino, E., et al.: A viscous fluid model for multimodal non-rigid image registration using mutual information. Med Image Anal 7(4), 565-575 (2003)

4. Davatzikos, C.: Spatial normalization of 3D brain images using deformable models. J. Comp. Assisted Tomography 20(4), 656-665 (1996)

5. Davies, R.H., et al.: A minimum description length approach to statistical shape modeling. IEEE Trans Med Imaging 21(5), 525-537 (2002)

6. Durrleman, S., et al.: Sparse approximation of currents for statistics on curves and surfaces. In: Med. Image Comp. Comput.-Assist. Intervention, Proceedings. vol. 11, pp. 390-398 (2008)

7. Fischl, B., et al.: High-resolution inter-subject averaging and a coordinate system for the cortical surface. Human Brain Mapping 8(4), 272-284 (1999)

8. Frisoni, G.B., et al.: The clinical use of structural MRI in Alzheimer disease. Nature Reviews Neurology 6(2), 67-77 (2010)

9. Hermosillo, G.: Ph.D. thesis, Université de Nice (INRIAROBOTVIS), Sophia Antipolis, France (2002)

10. Holland, D., et al.: Subregional neuroanatomical change as a biomarker for Alzheimer's disease. In: Proc. Natl. Acad. Sci. USA (2009)

11. Kim, B., et al.: Mutual information for automated unwarping of rat brain autora-diographs. Neuroimage 5(1), 31-40 (1997)

<table><caption>**Table 2.** FDR corrected *p*-values on hippocampal surfaces.</caption><thead><tr><th></th><th>Our Method</th><th>SPHARM [18]</th></tr></thead><tbody><tr><td>AD-CTL</td><td>0.0485</td><td>0.0459</td></tr><tr><td>MCI-AD</td><td>0.0215</td><td>0.0047</td></tr><tr><td>MCI-CTL</td><td>0.0483</td><td>0.0288</td></tr></tbody></table>
12. Lui, L.M., et al.: Computation of curvatures using conformal parameterization. Communications in Information and Systems 8(1), 1–16 (2008)

13. Meyer, C.R., et al.: Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using affine and thin-plate spline warped geometric deformations. Medical Image Analysis 1(3), 195–206 (1997)

14. Morra, J.H., et al.: Automated mapping of hippocampal atrophy in 1-year repeat MRI data from 490 subjects with Alzheimer’s disease, mild cognitive impairment, and elderly controls. Neuroimage 45(1 Suppl), S3–15 (2009)

15. Pitiot, A., et al.: Learning object correspondences with the observed transport shape measure. In: Proceedings of Information Processing in Medical Imaging (IPMI 2003). pp. 25–37 (2003)

16. Qiu, A., et al.: Regional shape abnormalities in mild cognitive impairment and Alzheimer’s disease. Neuroimage 45(3), 656–661 (2009)

17. Rueckert, D., et al.: Nonrigid registration using free-form deformations: application to breast MR images. IEEE Trans Med Imaging 18(8), 712–721 (1999)

18. Styner, M., et al: Framework for the statistical shape analysis of brain structures using SPHARM-PDM. The Insight Journal pp. 242–250 (2006)

19. Thompson, P.M., et al.: Mapping hippocampal and ventricular change in Alzheimer disease. Neuroimage 22(4), 1754–1766 (2004)

20. Vaillant, M., et al.: Surface matching via currents. In: Proceedings of Information Processing in Medical Imaging (IPMI 2005). vol. 19, pp. 381–392 (2005)

21. Wang, Y., et al.: Automated surface matching using mutual information applied to Riemann surface structures. In: Med. Image Comp. Comput.-Assist. Intervention, Proceedings, Part II. pp. 666–674 (2005)

22. Wang, Y., et al.: Mutual information-based 3D surface matching with applications to face recognition and brain mapping. In: Proc. Intl Conf. Computer Vision. pp. 527–534 (2005)

23. Wang, Y., et al.: Multivariate tensor-based brain anatomical surface morphometry via holomorphic one-forms. In: Med. Image Comp. Comput.-Assist. Intervention, Proceedings. vol. 12, pp. 337–344 (2009)

24. Wang, Y.: Multivariate tensor-based subcortical morphometry system (2011), http://gsl.lab.asu.edu/conformal.htm

25. West, J., et al.: Comparison and evaluation of retrospective intermodality brain image registration techniques. J. Comp. Assisted Tomography 21(4), 554–566 (1997)

26. Yushkevich, P.A.: Continuous medial representation of brain structures using the biharmonic PDE. Neuroimage 45(1 Suppl), S99–110 (2009)

27. Zhong, J., et al.: Multi-manifold diffeomorphic metric mapping for aligning cortical hemispheric surfaces. Neuroimage 49(1), 355–365 (2010)

<table><caption>**Table 3.** FDR corrected *p*-values on hippocampal surfaces (with failed subjects excluded).</caption><thead><tr><th></th><th>Our Method</th><th>SPHARM [18]</th></tr></thead><tbody><tr><td>AD-CTL</td><td>0.0486</td><td>0.0459</td></tr><tr><td>MCI-AD</td><td>0.0195</td><td>0.0047</td></tr><tr><td>MCI-CTL</td><td>0.0483</td><td>0.0288</td></tr></tbody></table>
<figure><img src="image_64.png" /><figcaption>**Fig. 3.** Comparison of surface fluid registration and the previous method on map of local shape differences (*p*-values), based on the mutivariate TBM method with hippocampal surfaces in 40 AD, 40 MCI, and 40 control subjects which were automatically segmented [14]. (a), (c), (e) are our results on group difference between AD and control, MCI and AD, MCI and control, respectively. Similarly, (b), (d), (f) are results of a prior work [23] on AD and control, MCI and AD, MCI and control, respectively. The *p*-map color scale is the same as Figure 2.</figcaption></figure>
<figure><img src="image_65.png" /><figcaption>**Fig. 4.** Comparison of surface fluid registration and SPHARM on map of local shape differences (*p*-values), based on the mutivariate TBM method with hippocampal surfaces from ADNI baseline data, which were automatically segmented by FIRST. (a), (c), (e) are our results on group difference between AD and control, MCI and AD, MCI and control, respectively, in 231 control, 199 AD, and 407 MCI subjects. Similarly, (b), (d), (f) are results of SPHARM on AD and control, MCI and AD, MCI and control, respectively, in 226 control, 182 AD, and 393 MCI subjects. The *p*-map color scale is the same as Figure 2.</figcaption></figure>
# A New Riemannian Setting for Surface Registration

Martin Bauer¹, Martins Bruveris²

¹ Fakultät für Mathematik, Universität Wien, Nordbergstrasse 15, A-1090 Wien

² Dep. of Mathematics, Imperial College, London SW7 2AZ, UK

**Abstract.** We present a new approach for matching regular surfaces in a Riemannian setting. We use a Sobolev type metric on deformation vector fields which form the tangent bundle to the space of surfaces. In this article we compare our approach with the diffeomorphic matching framework. In the latter approach a deformation is prescribed on the ambient space, which then drags along an embedded surface. In contrast our metric is defined directly on the deformation vector field and can therefore be called an *inner metric*. We also show how to discretize the corresponding geodesic equation and compute the gradient of the cost functional using finite elements.

**Keywords:** Registration, Surface Matching, LDDMM, Computational Anatomy, Geodesic Shooting, Adjoint Equations

## 1 Introduction

The field of computational anatomy concerns itself with the study and classification of the variability of biological shapes, including their statistical variance. Since the space of all shapes is inherently nonlinear, the usual methods of linear statistics cannot be applied. In particular, the addition of two surfaces cannot be meaningfully defined. One way to overcome this difficulty is to introduce a Riemannian structure on the space of shapes, which locally linearizes the space and allows the development of statistical methods that are analogous to the linear case. This approach was taken, e.g., in [7]. In the Riemannian setting, the average of two shapes may be defined as the middle point of a geodesic joining these two shapes. In a similar way one may define the corresponding geodesic mean of a collection of *n* shapes.

One class of shapes which are of interest in computational anatomy consists the surfaces embedded in $\mathbb{R}^3$. The cortical surface, the surfaces of hippocampi, thalami, and nasal cavities are all examples of shapes which are represented as two dimensional surfaces in $\mathbb{R}^3$. This is also an example, where the Riemannian setting may be applied to study collections of shapes.

Throughout the last decade various (Riemannian) metrics have been proposed. They include a metric that preserves local rigidity [11], a generalization of the elastic metric for curves to higher dimensions [18, 13], a metric inspired by a continuum mechanics which is defined in the interior of a two dimensional
shape [23, 8] and a metric based on the square-root representation of surfaces [12]. Other approaches include the representation of surfaces via densities [20] and metrics defined on surfaces via level sets [15, 19].

Another method for comparing anatomical shapes in the Riemannian setting is the method of large deformation diffeomorphic metric matching (LDDMM), based on the deformable template paradigm of Grenander [10]. In this setting, a template shape is matched to a target shape by finding a transformation in a suitable group of deformations of the ambient space that transforms the template into the target. This approach has been systematically developed in [3, 4, 9, 21, 22] and applied to various problems in computational anatomy. Registering two surfaces in this framework involves finding a diffeomorphism of the whole space, which transforms one surface to the other. Because of its widespread use in the field of computational anatomy for registering volumetric images we will use LDDMM as a reference to highlight the features of our proposed framework.

In this paper we propose a different way of defining a Riemannian structure on the space of surfaces, which also provides the full range of tools for nonlinear statistics. We use a Sobolev type norm to enforce regularity of the deformation vector field and measure the cost of the deformation. Our approach to characterizing a deformation is intrinsic to the surface, rather than resulting from a transformation of the surface induced by a deformation of the ambient space in which the image is embedded. For this reason, the Riemannian metrics used here are called *inner metrics* as opposed to the *outer metrics* used in LDDMM, where deformations are imposed via the ambient space. Other examples for inner metrics can be found in [18, 12, 23, 8, 1].

Inner metrics of Sobolev type on planar curves were introduced and studied previously in [17, 24]. Recently they were generalized to surfaces and higher dimensional hypersurfaces in Euclidean space in [2]. The numerical implementation of matching with these metrics differs from LDDMM, because the metric on the tangent space at each surface depends nonlinearly on the surface. This means the metric will change adaptively as one moves around in shape space. This adaptive property is in marked contrast to LDDMM, where the metric is defined on the diffeomorphisms and projected down to the shapes, so it doesn't depend on the particular shape.

The outline of the paper is as follows. In Sec. 2 we review the registration problem for surfaces and recall how it is solved using outer metrics in LDDMM. Then we present the approach via inner metrics of Sobolev type and point out the differences between the two methods. For definiteness, we will concentrate our attention on the Sobolev metric of order one. In Sec. 3 we discuss how to discretize and implement the geodesic equations for this metric and how to solve the registration problem via geodesic shooting. Finally, in Sec. 4 we show how this metric performs in some examples using synthetic data.
## 2 The Mathematical Formulation

We are dealing with the registration of parametrized regular surfaces. Such a surface is given by a smooth function $q: M \to \mathbb{R}^3$ from a model surface $M$ into the Euclidean space $\mathbb{R}^3$. We will consider different choices of the model surface $M$ in this paper: the plane sheet $M = [0, 1] \times [0, 1]$, the cylinder $M = S^1 \times [0, 1]$ and the torus $M = S^1 \times S^1$. Another interesting choice would be the sphere $M = S^2$, which however is not considered in this paper. The metric can be defined in the same way as for the other topologies, however the numerical treatment is more challenging, because the sphere cannot be covered by a single global coordinate chart. We require the parametrization of the surface $q$ to be regular in the following sense: at each point $x \in M$ the partial derivatives $\frac{\partial q}{\partial x^1}$, $\frac{\partial q}{\partial x^2}$ are required to be linearly independent. We will denote the space of all such surfaces by $\mathcal{S}$.

### 2.1 Registration with LDDMM

In the LDDMM framework, the registration of a template surface $q_0$ to a target surface $q_{\text{targ}}$ involves finding a curve $\varphi_t$ of diffeomorphisms of the ambient space $\mathbb{R}^3$, whose deformation carries the template surface to the target surface. Mathematically, one constructs these deformations using time dependent vector fields $v_t(y)$, which generate $\varphi_t$ as their flow, i.e.

$$ \partial_t \varphi_t = v_t(\varphi_t) . \qquad (1) $$

The registration problem consists of finding a vector field, which minimizes the following sum of a kinetic energy and a matching term

$$ E(v_t) = \frac{1}{2} \int_0^1 \|v_t\|_V^2 dt + \frac{1}{2\sigma^2} d(\varphi_1(q_0), q_{\text{targ}}) . \qquad (2) $$

The kinetic energy is usually measured using a norm $\|\cdot\|_V$ defined on a reproducing kernel Hilbert space $V$ of vector fields on $\mathbb{R}^3$ with kernel $K$. The norm is then given by $\|u\|_V^2 = \int u \cdot K^{-1} \star u dx$. We will discuss possible choices of the matching term in Section 2.5.

It is possible to reduce the complexity of the problem, since one can show that minimizing vector fields $v_t$ have to obey an evolution equation. This enables us to describe the whole vector field $v_t$ by knowing only its value $v_0$ at time $t = 0$. The equations

$$ \partial_t q_t(x) = v_t(q_t(x)) \qquad (3a) $$

$$ \partial_t p_t(x) = (Dv_t(q_t(x)))^T p_t(x) \qquad (3b) $$

$$ v_t(y) = \int_M p_t(x) K(y - q_t(x)) dx \qquad (3c) $$

are given in terms of a momentum $p_t$, which lives on the surface. The momentum is convolved with the kernel in (3c) to reconstruct the minimizing vector field,
which defines the deformation of the ambient and drags the surface along in (3a). Details about LDDMM and the reformulation as evolution equations can be found in [3, 25].

## 2.2 Registration with Inner Metrics

We propose to use a different approach, described from a mathematical point of view in [1, 2]. In this approach we describe the deformation of the surface directly, without assuming an underlying deformation of the whole space. In our approach we will replace (3a) by

$$ \partial_t q_t(x) = u_t(x), \qquad (4) $$

where $u_t(x) \in C^\infty(M, \mathbb{R}^3)$ is a time dependent vector field, defined only on the surface. Note the difference between the vector field $v_t(y)$ which is defined on $\mathbb{R}^3$ and $u_t(x)$, which is defined on the model space $M$, c.f. Fig. 1.

<figure><img src="image_66.png" /><figcaption>**Fig. 1.** For inner metrics the vector field governing the deformation is defined directly on the surface (left picture). In contrast the LDDMM approach defines a deformation vector field on all of $\mathbb{R}^3$ (right picture). The latter vector field deforms the ambient space and in the process induces a deformation of the surface.</figcaption></figure>

In the new framework the registration problem still consists of minimizing the following sum of a kinetic energy term and a matching functional

$$ E(u_t) = \frac{1}{2} \int_0^1 \langle u_t, u_t \rangle_{q_t} dt + \frac{1}{2\sigma^2} d(q_1, q_{\text{targ}}) . \qquad (5) $$

The kinetic energy is measured via an inner product $\langle \cdot, \cdot \rangle_{q_t}$ on the space of vector fields along the surface. This inner product can and usually will depend nonlinearly on the surface $q_t$ itself. This is another difference with the LDDMM framework, where the kinetic energy of the vector fields didn't depend on the surfaces that were matched. There is a whole variety of inner metrics, one can choose from. We will concentrate in this paper on an $H^1$-type metric, which will be discussed in Sec.2.4. Again we postpone the discussion of the matching term to Sec. 2.5.
## 2.3 Geometry

Both the LDDMM approach and the inner metrics can be seen as a special case of a constructions in the general framework of Riemannian geometry. In the case of inner metrics the collection of inner products $\langle \cdot, \cdot \rangle_q$ defines a Riemannian metric on the space $\mathcal{S}$ of all surfaces. In Riemannian geometry, curves, which minimize the energy $\int_0^1 \langle u_t, u_t \rangle_{q_t} dt$ for fixed endpoints $q_0, q_1$ are called geodesics. We can see that minimizers of (5) have to be geodesics in the space $\mathcal{S}$.

In LDDMM the inner product on the space of vector fields $V$ also defines a Riemannian metric, this time on the group of diffeomorphisms. Therefore the minima of the registration problem (2) generates geodesics in the diffeomorphism group.

Why is it advantageous to work in a Riemannian setting? In this setting, the minima of the matching energy are geodesics, so one may describe the nonlinear space of shapes $\mathcal{S}$ in terms of the initial velocity or momentum of a geodesic, which is an element in a linear vector space. This is possible, because geodesics obey an evolution equation like (3). Using the initial velocity, which encodes the whole solution of the registration problem, we are able to view the space of surfaces from the template surface $q_0$ as a linear space. This enables us to use statistics, compute average surfaces and measure distances.

## 2.4 $H^1$-type Metric on Surfaces

We will consider a metric on the space of surfaces, which is the analogue of the $H^1$-norm $\|f\|_{H^1}^2 = \int_{\mathbb{R}^2} |f(x)|^2 + |\nabla f(x)|^2 dx$ for functions on $\mathbb{R}^2$. We will replace functions on $\mathbb{R}^2$ by vector fields living on the curved surface $q$ and adjust the definition of the $H^1$-norm to take into account the curved nature of the surface $q$. Since $q$ is a surface in $\mathbb{R}^3$, we can measure angles and distances of vectors tangent to the surface, using the Euclidean inner product on $\mathbb{R}^3$. At each point $q(x)$ of the surface we also have a canonical basis for the plane tangent to $q$, given by the vectors $\frac{\partial q}{\partial x^1}$, $\frac{\partial q}{\partial x^2}$. We denote the inner product induced on the tangent plane to the surface by $g$. This inner product has the following coordinate matrix with respect to the basis $\frac{\partial q}{\partial x^1}$, $\frac{\partial q}{\partial x^2}$:

$$ (g_{ij}) = \left( \sum_{k=1}^{3} \frac{\partial q^k}{\partial x^i} \frac{\partial q^k}{\partial x^j} \right)_{i,j=1}^{2}, \quad (g^{ij}) = (g_{ij})^{-1}, \quad \mathrm{vol}(g) = \sqrt{\det(g_{ij})} . \quad (6) $$

We denote by $(g^{ij})$ the inverse matrix of $(g_{ij})$ and by $\mathrm{vol}(g)$ the volume density of the surface $q$ with respect to the metric $g$.
In this paper we will use the $H^1$-type inner metric on surfaces defined by the expression

$$ \langle u, v \rangle_q := \sum_{k=1}^{3} \int_M u^k v^k + \alpha^2 g (\mathrm{grad}^g(u^k), \mathrm{grad}^g(v^k)) \mathrm{vol}(g) \, dx \quad (7) $$

$$ = \sum_{k=1}^{3} \int_M u^k v^k + \alpha^2 g^{ij} \frac{\partial u^k}{\partial x^i} \frac{\partial v^k}{\partial x^j} \mathrm{vol}(g) \, dx . \quad (8) $$

One reason to use this generalization of the $H^1$-metric is, that this metric is invariant under reparametrizations of the surface and only depends on the image $q(M)$ as a subset of $\mathbb{R}^3$, in a similar way as the length of a curve in two dimensions only depends on the image of the curve, and not on a particular parametrization. This is necessary, if one wants to match unparametrized surfaces. For this task it is possible to use the same framework with this metric, only the matching term has to be chosen to be invariant under reparametrizations.

The constant $\alpha$, which appears in the metric, is a parameter, which has to be chosen for each problem. It represents the characteristic length scale, at which deformations take place. Another interpretation of $\alpha$ is the scale across which the momentum is smoothed, when passing from momenta to velocities. It can be compared with the kernel size in LDDMM.

Other choices for the metric are possible. One could use Sobolev type metrics of higher order as in [2] or multiply the components of the metric with a function depending on geometric quantities of the surface, like area, mean or Gaussian curvature as was done in [1].

## 2.5 The Matching Functional

There are different possible choices for the matching term. In this paper we will use the squared $L^2$-distance

$$ d(q_0, q_1) = \int_M |q_0(x) - q_1(x)|^2 dx. \quad (9) $$

Since we are dealing with parametrized surfaces this is a natural choice for the matching functional.

When matching unparametrized surfaces, natural choices of the matching functional would include currents, see [9], or one could use the reparametrization framework of [5].

## 3 Discretization

In this section we will describe how to discretize the variational problem (5) and compute the optimal path between two surfaces. Starting with an initial guess for the velocity $u_0^1$, we will use a gradient descent scheme

$$ u_0^{i+1} = u_0^i - \varepsilon \nabla_{u_0} E(u_0^i) $$
to converge towards the initial velocity of the optimal geodesic. The discretization thus consists of two parts:

- compute the geodesic, given the initial velocity to evaluate $E(u_0^i)$

- compute the gradient $\nabla_{u_0} E(u_0^i)$ to update the initial velocity.

We show how to discretize the geodesic equation in Sec. 3.1 and how to compute the gradient in Sec. 3.2.

## 3.1 The Geodesic Equation

We discretize the time-evolution of the surface $q(t)$ using the explicit Euler method

$$q_{i+1} = q_i + \Delta t u_i, \quad (10)$$

where $q_i = q(i\Delta t)$ is the discretized version of the curve and $\Delta t = 1/N$ is the time step, if we divide the interval $[0, 1]$ into $N$ parts. To compute $u_i$ we note that a geodesic is a critical point of the energy $E(u_i) = \frac{\Delta t}{2} \sum_{i=0}^{N-1} \langle u_i, u_i \rangle_{q_i}$, i.e. $\nabla_{u_i} E(u_i) = 0$. Following [6] we introduce the Lagrangian multiplier $p_i$ in the discrete variational principle

$$E(u_0, \dots, u_{N-1}) = \sum_{i=0}^{N-1} \frac{\Delta t}{2} \langle u_i, u_i \rangle_{q_i} + \langle p_i, q_{i+1} - q_i - \Delta t u_i \rangle_{L^2}. \quad (11)$$

and take variations. From variations in $u_i$ we see that $p_i$ is the momentum dual to the velocity $u_i$ in the sense that $\langle u_i, \delta u_i \rangle_{q_i} = \langle p_i, \delta u_i \rangle_{L^2}$ and we obtain the evolution equation for $u_i$ in the form

$$\langle u_{i+1}, \delta q_{i+1} \rangle_{q_{i+1}} = \langle u_i, \delta q_{i+1} \rangle_{q_i} + \Delta t \left\langle \frac{\delta \ell}{\delta q}(u_{i+1}, u_{i+1}; q_{i+1}), \delta q_{i+1} \right\rangle_{L^2}, \quad (12)$$

with $\delta q_{i+1}$ an arbitrary variation. Here we use the notation

$$\ell(u, v; q) = \frac{1}{2} \langle u, v \rangle_q. \quad (13)$$

We denote by $\frac{\delta \ell}{\delta q}$ the variational derivative of $\ell(u, v; q)$ with respect to the variable $q$, defined via

$$\left\langle \frac{\delta \ell}{\delta q}(u, v; q), \delta q \right\rangle_{L^2} = \lim_{h \to 0} \frac{\ell(u, v; q + h\delta q) - \ell(u, v; q)}{h}. \quad (14)$$

Equation (12) is an implicit time step for $u_{i+1}$, since $u_{i+1}$ appears on the right hand side in a quadratic term. To make computations easier and avoid having to solve a nonlinear equation, we changed the right hand side to an explicit Euler time step

$$\langle u_{i+1}, \delta q \rangle_{q_{i+1}} = \langle u_i, \delta q \rangle_{q_i} + \Delta t \left\langle \frac{\delta \ell}{\delta q}(u_i, u_i; q_i), \delta q \right\rangle_{L^2}. \quad (15)$$
## 3.2 Computing the Gradient

Given (10), (15) for the evolution of a geodesic we again use method of adjoint equations from [6] to compute the gradient of the energy with respect to the initial velocity. The resulting equations for the variables $\hat{u}_i, \hat{v}_i$ have to be integrated backwards in time

$$ 
\begin{align} 
\langle \hat{u}_i, \delta u_i \rangle_{q_i} &= \langle \hat{u}_{i+1}, \delta u_i \rangle_{q_i} + \Delta t \langle \hat{v}_{i+1}, \delta u_i \rangle_{q_{i+1}} + 2 \Delta t \langle \frac{\delta \ell}{\delta q}(u_i, \delta u_i; q_i), \hat{u}_{i+1} \rangle \nonumber \\
\langle \hat{v}_i, \delta q_i \rangle_{q_i} &= \langle \hat{v}_{i+1}, \delta q_i \rangle_{q_{i+1}} + 2 \left\langle \frac{\delta \ell}{\delta q}(\hat{u}_{i+1} - \hat{u}_i, u_i; q_i), \delta q_i \right\rangle \tag{16} \\
&\quad + \Delta t \left\langle \frac{\delta^2 \ell}{\delta q^2}(u_i, u_i; q_i), (\hat{u}_{i+1}, \delta q_i) \right\rangle \nonumber
\end{align} 
$$

with the initial conditions

$$ \hat{u}_N = 0 \qquad \langle \hat{v}_N, \delta q_N \rangle_{q_N} = -\frac{1}{\sigma^2} \langle q_N - q_{\text{targ}}, \delta q_N \rangle \quad (17) $$

at time $t = 1$. The gradient is then given by

$$ \nabla_{u_0} E(u_0) = u_0 - \hat{u}_0 . \quad (18) $$

## 4 Numerical Experiments

We implemented the geodesic and adjoint equations (10), (15) and (16) in Python using the finite element library FEniCS [14]. All model manifolds ($[0, 1] \times [0, 1]$, $S^1 \times [0, 1]$, $S^1 \times S^1$) were modelled on the rectangle $[0, 1] \times [0, 1]$ with periodic boundary conditions prescribed where necessary. The domain was subdivided into a regular triangular mesh, on which Lagrangian finite elements of order 1 were defined.

In the first example we apply our method to compute the geodesic path between two shapes, which includes both large and small deformations. The template shape is a straight cylinder of height 1 and radius 0.25, which is discretized using a regular triangular mesh of $2 \times 30 \times 30$ elements. The target shape is a cylinder, which is bent by 90° and has 5 small ripples added to it along the vertical axis. Compared to the bending the ripples constitute a small and local deformation of the shape. The target shape is discretized in the same way as the template. We use $\alpha = 0.6$ as the length scale parameter and 10 time steps for the time integration. The gradient descent takes 80 steps to converge to an $L^2$-error of 0.008. We can see in Fig. 2 that both the large and the small deformations are captured by the geodesic.

In the second example we want to illustrate the curved nature of shape space. To do so we pick three asymmetric tori, lying in different positions in space. Each two tori differ by a composition of two rigid rotations. We compute the geodesics between each pair of tori to measure the angles and side lengths of the triangle with the tori as vertices and the geodesics as edges. By comparing the sum of the
<figure><img src="image_67.png" /><figcaption>**Fig. 2.** Samples are shown from a geodesic in the space of surfaces between a straight cylinder and a bent cylinder with ripples, at time points $t = 0, 0.3, 0.5, 0.8, 1$. The color encodes the Euclidean length of the deformation vector field at each point of the surface.</figcaption></figure>

angles with $\pi$ one can estimate, whether the curvature of shape space along the plane containing the triangle is positive or negative. We measured $\alpha = 33.766^{\circ}$, $\beta = 34.802^{\circ}$ and $\gamma = 34.675^{\circ}$. The sum $\alpha + \beta + \gamma = 103.243^{\circ}$ is smaller than $180^{\circ}$, which indicates that the space is negatively curved in this area (c.f. [16, Sec. 5.4]). In negatively curved spaces geodesics tend to be attracted towards a common point. In this example the geodesics are attracted towards the surface, which is degenerated to a point. We can see in Fig. 3 that the midpoints of the geodesics between the vertices are slightly shrunk. This is another indication for the negatively curved nature of the space.

In the third example we show, that our framework is indeed capable to do nonlinear statistics on shape space. We generate five sample shapes and compute the mean shape between them. The five shapes are cylindrical vases with an open top and bottom, discretized again using a triangular mesh of $2 \times 30 \times 30$ elements. As the initial guess for the mean we use a straight cylinder. First we register this initial shape to the five target shapes, compute the average of the initial velocities and then shoot with this average velocity to obtain a next guess for the mean shape. We iterate this procedure until the average velocity is close to zero. This method of computing the Karcher mean was proposed in [7]. After four iterations we obtained an average velocity with norm 0.006. As can be seen in Fig. 4 the average shape indeed combines the characteristics of the five shapes.

# 5 Conclusions

In this paper we propose a metric to match regular surfaces in a Riemannian setting. Although this metric has been studied from a mathematical point of view in several papers, including [2, 17], so far it hasn't been applied to problems in computational anatomy. The aim of this work is to argue, that this is a promising approach, which is worthwhile to be studied further.
<figure><img src="image_68.png" /><figcaption>**Fig. 3.** This figure shows a geodesic triangle in the space of surfaces with asymmetric tori as vertices. The tori along the edges are the middle points of the geodesics connecting the vertices. One can see that the shapes tend to shrink along the geodesics before expanding again towards the ends. This effect implies negative curvature in this region of shape space.</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>**Fig. 4.** In this figure we show the Karcher mean of five vase-shaped objects. The mean shape, which is displayed in the center of the figure is computed using an iterated shooting method. The colored regions on the averaged shapes encode the Euclidean length of the initial velocity of the geodesic, which connects each shape to the mean. The color of the mean was chosen for artistic purposes only.</figcaption></figure>
In contrast to LDDMM, where surfaces are deformed via a deformation of the ambient space, in this approach the deformation is prescribed directly on the surface, while the ambient space stays constant. Because of this we call this approach matching with *inner metrics* as opposed to LDDMM, which can be described as matching with *outer metrics*. Other inner metrics, which have been proposed in the literature include [17, 12, 22, 7, 1] We show how to discretize the geodesic equation and how to compute the gradient of the matching functional with respect to the initial velocity. In the last part of the paper we present numerical results on synthetic data of different topologies, which demonstrate the versatility and applicability of our method.

At the present we applied this metric to match parametrized surfaces, which is an unwelcome restriction in applications. This is not a restriction of the framework itself, but only of the matching functional. By choosing a matching functional, which is independent of the parametrization of the surface, one can apply the same framework to unparametrized surfaces. In future work we plan to implement this capability and test the method on real anatomical data.

## Acknowledgments

We are enormously grateful to our friends and colleagues Colin Cotter, Phillip Harms, Darryl Holm, Peter Michor and David Mumford for their thoughtful comments and encouragement in the course of this work.

We also acknowledge partial support by the Royal Society of London Wolfson Award, the European Research Council Advanced Grant, the Imperial College London SIF Programme and the Austrian Science Fund (FWF).

## References

1. Bauer, M., Harms, P., Michor, P.W.: Almost local metrics on shape space of hypersurfaces in n-space. (Submitted. Preprint available at arXiv:1001.0717) (2010)

2. Bauer, M., Harms, P., Michor, P.W.: Sobolev metrics on shape space of surfaces in n-space. (To appear in Journal of Geometric Mechanics.) (2011)

3. Beg, M.F., Miller, M.I., Trouvé, A., Younes, L.: Computing large deformation metric mappings via geodesic flows of diffeomorphisms. Int. J. Comput. Vision 61(2), 139–157 (2005)

4. Bruveris, M., Gay-Balmaz, F., Holm, D., Ratiu, T.: The momentum map representation of images. Journal of Nonlinear Science 21, 115–150 (2011)

5. Cotter, C.J., Holm, D.D.: Geodesic boundary value problems with symmetry. J. Geom. Mech. 2(1), 51–68 (2010)

6. Dimet, F.X.L., Ngodock, H.E., Navon, I.M.: Sensitivity analysis in variational data assimilation. J. Meteorol. Soc. Japan pp. 145–155 (1997)

7. Fletcher, P., Venkatasubramanian, S., Joshi, S.: Robust statistics on riemannian manifolds via the geometric median. In: Computer Vision and Pattern Recognition, 2008. CVPR 2008. IEEE Conference on. pp. 1–8 (June 2008)

8. Fuchs, M., Jüttler, B., Scherzer, O., Yang, H.: Shape metrics based on elastic deformations. J. Math. Imaging Vis. 35, 86–102 (September 2009)
9. Glaunès, J., Qiu, A., Miller, M.I., Younes, L.: Large deformation diffeomorphic metric curve mapping. Int. J. Comput. Vision 80(3) (2008)

10. Grenander, U.: General Pattern Theory. Oxford Science Publications (1993)

11. Kilian, M., Mitra, N.J., Pottmann, H.: Geometric modeling in shape space. ACM Trans. Graph. 26 (July 2007)

12. Kurtek, S., Klassen, E., Ding, Z., Srivastava, A.: A novel riemannian framework for shape analysis of 3d objects. Computer Vision and Pattern Recognition, IEEE Computer Society Conference on 0, 1625–1632 (2010)

13. Liu, X., Shi, Y., Dinov, I., Mio, W.: A computational model of multidimensional shape. Int. J. Comput. Vision 89, 69–83 (August 2010)

14. Logg, A., Wells, G.N.: Dolfin: Automated finite element computing. ACM Trans. Math. Softw. 37, 20:1–20:28 (April 2010),

15. Bertalmio, M., Cheng, L-T., Osher, S., Sapiro, G.: Variational problems and par-
tial differential equations on implicit surfaces. Journal of Computational Physics
174(2), 759–780 (2001)

16. Michor, P.W., Mumford, D.: Riemannian geometries on spaces of plane curves. J. Eur. Math. Soc. (JEMS) 8 (2006), 1-48 (2006),

17. Michor, P.W., Mumford, D.: An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach. Appl. Comput. Harmon. Anal. 23(1), 74–113 (2007)

18. Mio, W., Srivastava, A., Joshi, S.: On shape of plane elastic curves. International Journal of Computer Vision 73, 307–324 (2007)

19. Mmoli, F., Sapiro, G., Osher, S.: Solving variational problems and partial differ-
ential equations mapping into general target manifolds. Journal of Computational
Physics 195(1), 263 – 292 (2004)

20. Peter, A., Rangarajan, A.: Maximum likelihood wavelet density estimation with applications to image and shape matching. Image Processing, IEEE Transactions on 17(4), 458–468 (april 2008)

21. Risser, L., Vialard, F.X., Murgasova, M., Holm, D., Rueckert, D.: Large deforma-
tion diffeomorphic registration using fine and coarse strategies. In: Fischer, B.,
Dawant, B., Lorenz, C. (eds.) Biomedical Image Registration, Lecture Notes in
Computer Science, vol. 6204, pp. 186–197. Springer Berlin / Heidelberg (2010)

22. Vialard, F.X., Risser, L., Rueckert, D., Cotter, C.J.: Diffeomorphic 3D image registration via geodesic shooting using an efficient adjoint calculation. (Submitted.) (2011)

23. Wirth, B., Bar, L., Rumpf, M., Sapiro, G.: A continuum mechanical approach to geodesics in shape space. IJCV 93(3), 293–318 (2011)

24. Younes, L., Michor, P.W., Shah, J., Mumford, D.: A metric on shape space with explicit geodesics. Atti Accad. Naz. Lincei Cl. Sci. Fis. Mat. Natur. Rend. Lincei (9) Mat. Appl. 19(1), 25–57 (2008)

25. Younes, L., Arrate, F., Miller, M.I.: Evolutions equations in computational anatomy. NeuroImage 45(1, Supplement 1), S40 – S50 (2009), mathematics in Brain Imaging
Author Index

Allassonière, Stéphanie, 111

Barysheva, Marina, 135
Bauer, Martin, 182
Bellani, Marcella, 148
Brambilla, Paolo, 148
Bruveris, Martins, 182

Castellani, Umberto, 148

de Bruijne, Marleen, 123
De Zubicaray, Greig I., 135
Dirksen, Asger, 123

Feng, Albert, 135
Feragen, Aasa, 123
Fletcher, P. Thomas, 75

Gilmore, John H., 87
Giraud, Christophe, 111
Gorbunova, Vladlena, 123
Gu, Xianfeng David, 38
Gunther, Andreas, 1

Hanaoka, Shouhei, 159
Hayashi, Naoto, 159
Ho, Jeffrey, 99

J. Sahn, David, 15
Jolivet, Pierre, 111
Joshi, Anand A., 135

Lamecker, Hans, 1
Lauze, François, 123
Lo, Pechin, 123
Lorenzi, Marco, 64

Marron, J.S., 87
Marsland, Stephen, 52
Martin, Nicholas G., 135
Masutani, Yoshitaka, 159
McLachlan, Robert, 52
McMahon, Katie L., 135
Modin, Klas, 52
Murino, Vittorio, 148

Nemoto, Mitsutaka, 159
Nielsen, Mads, 123
Nomura, Yukihiro, 159

Ohtomo, Kuni, 159

Pennec, Xavier, 64
Perlmutter, Matthew, 52
Prasad, Gautam, 135

Rambaldelli, Gianluca, 148
Reinhardt, Joseph M., 123

Seo, Dohyung, 99
Shi, Jie, 170
Shi, Rui, 38
Song, Xubo, 15
Srivastava, Anuj, 25
Styner, Martin, 87

Tansella, Michele, 148
Terzopoulos, Demetri, 135
Thompson, Paul M., 135, 170
Toga, Arthur W., 135

Ulaş, Aydın, 148

Vemuri, Baba, 99

Wang, Yalin, 38, 170
Wang, Yuting, 170
Weiser, Martin, 1
Wright, Margaret J., 135

Xie, Qian, 25

Yoshikawa, Takeharu, 159
Yoshioka, Naoki, 159
Yuan, Ying, 87

Zeng, Wei, 38
Zhang, Zhijun, 15
Zhu, Hongtu, 87