connections of a Lie group are:

$$ \tilde{\nabla}_{\tilde{\mathbf{X}}} \tilde{\mathbf{Y}} = 0 \qquad \text{Left} $$

$$ \tilde{\nabla}_{\tilde{\mathbf{X}}} \tilde{\mathbf{Y}} = \frac{1}{2} \left[ \tilde{\mathbf{X}}, \tilde{\mathbf{Y}} \right] \qquad \text{Symmetric} $$

$$ \tilde{\nabla}_{\tilde{\mathbf{X}}} \tilde{\mathbf{Y}} = \left[ \tilde{\mathbf{X}}, \tilde{\mathbf{Y}} \right] \qquad \text{Right} $$

## 2.2 Parallel Transport on Cartan connections

Once described the conditions for the one-parameters subgroups to be geodesics, it is natural to ask how to parallel transport along these paths, and each Cartan connection lead to a specific parallel transport method.

For the left Cartan connection, the unique fields that are covariantly constant are the *left-invariant* vector fields, and the parallel transport is induced by the left multiplication, i.e. $\Pi^L : T_pG \to T_qG$ is defined as

$$ \Pi^L(\mathbf{X}) = DL_{qp^{-1}}\mathbf{X} \qquad (4) $$

Conversely, the *right-invariant* vector fields are covariantly constant with respect to the right invariant connection. As above, the parallel transport is given by the differential of the right translation $\Pi^R(\mathbf{X}) = DR_{p^{-1}q}\mathbf{X}$. Finally, for the symmetric Cartan connection the parallel transport is given by the combination of the left and right transports. In fact it can be shown [5] that the parallel transport of $\mathbf{X}$ along the curve $\exp(t\mathbf{Y})$ is

$$ \Pi^S(\mathbf{X}) = DL_{\exp(\frac{1}{2}\mathbf{Y})} DR_{\exp(\frac{1}{2}\mathbf{Y})} \mathbf{X} \qquad (5) $$

## 3 Application to image registration

The Lie group theory is of relevant interest in the image registration context. For instance, the Lie group exponential has been already used for the diffeomorphic registration parametrized by stationary velocity fields. Of course, when moving to the infinite dimensional group of the diffeomorphisms, some caution is required for the generalization of the standard Lie theory and further research is still needed in order to quantify the impact of using such mathematical framework. However, the effectiveness of the SVF parametrization in terms of registration accuracy and computational efficiency encourage the adoption of the SVF as a valid instrument for the computational anatomy.

Given that one-parameter subgroups are geodesics for all the Cartan connections, we can implement the associated parallel transport for the SVF. Given the left and right actions on $\text{diff}(M)$,

$$ L_f g = f \circ g \qquad R_f g = g \circ f $$
we have

$$DL_f \simeq Df \cdot g \qquad DR_f \simeq g \circ f$$

We can therefore provide an explicit closed form formula for the parallel transport with respect to the canonical Cartan connections. In particular, if $X$ is a vector to be transported, and $\exp(tY)$ is the one-parameter subgroup we have:

$$\Pi_{\mathbf{Y}}^{L}(\mathbf{X}) = D\exp(\mathbf{Y}) \cdot \mathbf{X} \qquad (6)$$

$$\Pi_{\mathbf{Y}}^{R}(\mathbf{X}) = \mathbf{X} \circ \exp(\mathbf{Y}) \qquad (7)$$

$$\Pi_{\mathbf{Y}}^{S}(\mathbf{X}) = D\exp(\frac{\mathbf{Y}}{2}) \cdot (\mathbf{X} \circ \exp(\frac{\mathbf{Y}}{2})) \qquad (8)$$

**Remark 1.** The geodesics given by the Cartan connection are intrinsically different from the metric Riemannian ones, in the sense that the underlying connection is different from the Levi-Civita one. In particular the geodesics are not related to a positive definite from in $T\mathcal{G} \times T\mathcal{G}$. As consequence the space is not metrically complete, i.e. not all the elements of the space $\mathcal{G}$ might be reached by the one-parameter subgroups. The effect of such geometric property in the image registration context requires further investigation, in order to characterize the transformations that cannot be parametrized by SVF. However, we observe that in the image registration we are not interested in recovering “all” the possible diffeomorphisms, but only those which lead to admissible anatomical transformations.

**Remark 2.** From the computational point of view, we notice that among the three transport methods, $\Pi^R$ requires the simple resampling of the velocity field, while both $\Pi^L$ and $\Pi^S$ involve the computation of the Jacobian Matrix. The involvement of high order terms can raise accuracy problems, especially in case of noisy data and numerical approximations. We can alleviate the computational inaccuracy by taking advantage of the scaling properties of the one-parameter subgroups. Rather than directly compute the Jacobian $D\exp(\mathbf{Y})$, from the property

$$\begin{aligned} \exp(\mathbf{Y}) &= \exp(\frac{\mathbf{Y}}{2}) \circ \exp(\frac{\mathbf{Y}}{2}) = \\ &= \exp(\frac{\mathbf{Y}}{n}) \circ \dots \circ \exp(\frac{\mathbf{Y}}{n}). \end{aligned}$$

we can derive an iterative scheme for the Jacobian computation.

In fact, given a suitable first approximation $D\exp(\mathbf{Y})^{[0]} \simeq D\frac{\mathbf{Y}}{n}$ for an opportune scaling factor $n$, we have the iterative formula

$$D\exp(\mathbf{Y})^{[N+1]} = D\exp(\mathbf{Y})^{[N]}|_{\exp(\frac{\mathbf{Y}}{n})} \cdot D\exp(\frac{\mathbf{Y}}{n}) \qquad (9)$$

Thanks to the iterative scheme (9) the Jacobian is updated for a sufficient number of small steps along the one-parameter subgroup. Thus, the scheme avoids the computation of high order quantities on the final deformation field, that could
introduce biases due to the discretization inaccuracy. In fact, the derivatives here
are more robustly computed only for an opportunely scaled velocity field, and
the iterative formula evaluates the final Jacobian by successive resampling and
multiplications. Although the resampling scheme has an important impact on
the final computational accuracy, in the following it will be performed by simple
scalar interpolation.

4 Transport of longitudinal atrophy

**4.1 Synthetic experiment on a simplified geometry**

A synthetic progression of longitudinal atrophy was simulated on a simplified
geometry, represented by a 3D gray matter sphere S₀ enclosing a white/black
matter region. The atrophy was simulated by decreasing the gray layer thickness
on four subsequent time points to generate the sequence Sᵢ, i = 1 – 4 (Figure 2).
The longitudinal trajectories of deformation fields exp(Xᵢ) were then evaluated
by registering the images to the baseline with the Log-Demons algorithm [15].
The sequence of deformations exp(Xᵢ) was then transported on a target el-
lipsoidal geometry E₀ along the inter-subject deformation exp(Y) such that
exp(Y) * S₀ = E₀. The transport methods that we tested were:

- $\Pi^R$,

- $\Pi^L$ and $\Pi^S$ with the iterative scheme,

- the conjugate action $\text{Conj}(\exp(\mathbf{X}^i)) = \exp(\mathbf{Y})\exp(\mathbf{X}^i)\exp(\mathbf{Y})^{-1}$.

Moreover, the velocity fields $\mathbf{X}^i$ were transported with the Schild's Ladder, which operates along the “diagonal” inter-subject deformations $\exp(\mathbf{Y}^i)$ such that $\exp(\mathbf{Y}^i) * S_i = E_0$ (Figure 1).

The methods were quantitatively assessed by evaluating the features of in-
terest in the ellipsoid gray layer: the average $L^2$ Norm of the transported sta-
tionary velocity field and the Jacobian determinant, log-Jacobian determinant
and Elastic energy of the associated deformation fields. Since we are interested
in preserving the interesting features of the transported trajectories, the trans-
ported quantities were compared to the original values in the reference sphere
space. Moreover, the stability of the methods was tested by checking the scalar
spatial maps associated to the features.

**Results** Table 1 shows the accuracy of the transport methods in the preser-
vation of the measure of changes in the gray matter layer. Among the different
methods, the transport $\Pi^R$ was the most accurate in preserving the average
measures, while the Schild's Ladder performed better on the Log-Jacobian.

From the inspection of the related scalar log-Jacobian maps (Figure 2), the
transport $\Pi^L$ is the less stable and leads to noisy maps. Moreover, we notice
that the areas of expansions does not fit the boundary of the ellipsoid. On the
<figure><img src="image_28.png" /><figcaption>**Fig. 1.** Synthetic example: Intra and inter-subject variations from the sphere source space to the ellipsoid target space with related deformations.</figcaption></figure>

other hand, the transport $\Pi^R$ leads to smooth maps of changes, consistent with the target geometry, while the transport $\Pi^S$ lies “in between”, as one could reasonably expect. The Schild’s Ladder lead to smooth maps as well, although the inner spherical shape seems corrupted for higher deformations. This could explain the lower performance on the quantitative measurements for the time points 3 and 4. Finally, the log-Jacobian maps associated to the conjugate actions are smooth but fail to preserve the target ellipsoidal geometry, especially for the higher deformations.

<table><caption>**Table 1.** Average measures of changes on the gray matter layer. Top-row (Source Space): changes measured on the reference sphere at each time point 1–4. Bottom-rows: changes measured from the transported longitudinal deformations on the ellipsoid. For the conjugate action it was not possible to compute the $L^2$ Norm of the associated stationary velocity field, since it acts on deformation fields.</caption><thead><tr><th rowspan="2"></th><th colspan="4">L<sup>2</sup> Norm</th><th colspan="4">Log Jacobian</th><th colspan="4">Jacobian</th><th colspan="4">Elastic energy</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>Source Space</td><td>2.97</td><td>9.85</td><td>22.68</td><td>44.62</td><td>-4.77</td><td>-9.54</td><td>-14.76</td><td>-19.14</td><td>0.68</td><td>0.47</td><td>0.35</td><td>0.37</td><td>3.47</td><td>3.93</td><td>4.5</td><td>5.23</td></tr><tr><td>Π<sup>L</sup></td><td>3.02</td><td>9.57</td><td>22.14</td><td>42.32</td><td>-5</td><td>-9.82</td><td>-14.88</td><td>-20.43</td><td><strong>0.69</strong></td><td>0.51</td><td>0.43</td><td>0.45</td><td>3.51</td><td>4.01</td><td>4.67</td><td>5.53</td></tr><tr><td>Π<sup>R</sup></td><td><strong>2.94</strong></td><td><strong>10</strong></td><td><strong>22.81</strong></td><td><strong>44.58</strong></td><td>-4.70</td><td><strong>-9.36</strong></td><td>-14.51</td><td>-19.18</td><td><strong>0.69</strong></td><td><strong>0.49</strong></td><td><strong>0.36</strong></td><td><strong>0.37</strong></td><td><strong>3.49</strong></td><td><strong>3.9</strong></td><td><strong>4.44</strong></td><td><strong>5.15</strong></td></tr><tr><td>Π<sup>S</sup></td><td>3.3</td><td>11.17</td><td>25.7</td><td>50.37</td><td>-5.74</td><td>-11.2</td><td>-17.13</td><td>-23.65</td><td><strong>0.67</strong></td><td>0.50</td><td>0.42</td><td>0.48</td><td>3.58</td><td>4.2</td><td>4.99</td><td>6.05</td></tr><tr><td>Schild's Ladder</td><td>3.65</td><td>10.74</td><td>24.3</td><td>51.49</td><td><strong>-4.83</strong></td><td>-9.86</td><td><strong>-14.65</strong></td><td><strong>-19.11</strong></td><td>0.71</td><td>0.51</td><td>0.45</td><td>0.49</td><td>3.57</td><td>4.14</td><td>4.84</td><td>6.21</td></tr><tr><td>Conjugate</td><td>/</td><td>/</td><td>/</td><td>/</td><td>-2.6</td><td>-5.5</td><td>-9.18</td><td>-13.93</td><td>0.8</td><td>0.63</td><td>0.47</td><td>0.32</td><td>3.43</td><td>3.83</td><td>4.36</td><td>5.04</td></tr></tbody></table>
# 5 Conclusions

The study shows how the straightforward application of the Lie group theory to the diffeomorphic registration can lead to simple and efficient solutions for the transport of deformations. In particular the one-parameters subgroups are the geodesics with respect to the Cartan connections, and this mathematical setting leads to a closed form solution for the parallel transport. The geodesic of the Cartan connections generally differ from those of the Riemannian framework like the LDDMM, in the sense that they are not defined from a metric on the tangent space, and consequently the parallel transport is not related to the preservation of metric properties.

The present study highlights the trade-off between the choice of proper mathematical constructions and the related numerical implementation. In fact, among the parallel transports from the Cartan connections, the right one showed greater accuracy and smoothness, due to the simple computational requirements. However, the transport $\Pi^R$ operates according to a specific geometry corresponding to the right Cartan connection. In this case, we are working in a zero-curvature space with torsion, while from the theoretical point of view it might be preferable to work with respect to a symmetric connection which leads to torsion-free spaces. At this purpose, further studies are required in order to clarify the effects in the image registration context of imposing a specific connection, and for defining more robust numerical schemes for the computation of high order quantities.

*Acknowledgements* This work was partially supported by the ANR program *Karame-tria* ANR-09-BLAN-0332 from the French Agence Nationale pour la Recherche.

# References

1. V. Arsigny, O. Commowick, X. Pennec, and N. Ayache. A log-euclidean framework for statistics on diffeomorphisms. In R. Larsen, M. Nielsen, and J. Sporring, editors, *Medical Image Computing and Computer-Assisted Intervention - MICCAI*, volume 4190, page 924. Springer, Heidelberg, Sep 2006.

2. M.N. Bossa, E. Zacur, and S. Olmos. On changing coordinate systems for longitudinal tensor-based morphometry. *Spatio Temporal Image Analysis Workshop (STIA), MICCAI 2010, 2010.*

3. S. Durrleman, P. Fillard, X. Pennec, A. Trouve, and N. Ayache. Registration, atlas estimation and variability analysis of white matter fiber bundles modeled as currents. *NeuroImage*, 2011.

4. S. Durrleman, X. Pennec, A. Trouv, G. Gerig, and N. Ayache. Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets. In G. Yang, D. Hawkes, D. Rueckert, J Noble, and C. Taylor, editors, *Medical Image Computing and Computer-Assisted Intervention - MICCAI*, volume 5761. Springer, Heidelberg, Sep 2009.

5. S. Helgason. *Differential geometry, Lie groups, and symmetric spaces*. American Mathematical Soc, 1978.

6. M. Lorenzi, N. Ayache, and X. Pennec. Schilds ladder for the parallel transport of deformations in time series of images. In G. Székely and H. Hahn, editors, *Information Processing in Medical Imaging - IPMI*. Springer, Heidelberg, Jul 2011.
7. M. Lorenzi, G.B Frisoni, N. Ayache, and X. Pennec. Mapping the effects of aβ₁⁻₄₂ levels on the longitudinal changes in healthy aging: hierarchical modeling based on stationary velocity fields. In G. Fichtinger, A. Martel, and T. Peters, editors, *Medical Image Computing and Computer-Assisted Intervention - MICCAI*. Springer, Heidelberg, Sep 2011.

8. T. Mansi, X. Pennec, M. Sermesant, H. Delingette, and N. Ayache. Logdemons revisited: Consistent regularisation and incompressibility constraint for soft tissue tracking in medical images. In T. Jiang, N. Navab, J.P.W. Pluim, and M. Viergever, editors, *Medical Image Computing and Computer-Assisted Intervention - MICCAI*. Springer, Heidelberg, Sep 2010.

9. T. Mansi, I. Voigt, B. Leonardi, X. Pennec, S. Durrleman, M. Sermesant, H. Delingette, A.M. Taylor, Y. Boudjemline, G. Pongiglione, and N. Ayache. A statistical model for quantification and prediction of cardiac remodelling: Application to tetralogy of Fallot. *IEEE Transactions on Medical Images*, 2011.

10. M. Miller, A. Trouvé, and L. Younes. On the metrics and Euler-Lagrange equations of computational anatomy. *Annu Rev Biomed Eng*, 4(1):375–405, 2002.

11. M.M. Postnikov. *Geometry VI*. Springer, 2001.

12. A. Rao, R. Chandrashekara, G. Sanchez-Hortiz, R. Mohiaddin, P. aljabar, J. Hajnal, B. Puri, and D. Rueckert. Spatial trasformation of motion and deformation fields using nonrigid registration. *IEEE Transactions on Medical Imaging*, 23(9), 2004.

13. C. Seiler, X. Pennec, and M. Reyes. Geometry-Aware Multiscale Image Registration Via OBBTree-Based Polyaffine Log-Demons. In G. Fichtinger, A. Martel, and T. Peters, editors, *Medical Image Computing and Computer-Assisted Intervention - MICCAI*. Springer, Heidelberg, Sep 2011.

14. P. Thompson, K.M. Ayashi, G. Zubicaray, A.L. Janke, S.E. Rose, J. Semple, D. Herman, M.S. Hong, S.S. Dittmer, D.M. Dodrell, and A.W. Toga. Dynamics of gray matter loss in alzheimer's disease. *The Journal of Neuroscience*, 23:994–1005, 2003.

15. T. Vercauteren, X. Pennec, A. Merchant, and N. Ayache. Symmetric Log-domain diffeomorphic registration: A Demons-based approach. In D. Metaxas, L. Axel, G. Fichtinger, and G. Szekely, editors, *Medical Image Computing and Computer-Assisted Intervention - MICCAI*, volume 5241, pages 754–761. Springer, Heidelberg, Sep 2008.

16. L. Younes, A. Qiu, R. Winslow, and M. Miller. Transport of relational structures in groups of diffeomorphisms. *J Math Imaging Vis*, 32(1):41–56, Sep 2008.
<figure><img src="image_29.png" /><figcaption>**Fig. 2.** Top row: Spherical source and ellipsoidal target geometrical references. From top to bottom: Longitudinal atrophy sequence in the spherical space, associated log-Jacobian determinant scalar maps, and log-Jacobian determinant maps associated to the different methods of transport.</figcaption></figure>
# Geodesic Regression on Riemannian Manifolds

P. Thomas Fletcher

School of Computing, University of Utah

**Abstract.** This paper introduces a regression method for modeling the relationship between a manifold-valued random variable and a real-valued independent parameter. The principle is to fit a geodesic curve, parameterized by the independent parameter, that best fits the data. Error in the model is evaluated as the sum-of-squared geodesic distances from the model to the data, and this provides an intrinsic least squares criterion. Geodesic regression is, in some sense, the simplest parametric model that one could choose, and it provides a direct generalization of linear regression to the manifold setting. A hypothesis test for determining the significance of the estimated trend is also developed. While the method can be generally applied to data on any manifold, specific examples are given for a set of synthetically generated rotation data and an application to analyzing shape changes in the corpus callosum due to age.

## 1 Introduction

Regression analysis is a fundamental statistical tool for determining how a measured variable is related to one or more potential explanatory variables. The most widely used regression model is linear regression, due to its simplicity, ease of interpretation, and ability to model many phenomena. However, if the response variable takes values on a nonlinear manifold, a linear model is not applicable. Such manifold-valued measurements arise in many applications, including those involving directional data, transformations, tensors, and shape. For example, in biology and medicine it is often critical to understand processes that change the shape of anatomy. The difficulty is that shape variability is inherently high-dimensional and nonlinear. An effective approach to capturing this variability has been to parameterize shape as a manifold, or shape space.

Several works have studied the regression problem on manifolds. Jupp and Kent [6] propose an unrolling method on shape spaces. Regression analysis on the group of diffeomorphisms has been proposed as growth models by Miller [10], nonparametric regression by Davis, et al. [3], and second order splines by Trouvé and Vialard [12]. Finally, Shi, et al. [11] proposed a semiparametric model with multiple covariates for manifold response data. None of these methods provide a direct generalization of linear regression to manifolds. The purpose of this paper is to develop such a generalization, called geodesic regression, which models the relationship between an independent scalar variable with a dependent manifold-valued random variable as a geodesic curve. Like linear regression, the advantages of this model are its simplicity and ease of interpretation. As will be shown, the geodesic regression model also leads to a straightforward generalization of the $R^2$ statistic and a hypothesis test for significance of the estimated geodesic trend.
## 2 Multiple Linear Regression

Before formulating geodesic regression on general manifolds, we begin by reviewing multiple linear regression in $\mathbb{R}^n$. Here we are interested in the relationship between a non-random *independent* variable $X \in \mathbb{R}$ and a random *dependent* variable $Y$ taking values in $\mathbb{R}^n$. A multiple linear model of this relationship is given by

$$Y = \alpha + X\beta + \epsilon, \quad (1)$$

where $\alpha \in \mathbb{R}^n$ is an unobservable *intercept* parameter, $\beta \in \mathbb{R}^n$ is an unobservable *slope* parameter, and $\epsilon$ is an $\mathbb{R}^n$-valued, unobservable random variable representing the error. Geometrically, this is the equation of a one-dimensional line through $\mathbb{R}^n$ (plus noise), parameterized by the scalar variable $X$. For the purposes of generalizing to the manifold case, it is useful to think of $\alpha$ as the starting point of the line and $\beta$ as a velocity vector.

Given realizations of the above model, i.e., data $(x_i, y_i) \in \mathbb{R} \times \mathbb{R}^n$, for $i = 1, \dots, N$, the least squares estimates, $\hat{\alpha}, \hat{\beta}$, for the intercept and slope are computed by solving the minimization problem

$$(\hat{\alpha}, \hat{\beta}) = \arg \min_{(\alpha, \beta)} \sum_{i=1}^{N} \|y_i - \alpha - x_i\beta\|^2. \quad (2)$$

This equation can be solved analytically, yielding

$$\hat{\beta} = \frac{\frac{1}{N} \sum x_i y_i - \bar{x} \bar{y}}{\sum x_i^2 - \bar{x}^2}, \\ \hat{\alpha} = \bar{y} - \bar{x} \hat{\beta},$$

where $\bar{x}$ and $\bar{y}$ are the sample means of the $x_i$ and $y_i$, respectively. If the errors in the model are drawn from distributions with zero mean and finite variance, then these estimators are unbiased and consistent.

<figure><img src="image_30.png" /><figcaption>Fig. 1. Schematic of the geodesic regression model.</figcaption></figure>
## 3 Geodesic Regression

Let $y_1, \dots, y_N$ be points on a smooth Riemannian manifold $M$, with associated scalar values $x_1, \dots, x_N \in \mathbb{R}$. The goal of geodesic regression is to find a geodesic curve $\gamma$ on $M$ that best models the relationship between the $x_i$ and the $y_i$. Just as in linear regression, the speed of the geodesic will be proportional to the independent parameter corresponding to the $x_i$. Estimation will be set up as a least-squares problem, where we want to minimize the sum-of-squared Riemannian distances between the model and the data. A schematic of the geodesic regression model is shown in Figure 1.

Before formulating the model, we review a few basic concepts of Riemannian geometry. We will write an element of the tangent bundle as the pair $(p, v) \in TM$, where $p$ is a point in $M$ and $v \in T_pM$ is a tangent vector at $p$. Recall that for any $(p, v) \in TM$ there is a unique geodesic curve $\gamma$, with initial conditions $\gamma(0) = p$ and $\gamma'(0) = v$. This geodesic is only guaranteed to exist locally. When $\gamma$ is defined over the interval $[0, 1]$, the exponential map at $p$ is defined as $\text{Exp}_p(v) = \gamma(1)$. In other words, the exponential map takes a position and velocity as input and returns the point at time 1 along the geodesic with these initial conditions. The exponential map is locally diffeomorphic onto a neighborhood of $p$. Let $V(p)$ be the largest such neighborhood. Then within $V(p)$ the exponential map has an inverse, the Riemannian log map, $\text{Log}_p : V(p) \to T_pM$. For any point $q \in V(p)$ the Riemannian distance function is given by $d(p, q) = \|\text{Log}_p(q)\|$. It will be convenient to include the point $p$ as a parameter in the exponential and log maps, i.e., define $\text{Exp}(p, v) = \text{Exp}_p(v)$ and $\text{Log}(p, q) = \text{Log}_p(q)$.

Notice that the tangent bundle $TM$ serves as a convenient parameterization of the set of possible geodesics on $M$. An element $(p, v) \in TM$ provides an intercept $p$ and a slope $v$, analogous to the $\alpha$ and $\beta$ parameters in the multiple linear regression model (1). In fact, $\beta$ is a vector in the tangent space $T_\alpha \mathbb{R}^n \cong \mathbb{R}^n$, and thus $(\alpha, \beta)$ is an element of the tangent bundle $T\mathbb{R}^n$. Now consider an $M$-valued random variable $Y$ and a non-random variable $X \in \mathbb{R}$. The generalization of the multiple linear model to the manifold setting is the geodesic model,

$$ Y = \text{Exp}(\text{Exp}(p, Xv), \epsilon), \qquad (3) $$

where $\epsilon$ is a random variable taking values in the tangent space at $\text{Exp}(p, Xv)$. Notice that for Euclidean space, the exponential map is simply addition, i.e., $\text{Exp}(p, v) = p + v$. Thus, the geodesic model coincides with (1) when $M = \mathbb{R}^n$.

### 3.1 Least Squares Estimation

Consider a realization of the model (3): $(x_i, y_i) \in \mathbb{R} \times M$, for $i = 1, \dots, N$. Given this data, we wish to find estimates of the parameters $(p, v) \in TM$. First, define the sum-of-squared error of the data from the geodesic given by $(p, v)$ as

$$ E(p, v) = \frac{1}{2} \sum_{i=1}^{N} d(\text{Exp}(p, x_i v), y_i)^2. \qquad (4) $$
Following the ordinary least squares minimization problem given by (2), we formulate a least squares estimator of the geodesic model as a minimizer of the above sum-of-squares energy, i.e.,

$$ (\hat{p}, \hat{v}) = \arg \min_{(p,v)} E(p, v). \qquad (5) $$

Again, notice that this problem coincides with the ordinary least squares problem when $M = \mathbb{R}^n$.

Unlike the linear setting, the least squares problem in (5) for a general manifold $M$ will typically not yield an analytic solution. Instead we derive a gradient descent algorithm. Computation of the gradient of (4) will require two parts: the derivative of the Riemannian distance function and the derivative of the exponential map. Fixing a point $p \in M$, the gradient of the squared distance function is $\nabla_x d(p, x)^2 = -2 \text{Log}_x(p)$ for $x \in V(p)$.

<figure><img src="image_31.png" /><figcaption>Fig. 2. Jacobi fields as derivatives of the exponential map.</figcaption></figure>

The derivative of the exponential map $\text{Exp}(p, v)$ can be separated into a derivative with respect to the initial point $p$ and a derivative with respect to the initial velocity $v$. To do this, first consider a variation of geodesics given by $c_1(s, t) = \text{Exp}(\text{Exp}(p, su_1), tv(s))$, where $u_1 \in T_pM$ defines a variation of the initial point along the geodesic $\eta(s) = \text{Exp}(p, su_1)$. Here we have also extended $v \in T_pM$ to a vector field $v(s)$ along $\eta$ via parallel translation. This variation is illustrated on the left side of Figure 2. Next consider a variation of geodesics $c_2(s, t) = \text{Exp}(p, su_2 + tv)$, where $u_2 \in T_pM$. (Technically, $u_2$ is a tangent to the tangent space, i.e., an element of $T_v(T_pM)$, but there is a natural isomorphism $T_v(T_pM) \cong T_pM$.) The variation $c_2$ produces a “fan” of geodesics as seen on the right side of Figure 2.

Now the derivatives of $\text{Exp}(p, v)$ with respect to $p$ and $v$ are given by

$$ d_p \text{Exp}(p, v) \cdot u_1 = \left. \frac{d}{ds} c_1(s, t) \right|_{s=0} = J_1(1) $$

$$ d_v \text{Exp}(p, v) \cdot u_2 = \left. \frac{d}{ds} c_2(s, t) \right|_{s=0} = J_2(1), $$
where $J_i(t)$ are Jacobi fields along the geodesic $\gamma(t) = \text{Exp}(p, tv)$. Jacobi fields are solutions to the second order equation

$$ \frac{D^2}{dt^2} J(t) + R(J(t), \gamma'(t)) \gamma'(t) = 0, \qquad (6) $$

where $R$ is the Riemannian curvature tensor. For more details on the derivation of the Jacobi field equation and the curvature tensor, see for instance [2]. The initial conditions for the two Jacobi fields above are $J_1(0) = u_1$, $J'_1(0) = 0$ and $J_2(0) = 0$, $J'_2(0) = u_2$, respectively. If we decompose the Jacobi field into a component tangential to $\gamma$ and a component orthogonal, i.e., $J = J^\top + J^\perp$, the tangential component is linear: $J^\top(t) = u_1^\top + tu_2^\top$. Therefore, the only challenge is to solve for the orthogonal component.

Finally, the gradient of the sum-of-squares energy in (4) is given by

$$ \begin{aligned} \nabla_p E(p, v) &= -\sum_{i=1}^{N} d_p \text{Exp}(p, x_i v)^\dagger \text{Log}(\text{Exp}(p, x_i v), y_i), \\ \nabla_v E(p, v) &= -\sum_{i=1}^{N} x_i d_v \text{Exp}(p, x_i v)^\dagger \text{Log}(\text{Exp}(p, x_i v), y_i), \end{aligned} $$

where we have taken the adjoint of the exponential map derivative, e.g., defined by $\langle d_p \text{Exp}(p, v)u, w \rangle = \langle u, d_p \text{Exp}(p, v)^\dagger w \rangle$. As we will see in the next section, formulas for Jacobi fields and their respective adjoint operators can often be derived analytically for many useful manifolds.

## 3.2 $\mathbb{R}^2$ Statistics and Hypothesis Testing

In regression analysis the most basic question one would like to answer is whether the relationship between the independent and dependent variables is significant. A common way to test this is to see if the amount of variance explained by the model is high. For geodesic regression we will measure the amount of explained variance using a generalization of the $\mathbb{R}^2$ statistic, or coefficient of determination, to the manifold setting. To do this, we first define predicted values of $y_i$ and the errors $\epsilon_i$ as

$$ \begin{aligned} \hat{y}_i &= \text{Exp}(\hat{p}, x_i \hat{v}), \\ \hat{\epsilon}_i &= \text{Log}(\hat{y}_i, y_i), \end{aligned} $$

where $(\hat{p}, \hat{v})$ are the least squares estimates of the geodesic parameters defined above. Note that the $\hat{y}_i$ are points along the estimated geodesic that are the best predictions of the $y_i$ given only the $x_i$. The $\hat{\epsilon}_i$ are the residuals from the model predictions to the true data.

Now to define the total variance of data, $y_1, \dots, y_N \in M$, we use the Fréchet variance, intrinsically defined by

$$ \text{var}(y_i) = \min_{y \in M} \frac{1}{N} \sum_{i=1}^{N} d(y, y_i)^2. $$
The unexplained variance is the variance of the residuals, $\text{var}(\hat{\epsilon}_i) = \frac{1}{N} \sum \|\hat{\epsilon}_i\|^2$. From the definition of the residuals, it can be seen that the unexplained variance is the mean squared distance of the data to the model, i.e., $\text{var}(\hat{\epsilon}_i) = \frac{1}{N} \sum d(\hat{y}_i, y_i)^2$. Using these two variance definitions, the generalization of the $R^2$ statistic is then given by

$$R^2 = 1 - \frac{\text{unexplained variance}}{\text{total variance}} = 1 - \frac{\text{var}(\hat{\epsilon}_i)}{\text{var}(y_i)} \quad (7)$$

Fréchet variance coincides with the standard definition of variance when $M = \mathbb{R}^n$. Therefore, it follows that the definition of $R^2$ in (7) coincides with the $R^2$ for linear regression when $M = \mathbb{R}^n$. Also, because Fréchet variance is always nonnegative, we see that $R^2 \le 1$, and that $R^2 = 1$ if and only if the residuals to the model are exactly zero, i.e., the model perfectly fits the data. Finally, it is clear that the residual variance is always smaller than the total variance, i.e., $\text{var}(\hat{\epsilon}_i) \le \text{var}(y_i)$. This is because we could always choose $\hat{p}$ to be the Fréchet mean and $v = 0$ to achieve $\text{var}(\hat{\epsilon}_i) = \text{var}(y_i)$. Therefore, $R^2 \ge 0$, and it must lie in the interval $[0, 1]$, as is the case for linear models.

We now describe a permutation test for testing the significance of the estimated slope term, $\hat{v}$. Notice that if we constrain $v$ to be zero in (5), then the resulting least squares estimate of the intercept, $\hat{p}$, will be the Fréchet mean of the $y_i$. The desired hypothesis test is whether the fraction of unexplained variance is significantly decreased by also estimating $v$. The null hypothesis is $H_0: R^2 = 0$, which is the case if the unexplained variance in the geodesic model is equal to the total variance. Under the null hypothesis, there is no relationship between the $X$ variable and the $Y$ variable. Therefore, the $x_i$ are exchangeable under the null hypothesis, and a permutation test may randomly reorder the $x_i$ data, keeping the $y_i$ fixed. Estimating the geodesic regression parameters for each random permutation of the $x_i$, we can calculate a sequence of $R^2$ values, $R_1^2, \dots, R_m^2$, which approximate the sampling distribution of the $R^2$ statistic under the null hypothesis. Computing the fraction of the $R_k^2$ that are greater than the $R^2$ estimated from the unpermuted data gives us a $p$-value.

# 4 Results

## 4.1 Regression of 3D Rotations

**Overview of Unit Quaternions** We represent 3D rotations as the unit quaternions, $\mathbb{Q}_1$. A quaternion is denoted as $q = (a, v)$, where $a$ is the “real” component and $v = bi + cj + dk$. Geodesics in the rotation group are given simply by constant speed rotations about a fixed axis. Let $e = (1, 0)$ be the identity quaternion. The tangent space $T_e\mathbb{Q}_1$ is the vector space of quaternions of the form $(0, v)$. The tangent space at an arbitrary point $q \in \mathbb{Q}_1$ is given by right multiplication of $T_e\mathbb{Q}_1$ by $q$. The Riemannian exponential map is $\text{Exp}_q((0, v) \cdot q) = (\cos(\theta/2), 2v \cdot \sin(\theta/2)/\theta) \cdot q$, where $\theta = 2\|v\|$. The log map is given by $\text{Log}_q((a, v) \cdot q) = (0, \theta v / \|v\|) \cdot q$, where $\theta = \arccos(a)$.

Being a unit sphere, $\mathbb{Q}_1$ has constant sectional curvature $K = 1$. In this case the orthogonal component of the Jacobi field equation (6) along a geodesic $\gamma(t)$
<figure><img src="image_32.png" /><figcaption>Fig. 3. Results for simulated rotation data: MSE of the geodesic regression estimates for the intercept (left) and slope (right) as a function of sample size.</figcaption></figure>

has the analytic solution

$$J(t)^{\perp} = u_1(t) \cos(Lt) + u_2(t) \sin(Lt),$$

where $u_1, u_2$ are parallel vector fields along $\gamma$, with initial conditions $u_1(0) = J(0)^{\perp}$ and $u_2(0) = J'(0)^{\perp}$, and $L = \|\gamma'\|$. While the Jacobi field equation gives us the differential of the exponential map, we really need the adjoint of this operator for geodesic regression. However, from the above equation it is clear that $d_p \text{Exp}$ and $d_v \text{Exp}$ are both self-adjoint operators. That is, the above Jacobi field equation provides us both the differential and its adjoint.

**Geodesic Regression of Simulated Rotation Data** To test the geodesic regression least squares estimation on $\mathbb{Q}_1$, synthetic rotation data was simulated according to the geodesic model (3). The intercept was the identity rotation: $p = (1, 0, 0, 0)$, and the slope was a rotation about the z-axis: $v = (0, 0, 0, \pi/4)$. The $x_i$ data were drawn from a uniform distribution on $[0, 1]$. The errors in the model were generated from an isotropic Gaussian distribution in the tangent space, with $\sigma = \pi/8$. The resulting data $(x_i, y_i)$ were used to compute estimates of the parameters $(\hat{p}, \hat{v})$. This experiment was repeated 1,000 times each for sample sizes $N = 2^k$, $k = 1, \dots, 8$. We would expect that as the sample size increases, the mean squared error (MSE) in the estimates $(\hat{p}, \hat{v})$, relative to the true parameters, would approach zero. The MSE is defined as

$$\text{MSE}(\hat{p}) = \frac{1}{M} \sum_{i=1}^{M} d(\hat{p}_i, p)^2, \quad \text{MSE}(\hat{v}) = \frac{1}{M} \sum_{i=1}^{M} \|\hat{v}_i \cdot (\hat{p}_i^{-1} p) - v\|^2,$$

where $M = 1,000$ is the number of repeated trials, and $(\hat{p}_i, \hat{v}_i)$ is the estimate from the $i$th trial. Notice the multiplication by $(\hat{p}_i^{-1} p)$ in the second equation
is a right-translation of $\tilde{v}_i$ to the tangent space of $p$. Figure 3 shows plots of the resulting MSE for the slope and intercept estimates. As expected, the MSE approaches zero as sample size increases, indicating at least empirically that the least squares estimates are consistent.

## 4.2 Regression in Shape Spaces

One area of medical image analysis and computer vision that finds the most widespread use of Riemannian geometry is the analysis of shape. Dating back to the groundbreaking work of Kendall [7] and Bookstein [1], modern shape analysis is concerned with the geometry of objects that is invariant to rotation, translation, and scale. This typically results in representing an object's shape as a point in a nonlinear Riemannian manifold, or *shape space*. Recently, there has been a great amount of interest in Riemannian shape analysis, and several shape spaces for 2D and 3D objects have been proposed [5, 8, 9, 13]. We choose here to use Kendall's shape space, but geodesic regression is applicable to other shape spaces as well. It could also be applied to spaces of diffeomorphisms, using the Jacobi field calculations given by Younes [14].

**Overview of Kendall's Shape Space** We begin with derivations of the necessary computations for geodesic regression on Kendall's shape space. A configuration of $k$ points in the 2D plane is considered as a complex $k$-vector, $z \in \mathbb{C}^k$. Removing translation, by requiring the centroid to be zero, projects this point to the linear complex subspace $V = \{z \in \mathbb{C}^k : \sum z_i = 0\}$, which is equivalent to the space $\mathbb{C}^{k-1}$. Next, points in this subspace are deemed equivalent if they are a rotation and scaling of each other, which can be represented as multiplication by a complex number, $\rho e^{i\theta}$, where $\rho$ is the scaling factor and $\theta$ is the rotation angle. The set of such equivalence classes forms the complex projective space, $\mathbb{CP}^{k-2}$. As Kendall points out, there is no unique way to identify a shape with a specific point in complex projective space. However, if we consider that the geodesic regression problem only requires computation of exponential/log maps and Jacobi fields, we can formulate these computations without making an explicit identification of shapes with points in $\mathbb{CP}^{k-2}$.

Thus, we think of a centered shape $x \in V$ as representing the complex line $L_x = \{z \cdot x : z \in \mathbb{C}\setminus\{0\}\}$, i.e., $L_x$ consists of all point configurations with the same shape as $x$. A tangent vector at $L_x \in V$ is a complex vector, $v \in V$, such that $\langle x, v \rangle = 0$. The exponential map is given by rotating (within $V$) the complex line $L_x$ by the initial velocity $v$, that is,

$$ \mathrm{Exp}_x(v) = \cos \theta \cdot x + \frac{\|x\| \sin \theta}{\theta} \cdot v, \quad \theta = \|v\|. \quad (8) $$

Likewise, the log map between two shapes $x, y \in V$ is given by finding the initial velocity of the rotation between the two complex lines $L_x$ and $L_y$. Let $\pi_x(y) = x \cdot \langle x, y \rangle / \|x\|^2$ denote the projection of the vector $y$ onto $x$. Then the log map is given by

$$ \mathrm{Log}_x(y) = \frac{\theta \cdot (y - \pi_x(y))}{\|y - \pi_x(y)\|}, \quad \theta = \arccos \frac{|\langle x, y \rangle|}{\|x\|\|y\|}. \quad (9) $$
<figure><img src="image_33.png" /><figcaption>Fig. 4. Corpus callosum segmentation and boundary point model for one subject.</figcaption></figure>

Notice that we never explicitly project a shape onto $CP^{k-2}$. This has the effect that shapes computed via the exponential map (8) will have the same orientation and scale as the base point $x$. Also, tangent vectors computed via the log map (9) are valid only at the particular representation $x$ (and not at a rotated or scaled version of $x$). This works nicely for our purposes and implies that shapes along the estimated geodesic will have the same orientation and scale as the intercept shape, $\hat{p}$.

The sectional curvature of $CP^{k-2}$ can be computed as follows. Let $u, w$ be orthonormal vectors at a point $p \in CP^{k-2}$. These vectors may be thought of as vectors in $\mathbb{R}^{2k-2}$. Writing the vector $w$ as $w = (w_1, \dots, w_{2k-2})$, define the operator

$$j(w) = (-w_k, \dots, -w_{2k-2}, w_1, \dots, w_{k-1}).$$

(This is just multiplication by $i = \sqrt{-1}$ if we take $w$ as a complex vector with the $k-1$ real coordinates listed first.) Using this operator, the sectional curvature is given by

$$K(u, w) = 1 + 3\langle u, j(w) \rangle^2.$$

When $k=3$, $CP^1$ is the space of triangle shapes and is isomorphic to the sphere, $S^2$, and thus has constant sectional curvature, $K=1$. For $k>3$, $CP^{k-2}$ has sectional curvature in the interval $K \in [1, 4]$. Furthermore, let $u \in T_pCP^{k-2}$ be any unit length vector. If we decompose the tangent space into an orthonormal basis $e_1, \dots, e_{2k-2}$, such that $e_1 = j(u)$, then we have $K(u, e_1) = 4$ and $K(u, e_i) = 1$ for $i > 1$. This leads to the following procedure for computing the Jacobi field equation on $CP^{k-2}$ along a geodesic $\gamma$. Given initial conditions for $J(0)^\perp$ and $J'(0)^\perp$, decompose $J(0)^\perp = u_1 + w_1$, so that $u_1$ is orthogonal to $j(\gamma')$ and $w_1$ is tangential to $j(\gamma')$. Do the same for $J'(0)^\perp = u_2 + w_2$. As before, extend these vectors to parallel fields, $u_i(t), w_i(t)$, along $\gamma$. Then the orthogonal component of the Jacobi field along $\gamma$ is given by

$$J(t)^\perp = u_1(t) \cos(Lt) + u_2(t) \sin(Lt) + w_1(t) \frac{\cos(2Lt)}{2} + w_2(t) \frac{\sin(2Lt)}{2}.$$

As was the case for rotations, both $d_p$ Exp and $d_v$ Exp are self-adjoint operators.

**Application to Corpus Callosum Aging** The corpus callosum is the major white matter bundle connecting the two hemispheres of the brain. A midsagittal
<figure><img src="image_34.png" /><figcaption>Fig. 5. The input corpus callosum shape data and corresponding subject ages in years.</figcaption></figure>

slice from a magnetic resonance image (MRI) with segmented corpus callosum is shown in Figure 4. Several studies have shown that the volume of the corpus callosum decreases with normal aging [4]. However, less is known about how the *shape* of the corpus callosum changes with age. Understanding shape changes may provide a deeper understanding of the anatomical and biological processes underlying aging. For example, does the corpus callosum shrink uniformly in size, or do certain regions deteriorate faster than others? This type of question can be answered by geodesic regression in shape spaces.

To understand age-related changes in the shape of the corpus callosum, geodesic regression was applied to corpus callosum shape data derived from the OASIS brain database (www.oasis-brains.org). The data consisted of MRI from 32 subjects with ages ranging from 19-90 years old. The corpus callosum was segmented in a midsagittal slice using the ITK SNAP program (www.itksnap.org). These boundaries of these segmentations were sampled with 128 points using ShapeWorks (www.sci.utah.edu/software.html). This algorithm generates a sampling of a set of shape boundaries while enforcing correspondences between different point models within the population. An example of a segmented corpus callosum and the resulting boundary point model is shown in Figure 4. The entire collection of input shapes and their ages is shown in Figure 5 (boundary points have been connected into a boundary curve for visualization purposes). Each of these preprocessing steps were done without consideration of the subject age, to avoid any bias in the data generation.

Geodesic regression was applied to the data $(x_i, y_i)$, where $x_i$ was the $i$th subject's age, and $y_i$ was the $i$th subject's corpus callosum, generated as above and represented as a point in Kendall's shape space. First, the average age of the group, $\bar{x}$, was subtracted from each $x_i$, which was done to make the intercept
<figure><img src="image_35.png" /><figcaption>**Fig. 6.** Geodesic regression of the corpus callosum. The estimated geodesic is shown as a sequence of shapes from age 19 (blue) to age 90 (red).</figcaption></figure>

term correspond to the shape at the mean age, rather than the shape at age zero, which would be far outside the data range. Least squares estimates ($\hat{p}, \hat{v}$) were generated according to (5), and using the above calculations for $\mathbb{C}P^{k-2}$. The resulting estimated geodesic is shown in Figure 6 as a sequence of shapes: $\hat{\gamma}(t_k) = \text{Exp}(\hat{p}, (t_k - \bar{x})\hat{v})$, for $t_k = 19, 36, 54, 72, 90$. The shape trend shows a very clear thinning of the corpus callosum, with the largest effects in the posterior part of the body and in the genu (anterior end).

Finally, the statistical significance of the estimated trend was tested using the permutation test described in Section 3.2, using 10,000 permutations. The $p$-value for the significance of the slope estimate, $\hat{v}$, was $p = 0.009$. The coefficient of determination (for the unpermuted data) was $R^2 = 0.12$. The low $R^2$ value must be interpreted carefully. It says that age only describes a small fraction of the shape variability in the corpus callosum. This is not surprising: we would expect the intersubject variability in corpus callosum shape to be difficult to fully describe with a single variable (age). However, this does not mean that the age effects are not important. In fact, the low $p$-value says that the estimated age changes are highly unlikely to have been found by random chance.

# 5 Conclusion

We introduced a geodesic regression analysis method for Riemannian manifolds. The geodesic regression model is the natural generalization of linear regression and is parameterized by an intercept and slope term. We also developed a generalization of the $R^2$ statistic and a permutation test for the significance of the estimated geodesic trend. There are several avenues for future work. First, the hypothesis test presented here could be extended to test for group differences, for example, to test if age-related anatomical changes are different in a disease population compared to controls. Second, theoretical properties of geodesic regression, such as unbiasedness and consistency, would be of interest. Finally, regression diagnostics and model selection procedures need to be developed to assess the appropriateness of a geodesic model for a particular data set.

**Acknowledgments:** This work was supported by NSF CAREER Grant 1054057.
References

1. Bookstein, F.L.: Size and shape spaces for landmark data in two dimensions (with discussion). Statistical Science 1(2), 181–242 (1986)

2. do Carmo, M.: Riemannian geometry. Birkhäuser (1992)

3. Davis, B., Fletcher, P.T., Bullitt, E., Joshi, S.: Population shape regression from random design data. In: Proceedings of IEEE International Conference on Computer Vision (2007)

4. Driesen, N., Raz, N.: The influence of sex, age, and handedness on corpus callosum morphology: A meta-analysis. Psychobiology 23(3), 240–247 (1995)

5. Fletcher, P.T., Lu, C., Joshi, S.: Statistics of shape via principal geodesic analysis on Lie groups. In: IEEE CVPR. pp. 95–101 (2003)

6. Jupp, P.E., Kent, J.T.: Fitting smooth paths to spherical data. Applied Statistics 36(1), 34–46 (1987)

7. Kendall, D.G.: Shape manifolds, Procrustean metrics, and complex projective spaces. Bulletin of the London Mathematical Society 16, 18–121 (1984)

8. Klassen, E., Srivastava, A., Mio, W., Joshi, S.: Analysis of planar shapes using geodesic paths on shape spaces. IEEE PAMI 26(3), 372–383 (2004)

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

10. Miller, M.: Computational anatomy: shape, growth, and atrophy comparison via diffeomorphisms. NeuroImage 23, S19–S33 (2004)

11. Shi, X., Styner, M., Lieberman, J., Ibrahim, J., Lin, W., Zhu, H.: Intrinsic regression models for manifold-valued data. J Am Stat Assoc 5762, 192–199 (2009)

12. Trouvé, A., Vialard, F.X.: A second-order model for time-dependent data interpolation: Splines on shape spaces. In: MICCAI STIA Workshop (2010)

13. Younes, L.: Computable elastic distances between shapes. SIAM J. Appl. Math. 58, 565–586 (1998)

14. Younes, L.: Jacobi fields in groups of diffeomorphisms and applications. Quarterly of Applied Mathematics (2006)
Varying Coefficient Models for Modeling
Diffusion Tensors Along White Matter Bundles

Ying Yuan¹, Hongtu Zhu¹, J. S. Marron², John H. Gilmore³, and Martin
Styner⁴

¹ Dept Biostatistics, UNC Chapel Hill

² Dept Operations Research and Statistics, UNC Chapel Hill

³ Dept Psychiatry, UNC Chapel Hill

⁴ Dept Psychiatry and Computer Science, UNC Chapel Hill

**Abstract.** This paper develops a functional data analysis framework to model diffusion tensors along fiber bundles as functional responses with a set of covariates of interest, such as age, diagnostic status and gender. This framework has a wide range of clinical applications including the characterization of normal brain development, the neural bases of neuropsychiatric disorders, and the joint effects of environmental and genetic factors on white matter fiber bundles. A challenging statistical issue is how to appropriately handle diffusion tensors along fiber bundles as functional data in a Riemannian manifold. We propose a statistical model with varying coefficient functions, called VCTF to characterize the dynamic association between functional SPD matrix-valued responses and covariates. We calculate a weighted least squares estimation of the varying coefficient functions under the Log-Euclidean metric in the space of SPD matrices. We also develop a global test statistic to test specific hypotheses about these coefficient functions. Simulated data are further used to examine the finite sample performance of VCTF. We apply our VCTF to study potential gender differences and find statistically significant aspect of the development of diffusion tensors along the right internal capsule tract in a clinical study of neurodevelopment.

# 1 Introduction

Diffusion Tensor Imaging (DTI), which can track the effective diffusion of water in the human brain *in vivo*, has been widely used to map the microstructure and organization of fiber tracts and to assess the integrity of anatomical connectivity in white matter [1]. In DTI, the degree of diffusivity and the directional dependence of water diffusion in each voxel can be quantified by a 3 × 3 symmetric positive definite (SPD) matrix, called a diffusion tensor (DT), and its three eigenvalue-eigenvector pairs {($\lambda_k$, $v_k$) : k = 1, 2, 3} with $\lambda_1 \ge \lambda_2 \ge \lambda_3$. Fiber tracts in white matter can be constructed by consecutively connecting the principal directions ($v_1$) of DTs in adjacent voxels [2]. Therefore, DTs and tensor-derived quantities (e.g., fractional anisotropy (FA)) are distributed along these white matter fiber tracts for each subject. As an illustration, Figure 1 (a)
presents the right internal capsule tract and Figure 1 (b) presents DTs along
this tract obtained from 10 subject's, in which each DT is geometrically repre-
sented by an ellipsoid. In this representation, the lengths of the semiaxes of the
ellipsoid equal the square root of the eigenvalues of a DT, while the eigenvectors
define the direction of the three axes. Mathematically, these diffusion tensors
along the fiber tract are functionals of SPD matrices. Our research of interest is
to statistically model SPD functionals as responses with covariates of interest,
such as age and gender, across multiple subjects.

<figure><img src="image_36.png" /><figcaption>**Fig. 1.** (a) The right internal capsule tract. (b) The ellipsoidal representation of full tensors on the fiber tract from 10 selected subjects, colored with FA values.</figcaption></figure>

Statistical approaches have been developed for the statistical analysis of tensor-derived quantities along fiber tracts. A tract-based spatial statistics framework was developed to construct local diffusion properties along a white matter skeleton and then perform pointwise hypothesis tests at each grid point of the skeleton [3]. A model-based framework was developed to construct the medial manifolds of fiber tracts and then to test pointwise hypotheses based on diffusion properties along the medial manifolds [4]. However, since these two methods ignore the functional nature of diffusion properties along fiber tracts, they can suffer from low statistical power in detecting interesting features and in exploring variability in tract-based diffusion properties. A functional data analysis framework was used to compare a univariate diffusion property along fiber tracts across two (or more) populations for a single hypothesis test per tract by using functional principal component analysis and the Hotelling T² statistic [5]. Their method has two major limitations including only consideration of a univariate diffusion property and the lack of control for other covariates of interest, such as age. To address these two limitations, a functional regression framework was proposed to analyze multiple diffusion properties along fiber tracts as functional responses with a set of covariates of interest, such as age, diagnostic status and gender [6]. An alternative approach, called generalized functional linear models, was developed with a scalar outcome (e.g., diagnostic group) as responses and fiber bundle diffusion properties as varying covariate functions (or functional predictors) [7].
The calculated diffusion properties, which are nonlinear and linear functions of the estimated three eigenvalues of DT containing inherent bias, may be sub-stantially different from the true diffusion properties [8]. Numerical simulations have shown that estimates of the largest eigenvalue in a DT usually overesti-mate the true value of $\lambda_1$ and that estimates of the smallest eigenvalue usually underestimate $\lambda_3$. These differences between the estimated and true eigenval-ues subsequently bias the estimation of diffusion properties that are calculated from the values of these estimated eigenvalues. The sorting bias is pronounced in three types of degenerate DT including isotropic ($\lambda_1 = \lambda_2 = \lambda_3$), oblate ($\lambda_1 = \lambda_2 > \lambda_3$), or prolate ($\lambda_1 > \lambda_2 = \lambda_3$). Previous studies have shown that a major portion of DTs along fiber tracts are prolate tensors [9], and thus directly comparing these biased diffusion properties along fiber tracts can create ‘statis-tical artifacts’ including biased parameter estimates and incorrect test statistics and *p*-values for hypotheses of interest as shown in Section 3.

To avoid these statistical artifacts, it is important to directly analyze esti-mated DTs along fiber tracts. There are several advantages of comparing the es-timated DTs along fiber tracts with covariates. The first one is that the standard weighted least squared estimates of true DTs are almost unbiased [8]. Moreover, as shown in Section 3, directly modeling DTs along fiber tracts as a smooth SPD process allows us to incorporate smoothness constraint to further reduce noise in the estimated DTs along fiber tracts, which subsequently leads to reduced noise in estimated diffusion properties along the fiber tracts. Furthermore, the use of scalar diffusion properties ignores the direction information of DT, and thus it can lose the statistical power in detecting the differences in DT oriented in different directions.

There is a growing interest in the DTI literature in developing statistical methods for direct analysis of DTs in the space of SPD matrices. [10] proposed several parametric models for SPD matrices and derived the distributions of several test statistics for comparing differences between the means of the two (or multiple) groups of SPD matrices. [11] developed a nonparametric estimator for the common density function of a random sample of positive definite ma-trices. [12] developed a semi-parametric regression model with SPD matrices as responses in a Riemannian manifold and covariates in a Euclidean space. [13] and [14] proposed tensor splines and local constant regressions for interpolating DTI tensor fields based on the Riemannian metric.

In this paper, we propose a varying coefficient model (VCTF) to use varying coefficient functions to characterize the association between fiber bundle diffu-sion tensors and a set of covariates. This model is different from that in [12] because the former is applicable to DT-valued functional data and consider the within-subject correlations while the latter is only applicable to data with a DT response for each subject. Since the space of SPD matrices is a Riemanian man-ifold, to the best of our knowledge, our VCTF is the first paper for developing a functional data analysis framework for modeling functional manifold-valued responses with covariates in Euclidean space. To account for the curved na-ture of the SPD space, we employ the Log-Euclidean metric in [15] and then
use a weighted least squares estimation method based on the geodesic distance under the Log-Euclidean metric to estimate the varying coefficient functions. Furthermore, we develop global test statistics to test hypotheses on the varying coefficient functions and use a resampling method for approximating its p-value.

The rest of the paper is organized as follows. Section 2 presents VCTF and related statistical inference. Section 3 examines the finite sample performance of VCTF via simulation studies. Section 4 illustrates an application of VCTF in a clinical study of neurodevelopment. Section 5 presents concluding remarks.

## 2 Methodologies

In this section, we present our VCTF for the statistical analysis of DTs along fiber tracts as functional responses with a set of covariates. To compare DTs in populations of DTIs, we use the DTI atlas building followed by atlas fiberography and fiber parametrization as described in [5] to extract DTI fibers and establish DTI fiber correspondence across all DTI fiber correspondence across all DTI datasets from different subjects. For the sake of simplicity, we do not include these image processing steps here, which have been discussed in details in [5].

**Varying Coefficient Model for Functional SPD data** Let $Sym^+(3)$ and $Sym(3)$ be, respectively, the set of $3 \times 3$ SPD matrices and the set of $3 \times 3$ symmetric matrices with real entries. Let $vecs(C) = (c_{1,1}, c_{2,1}, c_{2,2}, \dots, c_{m_{1,1}}, \dots, c_{m_{1,m_1}})^T$ for any $m_1 \times m_1$ symmetric matrix $C = (c_{k,l})$. Let $Ivecs(\cdot)$ be the inverse operator of $vecs(\cdot)$ and $(a_l)$ be a $q \times 1$ vector with the $l$-th element $a_l$. Let $C \otimes D$ denote the Kronecker product of two matrices $C$ and $D$.

Let $x \in [0, L_0]$ be the arc length of any point on a specific fiber bundle relative to a fixed end point of the fiber bundle, where $L_0$ is the longest arc length on the fiber bundle. For the $i$-th subject, we measure a diffusion tensor, denoted by $S_i(x_j) \in Sym^+(3)$, at the arc length $x_j \in [0, L_0]$ for the $j$-th location grid point on the fiber bundle for $j = 1, \dots, n_G$ and $i = 1, \dots, n$, where $n_G$ and $n$ denote the numbers of grid points and subjects, respectively. We consider a varying coefficient model given as follows:

$$ \log(S_i(x)) = Ivecs((\mathbf{z}_i^T \beta_l(x))) + U_i(x) + \mathcal{E}_i(x) \quad \text{for } i = 1, \dots, n, \quad (1) $$

where $log(\cdot)$ denotes the matrix logarithm, $\mathcal{E}_i(x) \in Sym(3)$ is a $3 \times 3$ symmetric matrix of measurement errors, and $U_i(x) \in Sym(3)$ characterizes both individual matrix variations from $Ivecs((\mathbf{z}_i^T \beta_l(x)))$ and the correlation structure between $log(S_i(x))$ and $log(S_i(x'))$ for different $x$ and $x'$. Moreover, $\mathbf{z}_i$ and $\beta_l(x) = (\beta_{1l}(x), \dots, \beta_{rl}(x))^T$ are, respectively, a $r \times 1$ vector of covariates of interest with $z_{i,1} = 1$ and its associated vector of varying coefficient functions of $x$ for $l = 1, \dots, 6$. Model (1) can be regarded as a generalization of varying coefficient models, which have been widely studied and developed for longitudinal, time series, and functional data.
Let SP($\mu$, $\Sigma$) denote a stochastic process with mean $\mu(x)$ and covariance matrix function $\Sigma(x, x')$ for any $x, x' \in [0, L_0]$. It is also assumed that vecs($\mathcal{E}_i(x)$) and vecs($\mathcal{U}_i(x)$) are independent and respectively, independent and identical copies of SP($0$, $\Sigma_\mathcal{E}$) and SP($0$, $\Sigma_\mathcal{U}$). Moreover, vecs($\mathcal{E}_i(x)$) and vecs($\mathcal{E}_i(x')$) for $x \neq x'$ are assumed to be independent and thus $\Sigma_\mathcal{E}(x, x')$ takes the form of $\Sigma_\mathcal{E}(x)\mathbf{1}_{x=x'}$.

**Weighted Least Squares Estimation** To estimate the coefficient functions in $\beta(x) = (\beta_1^T(x), \dots, \beta_q^T(x))^T$, we develop a weighted least squares estimation method based on an adaptive local polynomial kernel (LPK) smoothing technique [16] and the geodesic distance under the Log-Euclidean metric (see [15] for details). Specifically, using Taylor's expansion, we can expand $\beta_l(x_j)$ at $x$ to obtain $\beta_l(x_j) \approx \beta_l(x) + \hat{\beta}_l(x)(x_j - x)$,

For a fixed bandwidth $h^{(1)}$, we first calculate a weighted least squares estimate of $\beta_l(x)$, denoted by $\hat{\beta}_l(x)$, which minimizes an objective function given by $\sum_{i=1}^n \sum_{j=1}^{n_G} K_{h^{(1)}}(x_j - x)\text{tr}\{[\log(S_i(x_j)) - \text{Ivecs}((\mathbf{z}_i^T(\hat{\beta}_l(x) + \hat{\beta}_l(x)(x_j - x))))]^\top 2\}$, where $K_{h^{(1)}}(\cdot) = K(\cdot/h^{(1)})/h^{(1)}$ is a rescaled kernel function, $K(\cdot)$ be a kernel function, such as the Gaussian and uniform kernels [16] and $\mathbf{a}^{\otimes 2} = \mathbf{a}\mathbf{a}^T$ for any vector or any matrix **a**. Then with some calculation, we can have

$$ \hat{\beta}_l(x) = (\hat{\beta}_{1k}(x), \dots, \hat{\beta}_{rl}(x))^T = [I_r \otimes (1,0)]A_l(x), \quad (2) $$

where $I_r$ is an $r \times r$ identity matrix, $A_l(x) = \Sigma(h^{(1)}, x)^{-1} \sum_{i=1}^n \sum_{j=1}^{n_G} K_{h^{(1)}}(x_j - x)[\mathbf{z}_i^T \mathbf{y}_{h^{(1)}}(x_j - x)](\log(S_i(x_j)))_l$ and $\Sigma(h^{(1)}, x) = \sum_{i=1}^n \sum_{j=1}^{n_G} K_{h^{(1)}}(x_j - x)[\mathbf{z}_i^{\otimes 2} \otimes \mathbf{y}_{h^{(1)}}(x_j - x)]^{\otimes 2}$ with $\mathbf{y}_{h^{(1)}}(x_j - x) = (1, (x_j - x)/h^{(1)})^T$.

We pool the data from all $n$ subjects and select an estimated bandwidth $h^{(1)}$, denoted by $\hat{h}_e^{(1)}$ by minimizing the cross-validation score given by $(nn_G)^{-1} \sum_{i=1}^n \sum_{j=1}^{n_G} \text{tr}\{[\log(S_i(x_j)) - \text{Ivecs}((\mathbf{z}_i^T \hat{\beta}_l(x_j, h^{(1)})(-i)))]^\top 2\}$, where $\hat{\beta}_l(x, h^{(1)})(-i)$ is the weighted least squares estimator of $\beta_l(x)$ for the bandwidth $h^{(1)}$ based on observed data with the observations from the $i$-th subject excluded. Finally, by substituting $\hat{h}_e^{(1)}$ into (2), we can obtain an estimate of $\beta_l(x)$, denoted by $\hat{\beta}_{l,e}(x)$. Combining all $\hat{\beta}_{l,e}(x)$ leads to $\hat{\beta}_e(x) = [\hat{\beta}_{1,e}(x), \dots, \hat{\beta}_{q,e}(x)]$.

**Smoothing Individual Functions and Estimating Covariance Matrices**

To simultaneously construct the individual function $\mathcal{U}_i(x)$, we also employ the local polynomial kernel smoothing technique. Specifically, using Taylor's expansion, we can expand $\mathcal{U}_i(x_j)$ at $x$ to obtain $\mathcal{U}_i(x_j) \approx \mathcal{U}_i(x) + \hat{\mathcal{U}}_i(x)(x_j - x)$. For each fixed $x$ and each bandwidth $h^{(2)}$, we calculate the weighted least square estimator of $\mathcal{U}_i(x)$, denoted as $\hat{\mathcal{U}}_i(x)$, by minimizing an objective function given by $\sum_{j=1}^{n_G} K_{h^{(2)}}(x_j - x)\text{tr}\{[\log(S_i(x_j)) - \text{Ivecs}((\mathbf{z}_i^T \hat{\beta}_{l,e}(x_j))) - \mathcal{U}_i(x) + \hat{\mathcal{U}}_i(x)(x_j - x)]^\top 2\}$. With some calculation, it can be shown that the weighted least square estimate of $\mathcal{U}_i(x)$, denoted by $\hat{\mathcal{U}}_i(x)$, given by

$$ \text{vecs}(\hat{\mathcal{U}}_i(x))^T = \sum_{j=1}^{n_G} \tilde{K}_{h^{(2)}}^0(x_j - x, x) (\text{vecs}(\log(S_i(x_j))) - (\mathbf{z}_i^T \hat{\beta}_{l,e}(x_j)))^T, \quad (3) $$
where $\tilde{K}_{h(2)}^0(x_j - x, x) = \Sigma_1(h^{(2)}, x)^{-1} K_{h(2)}(x_j - x) y_{h(2)}(x_j - x)$ is the empirical equivalent kernel and $\Sigma_1(h^{(2)}, x) = \sum_{j=1}^{n_G} K_{h(2)}(x_j - x) y_{h(2)}(x_j - x)^{\otimes 2}$ with $y_{h(2)}(x_j - x) = (1, (x_j - x)/h^{(2)})^T$

Let $R_i$ be a matrix with the j-th row vecs($\log(S_i(x_j))$) - $(\mathbf{z}_i^T \hat{\beta}_{l,e}(x_j))^T$ and $S$ be a $n_G \times n_G$ smoothing matrix with the $(i, j)$-th element $\tilde{K}_{h(2)}^0(x_j - x_i, x_i)$. We pool the data from all $n$ subjects and select an estimated bandwidth of $h^{(2)}$, denoted as $\hat{h}_e^{(2)}$, by minimizing the generalized cross-validation score given by $n^{-1} \frac{\sum_{i=1}^n \text{tr}\{(R_i - SR_i)^{\otimes 2}\}}{(1-n^{-1}\text{tr}(S))^2}$. Based on $\hat{h}_e^{(2)}$, we can use (3) to estimate vecs($\mathcal{U}_i(x)$) and $\mathcal{U}_i(x)$, denoted by vecs($\hat{\mathcal{U}}_{i,e}(x)$) and $\hat{\mathcal{U}}_{i,e}(x)$, respectively, for all $i$. After obtaining $\hat{\mathcal{U}}_{i,e}(x)$, we can estimate the mean function $\mathcal{U}(x)$ and the covariance function $\Sigma_{\mathcal{U}}(x, x')$ by using their empirical counterparts.

We construct a nonparametric estimator of the covariance matrix $\Sigma_{\mathcal{E}}(x, x)$ as follows. Let $\hat{\mathcal{E}}_i(x_j) = \log(S_i(x_j)) - \text{Ivecs}((\mathbf{z}_i^T \hat{\beta}_{l,e}(x_j))) - \hat{\mathcal{U}}_{i,e}(x_j)$ be estimated residuals for $i = 1, \dots, n$ and $j = 1, \dots, n_G$. We consider the kernel estimate of $\Sigma_{\mathcal{E}}(x, x)$ given by

$$ \hat{\Sigma}_{\mathcal{E}}(x, x) = (n-q)^{-1} \sum_{i=1}^{n} \sum_{j=1}^{n_G} \frac{K_{h^{(3)}}(x_j - x) \text{vecs}(\hat{\mathcal{E}}_i(x_j))^{\otimes 2}}{\sum_{j=1}^{n_G} K_{h^{(3)}}(x_j - x)}. \quad (4) $$

Let $\tilde{\Sigma}_{\mathcal{E}}(x_j, x_j) = (n-q)^{-1} \sum_{i=1}^n \text{vecs}(\hat{\mathcal{E}}_i(x_j))^{\otimes 2}$. To select an estimated bandwidth $h^{(3)}$, denoted by $\hat{h}_e^{(3)}$, we minimize the cross-validation score given by $(nn_G)^{-1} \sum_{i=1}^n \sum_{j=1}^{n_G} \text{tr}\{[\text{vecs}(\hat{\mathcal{E}}_i(x_j))^{\otimes 2} - \hat{\Sigma}_{\mathcal{E}}(x_j, x_j, h^{(3)})^{(-i)}]^{\otimes 2} \hat{\Sigma}_{\mathcal{E}}(x_j, x_j)^{-1}]^{\otimes 2}\}$, where $\hat{\Sigma}_{\mathcal{E}}(x, x, h^{(3)})^{(-i)}$ is the weighted least squares estimator of $\hat{\Sigma}_{\mathcal{E}}(x, x)$ based on observed data with the observations from the $i$-th subject excluded. Based on $\hat{h}_e^{(3)}$, we can use (4) to estimate $\Sigma_{\mathcal{E}}(x, x)$, denoted by $\hat{\Sigma}_{\mathcal{E},e}(x, x)$.

**Hypothesis Test** In neuroimaging studies, some scientific questions require the comparison of fiber bundle diffusion tensors along fiber bundles across two (or more) diagnostic groups and the assessment of the development of fiber bundle diffusion properties along time. Such questions can often be formulated as linear hypotheses of $\beta(x)$ as follows:

$$ H_0: \mathbf{R}\beta(x) = \mathbf{b}_0(x) \quad \text{for all } x \text{ vs.} \quad H_1: \mathbf{R}\beta(x) \neq \mathbf{b}_0(x), \quad (5) $$

where $\mathbf{R}$ is a $t \times 6r$ matrix of full row rank and $\mathbf{b}_0(x)$ is a given $t \times 1$ vector of functions.

We propose both local and global test statistics. The local test statistic can identify the exact location of significant grid point on a specific tract. At a given grid point $x_j$ on a specific tract, we test the local null hypothesis $H_0(x_j): \mathbf{R}\beta(x_j) = \mathbf{b}_0(x_j)$ against $H_1(x_j): \mathbf{R}\beta(x_j) \neq \mathbf{b}_0(x_j)$. We use a local test statistic $T_n(x_j)$ defined by

$$ T_n(x_j) = n\mathbf{d}(x_j)^T \{\mathbf{R}(\hat{\Sigma}_{\mathcal{U}}(x_j, x_j) \otimes \hat{\Omega}_Z^{-1})\mathbf{R}^T\}^{-1}\mathbf{d}(x_j), \quad (6) $$
where $\hat{\Omega}_Z = n^{-1} \sum_{i=1}^n z_i^{\otimes 2}$ and $d(x) = R(\hat{\beta}_o(x)^T - \text{bias}(\hat{\beta}_o(x)^T)) - b_0(x)$. Following [17], a smaller bandwidth leads to a small value of $\text{bias}(\hat{\beta}_o(x))$. Moreover, according to our simulation studies below, we have found that the effect of dropping bias($\hat{\beta}_o(x)$) is negligible and therefore, we drop it from now on.

We test the null hypothesis $H_0: R\beta(x) = b_0(x)$ for all $x$ using a global test statistic $T_n$ defined by

$$T_n = n \int_0^{L_0} d(x)^T [R(\hat{\Sigma}_U(x, x) \otimes \hat{\Omega}_X^{-1}) R^T]^{-1} d(x) dx. \quad (7)$$

It follows from Theorem 1 [18] along with the continuous mapping theorem that as both $n$ and $n_G$ converges to infinity, we have $T_n$ converges to some distribution (weighted $\chi^2$). Based on this result, we develop a wild bootstrap method to approximate the p-value of $T_n$.

## 3 Simulation Studies

We conducted two sets of Monte Carlo simulations. The first set of simulations was to evaluate the Type I and II error rates of the global test statistic $T_n$. The second set was to compare the power in detecting the group effect using either whole diffusion tensor or the diffusion properties.

**Simulation 1** In the first set of simulations, we evaluated the Type I and II error rates by simulating diffusion tensors along the right internal capsule tract (Figure 1 (a)) according to $S_i(x) = \exp(\text{Ivecs}((\beta_{1l}(x) + \beta_{2l}(x) \times G_i + \beta_{3l}(x) \times Gage_i)) + U_i(x) + \mathcal{E}_i(x))$, where $Gage_i$ and $G_i$, respectively, denote the gestational age at the scan time and gender of the $i$-th infant, $\text{Ivecs}(U_i(x))$ is a Gaussian process with zero mean and covariance matrix $\Sigma_U(x, x')$ and $\text{Ivecs}(\mathcal{E}_i(x))$ is a Gaussian random vector with zero mean and covariance matrix $\Sigma_{\mathcal{E}}(x, x) \mathbf{1}(x = x')$. To mimic imaging data, we used the diffusion tensors along the right internal capsule tract from all the 96 infants in our clinical data to estimate $\hat{\beta}(x)$ of $\beta(x)$ via equation (2), $\hat{U}(x)$ of $U(x)$ via equation (3), and $\hat{\mathcal{E}}(x)$ of $\mathcal{E}(x)$ via $\hat{\mathcal{E}}(x) = \log(S_i(x)) - \text{Ivecs}((\hat{\beta}_l(x)^T \mathbf{z})) - \hat{U}(x)$. We fixed all the parameters at their values obtained from our clinical data, except that we assumed $(\beta_{31}(x), \dots, \beta_{36}(x)) = c(\hat{\beta}_{31}(x), \dots, \hat{\beta}_{36}(x))$, where $c$ is a scalar specified below and $(\hat{\beta}_{31}(x), \dots, \hat{\beta}_{36}(x))$ were estimators obtained from our clinical data. Figure 2 displays the estimated diffusion tensors using VCTF method when $c = 1$. Note that the method did an excellent job of recovering ground truth.

In neuroimaging studies, some scientific questions require the assessment of the development of fiber bundles diffusion tensors across time. In this simulation study, the questions were formulated as the hypotheses test $H_0: \beta_{31}(x) = \dots = \beta_{36}(x) = 0$ for all $x$ along the right internal capsule tract against $H_1: \beta_{31}(x) \neq 0$ for at least one $x$ on the tract for some $l = 1, \dots, 6$. We first assumed $c = 0$ to assess the Type I error rates for the global test statistic $T_n$, and then we assumed $c = .2, .4, .6$, and $.8$ to examine the Type II error rates for $T_n$ at different effect
<figure><img src="image_37.png" /><figcaption>**Fig. 2.** Ellipsoidal representations of the true (a), simulated (b) and estimated (c) diffusion tensors along the right internal capsule tract, colored with FA values.</figcaption></figure>

sizes. To evaluate the Type I and II error rates at different sample sizes, we let $n = 96$. The values of gender and gestational age were set the same as all the 96 infants in our clinical study. Note that the number of grid points on the right internal capsule equals $n_G = 112$ for both cases.

We applied the VCTF procedure to the simulated diffusion tensors. Particularly, we approximated the *p*-value of $T_n$ using the wild bootstrap method described in the hypothesis test section. For each simulation, the significance levels were set at $\alpha = .05$ and .01, and 100 replications were used to estimate the rejection rates. For a fixed $\alpha$, if the Type I rejection rate is smaller than $\alpha$, then the test is conservative, whereas if the Type I rejection rate is greater than $\alpha$, then the test is anticonservative, or liberal.

Figure 3(a) displays the rejection rates for $T_n$ based on the resampling method for sample size 96) and all effect sizes ($c = 0, .2, .4, .6$, or $.8$) at both significance levels ($\alpha = .01$ or .05) using full diffusion tensors. The statistical power for rejecting the null hypothesis increases with the effect size and the significance level, which is consistent with our expectation.

**Simulation 2** In the second set of simulations, the diffusion tensors along the fiber tract were generated as in the first set of simulation. For each simulated diffusion tensor, we calculated its FA and MD values. Then we applied the VCTF procedure for multiple measures to the simulated values of FA, MD, joint FA and MD, respectively, and then tested the significance of the gestational age effect. It is observed from Fig. 3 (a)-(c) that the statistical power is much higher when we use the whole tensor instead of FA and MD. The power is a little higher at small effect size when we use the joint FA and MD values instead of full diffusion tensors while at higher effect size, they almost have the same power (Fig. 3(a) and (d)). However, the Type I error is greater than the .05 significance level when we use the joint FA and MD values (Fig. 3(d)). It means that the test is liberal. As mentioned in Section 1, this is because there is some inherent bias in these calculated diffusion properties.
<figure><img src="image_38.png" /><figcaption>**Fig. 3.** Simulation study: Type I and Type II error rates. Rejection rates of $T_n$ based on the resampling method are calculated at five different values of $c$ for sample sizes of 96(solid lines) subjects at the (a) .05 and (b) .01 significance levels using diffusion tensor, FA values, MD values, joint values of FA and MD.</figcaption></figure>

## 4 A Real Example

We investigate early brain development by using DTI and our VCTF. We consider 96 healthy infants (36 males and 60 females) whose mean gestational age was 245.6 days with SD: 18.5 days (range: 192-270 days). A 3T Allegra head only MR system was used to acquire all the images. The system was equipped with a maximal gradient strength of 40 mT/m and a maximal slew rate of 400 mT/(m·msec). The DTI images were obtained by using a single shot EPI DTI sequence (TR/TE=5400/73 msec) with eddy current compensation. The six non-collinear directions at the b-value of 1000 s/mm² with a reference scan (b=0) was applied. To improve the signal-to-noise ratio of the images, a total of five scans were acquired and averaged. A weighted least square estimation method [1] was used to construct the diffusion tensors. Then DTI atlas building followed by atlas based tractography procedure was employed to process all 96 DTI datasets. A nonlinear fluid deformation based highdimensional, unbiased atlas computation method is used to carry out a large deformation non-linear registration [19]. We chose the right internal capsule tract to illustrate the applicability of our method. Diffusion tensors were extracted along this fiber tract for all the 96 infants [5].

In this study, we have two specific aims. The first one is to compare diffusion tensors along the selected fiber bundle across the male and female groups and thus illuminate the gender effect on the development of these fiber bundle diffusion tensors. The second one is to delineate the development of fiber bundle diffusion tensor across the gestational age effect. To statistically test the effects, we applied our VCTF to diffusion tensors along the fiber tract. For the selected tract, we fitted the VCTF model (1) to the diffusion tensors from all 96 subjects, in which $z = (1, gender, Gage)^T$. Then, we used equation (2) to estimate the functional coefficients $\beta(x)$. For the hypothesis testing, we constructed the global test statistic $T_n$ via equation (7) to test the gender and age effects for the diffusion tensors. The $p$ value of $T_n$ was approximated using the resampling
method with $G = 1,000$ replications. Furthermore, we fitted VCTF model to FA, MD, joint FA and MD, each of three eigenvalues and three eigenvalues together values along the fiber tract.

In this study, we statistically test the effects of gender and gestational age on the diffusion tensors along the right internal capsule tract. To test the gender effect, we calculated the local test statistics $T_n(x_j)$ and their corresponding $p$ values across all grid points on the right internal capsule tract. It is observed from Figure 4 (a) that most grid points do not have significant $-\log_{10}(p)$ values, which are less than 1.3. Then, we also computed the $p$-value of the global test statistic $T_n$, $p = .335$ indicating no gender effect. It is observed from Figure 4 (b) that the $-\log_{10}(p)$ values of $T_n(x_j)$ for testing the gestational age effect at some grid points of the right tail are greater than 1.5 while a very high significant gestational age effect was found with the $p$-value of the global test statistic, $p = .01$. This indicates that diffusion tensors along the right internal capsule tract do not differ significantly between male and female groups but are significantly associated with the gestational age. We picked a grid point with the significant $p$ value of $T_n(x)$ and observed that the diffusion tensors become more spherical with the gestational age (Figure 5). This indicates a decreasing pattern of diffusion.

<figure><img src="image_39.png" /><figcaption>**Fig. 4.** The $-\log_{10}(p)$ values of test statistics $T_n(x_j)$ for testing gender (a) or gestational age (b) effect of diffusion tensors on the right internal capsule tract.</figcaption></figure>

# 5 Discussion

We have developed VCTF methods for diffusion tensors along fiber tracts in the Riemannian manifold of SPD matrices under the log-Euclidean metric. From the application end, VCTF is demonstrated in a clinical study of neurodevelopment for revealing the complex inhomogeneous spatiotemporal maturation patterns as the apparent changes in fiber bundle diffusion tensors. We have shown that this novel statistical tool leads to new findings in our clinical applications.

Another commonly used metric on the Riemannian manifold of SPD matrices is the Riemannian metric. In contrast, some operations, e.g. average or interpolation of a set of tensors under the Log-Euclidean and Riemannian metrics, are
<figure><img src="image_40.png" /><figcaption>**Fig. 5.** The ellipsoidal representations of (a) Raw diffusion tensor of all 96 infants and (b) smoothed diffusion tensors changing with the gestational age at one point on the right internal capsule tract with significant gestational age effect.</figcaption></figure>

theoretically and practically very similar [15]. Moreover, some statistical meth-
ods based on the two metrics have very similar results. However, the Riemannian
metric is affine invariant. Affine invariance is a desirable feature for imaging pro-
cessing, e.g. segmentation. In this scenario, as shown in [13], the method using
the Riemannian metric outperformed the method using the Log-Euclidean met-
ric. So it is interesting to develop VCTF under the Riemannian metric and then
compare the statistical powers of detecting group differences under these two
different metrics.

References

1. P. J. Basser, J. Mattiello, and D. LeBihan. Estimation of the effective self- diffusion tensor from the nmr spin echo. *Journal of Magnetic Resonance Ser.* B, 103:247–254, 1994.

2. P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi. In vivo fiber tractography using dt-mri data. *Magnetic Resonance in Medicine*, 44:625–632, 2000.

3. S. M. Smith, Jenkinson M., H. Johansen-Berg, D. Rueckert, T. E. Nichols, C. E. Mackay, K. E. Watkins, O. Ciccarelli, M.Z. Cader, P.M. Matthews, and T. E. Behrens. Tract-based spatial statistics: voxelwise analysis of multi-subject diffusion data. *NeuroImage*, 31:1487–1505, 2006.

4. P. A. Yushkevich, H. Zhang, T.J. Simon, and J. C. Gee. Structure-specific statis-
tical mapping of white matter tracts. *Neuroimage*, 41:448–461, 2008.

5. C. B. Goodlett, P. T. Fletcher, J. H. Gilmore, and G. Gerig. Group analysis of dti fiber tract statistics with application to neurodevelopment. *NeuroImage*, 45:S133–S142, 2009.

6. H. T. Zhu, M. Styner, N. S. Tang, Z. X. Liu, W. L. Lin, and J.H. Gilmore. Frats: Functional regression analysis of dti tract statistics. *IEEE Transactions on Medical Imaging*, 29:1039–1049, 2010.
7. J. Goldsmith, J. Feder, C. M. Crainiceanu, B. Caffo, and D. Reich. Penalized functional regression. Technical report, Working Paper 2010-204, Johns Hopkins University Dept. of Biostatistics., 2010.

8. H. T. Zhu, H. P. Zhang, J. G. Ibrahim, and B. G. Peterson. Statistical analysis of diffusion tensors in diffusion-weighted magnetic resonance image data (with discussion). *Journal of the American Statistical Association*, 102:1085–1102, 2007.

9. H. T. Zhu, D. Xu, R. Amir, X. Hao, H. Zhang, K. Alayar, B. Ravi, and B. Peterson. A statistical framework for the classification of tensor morphologies in diffusion tensor images. *Magnetic Resonance Imaging*, 24:569–582, 2006.

10. A. Schwartzman, W. Mascarenhas, and J. E. Taylor. Inference for eigenvalues and eigenvectors of gaussian symmetric matrices. *Ann. Statist.*, 36:1423–1431, 2008.

11. P. T. Kim and D. S. Richards. *Deconvolution density estimation on spaces of positive definite symmetric matrices*. IMS Lecture Notes Monograph Series. A Festschrift of Tom Hettmansperger, 2010.

12. H. T. Zhu, Y. S. Cheng, J. G. Ibrahim, Y. M. Li, C. Hall, and W. L. Lin. Intrinsic regression models for positive definitive matrices with applications in diffusion tensor images. *Journal of the American Statistical Association*, 104:1203–1212, 2009.

13. A. Barmpoutis, B. C. Vemuri, T. M. Shepherd, and J. R. Forder. Tensor splines for interpolation and approximation of dt-mri with applications to segmentation of isolated rat hippocampi. *IEEE Transactions on Medical Imaging*, 26:1537–1546, 2007.

14. B. C. Davis, E. Bullitt, P. T. Fletcher, and S. Joshi. Population shape regression from random design data. *IEEE 11th International Conference on Computer Vision*, pages 1–7, 2007.

15. V. Arsigny. Processing data in lie groups: An algebraic approach. application to non-linear registration and diffusion tensor mri. Ph.D. dissertation, 2006.

16. J. Fan and I. Gijbels. *Local Polynomial Modelling and Its Applications*. Chapman and Hall, London, 1996.

17. Jianqing Fan and Wenyang Zhang. Simultaneous confidence bands and hypothesis testing in varying-coefficient models. *Scand. J. Statist.*, 27(4):715–731, 2000.

18. H. T. Zhu, R. Z. Li, and L. L. Kong. Multivariate varying coefficient models for functional responses. Technical report, University of North Carolina at Chapel Hill, 2010.

19. S. Joshi, B. Davis, M. Jomier, and G. Gerig. Unbiased diffeomorphic atlas construction for computational anatomy. *Neuroimage*, 23:S151–S160, 2004.
# Matching and Classification of Images Using The Space of Image Graphs

Dohyung Seo, Jeffrey Ho, and Baba C. Vemuri

Dept of CISE, University of Florida, Gainsville FL, 32611, USA

**Abstract.** This paper proposes a geometric approach for comparing tensor-valued images (tensor fields) that is based on the idea of matching intrinsically low-dimensional shapes embedded in a higher-dimensional ambient space. More specifically, instead of regarding the tensor fields as tensor-valued functions defined on a given (image) domain, we consider their image graphs. These tensorial image graphs can naturally be regarded as submanifolds (shapes) in an ambient space that is the cartesian product of their domain and the space of tensors. With this viewpoint, comparisons between tensor fields can naturally be formulated as comparisons between their corresponding shapes, and an intrinsic comparison measure can be developed based on matching these low-dimensional shapes. The proposed approach offers great conceptual clarity and transparency, and thorny issues such as parametric invariance and symmetric registration can be handled effortlessly in this novel framework. Furthermore, we show that the resulting variational framework can be satisfactorily optimized using a gradient descent-based method, and the computed similarities can be used as the affinity measures in a supervised learning framework to yield competitive results on challenging classification problems. In particular, experimental results have shown that the proposed approach is capable of producing impressive results on several classification problems using the OASIS image database, which include classifying the MR brain images of Alzheimer's disease patients.

## 1 Introduction

A central application of computational anatomy is to quantify the difference between normal anatomy and pathology. In many application problems, this frequently requires comparisons between 3D volumetric images or shapes extracted from these images. In particular, the difference is often characterized via a similarity measure computed by first registering the two different images/shapes followed by an integration step that sums over the measured difference at corresponding pixels/voxels or surfels (surface elements). Due to its simplicity, this $L^2$-based approach is very popular in the literature, especially for comparing images and shapes. However, for many applications, tensor fields derived from the scalar-valued images usually contain more useful information that have frequently been under-utilized until now, and for these applications, comparisons between tensor fields are important and fundamental. The examples include the
metric tensor for shape comparison, deformation tensor for morphometry [1],
the diffusivity tensor [2] for diffusion MRI and many others. Therefore, there is
an increasing need for developing a sound and principled approach for comparing
tensor fields, and in this paper we will present one such approach.

The approach taken in this paper is entirely geometrical and it is based on
the fundamental idea of comparing the intrinsic shapes of the tensor fields. More
precisely, suppose $\mathbf{T}_i$, and $\mathbf{T}_j$ are $\mathbf{P}(n)$-valued tensor fields defined over an im-
age domain $\Omega$ considered as a subset in $\mathbb{R}^n$, where $\mathbf{P}(n)$ denotes the space of
symmetric positive-definite $n \times n$ matrices (Figure 1). The tensor fields can nat-
urally be regarded as submanifolds in the product space $\Omega \times \mathbf{P}(n)$. Once we have
equipped $\Omega \times \mathbf{P}(n)$ with a Riemannian metric, the geometries of $\mathbf{T}_i$, and $\mathbf{T}_j$ are
then naturally defined using the induced metrics. In other words, each tensor
field is considered as a kind of high-dimensional (tensorial) image graph [5]. For
a grey-scale image, its image graph is a submanifold in $\mathbb{R}^3$; however, for a tensor
field, its tensorial image graph is contained in a the (nonlinear) ambient space
that typically has dimension greater than three. With this particular viewpoint,
comparisons between tensor fields can be formulated intrinsically as compar-
isons between their corresponding tensorial image graphs, and the method of [4]
can be generalized directly to this higher-dimensional context by computing the
registration between two shapes embedded in the ambient space.

<figure><img src="image_41.png" /><figcaption>**Fig. 1.** **Left:** The tensorial image graph of a tensor field defined on $\mathbb{R}^3$. The ambient space is the cartesian product $\mathbb{R}^3 \times \mathbf{P}(3)$. **Right:** Comparing shapes in $\mathbb{R}^3 \times \mathbf{P}(3)$ can be realized via a registration map $X$ defined between the two shapes.</figcaption></figure>

The approach outlined above offers several important advantages over the rel-
atively straightforward $\mathbf{L}^2$-based approach. First, it is conceptually transparent
and by comparing the shapes (instead of tensor fields) directly on the tensorial
image graphs, the issue of parametrization invariance becomes irrelevant and the
complicated procedures of finding a parametrically invariant representation [2]
can be completely bypassed. Second, the issue of symmetric comparison and reg-
istration [3] is naturally incorporated through the use of intrinsic volume forms
defined on the shapes. Third, while the mathematics is more involved when com-
pared with the $\mathbf{L}^2$-based methods, it is still tractable and we have developed a
gradient descent-based method to efficiently optimize the objective function. Fi-
nally, the proposed approach is also flexible in that it permits different metrics to define different shapes for a given tensor field. For example, the ambient space $\Omega \times \mathbf{P}(n)$ can be equipped with several different metrics (for example, the usual Frobenius metric and affine-invariant metric on $\mathbf{P}(n)$), and our approach can easily accommodate these different choices.

For validation, we apply the proposed method to classify brain MRI from dif- ferent population groups contained in the OASIS database [10] which includes the challenging problem of classifying MR brain scans of patients suffering from the alzheimer's disease and that of normal healthy old adults. The MR brain images are first converted into their corresponding tensor representations, and pairwise similarities between the resulting tensor fields are computed using the proposed method. Nonlinear dimensionality reduction of the data is achieved using a diffusion map [11] with the pairwise similarities as the basic affinity measures, and the final classification step is carried out efficiently and accu- rately in a low-dimensional Euclidean space. In terms of better and stronger classification results, the experimental results reported in section 4 validate the two novel points advocated in this paper that the tensor field comparisons are fundamental and important in many applications, and a sound and principled approach to tensor field comparison can form the basis of algorithmic solutions to challenging classification problems.

## 2 Cost Function for Shapes Matching

Shapes can be represented as $R^n$-dimensional point clouds, triangular meshes, or parametric or implicit surfaces. In this paper, shapes are considered as para- metric surfaces. For instance, surfaces in 3-D are represented parameterically by : $f: \Omega \to R^3$, where $\Omega$ is a domain. If we represent images as image graphs in order to represent images as shapes, the desired parameterization is $f: (u, v) \to (x, y, I(x, y))$ for 2-D gray scale images or $f: (u, v) \to (x, y, P(3)(x, y))$ for 2-D 3-by-3 symmetric positive definite (SPD) tensor-valued images. More precisely, an image can be considered as a section of a fiber bundle [5]. Then for a pair of shapes $(S_1, S_2)$, a matching problem can be formulated as follows :

$$D((S_1, S_2), \gamma) = \min_{\gamma} \int Dist(S1, S2 \circ \gamma) \sqrt{\kappa} d\Omega, \quad (1)$$

where $S_1$ and $S_2$ share common domain $\Omega$. On this stage we assume that shapes are already globally registered. In Eq.(2), $Dist(S_1, S_2)$ is a point-wise distance defined in the ambient space between two shapes and $\sqrt{\kappa}d\Omega$ is volume form defined as $||df_p(u) \times df_p(v)||$ where $f: (u, v) \to p$ and $\sqrt{\kappa}$ is invariant under re-parametrization.

For shape comparison, $D((S_1, S_2), \gamma)$ is required to be invariant under re- parametrization, and symmetric between $S_1$ and $S_2$. Recently, [2] introduces the notion of $q$-map as a novel representation for shapes in $\mathbb{R}^3$ with the aim of achieving re-parametrization invariance of $D((S_1, S_2), \gamma)$. However if the ambient space is not Euclidean, e.g., $\mathbf{P}(3)$, it is not clear how to complete the centering
step in the construction of the $q$-map in [2]. A similar cost function as Eq.(2) was proposed earlier in [4] which is designed for matching 3-D face meshes. However, the proposed matching algorithm is not symmetric and nor is it designed for matching tensorial image graphs.

In this paper, we introduce a novel symmetric matching framework for low-dimensional shapes embedded in an ambient space $X$ that generally has dimension greater than three. Such shapes will be represented by their parameterizations with domains in $\mathbb{R}^2$ or $\mathbb{R}^3$. Let $S_1$ and $S_2$ be two such surfaces with parameterizations $f_1$ and $f_2$ defined on two domains $\Omega_1$ and $\Omega_2$ respectively. We will use the following cost function to define a similarity measure between these two shapes

$$E((S_1, S_2), \gamma) = \min_{\gamma} \int (Dist(f_1, f_2 \circ \gamma))^2 (\sqrt{\kappa_1} + \sqrt{\kappa_2(\Omega_2 \circ \gamma)} J_{\gamma}) d\Omega, \quad (2)$$

where $\sqrt{\kappa_1}d\Omega_1$ and $\sqrt{\kappa_2}d\Omega_2$ are the pull-backs of the volume form on $S_1, S_2$ under $f_1$ and $f_2$, respectively. $\Omega = \Omega_1$, and $J_{\gamma}$ is the determinant of the Jacobian of $\gamma$. The details of derivation of Eq.(2) is given in the appendix.

## 2.1 Tensorial Image Graphs

In this section, we apply our matching framework to tensor-valued images, specifically 3-by-3 3-D SPD tensor-valued images. For this purpose, each tensor-valued image is represented as a section of a fiber bundle [5] with the map $\mathbf{X}: \mathbb{R}^3 \to \mathbb{R}^3 \times \mathbf{P}(3)$ or $\mathbf{X}: (u, v, w) \to (x(u, v, w), y(u, v, w), z(u, v, w), I(x, y, z))$. In this map, $x = u, y = v, z = w$ and $I(x, y, z) \in \mathbf{P}(3)$ at each voxel. This parametrization was introduced by Gur et. al [5] to achieve anisotropic smoothing of 2-D DTI. In this application, the volume form is $\sqrt{\kappa}dudvdw$ where $\kappa$ is the determinant of

$$\mathbf{K} = \begin{pmatrix} < \mathbf{X}_u, \mathbf{X}_u > & < \mathbf{X}_u, \mathbf{X}_v > & < \mathbf{X}_u, \mathbf{X}_w > \\ < \mathbf{X}_u, \mathbf{X}_v > & < \mathbf{X}_v, \mathbf{X}_v > & < \mathbf{X}_v, \mathbf{X}_w > \\ < \mathbf{X}_u, \mathbf{X}_w > & < \mathbf{X}_v, \mathbf{X}_w > & < \mathbf{X}_w, \mathbf{X}_w > \end{pmatrix}, \quad (3)$$

or

$$\mathbf{K} = \begin{pmatrix} \lambda + \mathrm{Tr}((I^{-1}I_x)^2) & \mathrm{Tr}((I^{-1}I_x)(I^{-1}I_y)) & \mathrm{Tr}((I^{-1}I_x)(I^{-1}I_z)) \\ \mathrm{Tr}((I^{-1}I_x)(I^{-1}I_y)) & \lambda + \mathrm{Tr}((I^{-1}I_y)^2) & \mathrm{Tr}((I^{-1}I_y)(I^{-1}I_z)) \\ \mathrm{Tr}((I^{-1}I_x)(I^{-1}I_z)) & \mathrm{Tr}((I^{-1}I_y)(I^{-1}I_z)) & \lambda + \mathrm{Tr}((I^{-1}I_y)^2) \end{pmatrix}, \quad (4)$$

and the metric in the ambient space is given by:

$$\mathbf{G} = \begin{pmatrix} \lambda & 0 & 0 & 0 \\ 0 & \lambda & 0 & 0 \\ 0 & 0 & \lambda & 0 \\ 0 & 0 & 0 & \mathbf{P} \end{pmatrix}, \quad (5)$$

where $\mathbf{P}$ is the metric in $\mathbf{P}(3)$ space. We will call the shapes (submanifolds) represented by the parametrization $\mathbf{X}$ as the tensorial image graph of the tensor-valued images $I$. Representing the tensorial image graphs $S_1$ and $S_2$ by $(x_1, y_1, z_1, I_1)$
and $(x_2, y_2, z_2, I_2)$ respectively, we define $Dist((S_1, S_2), \gamma)$ in Eq.(2) by:

$$Dist((S_1, S_2), \gamma) = \lambda((x_1 - x_2)^2 + (y_1 - y_2)^2 + (z_1 - z_2)^2) + dist(I_1, I_2), \quad (6)$$

and

$$dist(I_1, I_2) = Trace(\log(I_1^{-1/2} R^\dagger (I_2 \circ \gamma) R I_1^{-1/2})^2). \quad (7)$$

Eq.(7) is voxel-wise Riemannian distance between tensors [6], with R the matrix used to re-orient the tensors (as they undergo a registration). When we choose $\Omega$ as the common domain, $(x_1, y_1, z_1) = (u, v, w)$ and $(x_2, y_2, z_2) = (u + U(u, v, w), v + V(u, v, w), w + W(u, v, w))$, where $U, V, W$ are the components of the deformation field. In Eq.(7), R is reoriented with respect to the deformation vector field $\gamma$ (or $U, V, W$) at each voxel. Reorientation of tensors must be carried out during optimization to and the reorientation transformation is derived from the deformation field [7, 8]. Also, recall that $Trace((\log(A)^2) = \sum_{i=1}^3 (\log(\Lambda_i))^2$, where $\Lambda_i$'s are the eigenvalues of A.

To smooth the vector field, $(U, V, W)$, we add the regularization term to Eq.(2):

$$E_{reg} = \int |\nabla U|^2 + |\nabla V|^2 + |\nabla W|^2 d\Omega, \quad (8)$$

and the final cost function, $E_{tot}((S_1, S_2), \gamma)$ is given as:

$$E_{tot}((S_1, S_2), \gamma) = (1 - \alpha)E((S_1, S_2), \gamma) + \alpha E_{reg}, \quad (9)$$

where $\alpha$ is a small positive scalar.

To efficiently solve the resulting optimization problem, we first discretize the cost function to get,

$$E_{tot}(\mathbf{U}, \mathbf{V}, \mathbf{W}) = \sum_{i \in \Omega} [(1 - \alpha)Dist(U_i, V_i, W_j)(\sqrt{\kappa_{1i}} + \sqrt{\kappa_{2i}}J_{\gamma i}) + \alpha \sum_{k=\{x,y,z\}} (U_{ki}^2 + V_{ki}^2 + W_{ki}^2)] \quad (10)$$

where $U_i, V_i, W_j, U_{ki}, V_{ki}, W_{ki}$ are the values of $U, V, W$ at the given discrete points. We optimize the above cost function with respect to $U, V, W$ using non-linear conjugate gradient (NCG) method following [9]. To simplify the steps, we evaluate the gradient vector and the Hessian matrix of Eq.(10) with fixed $(\sqrt{\kappa_{1i}} + \sqrt{\kappa_{2i}}J_{\gamma i})$ in each NCG iteration, and in turn, $(\sqrt{\kappa_{1i}} + \sqrt{\kappa_{2i}}J_{\gamma i})$ is computed using the most recently updated deformation vector field, $(U_i, V_i, W_i)$.

# 3 Application : Classification of MRI data

To apply our matching framework to classification of tensor-valued image data set, we first convert the input images into their associated tensor-valued images (tensor fields) and compute pairwise matching. In the second step, we use the L2 norms of the deformation vector fields, $d(S_i, S_j)$ between the two tensorial
image graphs $S_i$ and $S_j$ to build a data graph with Gaussian weights $e^{-d(S_i, S_j)/\epsilon}$ according to [11]. The corresponding Markov matrix is used to compute the diffusion maps [11], which provides a dimensionality reduction of the input image data. The nearest neighbor classifier is then used for classification using the diffusion distances in the low-dimensional feature space.

## 3.1 Data Preparation

To create 3-D tensor images from a 3D gray-scale images, we used the OASIS MRI database [10]. The images in the database are the MR human brain scans of subjects aged between 18-96 of 416. Each image has a resolution of 208×176×176 voxels. Our goal is to classify the MR image data into different age groups, and we have chosen the ventricle as the ROI (region of interest) as it captures the part of brain showing the most significant difference across ages. Fig.2 shows ROIs from sagittal and longitudinal MR slices acquired from brains of 18, 43 and 81 year old subjects respectively.

Our tensor-valued images are the first fundamental forms (metric tensors) of image graphs of the 3-D intensity values with map $f: R^3 \to R^3 \times R$,

$$ \left( \begin{array}{c} < \mathbf{f}_u, \mathbf{f}_u > < \mathbf{f}_u, \mathbf{f}_v > < \mathbf{f}_u, \mathbf{f}_w > \\ < \mathbf{f}_u, \mathbf{f}_v > < \mathbf{f}_v, \mathbf{f}_v > < \mathbf{f}_v, \mathbf{f}_w > \\ < \mathbf{f}_u, \mathbf{f}_w > < \mathbf{f}_v, \mathbf{f}_w > < \mathbf{f}_w, \mathbf{f}_w > \end{array} \right), \qquad (11) $$

and we consider these tensor fields as tensorial image graphs. Fig.3 shows the tensor-valued images from subjects in Fig.2 according to Eq.11.

## 3.2 Diffusion Map and Diffusion Distance for Classification

In this paper, we represent the SPD tensor-valued images as sections of a fiber bundle, therefore we need to find a meaningful geometric description of the space of sections for classification purposes and diffusion maps can generate efficient representation of desired geometric structures based on the diffusion processes.[11] Once we build a graph with Gaussian weights $e^{-d(S_i, S_j)/\epsilon}$ and construct the corresponding Markov matrix, $M$, then the family of diffusion maps {$\Psi_t$}$_{t \in N}$ and diffusion distance $D_t(S_i, S_j)$ are defined as follows :

$$ \Psi_t(\mathbf{S}_i) = \begin{pmatrix} \Lambda_1^t \Psi_1(S_i) \\ \Lambda_2^t \Psi_2(S_i) \\ \vdots \\ \Lambda_s^t \Psi_s(S_i) \end{pmatrix}, \qquad (12) $$

and

$$ D_t(S_i, S_j) = \left( \sum_{l=1}^{s} \Lambda_l^{2t} (\Psi_l(S_i) - \Psi_l(S_j))^2 \right)^{1/2}, \qquad (13) $$
where {$A_l$}$_{l \ge 0}$ and {$\Psi_l$}$_{l \ge 0}$ are eigenvalues and eigenvectors of $\mathbf{M}$ respectively such that $1 = A_0 > |A_1| \ge |A_2| \ge \dots$ and $\mathbf{M}\Psi_l = A_l\Psi_l$ (Ref.), and

$$s = \max\{l \in \mathcal{N} \text{ such that } |A_l|^t > \delta |A_1|^t\}. \qquad (14)$$

The diffusion map $\Psi_t$ embeds all data in the set, {$S_i$} into the Euclidean space $\mathbf{R}^s$ and $D_t$ reflects the connectivity in the graph of the data in the set, {$S_i$} defining the Euclidean distance in $\mathbf{R}^s$: points in the set {$S_i$} are closer if they are highly connected in the graph.

We set $d(S_i, S_j)$ as L2 norm of deformation vector fields, or $\int (|U|^2 + |V|^2 + |W|^2)d\Omega$, after matching and use the diffusion distance as feature of nearest neighbor classifier in the low dimensional space.

Fig.4 shows $\Psi_1$ vs. $\Psi_2$ plots for classifications between groups. The details of classification results are reported in next section.

# 4 Experimental Results

In this section, we report the experimental results on classifications of MR brain images from the OASIS data set. We divided the subjects into three groups : the young group with age below 40, the old group with age 60 or above and the middle-aged group between 41-60. In these experiments, we use a four-fold cross-validation and leave-one-out validation to determine the classification scores. In the four-fold cross-validation, the subgroups are randomly selected 50 times and the maximum, minimum and average classification scores together with the variances are reported in the first three columns in Table 1. We also test our method to classify Alzheimer's disease (AD). We take 70 subjects from old age group and in the subgroup, 35 of them are diagnosed as AD and rest of them are control, and we use four-fold and leave-one-out validation tests within the subgroup[13]. The classifier used in the reduced dimension is the nearest-neighbor classifier, and the diffusion distance is used as distance measure. The criterion for determining the dimension of the diffusion map in these experiments is given by $\delta = 0.07$ in Eq.(14) with $t = 1$. The metric used for the ambient space $\mathbb{R}^3 \times \mathbf{P}(3)$ is the product metric of the Euclidean metric on $\mathbb{R}^3$ and the affine-invariant metric on $\mathbf{P}(3)$. Furthermore, tensor reorientation [7],[8] is applied to reorient the tensors after transformation. We note that for three different sets of classification, the classification rates are uniformly high.

In the second set of experiments, we test the effects of using tensor field (compared with only scalar field) and the choice of different metric on $\mathbf{P}(3)$. First, we change the metric on $\mathbf{P}(3)$ from the affine-invariant metric to the Frobenius norm (i.e., $\mathbf{L}^2$-norm), and four-fold cross-validation results between young and old age group are reported in the first two columns in Table 2. In these two columns, we test the effect of tensor reorientation on classifications, and the results show that for Frobenius norm, the effect is generally small. In the third column, the images are represented by their image graphs and the shape comparisons are carried out in the ambient space $\mathbb{R}^3 \times \mathbb{R}$ instead of $\mathbb{R}^3 \times \mathbf{P}(3)$. In this experiment, distances between image graphs are the Riemannian distances
and $\lambda = 0.000001$ in Eq.(6) and $\alpha = 0.02$ in Eq.(9). We remark that the result clearly demonstrates that the classification result using only scalar-valued images (image graphs in $\mathbb{R}^3$) is inferior to the one using tensor fields and the affine-invariant metric on $\mathbf{P}(3)$ provide superior classification result compared with the Frobenius metric.

In Table 3, we show the comparisons between our method and several previously published classification results on the OASIS database. [12] uses the deformation tensor field (computed from registering the image to an atlas) as the main feature for each image. Submanifold of each age group is constructed from the training samples and the geodesic distances between subjects and the submanifolds are used as the main discriminative feature for classification. In [13], alternatively, histograms of deformation vector fields have been used as features, and the CAVIAR method proposed in [13] takes a adaboost-like approach to integrate the results from a collection of weak classifiers into a strong classification result. We remark that our method compared favorably with these methods in terms of classification rates, and in particular, for the more challenging problem of classifying brain images of Alzheimer's disease patients, our method demonstrates a small but real improvement over these two methods.

<table><caption>**Table 1.** Scores of leave-one-out and four-fold cross-validation test of four subgroups randomly selected 50 times. The metric used for the ambient space $\mathbb{R}^3 \times \mathbf{P}(3)$ is the product metric of the Euclidean metric on $\mathbb{R}^3$ and the affine-invariant metric on $\mathbf{P}(3)$.</caption><thead><tr><th></th><th>Old vs. Young</th><th>Old vs. Middle</th><th>Middle vs. Young</th><th>AD vs. Control</th></tr></thead><tbody><tr><th>Maximum</th><td>100%</td><td>100%</td><td>100%</td><td>100%</td></tr><tr><th>Minimum</th><td>97.72%</td><td>90.77%</td><td>87.27%</td><td>75.0%</td></tr><tr><th>Average</th><td>99.25%</td><td>97.6%</td><td>94.36%</td><td>94.87%</td></tr><tr><th>Standard deviation</th><td>0.8384%</td><td>2.46%</td><td>3.09%</td><td>5.55%</td></tr><tr><th>Leave-one-out</th><td>99.15%</td><td>98.46%</td><td>96.36%</td><td>95.32%</td></tr></tbody></table>

<table><caption>**Table 2.** Scores of leave-one-out and four-fold cross-validation test of four subgroups randomly selected 50 times between Young and Old samples using Frobenius metric on $\mathbf{P}(3)$ with and without tensor re-orientation. Last column gives the classification result of using only gray-scale images.</caption><thead><tr><th>Methods</th><th>Frobenius with reo.</th><th>Frobenius without reo.</th><th>Image Graph</th></tr></thead><tbody><tr><th>Maximum</th><td>100%</td><td>100%</td><td>98.86%</td></tr><tr><th>Minimum</th><td>88%</td><td>88.63%</td><td>88.63%</td></tr><tr><th>Average</th><td>94.96%</td><td>96.52%</td><td>94.29%</td></tr><tr><th>Standard deviation</th><td>2.57%</td><td>2.15%</td><td>2.29%</td></tr><tr><th>Leave-one-out validation</th><td>97.03%</td><td>98.01%</td><td>94.89%</td></tr></tbody></table>
<table><caption>**Table 3.** Comparison of classification scores with 4-fold validation between classification methods</caption><thead><tr><th></th><th>Old vs. Young</th><th>Old vs. Middle</th><th>Middle vs. Young</th><th>AD vs. Control</th></tr></thead><tbody><tr><td>Image Graphs</td><td>99.25%</td><td>97.6%</td><td>94.36%</td><td>94.87%</td></tr><tr><td>CAVIAR [13]</td><td>99.14%</td><td>98.36%</td><td>97.76%</td><td>88.0%</td></tr><tr><td>Adaboost [13]</td><td>98.75%</td><td>96.80%</td><td>96.0%</td><td>90.25%</td></tr><tr><td>Submanifold projection [12]</td><td>96.43%</td><td>90.23%</td><td>84.32%</td><td>88.57%</td></tr><tr><td>Nearest Neighbor in PCA [12]</td><td>92.43%</td><td>87.74%</td><td>78.42%</td><td>84.29%</td></tr></tbody></table>

# 5 Conclusion

We have proposed a novel geometric approach for comparing tensor-valued images (tensor fields) that is based on the simple idea of matching the low-dimensional tensorial image graphs formed by the tensor fields. Our framework provides a registration method that is both symmetric and invariant under different parametrization, and the resulting cost function can be satisfactorily optimized using a gradient descent-based method. We have reported four different classification experiments using the OASIS image database, and our method has produced results that are in par or exceeding the current state-of-the-art results. In particular, our experiments have shown that tensor fields do indeed contain subtle information that can be useful for challenging classification problems, and the experimental results have demonstrated that the proposed method, although more elaborate and involved compared with L²-based method, is able to access and utilize this information to obtain good classification results.

# References

1. N. Lepore, C. Brun, Y. Chou, M. Chiang, R. A. Dutton, K. M. Hayashi, E. Luders, O. L. Lopez, H. J. Aizenstein, A. W. Toga, J. T. Becker, P. M. Thompson, "Generalized Tensor-Based Morphometry of HIV/AIDS Using Multicariate Statistics on Deformation Tensors", IEEE Transactions on Medical Imaging, 2008.

2. S. Kurtek, E. Klassen, Z. Ding and A. Srivastava, "A novel riemannian framework for shape analysis of 3D objects", Proc. CVPR pp. 1625-1632, 2010.

3. H. D. Tagare, D. Groisser and O. Skrinjar, "Symmetric Non-rigid Registration: A Geometric Theory of Some Numerical Techniques", J Math Imaging Vis pp. 61-88, 2009.

4. N. Litke, M. Droske, M. Rumpf and P. Schröder, "An Image Processing Approach to Surface Matching", in "Symposium on Geometry Processing", pp. 207-216, 2005.

5. Y. Gur and N. Sochen, "Coordinate-based diffusion over the space of symmetric positive-definite matrices", in "Visualization and Processing of Tensor Fields - Advances and Perspectives", pp 325-340, Springer-Verlag, Berlin Heidelberg, 2009.

6. M. Moakher and P. G. Batchelor, "Symmetric Positive-Definite Matrices: From Geometry to Applications and Visualization", in "Visualization and Processing of Tensor Fields", pp285-298, Springer-Verlag, New York, 2006.
7. G. Cheng, B. C. Vemuri, M. Parekh, P. R. Carney, and T. H. Mareci, “Non-rigid Registration of HARDI Data Represented by a Field of Gaussian Mixtures”, In LNCS (Springer) Proceedings of MICCAI09: Int. Conf. on Medical Image Computing and Computer Assisted Intervention, 2009.

8. D. C. Alexander, C. Pierpaoli, P. J. Basser, and J. C. Gee, “Spatial Transformations of diffusion tensor magnetic resonance images”, IEEE Trans. Med. Imag. 20 11 pp. 1131-1139, 2001.

9. S.H. Lai and B. C. Vemuri, “Reliable and Efficient Computation of Optical Flow”, IJCV 29(2), pp. 87-105, 1998.

10. D. S. Marcus, T. H. Wang, J. parker, J. G. Csernansky, J. C. Moris, R. L. Buckner “Open Access Series of Imaging Studies (OASIS): Cross-Sectional MRI Data in Young, Middle Aged, Nondemented, and Demented Older Adults”. Journal of Cognitive Neuroscience (2007)

11. R. R. Coifman and S. Lafon, “Diffusion Maps”, Appl. Comput. Harmon. Anal. 21 pp. 5-30, 2006

12. Y. Xie, B. C. Vemuri, and J. Ho, “Statistical Analysis of Tensor Fields ”, In Int. Conf. on Medical Image Computing and Computer Assisted Intervention (MICCAI), 2010.

13. T. Chen, A. Rangarajan, and B. C. Vemuri, “CAVIAR: Classification via Aggregated Regression and Its Application in Classifying the OASIS Brain Database”, In IEEE International Symposium on Biomedical Imaging, pp. 1337-1340, 2010.

Appendix. Cost Function for Symmetric Matching

If we have two parameterized shapes, $S_1$ and $S_2$ with domain $\Omega_1$ and $\Omega_2$, then
the cost function for symmetric matching problem is formulated as following :

$$
E((S_1, S_2), \phi, \psi) = \min_{\phi} \int_{\Omega_1} \mathrm{Dist}(S_1, S_2 \circ \phi) \sqrt{\kappa_1} d\Omega_1 + \min_{\psi} \int_{\Omega_2} \mathrm{Dist}(S_2, S_1 \circ \psi) \sqrt{\kappa_2} d\Omega_2, \quad (A-1)
$$

where $\phi: \Omega_1 \to \Omega_2$ and $\psi: \Omega_2 \to \Omega_1$, and $\sqrt{\kappa} d\Omega$ is volume form. If The second
term in the left hand side of Eq(A-1) can be rewritten as following :

$$
\int_{\Omega_1} Dist(S_2 \circ \psi^{-1}, S_1) \sqrt{\kappa_2 (\Omega_2 \circ \psi^{-1})} J_{\psi^{-1}} d\Omega_1, \quad (A-2)
$$

and if we require that Eq(A-1) is symmetric matching, $\psi^{-1}$ should be $\phi$ and
$J_{\psi^{-1}} = J_{\phi}$ which is determinant of Jacobian such as

$$
Det \left( \begin{matrix} \frac{\partial u_2}{\partial u_1} & \frac{\partial u_2}{\partial v_1} & \frac{\partial u_2}{\partial w_1} \\ \frac{\partial v_2}{\partial u_1} & \frac{\partial v_2}{\partial v_1} & \frac{\partial v_2}{\partial w_1} \\ \frac{\partial w_2}{\partial u_1} & \frac{\partial w_2}{\partial v_1} & \frac{\partial w_2}{\partial w_1} \end{matrix} \right) \quad (A-3)
$$

And $\kappa_2$ is determinant of $\mathbf{K}_2$ given as following :

$$
\mathbf{K}_2 = \begin{pmatrix} <S_{2u_2}, S_{2u_2}> <S_{2u_2}, S_{2v_2}> <S_{2u_2}, S_{2w_2}> \\ <S_{2u_2}, S_{2v_2}> <S_{2v_2}, S_{2v_2}> <S_{2v_2}, S_{2w_2}> \\ <S_{2u_2}, S_{2w_2}> <S_{2v_2}, S_{2w_2}> <S_{2w_2}, S_{2w_2}> \end{pmatrix} \quad (A-4)
$$

$$
= (J_{\phi}^{-1})^2 \mathbf{K}'_2, \qquad (A-5)
$$
and

$$ K_2' = \begin{pmatrix} <S_{2u_1}, S_{2u_1}> <S_{2u_1}, S_{2v_1}> <S_{2u_1}, S_{2w_1}> \\ <S_{2u_1}, S_{2v_1}> <S_{2v_1}, S_{2v_1}> <S_{2v_1}, S_{2w_1}> \\ <S_{2u_1}, S_{2w_1}> <S_{2v_1}, S_{2w_1}> <S_{2w_1}, S_{2w_1}> \end{pmatrix} \quad (A-6) $$

Finally the cost function is given as following :

$$ E((S_1, S_2), \phi, \psi) = \min_{\phi} \int_{\Omega_1} Dist(S_1, S_2 \circ \phi)(\sqrt{\kappa_1} + \sqrt{\kappa_2'}) d\Omega_1, \quad (A-7) $$

where $\kappa_2'$ is the determinant of $K_2'$.

<figure><img src="image_42.png" /><figcaption>Fig. 2. Slices of 3-D MRI images. (a)-(c) : Cross-sectional images of ventricles of 18, 43, and 81 years old respectively. (d)-(f): Longitudinal images of ventricles in the same order.</figcaption></figure>
<figure><img src="image_43.png" /><figcaption>**Fig. 3.** Slices of 3-D tensor images of ventricles created by Eq.(11). (a)-(c) : 18, 43, and 81 years old respectively.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>**Fig. 4.** 2-D plots of diffusion maps: (a) young vs. old, (b) young vs. middle, and (c) middle vs. old. In each plot, x-axis and y-axis are $\Psi_1$ and $\Psi_2$, respectively</figcaption></figure>
# Detecting Long Distance Conditional Correlations Between Anatomical Regions Using Gaussian Graphical Models

Stéphanie Allassonnière, Pierre Jolivet and Christophe Giraud *

CMAP Ecole Polytechnique, Route de Saclay, 91128 Palaiseau cedex, France
http://www.cmap.polytechnique.fr

**Abstract.** The conditional correlation patterns of an anatomical shape may provide some important information on the structure of this shape. We propose to investigate these patterns by Gaussian Graphical Modelling. We design a model which takes into account both local and long-distance dependencies. We provide an algorithm which estimates sparse long-distance conditional correlations, highlighting the most significant ones. The selection procedure is based on a criterion which quantifies the quality of the conditional correlation graph in terms of prediction. The preliminary results on AD versus control population show noticeable differences.

**Keywords:** Gaussian Graphical Model, long-distance conditional correlation, sparse dependencies, deformation of Alzheimer's disease

## 1 Introduction

It is quite natural to think that any region of the brain depends on all the other regions at least via some detours. This is equivalent to say that these regions are all correlated to each other. Several studies [10, 12, 5] have tried to design correlation patterns by computing the correlation matrix given by correlations between any two regions (using a PCA decomposition for example). Then, they highlight the most significant or stable ones by thresholding. However, correlations describe the global statistical dependencies between variables, corresponding to both direct and indirect interactions. The *direct* relations between two of these regions are less numerous and harder to capture but they carry some interesting information as well. Focusing on these direct dependencies, called *conditional correlations*, enables to see which areas directly affect the behaviour of a given region and avoids the burden given by the indirect interactions.

We want to study the geometric conditional dependencies of a response given by a real valued signal carried by a discrete grid. A natural approach to estimate these conditional correlations is given by Gaussian Graphical Modelling (GGM)

<footnote>* The authors thank Johns Hopkins University for the data and Dominique for the cartoon. Support for the JHU group was provided in part by P41-RR015241 (JHU) and U24 RR021382 (MGH) to BIRN, National Center for Research Resources (NIH), P01-AG0568 and P01-AG03991 (WU), and the NSF via DMS-0456253.</footnote>
[11, 15, 7]. The response is modelled as a Gaussian random vector. The conditional correlations of the response variables are depicted by a graph. Each node represents a variable and an edge in the graph is set between two variables if they are conditionally dependent given all the remaining ones. This corresponds to a non-zero entry of the inverse covariance matrix called the precision matrix.

In the context of medical data, we face two problems with Gaussian Graphical Modelling. The first one is the high dimension of the data compared to the low sample size. The second problem is that in this framework, the neighbour points are very likely to be conditionally correlated since the medical phenomenon is usually continuous and these local dependencies can be predominant in the estimated graph possibly hiding any other links. We propose here an algorithm and a selection criterion which address both questions.

The high-dimension-low-sample-size paradigm is a common issue in medical imaging. This usually leads to non accurate and non stable statistical analysis. In our graphical framework, this is equivalent to a high dimensional graph and small number of measures of the random signal. To face this situation, sparse representation is well known to be a very powerful tool in Computer Vision and Pattern recognition [26]. Although we know that the underlying real graph is not sparse, it is interesting to perform a *sparse estimation of its structure* [13, 14]. Such an estimation provides both a stable estimate of the conditional dependencies and a selection of the edges with the highest conditional correlations. The procedure proposed by [15, 13] is based on the LASSO algorithm [21]. In addition, the elastic net algorithm [27] slightly relaxes this sparse constrain providing more accurate and robust-to-noise sparsity patterns [9]. This sparsity constraint is crucial as it prunes the graph keeping a small number of edges.

As noticed above, local dependencies of medical images are predominant while estimating the conditional correlations. We want to go beyond that and focus on long-distance dependencies. To show only long distance dependencies, the authors in [10] post-process their graph merging some nodes and edges when the nodes are closer than a fixed threshold. However, we want to introduce this knowledge as a prior to the estimation. We have brought in the algorithm a neighbourhood prior: knowing that the neighbours are linked, we only look for the other conditional dependencies which are modelled as sparse. This enables to highlight the long distance conditional correlations and to gain in term of statistical accuracy since the estimation of long-distance relations is less affected by the noise and the neighbours. The prior knowledge is modelled as a neighbouring graph provided by the user. Because the neighbouring graph is given by the user, it can take into account geographical closeness or some anatomical proximity given for example by fibres.

The estimation procedure depends on some free parameters. In order to find suitable values, one can either use a cross-validation [13] or a theoretical analysis based on a restrictive condition on the data [14]. Unfortunately these choices are not satisfactory in practice and the final choice is usually done by hand looking at the estimations. To avoid this subjective choice, Giraud [7] proposes, for the usual sparse GGM framework, a criterion which optimises the choice of the graph
depending on its capacity to predict one random variable given the others. We
generalise this criterion to the case of neighbourhood prior.

In Computational Anatomy, learning what characterises the population in
terms of dependencies as in [19] for Alzheimer's disease (AD) is an important as-
pect of the analysis of the pathology. Different but complementary methods have
been used to analyse populations for instance Deformation Based-Morphometry
(DBM) [3], Tensor Based-Morphometry [25] or Feature Based-Morphometry [22].
The general idea is to compute some statistics from a sample of template-to-
subject deformations. Statistical analysis can also be done directly on the de-
formations [3, 6, 2] or by modelling the generation of images [18, 1] showing the
admissible deformations in a population. Although the global deformations give
an interesting interpretation, we may want to know what parts of the shape are
deformed jointly. Our proposed algorithm based on LASSO or elastic net has
been tested on a training set of high dimension data representing some deforma-
tions of a template hippocampus towards 101 patients' hippocampi. The results
show some correlated regions when dealing with deformations which suggest that
this structure has regions which move together. Moreover, the population is di-
vided into two groups, controls and patients suffering from AD. The comparison
of the estimated graphs of the two groups show some noticeable differences which
may be a new characterisation of AD.

## 2 Gaussian Graphical Models

Let us consider *p* points on a given shape that will compose
the nodes of the graph. On these points, we observe *n* random
responses such as a discretisation (with *p* points) of a quan-
tification of template-to-subject deformations. The *p* nodes of
the graph are thus identified to *p* random variables denoted
(*X*₁, ..., *X*ₚ), this vector is assumed to be distributed as a mul-
tivariate Gaussian *N*ₚ(0, *Σ*) (data are re-centred if necessary).
The graph *G*<sub>*Σ*</sub> of conditional dependencies is defined as fol-
lows: there exists an edge between nodes *a* and *b* if and only if
the variables *X*<sub>*a*</sub> and *X*<sub>*b*</sub> are dependent given all the remain-
ing variables. This will be denoted *a* <sup>*G*<sub>*Σ*</sub></sup> *b*. To illustrate the
notion of conditional dependency, let us give a toy example
illustrated in Fig. 1. The traffic jam intensity and the number
of snowmen in town are correlated due to snowstorm. But
conditionally on the height of snow, the number of snowmen
is independent of the intensity of traffic jams. This corresponds to the graph in
Fig 1 with only two edges. In particular, there is no edge between the traffic jam
and snowmen variables.

**Fig. 1.** Conditionally on the height of snow, the number of snowmen is independent of the intensity of traffic jams. This is represented by a two edges graph.

<figure><img src="image_45.png" /><figcaption>Fig. 1. Conditionally on the height of snow, the number of snowmen is independent of the intensity of traffic jams. This is represented by a two edges graph.</figcaption></figure>

## 2.1 Estimation processes

The estimation is done by a regression as presented in [15]. Let $X_a$ be the $a^{th}$ node of the graph. The goal is to estimate the matrix $\theta$ such that: $X_a =$
$\sum_{b \neq a} \theta_{a,b} X_b + \varepsilon_a$ where $\varepsilon_a$ is assumed to follow a normal distribution with zero-mean and variance $1/(\Sigma^{-1})_{a,a}$. An important property of this modelling is that $\theta_{a,b} = -\frac{K_{a,b}}{K_{a,a}}$ where $K = \Sigma^{-1}$ is the precision matrix. Therefore $\theta_{a,b} = 0$ is equivalent to $a \not\sim b$. Let $\theta_{.,b}$ be the $p$-dimensional vector of $(\theta_{a,b})_{1 \leq a \leq p}$ and $X_{-a}$ the $(p-1)$-dimensional random vector $(X_1, \dots, X_{a-1}, X_{a+1}, \dots, X_p)$. The estimation of $\theta$ is done by minimising an energy (which is actually the negative penalised log likelihood of the model) given by two terms. The first one is the $l^2$ distance between the observations and their estimation (coming from the Gaussian distribution of the noise). The second term is a prior on the parameters $\theta$ to estimate. In the case we consider, the prior is a Laplacian prior with parameter $\lambda$. This corresponds to an $l^1$ penalty on the energy side. The energy to minimise is therefore the following:

$$ \forall b \in \{1, \dots, p\}, \quad \hat{\theta}_{.,b} = \operatorname*{argmin}_{\theta_{.,b}} \|X_b - X_{-b}\theta_{.,b}\|_2^2 + \lambda \|\theta_{.,b}\|_1. \quad (1) $$

This regression technique is known as the LASSO algorithm [21]. Introducing a penalty on the estimated parameter $\theta$ enables to estimate it. Indeed, $\theta$ may be of large dimension which prevents from using a small sample size $n$. Reducing the estimation to more likely or expected matrices, enables to handle this typical case in medical imaging.

The minimum argument of this energy is sparse thanks to the $l^1$ penalty. Besides, it provides a stable estimation in terms of prediction of a variable given the other ones. Including this constraint on the candidate matrices enables to really focus on the most important conditional correlations. Therefore predicting one variable from the others is not much dependent on the training set.

However, one may want to relax slightly this constraint by balancing the effect of this $l^1$ penalty by an $l^2$ one which has the opposite behaviour. Instead of looking for sparse matrix with big coefficients, it tends to spread to all elements in the matrix. This is known as the elastic net algorithm [27] which provides more accurate sparsity patterns [9]. The energy is:

$$ \forall b \in \{1, \dots, p\}, \quad \hat{\theta}_{.,b} = \operatorname*{argmin}_{\theta_{.,b}} \|X_b - X_{-b}\theta_{.,b}\|_2^2 + \lambda \|\theta_{.,b}\|_1 + \gamma \|\theta_{.,b}\|_2^2. \quad (2) $$

These two algorithms take each node independently one after the other. This has a drawback considering that we would like to have a symmetric relation between nodes. Minimising either energy above does not guarantee that $\hat{\theta}_{a,b} \neq 0$ at the same time as $\hat{\theta}_{b,a} \neq 0$. The construction of the graph is therefore not straightforward. Two natural choices appear to enforce the symmetry [15]. On the one hand, we can say that only one of $\hat{\theta}_{a,b}$ or $\hat{\theta}_{b,a}$ needs to be non-zero to set an edge in the graph between $a$ and $b$. This yields $a \not\sim b$ if and only if $\hat{\theta}_{a,b} \neq 0$ or $\hat{\theta}_{b,a} \neq 0$. On the other hand, we can prefer a harder decision process by saying that only one of the $\hat{\theta}_{a,b}$ or $\hat{\theta}_{b,a}$ needs to be zero to remove the edge between $a$ and $b$ in the graph. This leads to $a \not\sim b$ if and only if $\hat{\theta}_{a,b} \neq 0$ and $\hat{\theta}_{b,a} \neq 0$. Both will be tested in the sequel.
## 2.2 Penalised criterion to choose among a family of estimated graphs

In the two previous energies (1) and (2), there are two parameters which have to be chosen. This may be very difficult and data depending. A classical choice is to use cross-validation. Giraud in [7] suggests instead to set a criterion so that after computing a family of possible graphs with different parameters, it enables to select the best one with respect to this criterion which aims at answering the question of how good the network is to predict one variable from the others.

To define the criterion crit(.) we need to introduce a few notations. We associate to any graph G in $\hat{\mathcal{G}}$, the $p \times p$ matrix

$$ \hat{\theta}^G = \arg \min_{\theta \in \Theta_G} \sum_{1 \le i \le n, 1 \le a \le p} [X - X\theta]_{i,a}^2, $$

where $\Theta_G$ is the set of $p \times p$ matrices $\theta$ such that $\theta_{a,b}$ is non-zero if and only if there is an edge between a and b in G. Then, we define the criterion crit(G) by

$$ \mathrm{crit}(G) = \sum_{a=1}^{p} \left[ \|X_a - \sum_b X_b \hat{\theta}_{a,b}^G\|_n^2 \left( 1 + \frac{\mathrm{pen}[d_a(G)]}{n - d_a(G)} \right) \right], \quad (3) $$

where the penalty function is defined by

$$ \mathrm{pen}(d) = 1.1 \times \frac{n-d}{n-d-1} \mathrm{EDKhi} \left[ d+1, n-d-1, \left( \binom{p-1}{d} (d+1)^2 \right)^{-1} \right]. \quad (4) $$

The function EDKhi[d, N, .] is the inverse of the function

$$ \mathrm{DKhi}[d, N, .] : x \mapsto \mathbb{P} \left( F_{d+2,N} \ge \frac{x}{d+2} \right) - \frac{x}{d} \mathbb{P} \left( F_{d,N+2} \ge \frac{N+2}{Nd} x \right), $$

where $F_{d,N}$ denotes a Fisher random variable with $d$ and $N$ degrees of freedom. See [7] and [4] for more details and more explanations on this penalty. This Criterion (3) is implemented in a R library GGMselect [8].

# 3 Non local Gaussian Graphical Models

## 3.1 Introduction of $G_0$

As noticed above, the neighbouring points - neighbour nodes of the graph- are very likely to be conditionally correlated. We want to put more attention onto the other correlations -which will be called long-distance ones in the sequel since they do not affect the neighbours. To this purpose, we have introduced in the estimation a neighbouring graph $G_0$ which carries the neighbour nodes of all the graph nodes. We assume that there exist correlations between these points but we are not estimating them rather looking for the other ones. There are two reasons for that. On the one hand, we are not interested in the local correlations since they appear to be obvious. The long-distance ones however may reveal
some important behaviour which may be characteristic from the population we are studying. These non obvious relations between regions may show some non random effect on the shape. This will be illustrated in the experiments.

On the other hand, the value of the penalty is directly related to the complexity of the families of possible candidate graphs. Shrinking the families to those of graphs containing $G_0$ reduces their complexity and therefore the penalty. The effect it that more edges appear and therefore the long-distance ones.

The way we introduce $G_0$ in the estimation is as follows. The estimation of the conditional correlations is done in the orthogonal space of the neighbouring graph $G_0$. This leads to replace the random variable $X_a$ by

$$X_a - X_{m_a} (X_{m_a}^T X_{m_a})^{-1} X_{m_a}^T X_a, \qquad (5)$$

where $X^T$ is the transposition and $X_{m_a}$ is the matrix defined as follow: if we denote $X$ the $n \times p$ matrix of all the data and $m_a$ the list of neighbours of node $a$ in $G_0$, then $X_{m_a} = X(., m_a)$ is of dimension $n \times \text{card}(m_a)$.

This orthogonality constraint may have to be relaxed since the projection may lead to ill-conditioned matrices. Moreover, it will also enable to capture some edges with small projection onto the orthogonal of $G_0$ which may appear stronger since we know that the local dependencies are more likely to be the strongest and summarise the main information. This can be express in terms of introducing a small ridge (driven by a new parameter $\gamma_0$ chosen via the criterion). Denoting $Id$ the identity matrix, this yields

$$X_a - X_{m_a} (X_{m_a}^T X_{m_a} + \gamma_0 Id)^{-1} X_{m_a}^T X_a. \qquad (6)$$

## 3.2 Selection of graph: Change in the penalty

When we introduce an a priori graph $G_0$, we shall change the selection criterion to take into account this a priori knowledge. To fit with this situation, the Criterion (3) remains unchanged, but the penalty (4) is replaced by

$$\text{pen}_{G_0}(a, d) = 1.1 \times \frac{n-d}{n-d-1} \text{EDKhi} \left[ d+1, n-d-1, \left( \frac{p-d_0[a]-1}{d-d_0[a]} \right) (d-d_0[a]+1)^2 \right]^{-1} \qquad (7)$$

where $d_0[a]$ is the degree of the node $a$ in the graph $G_0$. This correction reflects the change of complexity induced by the graph $G_0$. It ensures a control of the prediction error similar to Theorem 1 in [7]. For further details, we refer to the Appendix. We have therefore modified GGMselect in two ways (which is included in the new update of the package). If a prior is given, we change the data so that the new random variable is given by eq. (5) or (6). Then we changed the penalty criterion for eq. (7).

## 4 Experiments

We have tested this model on meshed surfaces of hippocampi from Johns Hopkins University [16]. The data base contains $n = 101$ vectors of dimension $p = 713$
corresponding to the vertices of the graph. A template surface of a hippocam-
pus was mapped onto each subject surface using Large Diffeomorphic Defor-
mation Metric Mapping [17]. Prior rigid body registration has been done on
each subject. In this common coordinate system, the log of the absolute value of
the Jacobian determinant of this deformation was discretised onto the template
meshed surface. This leads to 101 vectors with real value signal which corre-
spond to scalar deformation fields over the surface. The subjects belong to three
sub-populations. The first 57 were control patients. The 32 next were suffering
from mild Alzheimer’s disease or semantic dementia (denoted mild AD) and the
12 last were in a late stage of the disease (denote late AD in the sequel). Some
examples of the training set are presented in Fig. 2.

<figure><img src="image_46.png" /><figcaption>**Fig. 2.** Examples of the training set. The colour depends on the intensity of the Jacobian of the deformation. Blue means a contraction and red dilatation. The intensity itself is not important but rather its relative value with respect to the others.</figcaption></figure>

4.1 Gaussian Graphical estimation without neighbour prior

The results of the graph estimations without any neighborhood prior are pre-
sented in the four graphs on the left hand side of Fig. 3. The four plots correspond
respectively to the solution of eq.(1) with both choice of symmetrisation (”or”
and ”and” resp.) and of eq.(2) again with both symmetrisation processes. For
each of these minimisation, we present the result with the lowest criterion. All
these graph present some similarities. Indeed, they all find the conditional corre-
lations between close nodes of the graph with respect to the euclidean distance.
This is what we expect since the deformations are smooth. Therefore, the defor-
mation at one point is close to the deformation on its neighbours leading to a
correlation between all nodes but only local conditional correlations. However,
one of the graphs (left one, with a slightly lower criterion) suggests that they
should be some conditional correlation not only between neighbour nodes. As we
can see, not all the part of the surface are linked through long edges. However,
this does not appear in the three other graphs which makes them uncertain.

4.2 Non local GGM estimation

To overcome the problem of the predominance of the conditional correlations
between close nodes, we have introduced a neighbouring graph $G_0$. Two examples
of $G_0$s are presented in Fig. 3 (right). The estimated graphs using our proposed
algorithm are presented in Fig. 4 where the top row corresponds to the first $G_0$ of Fig. 3 and the bottom row to the second one. All the eight estimations are presented: lasso (with the two possible symmetrisations) using the orthogonal projection and the ridge regularisation and next, in the same order, the four estimation with the elastic net energy. The orthogonal projection may appear too strong since it requires to invert the neighbouring matrix (cf eq. (5)). This yields no graph estimated with the lasso minimisation process with a strong symmetrisation constraint on the first line.

We can notice that they almost all carry some long-distance correlations. In particular, four of them show some common patterns pointing the same connected areas (row 1, columns 2,4,5,6). The estimation remains stable when adding or removing some edges in $G_0$. This is shown in the bottom row of Fig 4. The results are very similar. Noticing that the orthogonal of this new $G_0$ is smaller (the maximum degree of this $G_0$ is 7 whereas the first one was 3) the resulting graphs look alike with a little less edges. In this case, the estimation of the first graph is possible but presents very few connections. As well as for the first $G_0$, the lasso with a relaxed projection constrains provides very interesting graphs and the elastic net does not requires this relaxation. Since some graphs appear very dense, we removed the edges between close points w.r.t. the Euclidean distance. We did this pruning for the marked graph (see below for meaning of this mark) of Fig. 4. The result is presented in Fig. 5. Thanks to this representation, the long-distance conditional correlation appear clearer and can be easier interpreted. Indeed, this suggest that there is a strong conditional correlation between the bottom part of the head of the hippocampus and two regions of the body. This means that the hippocampus is subjected to deformations that are not random. When the bottom area of the head is deformed so are the two local parts of the body. Note that the conditional correlation can be positive or negative which correspond to same and opposite behaviour respectively.

Thanks to the criterion we can point the graph which, among the ones we estimate, has the best power for prediction purposes. This graph is marked with a red star. This graph looks indeed very interesting showing long-distance dependencies that are not trivial. The power of this criterion we propose is that

<figure><img src="image_47.png" /><figcaption>**Fig. 3.** Left: Estimated graphs using the lasso algorithm and both symmetrisations (called "or" and "and" in the sequel) and the elastic net algorithm with both symmetrisations as well ("or" and "and" resp.). The results show the predominance of the neighbours conditional correlations except for the first estimation where long-distance conditional correlation appear. Right: two examples of possible $G_0$s.</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>**Fig. 4.** Estimated graphs with the eight methods using the two neighbouring graphs presented in Fig. 3 on the top and bottom line resp. From left to right, each method with both symmetrisations ("or" and "and" resp.): lasso with the orthogonal projection, lasso with a relaxation on the projection, elastic net with projection, elastic net with relaxed projection. The projection constraint may be too strong preventing from computing the estimation (top left). In almost all graphs, long-distance dependencies appear and are stable with respect to the prior neighbours.</figcaption></figure>

it can quantify the predictive power of any other graph estimated by any other method. We not only provide some estimation technics which take into account a neighbourhood prior but also a way to compare them to other estimated graphs.

## 4.3 Clustering of the graph

In order to see the different areas which are conditionally correlated, we use the spectral clustering method which highlights these dependencies. Spectral clustering is a technique based on eigen-properties of a similarity matrix $Q$ that partitions data into disjoint clusters. We compute the $p$ eigenvectors of $Q$, $(V_i)_{1\le i\le p}$. Let the vectors $Y_j = (V_i^j)_{2\le i\le p}$ be the contatenation of the $j^{th}$ coordinates of each eigenvector. The clustering is done using a $k$-mean algorithm on these vectors. The choice of the similarity matrix is given as follows. If two nodes $i$ and $j$ are connected through the graph but not in the neighbouring graph $G_0$ then $Q(i,j) = 1/Z$. If this condition is not satisfied then $Q(i,j) = f(\|x_i - x_j\|)/Z$ where $Z$ is a renormalisation constant, $x_i$ is the pose of node $i$ and $f$ is the Gaussian density function with fixed variance (chosen with respect to the spreading of the data, here 50 whereas $\max(x_i - x_j)^2 = 7185$). $Q(i,i)$ is set so that the sum of the $i^{th}$ row is 0. The result using 5 clusters is presented in Fig. 5 where each node has a colour corresponding to its cluster. We can clearly see that the regions that are in the same cluster (same colour) correspond to areas that share a link in the pruned graph (Left of Fig. 5). This shows that several connected parts of the shape (red, black, green and dark blue) seem to have conditionally independent behaviour whereas some nodes (cyan) are clearly related. This result confirms the segmentation in [24].
<figure><img src="image_49.png" /><figcaption>**Fig. 5.** Best graph analysed. Left: the edges between close nodes have been removed for better visualisation. Right: k-mean clustering.</figcaption></figure>

## 4.4 Population different

Since we are provided with three sub-groups, it is interesting to see whether they carry differences with respect to conditional correlations. However, since the late AD group is pretty small, we only cluster it into two groups: controls versus AD (both mild and late together). The goal is to see if the disease changes the kind of conditional correlated deformations that are present in a population. This is actually what Fig 6 tends to show. The deformations of the hippocampus in a control population are not completely random since some important long distance conditional correlations appear. This reduces the safe deformations to a subset where these joint behaviour occurs. However, the AD population seems to have less long-distance conditional correlation. This would suggests that the disease affects the hippocampus by removing the direct dependencies of different regions. This effect cannot be seen when no neighbourhood obvious correlations are taken into account (graphs on the left of Fig. 6). However, when introducing the projection onto the orthogonal of a neighbouring graph (the first one presented above), the differences appear and are stable with respect to the estimation process (the criterions are of the same range for these four tests).

# 5 Conclusion

This paper presents a new way of analysing populations of shapes in terms of random graphs which carry sparse conditional correlations between areas of a shape. We have introduced a neighbourhood prior. It stabilises the estimation and highlights the long-distance conditional correlations which are the non obvious ones. This neighbouring graph is given by the user allowing for non trivial closeness notion as some anatomical ones. The results on the deformed hippocampi reveal some important conditional correlations between particular sub-regions which are stable along the estimation processes as well as with respect to a change in the neighbouring graph. Moreover, it emphasises differences between the control subjects which have more long-distance edges than the AD group. Provided with a larger database would help confirming this first trend. In
<figure><img src="image_50.png" /><figcaption>**Fig. 6.** Comparison of the estimated graphs from the two sub-groups. Left: estimations without a prior on the neighbours (lasso "and" and elastic net "or" resp.). Right: estimations using the first $G_0$ presented above (lasso "and" with orthogonal projection and the ridge relaxation, elastic net "or" with projection (when computable) and ridge relaxation resp.). The graph without the prior are very similar and does not enable any discrimination. Introducing the prior knowledge enables to catch the differences: the AD group present less long-distance conditional correlations.</figcaption></figure>

addition, we provide the user with a criterion which quantifies the quality of any graph (not only the one estimated with our algorithm). The estimation has been done here for anatomical shapes but is more general and it would be interested to test this on functional data to catch direct long-distance dependencies in the brain connectivity network and compare with for eg [23] and [20].

# 6 Appendix

The penalty pen$_{G_0}$ is a complexity penalty built from the theory of Giraud [7]. For a given integer $D$ and a graph $G_0$, we introduce the collection $\mathcal{M}_D(G_0)$ of graphs of degree less than $D$ and containing $G_0$. A careful inspection of the proof of Theorem 1 in [7] shows that for each node $a$ the penalty pen$_{G_0}(a, d)$ must ensures the control

$$ \sum_{m_a \in \mathcal{M}_0(a)} (|m_a| + 1) \text{DKhi} \left( |m_a| + 1, n - |m_a| - 1, \frac{(n - |m_a| - 1)\text{pen}(a, m_a)}{K(n - |m_a|)} \right) \le \log(n), $$

where $\mathcal{M}_0(a) = \{b : b \stackrel{G}{\sim} a\}, G \in \mathcal{M}_D(G_0)$. This control is achieved by eq (7).

# References

1. Allassonnière, S., Kuhn, E., Trouvé, A.: Bayesian deformable models building via stochastic approximation algorithm: A convergence study. Bernoulli J. 16(3), 641–678 (2010)

2. Allassonnière, S., Younes, L.: A stochastic algorithm for probabilistic independent component analysis. In revision for Annals of Applied Statistics

3. Ashburner, J., et al.: Identifying global anatomical differences: deformation-based morphometry. Human Brain Mapping 6(5-6), 348–357 (1998)
4. Baraud, Y., Giraud, C., Huet, S.: Gaussian model selection with an unknown variance. Annals of Statistics 37(2), 630–672 (2009)

5. Bickel, P., Levina, E.: Covariance regularization by thresholding. Annals of Statistics 36(6), 2577–2604 (2008)

6. Durrleman, S., et al.: Inferring brain variability from diffeomorphic deformations of currents: an integrative approach. Medical Image Analysis 12(5), 626–637 (2008)

7. Giraud, C.: Estimation of gaussian graphs by model selection. Electron. J. Stat. 2, 542–563 (2008)

8. Giraud, C., Huet, S., Verzelen, N.: Graph selection with ggmselect. arXiv:0907.0619

9. Jia, J., Yu, B.: On model selection consistency of the elastic net when p ≫ n (2008)

10. Kim, S.G., et al.: Structural connectivity via the tensor-based morphometry. Bayesian Meeting of Korean Society of Human Brain Mapping (2010)

11. Lauritzen, S.: Graphical models 17 (1996)

12. Lee, H., et al.: Discriminative persistent homology of brain networks (2011)

13. Lee, H., et al.: Sparse brain network using penalized linear regression. SPIE Medical Imaging (2011)

14. Lee, H.e.a.: Sparse brain network recovery under compressed sensing. IEEE Transactions on Medical Imaging (2011)

15. Meinshausen, N., Bühlmann, P.: High-dimensional graphs and variable selection with the lasso. Ann. Statist. 34(3), 1436–1462 (2006)

16. Miller, M., et., al.: Morphometry BIRN. collaborative computational anatomy: An MRI morphometry study of the human brain via diffeomorphic metric mapping. Human Brain Mapping 30(7), 2132–2141 (2009)

17. Miller, M.I., Trouvé, A., Younes, L.: On the metrics and Euler-Lagrange equations of computational anatomy. Annual Review of biomedical Engineering 4 (2002)

18. Sabuncu, M., Balci, S.K., Golland, P.: Discovering modes of an image population through mixture modeling. MICCAI LNCS(5242), 381–389 (2008)

19. Seeley, W., et al.: Neurodegenerative diseases target large-scale human brain networks 62(1), 42–52 (2009)

20. Smith, S., et al.: Network modelling methods for fMRI. NeuroImage 54(2), 875–891 (2011)

21. Tibshirani, R.: Regression shrinkage and selection via the lasso. J. Royal. Statist. Soc B 58 (1), 267–288 (1996)

22. Toewsa, M., et al.: Feature-based morphometry: Discovering group-related anatomical patterns. Neuroimage 49(3), 2318–2327 (2010)

23. Varoquaux, G., et al.: Brain covariance selection: better individual functional connectivity models using population prior. NIPS (2010)

24. Wang, L., et al.: Abnormalities of hippocampal surface structure in very mild dementia of the alzheimer type. Neuroimage 30:52-60 (2006)

25. Wang, Y., Chan, T.F., Toga, A.W., Thompson, P.M.: Shape analysis with multivariate tensor-based morphometry and holomorphic differentials (2009)

26. Wright, J., Ma, Y., Mairal, J., Sapiro, G., Huang, T., Yan, S.: Sparse representation for computer vision and pattern recognition. Proceedings of the IEEE (2010)

27. Zou, H., Hastie, T.: Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society, Series B (2005)
# An Airway Tree-shape Model for Geodesic Airway Branch Labeling

Aasa Feragen¹, Pechin Lo¹, Vladlena Gorbunova¹, Mads Nielsen¹, Asger Dirksen², Joseph M. Reinhardt³, François Lauze¹, Marleen de Bruijne¹ *

¹Department of Computer Science, University of Copenhagen, Denmark
{aasa, marleen}@diku.dk

²Department of Orthopedics and Internal Medicine, University of Copenhagen,
Denmark

³Department of Biomedical Engineering, University of Iowa, USA
http://image.diku.dk/aasa

**Abstract.** We present a mathematical airway tree-shape framework where airway trees are compared using geodesic distances. The framework consists of a rigorously defined shape space for treelike shapes, endowed with a metric such that the shape space is a geodesic metric space. This means that the distance between two tree-shapes can be realized as the length of the geodesic, or shortest deformation, connecting the two shapes. By computing geodesics between airway trees, as well as the corresponding airway deformation, we generate airway branch correspondences. Correspondences between an unlabeled airway tree and a set of labeled airway trees are combined with a voting scheme to perform automatic branch labeling of segmented airways from the challenging EXACT'09 test set. In spite of the varying quality of the data, we obtain robust labeling results.

**Keywords:** Airway branch registration, tree-shape model, airway shape model, tree metric, tree matching

## 1 Introduction

Medical imaging is an important diagnostic tool, and along with this tool comes a need for automatic analysis of medical images. Tree-structures are important in this context due to their roles as delivery systems for fluids and gases, which ties them directly and indirectly to a number of diseases. For instance, chronic obstructive pulmonary disease (COPD) is tied to properties of the airway, such as the airway wall thickness [11, 14].

In order to monitor progression of diseases and determine the range of normal variation in healthy anatomical trees, we need to compare measures between scans with varying characteristics. For instance, one needs to be able to compare measurements of airway dimensions, made at specified sites in the airway,

* This work is partially funded by the Lundbeck Foundation, the Danish Strategic Research Council (PSVT projects 09-065145 and 09-061346) and by AstraZeneca AB.
<figure><img src="image_51.png" /><figcaption>**Fig. 1.** ((a-b) Examples of segmented airways (CASE31 and CASE32 in the EXACT'09 dataset). The topological structures of the two trees are different, especially the circled "trifurcation-like" regions in the (image) left hand side, shown as trifurcations in fig. 3(a). c) By tracing points (illustrated by the circle, square and star) through the geodesic deformation between two tree-shapes, we obtain a registration of points and branches of the two endpoint tree-shapes. Branches which are collapsed during the tree deformation, are not traced further, as illustrated by the square.</figcaption></figure>

between several patients. In this paper, we provide a robust way to make such comparisons for airway trees by giving anatomical labels to branches in the airway based on the shape of the airway centerline tree. The technique uses a model of airway tree-shape which incorporates two notions: tree topology (i.e., parent/child connectivity) and geometry, defined by branch shapes (e.g., through landmark points).

Branch matching is, indeed, both a tree-topological and a geometric problem. Topological differences are particularly problematic for comparison and matching of anatomical trees such as the airways. These differences can come from noise in terms of spurious or missing branches, due to problems in the image recording and processing procedures. Topological differences can also come from anatomical variation between different patients, see figs. 1(a) and 1(b). Most current anatomical tree labeling methods focus, however, *either* on branch geometry or on tree topology. In this article we perform matching based on geodesics in an airway tree-space based on the tree-shape framework developed by Feragen et al. [3], where a varying tree topology becomes an integral part of the shape geometry. Not only is this framework geometrically very natural; it also handles tree-topological differences in a continuous, morphological way.

The airway tree model proposed in this paper consists of a shape space construction for treelike shapes, endowed with a geodesic metric. Any two treelike shapes are connected by a shortest possible deformation, or tree-space path (geodesic), whose length defines a distance between the shapes. Throughout the deformation, the initial tree changes its tree-topological structure to obtain an optimal match with the second tree. This makes the tree-shape framework well suited for branch registration in both inter- and intra-patient pairs of airways. The deformation, which is unique for sufficiently local data, induces a matching
of points along the tree-shape branches as illustrated in fig. 1(c). In a leave-one-
out fashion, we generate automatic branch labels on any given airway tree by
matching it with the other airway trees, which have been manually labeled by
an expert. The geodesic matching is combined with a branchwise vote among
anatomical labels induced by the matches, giving a robust automatic branch
labeling.

The main contributions of this paper are i) the adaptation of the tree-shape
framework of [3] to airway trees, giving a new version of the shape space; ii) a
thorough explanation of the underlying geometric ideas, making the tree-shape
model available for a broader community; iii) turning the computed geodesics
into an actual branch matching; and iv) the fusion of several branch matchings
through a voting scheme in order to obtain branch labelings.

The airway shape model has potential for applications beyond branch match-
ing and labeling. Labels or geodesic distances could be used to classify airways
into shape- and structure-dependent phenotypes, and the shape space frame-
work also opens for defining new biomarkers based on the whole airway shape
- topology and geometry combined. Moreover, the shape space framework used
here is very general, and can be transferred to other types of data, e.g., vascular
trees or medial axes, with little effort.

## 1.1 Related work.

Registration, branch matching and branch labeling in anatomical trees have been
studied in various ways for the past decade. Some of the most successful methods
are based on using association graphs [5,9,12]. Given two initial trees or graphs,
their association graph is a larger graph, which contains information from both
initial graphs. Branch matchings are induced by maximal cliques in the associ-
ation graph. The association graph and its maximal clique are predominantly
combinatorial constructions, although they can depend on geometric properties
of the initial trees. Separating the analysis of geometric and combinatorial prop-
erties like this is somewhat artificial, as the geometric and topological structures
play together in defining the efficiency of an anatomical tree as a space-filling
structure [8,15]. Moreover, finding the maximal clique is NP hard, making exact
computation intractable.

A more basic labeling method is given by van Ginneken et al. [13], where
an airway tree branch labeling is made recursively, starting at the trachea, as
part of the segmentation process. Labels are assigned using measures such as
radius, orientation and bifurcation angle. Such a labeling approach is likely to
be vulnerable to differences in topological structure, as is also noted by the
authors. A different approach is that of Kaftan et al. [4], who match tree *paths*
rather than branches. They avoid the difficulty with different tree-topological
structures, but also lose all information stored in the topological structure. Their
model does not seem to have applications beyond matching, and in particular
does not generate branch labels. Smeets et al. [10] match branches from lung
vessel trees using pairwise distances between nodes both in 3D Euclidean space
and along the tree to generate distance matrix "fingerprints", which are matched
in order to generate a matching. Bülow et al. [2] match airway tree branches without connectivity information, using only branch shape.

These methods all focus on one out of two properties of a tree-shape: tree topology [5, 9, 12], or branch-wise geometry [2, 4, 10]. However, the airway tree is both topology and branch geometry. This duality is precisely what makes matching and labeling difficult. Our airway tree-shape model is ideal for modeling airway trees because it considers topology and geometry simultaneously.

**The rest of the paper is organized as follows:** The airway tree-shape model and the tree-shape metric are presented in sec. 2; experimental details and results are presented in sec. 3, which are discussed in sec. 4 and concluded in sec. 5.

## 2 The tree-shape space

The proposed method uses airway centerlines as input. Each branch centerline is associated with an edge in a combinatorial airway tree, endowed with a hierarchical tree structure (parent/child) which describes branch connectivity. Shape information is represented by edgewise attributes, in this case a fixed number of equally spaced landmark points along each branch centerline (the equal spacing length varies from branch to branch).

We build a space of airway tree-shapes based on the tree-shape model defined by Feragen et al. [3]. Any tree-shape, such as an airway tree-shape, is represented as a pair $(T, f)$ of a combinatorial tree-structure $T$ with edgewise shape attributes $f$, see fig. 2. Here $T = (V, E, r)$ is a tree with vertices $V$, edges $E \subset V \times V$, and a root $r$. The shape attributes given by $n$ landmark points per edge are specified by a function $f: E \to \mathbb{R}^{3n}$.

$$ \mathcal{S} = \frac{\mathcal{A}(T)}{f} + \frac{(1, \mathcal{A}(V), \mathcal{A}(E), \mathcal{A}(r), -)}{f} $$

<caption>**Fig. 2.** A tree-shape is represented mathematically as a pair $(T, f)$ where $T$ is a binary, combinatorial tree and $f: E \to \mathbb{R}^{3n}$ assigns a shape-descriptor from $\mathbb{R}^{3n}$ to each edge.</caption>

### 2.1 Intuitive description of a geodesic tree-shape space

Our goal is to construct a continuous space of deformable trees, and this goal poses some constraints on the possible geometric structure of the space of tree-like shapes. To see this, first consider airway trees with a single fixed topological structure $T = (V, E, r)$. All such trees are described by a point in the same Euclidean product space $\prod_E \mathbb{R}^{3n}$, where each edge $e \in E$ has an associated attribute vector in $\mathbb{R}^{3n}$ describing its shape. Next, divide the set of all airway trees into classes with fixed topological structure. Each airway tree class lives
within its own Euclidean space $\prod_{E_i} \mathbb{R}^{3n}$, and the space of all airway tree-shapes
is a disjoint union $X = \bigsqcup (\prod_{E_i} \mathbb{R}^{3n})$ of shape spaces, one for each topological
structure. We consider each $\prod_{E_i} \mathbb{R}^{3n}$ as a component of the larger space of
all airway tree-shapes. We should tie these components together into a large
shape space that connects all airway trees. In order to understand how to tie
the components together, we consider how the tree-topological structure changes
throughout deformations.

<figure><img src="image_52.png" /><figcaption>**Fig. 3.** (a) This figure is best viewed in color. In the first generations of the airways, the branching structure is fairly pre-determined due to anatomy, and the branches have names. (b) The tree-shape $\mathcal{T}$ can be obtained as a limit of a sequence where a branch is disappearing – for instance, $\mathcal{T}$ can be obtained by shrinking a branch in either $\mathcal{T}_1$ or $\mathcal{T}_2$, or a number of other tree-shapes $\mathcal{T}_i$. This illustrates how trifurcations can be interpreted as pairs of bifurcations; their combinatorial representations are shown in the bottom row. This causes some problems in the representation space $X$ defined in sec. 2.2, since, for instance, the path from $\mathcal{T}_1$ to $\mathcal{T}$ to $\mathcal{T}_2$ is impossible in $X$. Passing from $\mathcal{T}_1$ or $\mathcal{T}_2$ to $\mathcal{T}$ is easy, but the two representations of $\mathcal{T}$ correspond to different points in $X$, and the path $\mathcal{T}_1 \to \mathcal{T} \to \mathcal{T}_2$ is not possible in $X$. This path is, however, possible in the quotient space $\tilde{X}$, where different representations of the same tree are glued together in one point. We define $\tilde{X}$ to be the space of treelike shapes.</figcaption></figure>

For example, in fig. 3(a), the LMB branch is shown as the parent of a tri-
furcation, with the LUL, LLB6 and LB6 branches emanating from it. In reality,
one will not find a trifurcation, but a pair of bifurcations, e.g., as seen in the
two airway subtrees $\mathcal{T}_1$ and $\mathcal{T}_2$ in fig. 3(b). A geodesic deformation from $\mathcal{T}_1$ to $\mathcal{T}_2$
should interchange the order of the LUL and LB6 branches by passing through
a tree of the type $\mathcal{T}$ shown in fig. 3(b). This tree can be obtained as a limit
of a sequence of trees with a fixed structure. For instance, such a sequence can
start at $\mathcal{T}_1$ or $\mathcal{T}_2$ and successively shrink the small branch until collapsing onto
one of the representations of $\mathcal{T}$. Similar sequences can be found for many differ-
ent topological structures, represented by tree-shapes $\mathcal{T}_i$, where slowly deleting
certain edges in $\mathcal{T}_i$ converges towards a tree of type $\mathcal{T}$. By following sequences
within different components of the tree-space, where all the tree-shapes in the $i^{th}$
component have the same structure as $\mathcal{T}_i$, we reach the same tree-shape $\mathcal{T}$. This must mean that in order to have a tree-space where the trees can continuously deform from one topological structure to another, as in fig. 3(b), the tree-space components of the $\mathcal{T}_i$ must all intersect along the component of $\mathcal{T}$. In the next section, we shall see how such a tree-space can be defined mathematically.

<figure><img src="image_53.png" /><figcaption>Fig. 4. Trees which are not binary, or which are smaller than $T$, are represented by cancelling extra edges, endowing them with the attribute zero (dotted lines). This leads to several representations of the same tree-shape. By identifying those different representations that represent the same shape, we construct our shape space as the quotient space $\bar{X}$.</figcaption></figure>

## 2.2 Mathematical definition of the tree-shape space

Having established an intuitive understanding of the geometry of the tree-space, we shall give a more concrete definition. Recall that any tree-shape is represented by a pair $(T, f)$ consisting of a combinatorial tree $T$ and edge shape attributes $f$. In order to study trees of different sizes and topologies in one unified setting, all shapes are represented by a fixed combinatorial tree $T = (V, E, r)$, which is large enough; e.g., if all the tree-shapes have depth $N$, then $T$ could be the full binary tree of depth $N$. Smaller trees are represented in the tree-space by endowing extra edges with zero attributes, see fig. 3(b). All trees are represented by a point in the Euclidean representation space

$$X = \prod_{e \in E} \mathbb{R}^{3n}. \qquad (1)$$

Assume, moreover, that $T$ is binary. All non-binary tree-shapes can be represented using a binary combinatorial tree by collapsing internal branches, as described in fig. 4. The space $X = \prod_{e \in E} \mathbb{R}^{3n}$ contains every tree-shape represented at least once. Some tree-shapes are even represented at several points in $X$, and as a consequence, some natural tree-shape deformations cannot be represented as paths in $X$, see fig. 3(b). This problem is solved by generating an equivalence relation $\sim$ on $X$ where different representations of the same shape are identified; now the space of tree-shapes is defined as the quotient space $\bar{X} = X/\sim$.

**Definition 1.** We say that two representations in $X$ define the "same tree-shape" when the following holds: Starting with the two tree-shape representatives $x_1 = (T, f_1)$ and $x_2 = (T, f_2)$, remove all branches with zero attribute, and consider the resulting (possibly no longer binary) tree-shape representations $(T_1, \tilde{f}_1)$ and $(T_2, \tilde{f}_2)$. Here, $T_1$ and $T_2$ are ordered, combinatorial, rooted trees. If, up to
a topology-preserving reordering of the branches, these two attributed trees are
exactly the same, then the tree-shapes are the same and the two representations
are defined to be equivalent. We write $x_1 \sim x_2$.

All equivalent points are identified by forming the quotient space $\bar{X}$ of the
equivalence relation [1], that is:

$$
\bar{X} = X / \sim = \{ \bar{x} | x \in X \}, \qquad (2)
$$

where

$$
\bar{x} = \{z \in X | z \sim x\}. \tag{3}
$$

The points in X̄ are equivalence classes x̄ in X – by abuse of notation we use x̄
to denote both points in X̄ and subsets of X. Geometrically, this construction
corresponds to gluing together those points in the representation space X that
represent the same tree-shape. As most of the points only have one represen-
tative, this creates a new space X̄ with self intersections at points with several
representatives. Many geometric properties of X are inherited by X̄. The quo-
tient space X̄ is the space of treelike shapes.

2.3 Metrics on the shape space

The Euclidean metric on X induces a standard quotient metric d on X̄ [1], called
Quotient Euclidean Distance (QED), defined as follows:

$$
d(\bar{x}, \bar{y}) = \inf_{k \in \mathbb{N}} \left\{ \sum_{i=1}^{k} \|a_i - b_i\| : a_1 \in \bar{x}, b_i \text{ equivalent to } a_{i+1}, b_k \in \bar{y} \right\}. \quad (4)
$$

Here, the norm $\| \cdot \|$ denotes the Euclidean norm on $X$, and $\bar{x}$ is the equivalence class of the point $x \in X$. Thus, when we write $a_1 \in \bar{x}$, this means that $a_1$ is a point belonging to the equivalence class $\bar{x}$ as a subset of $X$, or equivalently, $a_1 \sim x$.

One interpretation of eq. 4 is that a QED geodesic consists of a sequence
of k Euclidean lines, which are cut and concatenated whenever the geodesic
deformation switches between two representations bᵢ and aᵢ₊₁ of the same tree.
Typically, these identified points correspond to internal topological transitions in
the tree-shape structure. The infimum is taken over all possible concatenations
of lines for any k. E.g., the geodesic from T₁ to T₂ in fig. 3(b) consists of the
Euclidean line from T₁ to the first representative of T, concatenated with the
Euclidean line from the second representative of T, to T₂.

The QED metric is locally very well-behaved. In particular, geodesics between
data points and various forms of average shapes will exist and be unique [3].
This makes the tree-shape space and its geodesic deformations well suited for
registration. The tree-space construction described above is completely general,
and applies to any trees with continuous edge attributes, not only shapes.
## 2.4 Airway tree-shape space

In order to compute the geodesic connecting two 3D tree-shapes such as airways, we consider all possible branch orders on the trees in order to find an optimal branch alignment. The number of possible orders grows exponentially with the size of the trees, resulting in computational difficulties. In the case of airway trees, however, some branches are easy to identify, e.g., the main bronchi and some of the lobar bronchi seen in fig. 3(a). We significantly reduce computational complexity by identifying and fixing the branches $E_{fixed} = \{RMB, LMB, RUL, BronchInt, LB6, LLB6, LUL\}$, which are present and easy to identify in most data trees. In this way, we consider topological variation only in the lobar subtrees $T_{RUL}, T_{LUL}, T_{BronchInt}, T_{LLB6}, T_{LB6}$ following these, shown in colors in fig. 3(a). The airway shape space is then

$$ \left( \prod_{e \in E_{fixed}} \mathbb{R}^{3n} \right) \oplus \left( \prod_{e \in E_{RUL}} \mathbb{R}^{3n} \right) \oplus \cdots \oplus \left( \prod_{e \in E_{LLB6}} \mathbb{R}^{3n} \right) \quad (5) $$

where the first component is Euclidean and the others are quotient spaces. The symbol $\oplus$ denotes direct sum, or Cartesian product.

## 2.5 Computing geodesics

Computing the geodesics from eq. 4 is generally NP hard. The number $k$ of Euclidean concatenations that need to be checked will, in practice, be bounded for any given pair of trees, but it will grow exponentially with the number of edges in the trees. We make an approximation by bounding $k$ in eq. 4 for each lobar subtree. That is, we fix some number $K \in \mathbb{N}$, and approximate the distance in eq. 4 by computing

$$ d(\bar{x}, \bar{y}) = \inf_{k \le K} \left\{ \sum_{i=1}^{k} \|a_i - b_i\| : a_1 \in \bar{x}, b_i \text{ equivalent to } a_{i+1}, b_k \in \bar{y} \right\}. \quad (6) $$

This corresponds to assuming that the tree-shape represented throughout geodesic will undergo at most $K$ internal topological changes. In order to compute the approximated distance, a naive implementation is to list all the allowed combinations of topological changes, to compute the shortest version of a path going through each of those particular changes, and choose the shortest among the resulting paths. This is equivalent to Algorithm 1 from [3].

## 2.6 Branch label extraction

Using a geodesic deformation from a labeled airway tree to an unlabeled one, as seen in fig. 1(c), labels from a labeled tree can be propagated to the branches of an unlabeled tree.

For data with large variation in topology, we would normally need to choose a rather high bound $K$ on $k$ in eq. 6, in order to obtain the true geodesics between data points. This is punished by a large increase in computation time. We avoid
<figure><img src="image_54.png" /><figcaption>**Fig. 5.** The centerline trees of CASE26, CASE34 and CASE39, respectively. Note the large variation in size and topological structure; in particular CASE39 is missing the upper lobar trees.</figcaption></figure>

this problem by matching with multiple labeled trees, and thus extracting a number of candidate labels for each branch. This is followed by a voting scheme, where nearby trees, for which the approximation is good, are expected to win the vote. Given a training set of airway trees with manually assigned labels, a new airway tree is first matched with all the trees in the training set, and sets of labels are propagated to the branches of the unlabeled tree. Each branch is assigned a pool of labels, among which the majority vote is selected as a new, automatically extracted branch label.

## 3 Experimental validation

Experiments were made on airways from the EXACT'09 challenge test set [7], consisting of 20 CT scans from 15 different subjects, labeled as CASE21 - CASE40. The EXACT'09 challenge was a segmentation competition, and the CT scans come from a wide range of sources, states of breathing, and are processed using different methods [7]. As a result, the segmented airway trees display great variation in size, shape and noise level.

The airway trees were segmented from the CT scans using a voxel classification based airway tree segmentation algorithm by Lo et al. [6]. The centerlines were extracted from the segmented airway trees using a modified fast marching algorithm, directly giving a tree structure. Due to the segmentation method, the centerline trees are not connected, but have gaps at each bifurcation, as seen in fig. 5. This is not a problem, as the proposed matching method relies on the edge shapes, which are sufficiently well described in the disconnected model, along with parent/child connectivity information.

Leaves with segmented volume less than 10 mm³ were assumed to be noise and pruned away, and the centerlines were sampled with 6 equidistant landmark points along each edge. For each edge, the landmark points were translated so that the first landmark point was aligned with the origin. In this way, large differences in edges of low generation do not affect the whole subtree following them. To account for variation in size, each airway tree was normalized by the constant scaling factor 1/length(LMB), chosen since the LMB branch is present
<table><caption>**Table 1.** Results from the branch labeling; the percentage of correctly labeled branches among all branches labeled by the algorithm, and the number of branches correctly labeled by the algorithm. CASE39 is an outlier: both upper lobes were missing from the segmentation, making the algorithm fail. When the outlier is left out, we correctly label 83% of the branches on average.</caption><thead><tr><th>CASE</th><th>21</th><th>22</th><th>23</th><th>24</th><th>25</th><th>26</th><th>27</th><th>28</th><th>29</th><th>30</th></tr></thead><tbody><tr><th>% correct</th><td>75</td><td>88.2</td><td>92.9</td><td>80</td><td>77.8</td><td>86.7</td><td>88.9</td><td>94.4</td><td>66.7</td><td>89.5</td></tr><tr><th># correct</th><td>12</td><td>15</td><td>13</td><td>12</td><td>14</td><td>13</td><td>16</td><td>17</td><td>14</td><td>17</td></tr></tbody><thead><tr><th>CASE</th><th>31</th><th>32</th><th>33</th><th>34</th><th>35</th><th>36</th><th>37</th><th>38</th><th>39</th><th>40</th></tr></thead><tbody><tr><th>% correct</th><td>90</td><td>76.5</td><td>88.9</td><td>100</td><td>83.3</td><td>78.9</td><td>66.7</td><td>80</td><td>30</td><td>76.5</td></tr><tr><th># correct</th><td>18</td><td>13</td><td>16</td><td>13</td><td>15</td><td>15</td><td>12</td><td>8</td><td>4</td><td>13</td></tr></tbody></table>

and easy to measure in all segmentations. The trachea was left out in the experiment due to varying cut-off points. The next few branches (RMB, LMB, RUL, BronchInt, LB6, LLB6, LUL) were detected based on the orientation and extent of their subtrees. Branches from the first 6–7 generations were matched, whenever they were present in the segmentation.

## 3.1 Branch labeling

Each airway tree was aligned with all other airway trees in the dataset through computing geodesics with $K = 2$ for each lobar subtree (that is, permitting one structural transition in each of the five lobar subtrees). When computing the distances, not only the length of the geodesic from one tree to another, but also the geodesic from the first tree to the second were recorded. Using the geodesic we obtained a branch-wise matching by recording which edges are mapped where.

The reference airway trees from EXACT'09 were labeled manually by a trained image analyst. Up to 34 labels were assigned according to a standard nomenclature used in bronchoscopy. All cases were reviewed by a pulmonologist after labeling. For each pair of airway trees, a ground truth branch matching was induced from the anatomical labels. Several branches found by the registration were not labeled by the human experts, and were left out in evaluation.

The voting scheme was made based on the 20-airway dataset in a leave-one-out fashion, where each airway tree was matched with the remaining 19 and given a branch labeling based on the majority vote. In this way, labels down to the sixth generation of each airway tree were computed using the rest of the airway trees as training set. The best label for any given branch was obtained by voting among the labels. Branch labels with less than 55% voting consensus or less than 4 votes in total were discarded. The results are found in table 1; the overall success rate is 83%. We can obtain better scores by insisting on a higher majority vote; however, this implies labeling fewer branches.

# 4 Discussion

Tschirren et al. [12], van Ginneken et al. [13] and Bülow et al. [2] performed labeling on airway trees with success rates of 97%, 90% and 69%/40% (using
different features), respectively. Our success rate is 83%, which – taking the level of variation in the dataset into account – is high. Since the datasets used are all different, a direct comparison of percentages is not possible.

The large variance in our results is natural since the EXACT'09 data come from a wide range of sources, and are made with different subjects, scanners, scan- and reconstruction protocols etc, as described in [7, table 1]. This poses a great challenge: the structural differences between the segmented airway trees are very large, as illustrated in figs. 5 and 6, and the amount of noise is very different from tree to tree. We might have dealt with this problem by making a better approximation of the QED metric. Instead, based on the hypothesis that trees which are close together will be topologically similar, we choose to use a coarse but efficient approximation combined with a voting scheme.

<figure><img src="image_55.png" /><figcaption>**Fig. 6.** The number of detected branches in the EXACT'09 airway trees is very variable. Note also that this number is affected by the level of structural noise, and might not directly correspond to the number of anatomical branches.</figcaption></figure>

# 5 Conclusion

Based on a geometric model for shapes with a treelike structure, we have developed a technique for automatic inter- and intra-patient registration of airway tree centerlines. The method has been evaluated by performing an airway branch labeling on the EXACT'09 dataset. In spite of the variation in the dataset, our labeling results are good, illustrating the potential of the shape framework.

The shape space framework is very general, and the tree-space can be enriched with additional attributes if wanted. There are very few airway and vessel tree labeling algorithms available, most of which are ad hoc, specialized heuristics. In contrast, our proposed method takes a principled approach, which is easy to generalize to other biological tree-structures, e.g., vascular structures. This makes the proposed method both novel and important.

The airway shape model has potential for many applications beyond branch matching and labeling. The shape space has good properties for statistical analysis such as computation of means and modes of variation. Geodesic distances
could be used to classify airways into shape- and structure-dependent phenotypes, and the shape space framework also opens for defining new imaging biomarkers based on the whole airway shape – topology and geometry combined. These extensions are, however, by no means straightforward, as we are working in the non-smooth domain of tree-shape space. Further development of the shape space, numerical methods for tree-shape computations, and their applications are all topics of future work.

## References

1. M. R. Bridson and A. Haefliger. *Metric spaces of non-positive curvature*. Springer-Verlag, 1999.

2. T. Bülow, C. Lorenz, R. Wiemker, and J. Honko. Point based methods for automatic bronchial tree matching and labeling. In *SPIE Med. Im.*, volume 6143, pages 225–234, 2006.

3. A. Feragen, F. Lauze, P. Lo, M. de Bruijne, and M. Nielsen. Geometries in spaces of treelike shapes. In *ACCV*, pages 671–684, 2010.

4. J.N. Kaftan, A.P. Kiraly, D.P. Naidich, and C.L. Novak. A novel multi-purpose tree and path matching algorithm with application to airway trees. In *SPIE Med. Im.*, volume 6143, pages 215–224, 2006.

5. H. Kitaoka, Y. Park, J. Tschirren, J. M. Reinhardt, M. Sonka, G. McLennan, and E. A. Hoffman. Automated nomenclature labeling of the bronchial tree in 3D-CT lung images. In *MICCAI (2)*, pages 1–11, 2002.

6. P. Lo, J. Sporring, H. Ashraf, J.J.H. Pedersen, and M. de Bruijne. Vessel-guided airway tree segmentation: A voxel classification approach. *Medical Image Analysis*, 14(4):527–538, 2010.

7. P. Lo, B. van Ginneken, J.M. Reinhardt, and M. de Bruijne. Extraction of Airways from CT (EXACT09). In *The 2nd Intl WS Pulm. Im. Anal.*, pages 175–189, 2009.

8. B. Mauroy, M. Filoche, E. R. Weibel, and B. Sapoval. An optimal bronchial tree may be dangerous. *Nature*, 427:633–636, 2004.

9. J. H. Metzen, T. Kröger, A. Schenk, S. Zidowitz, H-O. Peitgen, and X. Jiang. Matching of anatomical tree structures for registration of medical images. *Im. Vis. Comp.*, 27:923–933, 2009.

10. D. Smeets, P. Bruyninckx, J. Keustermans, D. Vandermeulen, and P. Suetens. Robust matching of 3d lung vessel trees. In *MICCAI 2010, Proc 3rd Intl WS Pulm. Im. Anal.*, pages 61–70, 2010.

11. L. Sørensen, P. Lo, A. Dirksen, J. Petersen, and M. de Bruijne. Dissimilarity-based classification of anatomical tree structures. In *Information Processing in Medical Imaging*, LNCS, 2011.

12. J. Tschirren, G. McLennan, K. Palágyi, E. A. Hoffman, and M. Sonka. Matching and anatomical labeling of human airway tree. *IEEE TMI*, 24(12):1540–1547, 2005.

13. B. van Ginneken, W. Baggerman, and E. van Rikxoort. Robust segmentation and anatomical labeling of the airway tree from thoracic CT scans. In *MICCAI 2008*, pages 219–226, 2008.

14. G. R. Washko, T. Dransfield, R. S. J. Estepar, A. Diaz, S. Matsuoka, T. Yamashiro, H. Hatabu, E. K. Silverman, W. C. Bailey, and J. J. Reilly. Airway wall attenuation: a biomarker of airway disease in subjects with COPD. *J Appl Physiol*, 107(1):185–191, 2009.

15. E. R. Weibel. What makes a good lung? *Swiss Med. Weekly*, (139(27-28)):375–386, 2009.
Deformable Organisms and Error Learning for
Brain Segmentation

Gautam Prasad<fnref n="1" />,<fnref n="2" />, Anand A. Joshi<fnref n="3" />,

Albert Feng<fnref n="1" />, Marina Barysheva<fnref n="1" />, Katie L. McMahon<fnref n="4" />,

Greig I. de Zubicaray<fnref n="5" />, Nicholas G. Martin<fnref n="4" />, Margaret J. Wright<fnref n="4" />,

Arthur W. Toga<fnref n="1" />, Demetri Terzopoulos<fnref n="2" />, and Paul M. Thompson<fnref n="1" />

**Abstract.** Segmentation methods for medical images may not generalize well to different data sets or tasks, hampering their utility. We attempt to remedy these issues using deformable organisms to create an easily customizable segmentation plan. This plan is developed by borrowing ideas from artificial life to govern a set of deformable models that use control processes such as sensing, proactive planning, reactive behavior, and knowledge representation to segment an image. The image may have landmarks and features specific to that dataset; these may be easily incorporated into the plan. We validate this framework by creating a plan to locate the brain in 3D magnetic resonance images of the head (skull-stripping). This is important for surgical planning, understanding how diseases affect the brain, conducting longitudinal studies, registering brain data, and creating cortical surface models. Our plan dictates how deformable organisms find features in head images and cooperatively work to segment the brain. In addition, we use a method based on Adaboost to learn and correct errors in our segmentation. We tested our method on 630 T1-weighted images from healthy young adults, evaluating results using distance and overlap error metrics based on expert gold standard segmentations. We compare our segmentations with and without the error correction step; we also compare our results to three other widely used methods: BSE, BET, and the Hybrid Watershed algorithm. Our method had the least Hausdorff distance to expert segmentations on this dataset, but included slightly more non-brain voxels (false positives). Our framework captures diverse categories of information needed for skull-stripping, and produces competitive segmentations.

**Keywords:** deformable organisms, segmentation, MRI, Adaboost, Hausdorff, overlap, registration

<footnote n="1">¹Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA, USA</footnote>

<footnote n="2">²Department of Computer Science, UCLA, Los Angeles, CA, USA</footnote>

<footnote n="3">³Signal and Image Processing Institute, University of Southern California, Los Angeles, CA, USA</footnote>

<footnote n="4">⁴Queensland Institute of Medical Research, Brisbane, Australia</footnote>

<footnote n="5">⁵School of Psychology, University of Queensland, Brisbane, Australia</footnote>
# 1 Introduction

Deformable organisms label objects in images by integrating high level control mechanisms into a segmentation plan. They combine sensing, knowledge representation, reactive behavior, and proactive planning so our devised organisms may cooperatively segment an image. Deformable organisms were introduced into medical imaging by [13] who combined ideas from artificial life [29] and deformable models [14, 30]. Since their introduction, deformable organisms have been used for limb delineation [18], and segmentation of the spinal cord [16], vasculature [17], and corpus callosum in the brain [6]. [15] created a deformable organisms framework using the Insight Toolkit (ITK) [8], but we did not use it here, as its capabilities were too basic for our application. In contrast to several brain segmentation methods that work with low-level image processing and computer vision techniques, our deformable organisms can incorporate high-level knowledge and expectations regarding image data. In addition, almost every step of the plan presented here differs significantly from the one we previously presented [21].

Segmenting brain from non-brain tissues (such as the eyes, skull, scalp, and neck) in magnetic resonance (MR) images of the head is a vital pre-processing step for many types of image analysis. Accurate masks of the brain are helpful for longitudinal studies [22], for multi-subject analyses of brain structure and function [31], and as a pre-processing step prior to cortical surface modeling [32], surgical planning [4], and brain registration [34].

The process of segmenting brain versus non-brain tissue on MRI is commonly referred to as “skull-stripping” (although, strictly speaking, the skull generates almost no signal on T1-weighted MRI and the scalp and meninges are the main tissues removed). This has traditionally been done manually by trained experts, or by automated methods that are subsequently corrected by hand. Manually-created masks may also be used as gold standard delineations to validate performance of skull-stripping methods based on different principles. Though many approaches have been developed for this task, time consuming manual clean-up of these generated masks is almost always required. Most published methods do not perform well on all datasets, making improvements over existing methods critical.

There are a variety of existing skull-stripping methods. The Brain Extraction Tool (BET) [28] evolves a deformable model to find the boundary of the brain. It provides a robust way to find the boundary in unclear regions but does not incorporate prior knowledge of the brain's shape. The Brain Surface Extractor (BSE) [25] uses edge detection and morphological operations to find the brain/non-brain boundary. BSE quickly extracts the brain from an image but may include extra material in the mask, as it sometimes fails to remove connections between the brain and surrounding tissue. The Hybrid Watershed Algorithm (HWA) [24] uses the watershed algorithm to find the brain region, then fits a deformable model to the region, and finally deforms it based on a statistical atlas and geometric constraints. These methods have also been analyzed
in [1]. We chose these methods as they are the most widely used and are part of larger neuroimaging toolkits.

We create a deformable organism plan that governs a collection of organisms to segment different parts of the head and brain. The organisms evolve dynamically in the images and cooperatively compute an accurate and robust segmentation of the brain. We then use a learning method, based on an Adaboost wrapper, [33], to classify the error in our method. We evaluate the effectiveness of this additional error correction step in improving our segmentations. We test our method with 630 T1-weighted MR images from healthy young adults, aged 20-30 years. We compare our approach to three widely used methods and we validate our results using distance, overlap, and error metrics.

# 2 Methods

Our deformable organisms method aims to segment the brain in T1-weighted MR images of the head. We describe our deformable organism plan to segment the brain, a way to learn and correct errors in our method, validation metrics to compare our results to the gold standard and to other widely-used methods, and our experimental results.

## 2.1 Deformable Organisms

Deformable organisms are organized in five different layers that combine control mechanisms and different representations to segment an image. We adapt this general approach for segmenting the brain.

**Geometry and Physics** We represent each organism as a 3D triangulated mesh. These meshes are initialized on a standard brain template image. Our template was selected from the 40 images in the LONI Probabilistic Brain Atlas (LPBA40) [26], which have corresponding manual segmentations for 56 structures, and have manual delineations of the brain boundary. In the image we selected from this set, the voxels lying in each of our regions of interest are labeled. We fit our organisms to these labels to create a mesh using a marching cubes method [12] that goes through the image. The mesh is made up of polygons representing the border of the regions, which are then fused together.

These meshes deform to fit the 3D region that their corresponding organism is modeling. This iterative process moves each of the mesh's vertices along its normal direction with respect to the mesh surface. The surface is smoothed at every iteration using curvature weighted smoothing [2, 19]. This smoothing technique attenuates noise through a diffusion process as

$$ \frac{\partial \mathbf{S}}{\partial t} = \lambda L(\mathbf{S}), \qquad (1) $$

where S is the mesh, surface, or manifold and L is the Laplacian, which is equivalent to the total curvature of the surface. This Laplacian is linearly approximated
by

$$L(\mathbf{v}_i) = w_{ji} \sum \mathbf{v}_j - \mathbf{v}_{i_j \in N(i)}, \quad (2)$$

where $v_i$ is the vertex $i$ in the mesh, $N(i)$ are the neighbors of $i$, and $w_{ji}$ is a weight proportional to the curvature between vertices $j$ and $i$. We smooth the mesh to constrain its deformations to prevent intersections and artifacts from corrupting the boundary. On the boundary, we sample from the surrounding image.

Our prior based on an image from LPBA40 precludes us from using the Segmentation Validation Engine tool [27] to validate our method as it would lead to an unfair bias that would be favorable to our method; instead we use other metrics (below).

**Perception** Our organisms “sense” the encompassing image by sampling its intensities at vertices of the mesh. The vertices composing the mesh have real-valued coordinates, so we used nearest neighbor interpolation to find the intensities that correspond to them in the discrete grid of voxels in the image. The images may be any of the subject-derived volumes, which include the threshold image, 2-means classified image, 3-means classified image, or gradient image.

Paramount to the perception layer is the organisms' ability to sense each other's locations. We locate the voxels an organism resides on using a 3D rasterization algorithm [20, 5] that efficiently computes these values. These locations allow our organisms to dynamically change the way they deform based on their own positions and intensities of the original T1-weighted subject image.

**Motor Control** We move the vertices of the mesh along their normal direction with respect to the mesh surface by analyzing a set of intensities along this normal line. We describe the evolution of our mesh or surface $S(i, t)$ with respect to time $t$, where $i$ is a vertex or point on the surface, as

$$\frac{\partial S(i, t)}{\partial t} = F(\mathbf{P}, \mathbf{n}, \mathbf{I}_d, b(c, l)_j)\mathbf{n}, \quad (3)$$

with $F$ being the speed of evolution. $F$ samples a set of positions $P$ along point $i$'s normal and interpolates these values from any of the derived images $I_d$, where $d$ specifies the set of derived image. That set consists of the threshold image ($t$), 2-means image (2), 3-means image (3), and gradient image ($g$). The function $b(c)_j$ specifies any of a number of behaviors and decides how to move the point $i$ on the surface by analyzing the sampled intensities $l$ subject to a set of constraints, and weights the movement by the scalar $c$.

In practice, we evolve each vertex by the amount specified by $F$ along its normal and progress through time by iterating through all vertices in the mesh until there is no longer any significant movement.
**Behavior** Our behaviors are a higher level of abstraction to indicate how the organisms function and what information they need to find. Behaviors may prescribe a function for organisms to be attracted to or repelled from landmarks or help converge on the boundary of an object. The functions for our behaviors had specific tasks in mind in the context of skull-stripping but are general and simple enough for repeated use in any segmentation task.

We developed six behaviors that play an important role in almost every step of our skull-stripping plan.

1. We create a behavior that analyzes a binary image and locates a boundary in these images. It contracts if a vertex corresponds to an off value and expands if it corresponds to an on value, and may be described as

$$b(c, l)_1 = \begin{cases} -c & \text{if } l_i = 0 \\ c & \text{otherwise} \end{cases} \qquad (4)$$

In this case the set $l$ consists of a single value $l_i$, the value of the binary image that corresponds to the vertex $i$.

2. Our second behavior moves a vertex outwards if its corresponding intensity value is $q$ and may be described as

$$b(c, l)_2 = \begin{cases} c & \text{if } l_i = q \\ 0 & \text{otherwise} \end{cases} \qquad (5)$$

Its purpose is to expand into an area of an image with voxels having intensities $q$. In addition $q$ may be a set of labels that are appropriate for expansion.

3. Our third behavior is customized to move our skin mesh through the skull and skin into the cerebrospinal fluid (CSF). It contracts itself further if the boundary intersects another organism (in our plan, the other organisms are the eyes and we wish to deform through those areas) and will check the intensities ($l$) along the normal for those that correspond to CSF. If CSF markers $m$ are found (signified by very low intensity values and a specific label in our $k$-means images) we contract the mesh, and if they are deficient, we expand the mesh. This behavior may be represented as

$$b(c, l)_3 = \begin{cases} -c & \text{if } l_i \text{ intersects other organism} \\ -c & \text{if } m \in l \\ c & \text{if } m \notin l \\ 0 & \text{otherwise} \end{cases} \qquad (6)$$

In our framework, the intensities, $l$, are sampled from the $k$-means images and $m$ is the label corresponding to CSF. The sampled points are locations inside the surface with respect to $i$.

4. the fourth behavior we created was designed for the brain organism to locate the boundary of the brain. It contracts vertices away from other organisms, contracts if there are CSF marker values ($m$) present within the surface, expands if the label value at $i$ is not $q$, and expands if the gradient intensity
at $l_i$ is greater than or equal to a threshold $t$. The precedence of these constraints is ordered as

$$b(c, l)_4 = \begin{cases} -c & \text{if } l_i \text{ intersects other organism} \\ -c & \text{if } m \in l \\ c & \text{if } l_i \neq q \\ c & \text{if } l_i \leq t \\ 0 & \text{otherwise} \end{cases} \qquad (7)$$

The values at $l$ are sampled from either the 3D rasterization image of the other organisms' meshes, the 3-means classified image, or the gradient image.

5. Our final behavior is the same as behavior 2 but instead of expanding it contracts.

$$b(c, l)_5 = \begin{cases} -c & \text{if } l_i = q \\ 0 & \text{otherwise} \end{cases} \qquad (8)$$

**Cognition** We create a plan of different behaviors to perform a segmentation task. The plan may dynamically activate different behaviors depending on what features the organisms were able to find in the image. Our plan to skull-strip the brain is one such plan.

## 2.2 Skull-Stripping Plan

Our skull-stripping plan combines our image processing and deformable organisms to create objectives in the image to extract the locations of different regions, culminating in extracting the boundary of the brain. In what follows, we describe each step in detail and how it depends on previous knowledge obtained by organisms. This is just one plan and may be fashioned for any type of segmentation or specifics of the data. Table 1 summarizes the steps each organism takes during the segmentation.

1. We begin by registering the subjects T1-weighted MR image to the template we selected from the LPBA40. This registration step is important to transform subject images into a standard coordinate space as our organisms are tuned (iterations for deformations and labels for k-means classification) to images roughly corresponding to our template. Our template incorporates prior information and may be changed by users who need something closer to their own data. It provides initial locations and shapes for our skin, eye, and brain organisms.

We used an affine transformation for registration provided by FMRIB's Linear Image Registration Tool (FLIRT) [9]. FLIRT uses the correlation ratio [23] as the metric between the two images that takes the form

$$\frac{1}{\mathrm{Var}(Y)} \sum_{k} \frac{n_k}{N} \mathrm{Var}(Y_k). \qquad (9)$$

$Y$ represents one of the images, $\mathrm{Var}(Y)$ is the variance of $Y$, $Y_k$ is the k-th iso-set i.e. the set of intensities in $Y$ at positions where the other image $X$ has