Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2862393

Formats

Article sections

- Abstract
- I. Introduction
- II. Background - Demons Algorithm
- III. Spherical Demons
- IV. Experiments
- V. Discussion
- VI. Conclusion
- References

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2011 March 1.

Published in final edited form as:

Published online 2009 August 25. doi: 10.1109/TMI.2009.2030797

PMCID: PMC2862393

NIHMSID: NIHMS191511

B.T. Thomas Yeo,^{*} Mert R. Sabuncu,^{*} Tom Vercauteren, Nicholas Ayache, Bruce Fischl, and Polina Golland

B.T. Thomas Yeo, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, USA;

B.T. Thomas Yeo: ude.tim.liasc@samohty; Mert R. Sabuncu: ude.tim.liasc@ucnubasm; Tom Vercauteren: moc.hcetaekanuam@neretuacrev.mot; Nicholas Ayache: rf.airni.aihpos@ehcaya.salohcin; Bruce Fischl: ude.dravrah.hgm.rmn@lhcsif; Polina Golland: ude.tim.liasc@anilop

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

We present the Spherical Demons algorithm for registering two spherical images. By exploiting spherical vector spline interpolation theory, we show that a large class of regularizors for the modified Demons objective function can be efficiently approximated on the sphere using iterative smoothing. Based on one parameter subgroups of diffeomorphisms, the resulting registration is diffeomorphic and fast. The Spherical Demons algorithm can also be modified to register a given spherical image to a probabilistic atlas. We demonstrate two variants of the algorithm corresponding to warping the atlas or warping the subject. Registration of a cortical surface mesh to an atlas mesh, both with more than 160k nodes requires less than 5 minutes when warping the atlas and less than 3 minutes when warping the subject on a Xeon 3.2GHz single processor machine. This is comparable to the fastest non-diffeomorphic landmark-free surface registration algorithms. Furthermore, the accuracy of our method compares favorably to the popular FreeSurfer registration algorithm. We validate the technique in two different applications that use registration to transfer segmentation labels onto a new image: (1) parcellation of in-vivo cortical surfaces and (2) Brodmann area localization in ex-vivo cortical surfaces.

Motivated by many successful applications of the spherical representation of the cerebral cortex, this paper addresses the problem of registering two spherical images. Cortical folding patterns have been shown to correlate with both cytoarchitectural [25], [68] and functional regions [64], [27]. In group studies of cortical structure and function, determining corresponding folds across subjects is therefore important. There has been much effort focused on registering cortical surfaces in 3D [22], [23], [30], [58]. Since cortical areas – both structure and function – are arranged in a mosaic across the cortical surface, an alternative approach is to model the surface as a 2D closed manifold in 3D and to warp the underlying spherical coordinate system [27], [50], [59], [60], [64], [67]. Warping the spherical coordinate system establishes correspondences across the surfaces *without* actually deforming the surfaces in 3D.

There is frequently a need for invertible deformations that preserve the topology of structural or functional regions across subjects. Unfortunately, this causes many spherical warping algorithms to be computationally expensive. Previously demonstrated methods for cortical registration [27], [60], [67] rely on soft regularization constraints to encourage invertibility. These require unfolding the mesh triangles, or limit the size of optimization steps to achieve invertibility [27], [67]. Elegant regularization penalties that guarantee invertibility exist [5], [46] but they explicitly rely on special properties of the Euclidean image space that do not hold for the sphere.

An alternative approach to achieving invertibility is to work in the group of diffeomorphisms [4], [7], [9], [22], [31], [43], [66]. In this case, the underlying theory of flows of vector fields can be extended to manifolds [11], [44], [47]. The Large Deformation Diffeomorphic Metric Mapping (LDDMM) [7], [9], [22], [31], [43] is a popular framework that seeks a time-varying velocity field representation of a diffeomorphism. Because LDDMM optimizes over the entire path of the diffeomorphism, the resulting method is slow and memory intensive. By contrast, Ashburner [4] and Hernandez *et al.* [33] consider diffeomorphic transformations parameterized by a single stationary velocity field. While restricting the space of solutions reduces the memory needs relative to LDDMM, these algorithms still have to consider the entire trajectory of the deformation induced by the velocity field when computing the gradients of the objective function, leading to a long run time. We note that recent algorithmic advances [34], [43] promise to improve the speed and relieve the memory requirements of both LDDMM and the stationary velocity approach.

In this work, we adopt the approach of the Diffeomorphic Demons algorithm [66], demonstrated in the Euclidean image space, which constructs the deformation space that contains compositions of diffeomorphisms, each of which is parameterized by a stationary velocity field. Unlike the Euclidean Diffeomorphic Demons, the Spherical Demons algorithm utilizes velocity vectors tangent to the sphere and not arbitrary 3D vectors. This constraint need not be taken care of explicitly in our algorithm since we directly work with the tangent spaces. In each iteration, the method greedily seeks the locally optimal diffeomorphism to be composed with the current transformation. As a result, the approach is much faster than LDDMM [7], [9], [22], [31] and its simplifications [4], [33]. A drawback is that the path of deformation is no longer a geodesic in the group of diffeomorphisms.

Another challenge in registration is the tradeoff between the image similarity measure and the regularization in the objective function. Since most types of regularization favor smooth deformations, the gradient computation is complicated by the need to take into account the deformation in neighboring regions. For Euclidean images, the popular Demons algorithm [57] can be interpreted as optimizing an objective function with two regularization terms [14], [66]. The special form of the objective function facilitates a fast two-step optimization where the second step handles the warp regularization via a single convolution with a smoothing filter.

Using spherical vector spline interpolation theory [31] and other differential geometric tools, we show that the two-stage optimization procedure of Demons can be efficiently approximated on the sphere. We note that the problem is not trivial since tangent vectors at different points on the sphere are not directly comparable. We also emphasize that this decoupling of the image similarity and the warp regularization could also be accomplished with a different space of admissible warps, e.g., spherical thin plate splines [72].

Yet another reason why spherical image registration is slow is because of the difficulty in performing interpolation on a spherical grid, unlike a regular Euclidean grid. In this paper, we use existing methods for interpolation, requiring about one second to interpolate data from a spherical mesh of 160k vertices onto another spherical mesh of 160k vertices. Recent work on using different coordinate charts of the sphere [63] promises to further speed up our implementation of the Spherical Demons algorithm.

While most discussion in this paper concentrates on pairwise registration of spherical images, the proposed Spherical Demons algorithm can be modified to incorporate a probabilistic atlas. We derive and implement two variants of the algorithm for registration to an atlas corresponding to whether we warp the atlas or the subject. On a Xeon 3.2GHz single processor machine, registration of a cortical surface mesh to an atlas mesh, both with more than 160k nodes, requires less than 5 minutes when warping the atlas and less than 3 minutes when warping the subject. Note that the registration runtime reported includes registration components dealing with rotation, which takes up one quarter of the total runtime. The total runtime is comparable to other nonlinear landmark-free cortical surface registration algorithms whose runtime ranges from minutes [23], [60] to more than an hour [27], [67]. However, the other fast algorithms suffer from folding spherical triangles [60] and intersecting triangles in 3D [23] since only soft constraints are used. No runtime comparison can be made with spherical registration algorithm of the LDDMM type because to the best of our knowledge, no landmark-free LDDMM spherical registration algorithm that handles cortical surfaces has been developed yet.

Unlike landmark-based methods for surface registration [8], [22], [31], [50], [58], [64], we do not assume the existence of corresponding landmarks. Landmark-free methods have the advantage of allowing for a fully automatic processing and analysis of medical images. Unfortunately, landmark-free registration is also more challenging, because no information about correspondences are provided. The difficulty is exacerbated for the cerebral cortex since different sulci and gyri appear locally similar. Nevertheless, we demonstrate that our algorithm is accurate in both cortical parcellation and cytoarchitectonic localization applications.

The Spherical Demons algorithm for registering cortical surfaces presented here does not take into account the metric properties of the original cortical surface. FreeSurfer [27] uses a regularization that penalizes deformation of the spherical coordinate system based on the distortion computed on the original cortical surface. Thompson *et al.* [59] suggest the use of Christoffel symbols [39] to correct for the metric distortion of the initial spherical coordinate system during the registration process. However, it is unclear whether correcting for the metric properties of the cortex is important in practice, since we demonstrate that the accuracy of the Spherical Demons algorithm compares favorably to that of FreeSurfer. A possible reason is that we initialize the registration with a spherical parametrization that minimizes metric distortion between the spherical representation and the original cortical surface [27].

This paper is organized as follows. In the next section, we discuss the classical Demons algorithm [57] and its diffeomorphic variant [66]. In Section III, we present the extension of the Diffeomorphic Demons algorithm to the sphere. We conclude with experiments in Section IV and further discussion in Section V. The appendices provide technical and implementation details of the Spherical Demons algorithm and the extension to probabilistic atlases. This paper extends a previously presented conference article [69] and contains detailed derivations and discussions that were left out in the conference version. We note that our Spherical Demons code is freely available^{1}. To summarize,

- We demonstrate that the Demons algorithm can be efficiently extended to the sphere.
- We validate our algorithm by demonstrating an accuracy comparable to that of the popular FreeSurfer algorithm [27] in two different applications.

We choose to work with the modified Demons objective function [14], [66]. Let *F* be the fixed image, *M* be the moving image and Γ be the desired transformation that deforms the moving image *M* to match the fixed image *F*. Throughout this paper, we assume that *F* and *M* are scalar images, even though it is easy to extend the results to vector images [70]. We introduce a hidden transformation and seek

Data: A fixed image F and moving image M. |

Result: Transformation Γ so that M Γ is “close” to F. |

Set ^{0} = identity transformation (or some a-priori transformation, e.g., from a previous registration) |

repeat |

Step 1. Given ^{(t)}, |

Minimize the first two terms of Eq. (3) |

$${u}^{(t)}=\underset{u}{\text{argmin}}{\Vert {\sum}^{-1}(F-M\{{(t)}^{u\}}\Vert 2+\frac{1}{{\sigma}_{x}^{2}}\text{dist}({(t)}^{,}}^{}$$ (1) |

where u is any admissible transformation. Compute Γ^{(t)} = ^{(t)} u^{(t)}. |

Step 2. Given Γ^{(t)}, |

Minimize the last two terms of Eq. (3): |

$${(t+1)}^{=}$$ (2) |

until
convergence; |

$$\begin{array}{l}({,{\Gamma}^{})}^{=}\end{array}$$

(3)

In this case, the fixed image *F* and warped moving image *M* Γ are treated as *N* × 1 vectors. Typically, dist(, Γ) = ‖ − Γ‖^{2}, encouraging the resulting transformation Γ to be close to the hidden transformation and Reg() = ‖ ( − Id)‖^{2}, i.e., the regularization penalizes the gradient magnitude of the displacement field − Id of the hidden transformation . *σ _{x}* and

This formulation facilitates a two-step optimization procedure that alternately optimizes the first two (first and second) and last two (second and third) terms of Eq. (3). Starting from an initial displacement field ^{0}, the Demons algorithm iteratively seeks an update transformation to be composed with the current estimate, as summarized in Algorithm 1.

In the original Demons algorithm [57], the space of admissible warps includes all 3D displacement fields, and the composition operator corresponds to the addition of displacement fields. The resulting transformation might therefore be not invertible. In the Diffeormorphic Demons algorithm [66], the update *u* is a diffeormorphism from ^{3} to ^{3} parameterized by a stationary velocity field . Note that is a function that associates a tangent vector with each point in ^{3}. Under certain mild smoothness conditions, a stationary velocity field is related to a diffeomorphism through the exponential mapping *u* = exp(). In this case, the stationary ODE *x*(*t*)/*t* = (*x*(*t*)) with the initial condition *x*(0) ^{3} yields exp() as a solution at time 1, i.e., *x*(1) = exp()(*x*(0)) ^{3}. In this case, exp()(*x*(0)) maps point *x*(0) to point *x*(1).

The Demons algorithm and its variants are fast because for certain forms of dist(, Γ) and Reg(), Step 1 reduces to a non-linear least-squares problem that can be efficiently minimized via Gauss-Newton optimization and Step 2 can be solved by a single convolution of the displacement field Γ with a smoothing kernel. The proof of the reduction of Step 2 to a smoothing operation is illuminating and holds for dist(, Γ) = ‖ − Γ‖^{2} and any Sobolev norm Reg() = Σ_{i}
*σ _{i}*‖

In this section, we demonstrate suitable choices of dist(, Γ) and Reg() that lead to efficient optimization of the modified Demons objective function in Eq. (3) on the unit sphere *S*^{2}. We construct updates *u* as diffeomorphisms from *S*^{2} to *S*^{2} parameterized by a stationary velocity field . We emphasize that unlike Diffeomorphic Demons [66], is a tangent vector field on the sphere and not an arbitrary 3D vector field. A glossary of common terms used throughout the paper is found in Table I.

Suppose the transformations Γ and map a point *x _{n}*

Let *T*_{xn}
*S*^{2} be the tangent space at *x*_{n}. We define _{n}*T*_{xn}
*S*^{2} to be the tangent vector at *x _{n}* pointing along the great circle connecting

$${G}_{n}=\left(\begin{array}{ccc}0& -{x}_{n}(3)& {x}_{n}(2)\\ {x}_{n}(3)& 0& -{x}_{n}(1)\\ -{x}_{n}(2)& {x}_{n}(1)& 0\end{array}\right),$$

(4)

where *x _{n}*(

$${\overrightarrow{\Gamma}}_{n}=-{x}_{n}\times \left({x}_{n}\times \Gamma \left({x}_{n}\right)\right)=-{G}_{n}^{2}\Gamma \left({x}_{n}\right).$$

(5)

A more intuitive choice for the length of * _{n}* might be the geodesic distance between

Given *N* vertices
${\{{x}_{n}\}}_{n=1}^{N}$, the set of transformed points
${\{\Gamma ({x}_{n})\}}_{n=1}^{N}$ – or equivalently the tangent vectors
${\{{\overrightarrow{\Gamma}}_{n}\}}_{n=1}^{N}$ – together with a choice of an interpolation function define the transformation Γ everywhere on *S*^{2}. Similarly, we can define the transformation or the equivalent tangent vector field ϒ⃗ through a set of *N* tangent vectors
${\{{\stackrel{n}{\to}}_{\}}n=1N}_{}^{}$. We emphasize that these tangent vector fields are simply a convenient representation of the transformations and Γ and should not be confused with the stationary velocity field that will be used later on. We now set

$$\text{dist}(,\Gamma )=\sum _{n=1}^{N}{\Vert {\stackrel{n}{\to}}_{-}\Vert 2}^{},$$

(6)

which is well-defined since both _{n} and ϒ⃗_{n} belong to *T*_{xn}
*S*^{2} for each *n* = 1,…, *N*.

In this section, we show that the update for Step 1 of the Spherical Demons algorithm can be computed independently for each vertex. With our choice of dist(, Γ), step 1 of the algorithm becomes a minimization with respect to the velocity field
$\overrightarrow{\upsilon}{\{{n}_{{T}_{{x}_{n}}{S}^{2}}n=1N}_{}^{}$. By substituting *u* = exp() and
$\text{dist}(,\Gamma )={\sum}_{n=1}^{N}\left|\right|{\stackrel{n}{\to}}_{}-{\overrightarrow{\Gamma}}_{n}|{|}^{2}$ into Eq. (1), we obtain

$${\overrightarrow{\upsilon}}^{(t)}=\underset{\overrightarrow{\upsilon}}{\text{argmin}}\phantom{\rule{0.2em}{0ex}}f(\overrightarrow{\upsilon})$$

(7)

$$\begin{array}{l}=\underset{\overrightarrow{\upsilon}}{\text{argmin}}{\Vert {\sum}^{-1}(F-M\{{(t)}^{\text{exp}(\overrightarrow{\upsilon})\}}\Vert 2}^{\phantom{\rule{7em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}\text{dist}({(t)}^{,}}\end{array}$$

(8)

$$\begin{array}{l}=\underset{\overrightarrow{\upsilon}}{\text{argmin}}\sum _{n=1}^{N}\frac{1}{{\sigma}_{n}^{2}}{(F({x}_{n})-M\{{(t)}^{\text{exp}(\overrightarrow{\upsilon})\}({x}_{n})}2}^{}\phantom{\rule{4em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}\sum _{n=1}^{N}{\Vert {\stackrel{n}{\to}}_{+}^{{G}_{n}^{2}}2,}^{}\end{array}$$

(9)

where
${\sigma}_{n}^{2}$ is the *n*-th diagonal entry of Σ and denotes warp composition.

The cost function in Eq. (9) is a mapping from the tangent bundle *TS*^{2} to the real numbers . We can think of each tangent vector * _{n}* as a 3 × 1 vector in

Coordinate chart of the sphere *S*^{2}. The chart allows a reparameterization of the constrained optimization problem *f* in step 1 of the Spherical Demons algorithm into an unconstrained one.

It is a well-known fact in differential geometry that covering *S*^{2} requires at least two coordinate charts. Since the tools of differential geometry are coordinate-free [39], [44], our results are independent of the choice of the coordinate charts. Let *$\stackrel{\u20d7}{e}$ ^{n}*

$${\Psi}_{n}({x}^{\prime})=\frac{{x}_{n}+{E}_{n}{x}^{\prime}}{\Vert {x}_{n}+{E}_{n}{x}^{\prime}\Vert},\text{where}\phantom{\rule{0.3em}{0ex}}{E}_{n}=[{\overrightarrow{e}}^{n1}\phantom{\rule{0.3em}{0ex}}{\overrightarrow{e}}^{n2}].$$

(10)

As illustrated in Fig. 2, Ψ* _{n}*(0) =

$${\overrightarrow{\upsilon}}_{n}=D{\Psi}_{n}(0){\overrightarrow{z}}_{n}$$

(11)

$$=\frac{{I}_{3\times 3}-{\Psi}_{n}(0){\Psi}_{n}^{T}(0)}{\Vert {\Psi}_{n}(0)\Vert}{E}_{n}{\overrightarrow{z}}_{n}$$

(12)

$$=\frac{{I}_{3\times 3}-{x}_{n}{\mathit{\text{x}}}_{n}^{T}}{\Vert {x}_{n}\Vert}{E}_{n}{\overrightarrow{z}}_{n}$$

(13)

$$={E}_{n}{\overrightarrow{z}}_{n}=[{\overrightarrow{e}}^{n1}\phantom{\rule{0.3em}{0ex}}{\overrightarrow{e}}^{n2}]{\overrightarrow{z}}_{n}.$$

(14)

The above equation defines the mapping of a tangent vector *$\stackrel{\u20d7}{z}$ _{n}* at the origin of

From Eq. (14), we obtain exp() = exp({* _{n}*}) = exp({

$$\begin{array}{l}\{{\overrightarrow{z}}_{n}^{(t)}\}\\ =\underset{\{{\overrightarrow{z}}_{n}\}}{\text{argmin}}\sum _{n=1}^{N}\frac{1}{{\sigma}_{n}^{2}}{(F\left({x}_{n}\right)-M\{{(t)}^{\text{exp}\left(\{{E}_{n}{\overrightarrow{z}}_{n}\}\right)}\left({x}_{n}\right))2}^{}\phantom{\rule{4em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}\sum _{n=1}^{N}{\Vert {\stackrel{n}{\to}}_{+}^{{G}_{n}^{2}}2,}^{}\end{array}$$

(15)

$$\underset{=}{\Delta}\underset{\{{\overrightarrow{z}}_{n}\}}{\text{argmin}}\sum _{n=1}^{N}\frac{1}{{\sigma}_{n}^{2}}\phantom{\rule{0.2em}{0ex}}{f}_{n}^{2}(\overrightarrow{z})+\frac{1}{{\sigma}_{x}^{2}}\sum _{n=1}^{N}{\Vert {g}_{n}\Vert}^{2}(\overrightarrow{z})$$

(16)

This non-linear least-squares form can be optimized efficiently with the Gauss-Newton method, which requires finding the gradient of both terms with respect to {*$\stackrel{\u20d7}{z}$ _{n}*} at {

We let
${\overrightarrow{m}}_{n}^{T}$ be the 1 × 3 spatial gradient of the warped moving image *M* ^{(t)}(·) at *x _{n}* and note that
${\overrightarrow{m}}_{n}^{T}$ is tangent to the sphere at

$$\begin{array}{l}\frac{{\overrightarrow{z}}_{k}{[F\left({x}_{n}\right)-M\{{(t)}^{\text{exp}\left(\{{E}_{n}{\overrightarrow{z}}_{n}\}\right)}\left({x}_{n}\right)]\overrightarrow{z}=0}_{{=-\frac{{\overrightarrow{z}}_{k}M\{{(t)}^{\text{exp}\left(\{{E}_{n}{\overrightarrow{z}}_{n}\}\right)}\left({x}_{n}\right)|\overrightarrow{z}=0}{}}_{}}}{}\end{array}$$

(17)

$$\begin{array}{l}=-\frac{M\{{(t)}^{\text{exp}\left(\{{E}_{n}{\overrightarrow{z}}_{n}\}\right)}\left({x}_{n}\right)\text{exp}\left(\{{E}_{n}{\overrightarrow{z}}_{n}\}\right)\left({x}_{n}\right)\times \phantom{\rule{9em}{0ex}}{\times \phantom{\rule{0.5em}{0ex}}\left[\frac{\text{exp}\left(\{{E}_{n}{\overrightarrow{z}}_{n}\}\right)\left({x}_{n}\right){\overrightarrow{z}}_{k}}{]}\right|\overrightarrow{z}=0}_{}}{}\end{array}$$

(18)

$$={-\frac{M{(t)}^{\left({u}_{n}\right)}{u}_{n}|}{{u}_{n}={x}_{n}}{[\frac{\text{exp}\left(\{{E}_{n}{\overrightarrow{z}}_{n}\}\right)\left({x}_{n}\right){E}_{k}{\overrightarrow{z}}_{k}\frac{{E}_{k}{\overrightarrow{z}}_{k}{\overrightarrow{z}}_{k}}{]}}{|}\overrightarrow{z}=0}_{}}_{}$$

(19)

$$=-{\overrightarrow{m}}_{n}^{T}{E}_{n}\delta (k,n),$$

(20)

where *δ*(*k, n*) = 1 if *k* = *n* and 0 otherwise. Eq. (20) uses the fact that the differential of exp() at = 0 is the identity [47], i.e, [*D* exp(0)] = . In other words, a change in velocity * _{k}*. at vertex

Similarly, we define
${\mathit{\text{S}}}_{n}^{T}$ to be the 3 × 3 Jacobian of ^{(t)}(·) at *x _{n}*. The computation of
${\mathit{\text{S}}}_{n}^{T}$ is discussed in Appendix A-B. Differentiating the second term of the cost function

$$\begin{array}{l}\frac{{\overrightarrow{z}}_{k}{[{\stackrel{n}{\to}}_{+}^{{G}_{n}^{2}}}_{\phantom{\rule{12em}{0ex}}={G}_{n}^{2}{\mathit{\text{S}}}_{n}^{T}{E}_{n}\delta (k,n),}}{}\end{array}$$

(21)

where *G _{n}* is the skew-symmetric matrix defined in Eq. (4).

Once the derivatives are known, we can compute the corresponding gradients based on our choice of the metric of vector fields on *S*^{2}. In this work, we assume an *l*_{2} inner product, so that the inner product of vector fields is equal to the sum of the inner products of the individual vectors. The inner product of individual vectors is in turn specified by the choice of the Riemannian metric on *S*^{2}. Assuming the canonical metric, so that the inner product of two tangent vectors is the usual inner product in the Euclidean space [39], the gradients are equal to the transpose of the derivatives Eqs. (20), (21) (see Appendix A-C). We can then rewrite Eq. (15) as a linearized least-squares objective function:

$$\begin{array}{l}\{{\overrightarrow{z}}_{n}^{(t)}\}\\ \approx \underset{\{{\overrightarrow{z}}_{n}\}}{\text{argmin}}\sum _{n=1}^{N}\frac{1}{{\sigma}_{n}^{2}}{({f}_{n}\left(\overrightarrow{z}=0\right)+{{l}_{2}}_{{f}_{n}})2}^{}& \phantom{\rule{9em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}\sum _{n=1}^{N}{\Vert {g}_{n}(\overrightarrow{z}=0)+{{l}_{2}}_{{g}_{n}}\Vert 2}^{}\end{array}$$

(22)

$$\begin{array}{l}=\underset{\{{\overrightarrow{z}}_{n}\}}{\text{argmin}}\sum _{n=1}^{N}\frac{1}{{\sigma}_{n}^{2}}{(\left(F({x}_{n})-M{(t)}^{(})-{\overrightarrow{m}}_{n}^{T}{E}_{n}{\overrightarrow{z}}_{n}\right)2}^{}\phantom{\rule{9em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}\sum _{n=1}^{N}{\Vert {G}_{n}^{2}{\mathit{\text{S}}}_{n}^{T}{E}_{n}{\overrightarrow{z}}_{n}\Vert}^{2}\end{array}$$

(23)

$$\begin{array}{l}=\underset{\{{\overrightarrow{z}}_{n}\}}{\text{argmin}}\sum _{n=1}^{N}\Vert (\begin{array}{c}\frac{1}{{\sigma}_{n}}(F({x}_{n})-M{(t)}^{(})0\\ )\end{array}{\phantom{\rule{11em}{0ex}}+\phantom{\rule{0.2em}{0ex}}\left(\begin{array}{c}-\frac{1}{{\sigma}_{n}}{\overrightarrow{m}}_{n}^{T}\\ \frac{1}{{\sigma}_{x}}{G}_{n}^{2}{\mathit{\text{S}}}_{n}^{T}\end{array}\right){E}_{n}{\overrightarrow{z}}_{n}\Vert}^{2}.\end{array}$$

(24)

Because of the delta function *δ*(*k, n*) in the derivatives in Eqs. (20), (21), *$\stackrel{\u20d7}{z}$ _{n}* only appears in the

$$\begin{array}{l}{\overrightarrow{z}}_{n}^{(t)}=\frac{F({x}_{n})-M{(t)}^{(}{\sigma}_{n}^{2}({\mathit{\text{E}}}_{n}^{T}[\frac{1}{{\sigma}_{n}^{2}}{\overrightarrow{m}}_{n}{\overrightarrow{m}}_{n}^{T}+}{}{\phantom{\rule{8em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}{S}_{n}{({G}_{n}^{2})}^{T}{G}_{n}^{2}{\mathit{\text{S}}}_{n}^{T}]\phantom{\rule{0.2em}{0ex}}{E}_{n})}^{-1}{\mathit{\text{E}}}_{n}^{T}{\overrightarrow{m}}_{n}.\end{array}$$

(25)

For each vertex, we only need to perform matrix-vector multiplication of up to 3 × 3 matrices and matrix inversion of 2 × 2 matrices. This implies the update rule for *$\stackrel{\u20d7}{v}$ _{n}*:

$${\overrightarrow{\upsilon}}_{n}^{(t)}={E}_{n}{\overrightarrow{z}}_{n}^{(t)}$$

(26)

$$\begin{array}{l}=\frac{F({x}_{n})-M{(t)}^{(}{\sigma}_{n}^{2}{E}_{n}\phantom{\rule{0.2em}{0ex}}({\mathit{\text{E}}}_{n}^{T}\phantom{\rule{0.2em}{0ex}}[\frac{1}{{\sigma}_{n}^{2}}{\overrightarrow{m}}_{n}{\overrightarrow{m}}_{n}^{T}+}{}{\phantom{\rule{6em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}{S}_{n}{({G}_{n}^{2})}^{T}{G}_{n}^{2}{\mathit{\text{S}}}_{n}^{T}]\phantom{\rule{0.2em}{0ex}}{E}_{n})}^{-1}{\mathit{\text{E}}}_{n}^{T}{\overrightarrow{m}}_{n}.\end{array}$$

(27)

In practice, we use the Levenberg-Marquardt modification of Gauss-Newton optimization [49] to ensure matrix invertibility:

$$\begin{array}{l}{\overrightarrow{\upsilon}}_{n}^{(t)}=\frac{F({x}_{n})-M{(t)}^{(}{\sigma}_{n}^{2}{E}_{n}\phantom{\rule{0.2em}{0ex}}({\mathit{\text{E}}}_{n}^{T}\phantom{\rule{0.2em}{0ex}}[\frac{1}{{\sigma}_{n}^{2}}{\overrightarrow{m}}_{n}{\overrightarrow{m}}_{n}^{T}+}{}{\phantom{\rule{5em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}{S}_{n}{({G}_{n}^{2})}^{T}{G}_{n}^{2}{\mathit{\text{S}}}_{n}^{T}]\phantom{\rule{0.2em}{0ex}}{E}_{n}+\epsilon {I}_{2\times 2})}^{-1}{\mathit{\text{E}}}_{n}^{T}{\overrightarrow{m}}_{n}.\end{array}$$

(28)

where *ε* is a regularization constant. We note that in the classical Euclidean Demons [57], [14], the term
${\mathit{\text{E}}}_{n}^{T}{S}_{n}{\left({G}_{n}^{2}\right)}^{T}{G}_{n}^{2}{\mathit{\text{S}}}_{n}^{T}{E}_{n}$ turns out to be the identity, so it can also be seen as utilizing Levenberg-Marquardt optimization. Once again, we emphasize that a different choice of the coordinate charts will lead to the same update.

Given
${\{{\overrightarrow{\upsilon}}_{n}^{(t)}\}}_{n=1}^{N}$, we use “scaling and squaring” to compute exp(^{(t)}) [3], which is then composed with the current transformation estimate ^{(t)} to form Γ^{(t)} = ^{(t)} exp(^{(t)}). Appendix D discusses implementation details of extending the “scaling and squaring” procedure in Euclidean spaces to *S*^{2}.

We now define the Reg() term using the corresponding tangent vector field representation ϒ⃗. Following the work of [31], [61], we let *H* be the Hilbert space of square integrable vector fields on the sphere defined by the inner product:

$${{\overrightarrow{u}}_{1},{\overrightarrow{u}}_{2}H={\int}_{{S}^{2}}{{\overrightarrow{u}}_{1}(x),{\overrightarrow{u}}_{2}(x)R\phantom{\rule{0.2em}{0ex}}d{S}^{2}}_{},}_{}$$

(29)

where *$\stackrel{\u20d7}{u}$*_{1}, *$\stackrel{\u20d7}{u}$*_{2} *H* and ·,·_{R} refers to the canonical metric. Because vector fields from *H* are not necessarily smooth, we restrict the deformation ϒ⃗ to belong to the Hilbert space *V* *H* of vector fields obtained by the closure of the space of smooth vector fields on *S*^{2} via a choice of the so-called energetic inner product denoted by

$${\overrightarrow{u},\overrightarrow{\upsilon}V={L\overrightarrow{u},\overrightarrow{\upsilon}H,}_{}}_{}$$

(30)

where *L* could for example be the Laplacian operator on smooth vector fields on *S*^{2} [31], [61].

We define Reg() ‖ϒ⃗‖* _{V}*. With a proper choice of the energetic inner product (e.g., Laplacian), a smaller value of ‖ϒ⃗‖

With our choice of dist(, Γ) in Section III-A and Reg() in Section III-C, the optimization in Step 2 of the Spherical Demons algorithm

$${\stackrel{(t+1)}{\to}}^{=}$$

(31)

seeks a smooth vector field ϒ⃗ *V* that approximates the tangent vectors
${\{{\overrightarrow{\Gamma}}_{n}^{(t)}\}}_{n=1}^{N}$. This problem corresponds to the inexact vector spline interpolation problem solved in [31], motivating our use of tangent vectors in the definition of dist(, Γ) in Section III-A, instead of the more intuitive choice of geodesic distance.

We can express the tangent vectors * _{n}* and ϒ⃗

$$\stackrel{=}{^}$$

(32)

where *K* is a 2*N* × 2*N* matrix consisting of *N* × *N* blocks of 2 × 2 matrices: the (*i,j*) block corresponds to *k*(*x _{i},x_{j}*)

In [31], the spherical vector spline interpolation problem was applied to landmark matching on *S*^{2}, resulting in a reasonable sized linear system of equations. Solving the matrix inversion shown in Eq. (32) is unfortunately prohibitive for cortical surfaces with more than 100,000 vertices. If one chooses a relatively wide kernel *k*(*x _{i},x_{j}*), the system is not even sparse.

Inspired by the convolution method of optimizing Step 2 in the Demons algorithm [14], [57], [66] and the convolution-based fast fluid registration in the Euclidean space [12], we propose an iterative smoothing approximation to the solution of the spherical vector spline interpolation problem.

In each smoothing iteration, for each vertex *x _{i}*, tangent vectors of neighboring vertices

We note that the iterative smoothing approximation to spline interpolation is not exact because parallel transport is not transitive on *S*^{2} due to the non-flat curvature of *S*^{2} (unlike in Euclidean space), i.e., parallel transporting a tangent vector from point *a* to *b* to *c* results in a vector different from the result of parallel transporting a tangent vector from *a* to *c*. Furthermore, the approximation accuracy degrades as the distribution of points becomes less uniform. In Appendix B, we provide a theoretical bound on the approximation error and demonstrate empirically that iterative smoothing provides a good approximation of spherical vector spline interpolation for a relatively uniform distribution of points corresponding to those of the subdivided icosahedron meshes used in this work.

The Spherical Demons algorithm is summarized in Algorithm 2.

We run the Spherical Demons algorithm in a multi-scale fashion on a subdivided icosahedral mesh. We begin from a subdivided icosahedral mesh (ic4) that contains 2,562 vertices and work up to a subdivided icosahedral mesh (ic7) that contains 163,842 vertices, which is roughly equal to the number of vertices in the cortical meshes we work with. We perform 15 iterations of Step 1 and Step 2 at each level. Because of the fast convergence rate of the Gauss-Newton method, we find that 15 iterations are more than sufficient for our purposes. We also perform a rotational registration at the beginning of each multi-scale level via a sectioned search of the three Euler angles.

Empirically, we find the computation time of the Spherical Demons algorithm is roughly divided equally among the four components: registration by rotation, computing the Gauss-Newton update, performing “scaling and squaring” and smoothing the vector field.

Data: A fixed spherical image F and moving spherical image M. |

Result: Diffeomorphism Γ so that M Γ is “close” to F. |

Set ^{0} = identity transformation (or some a-priori transformation, e.g., from a previous registration) |

repeat |

Step 1. Given ^{(t)}, |

foreach
vertex n
do |

Compute ${\overrightarrow{\upsilon}}_{n}^{(t)}$ using Eq. (28). |

end |

Compute Γ^{(t)} = exp() using “scaling and squaring”. |

Step 2. Given Γ^{(t)}, |

foreach
vertex n
do |

Compute ${\stackrel{n}{\to}}_{}^{}$ using Eq. (48) implemented via iterative smoothing. |

end |

until
convergence; |

In practice, we work with spheres that are normalized to be of radius 100, because we find that at ic7, the average edge length of 1mm corresponds to that of the original cortical surface meshes. This allows for meaningful interpretation of distances on the sphere. This requires slight modification of the equations presented previously to keep track of the radius of the sphere.

The Spherical Demons algorithm presented here registers pairs of spherical images. To incorporate a probabilistic atlas defined by a mean image and a standard deviation image, we modify the Demons objective function in Eq. (3), as explained in Appendix C. This requires a choice of warping the subject or warping the atlas. We find that interpolating the atlas gives slightly better results, compared with interpolating the subject. However, interpolating the subject results in a runtime of under 3 minutes, while the runtime for interpolating the atlas is less than 5 minutes. In the next section, we report results for interpolating the atlas.

We use two sets of experiments to evaluate the performance of the Spherical Demons algorithm by comparing it to the widely used and freely available FreeSurfer [27] software. The FreeSurfer registration algorithm uses the same similarity measure as Demons, but explicitly penalizes for metric and areal distortion. As we will show, even though the Spherical Demons algorithm does not specifically take into account the original metric properties of the cortical surface, we still achieve comparable if not better registration accuracy than FreeSurfer. Furthermore, FreeSurfer runtime is more than an hour while Spherical Demons runtime is less than 5 minutes.

There are four parameters in the algorithm.
$1/{\sigma}_{x}^{2}$ and *ε* appear in Eq. (28). Larger values of
$1/{\sigma}_{x}^{2}$ and *ε* decrease the size of the update taken in Step 1 of the Spherical Demons algorithm. In the experiments that follow, we set
$1/{\sigma}_{x}^{2}=\epsilon $ and set their values such that the largest vector of the update velocity field is roughly two times the edge lengths of the mesh. The number of iterations *m* and the weight
$\text{exp}(-\frac{1}{2\mathit{\gamma}})$ determine the degree of smoothing. We set *γ* = 1 and explore a range of smoothing iterations *m* in the following experiments.

We validate Spherical Demons in the context of automatic cortical parcellation. Automatic labeling of cortical brain surfaces is important for identifying regions of interests for clinical, functional and structural studies [20], [52]. Recent efforts have ranged from the identification of sulcal/gyral ridge lines [56], [62] to the segmentation of sulcal/gyral basins [20], [28], [38], [41], [42], [51], [52], [67]. Similar to these prior studies, we are interested in parcellation of the entire cortical surface meshes, where each vertex is assigned a label.

We consider a set of 39 left and right cortical surface models extracted from in-vivo MRI [19]. Each surface is spherically parameterized and represented as a spherical image with geometric features at each vertex: mean curvature of the cortical surfaces, mean curvature of the inflated cortical surfaces and average convexity of the cortical surfaces, which roughly corresponds to sulcal depth [26]. These features are intrinsic and thus independent of the parameterization of the surface. The tools used for segmentation [19] and spherical parameterization [26] are freely available [29]. Both hemispheres of each subject were manually parcellated by a neuroanatomist into 35 labels, corresponding to the main sulci and gyri, enumerated in Table II.

We co-register all 39 spherical images of cortical geometry with Spherical Demons by iteratively building an atlas and registering the surfaces to the atlas. The atlas consists of the mean and variance of cortical geometry represented by the surface features described above. We then perform 4-fold cross-validation of the parcellation of the co-registered cortical surfaces. In each iteration of cross-validation, we leave out ten subjects and use the remainder of the subjects to train a classifier [20], [28] that predicts the labels based on location and geometric features. We then apply the classifier to the hold-out set of ten subjects. We perform each iteration with a different hold-out set, i.e., subjects 1-10, 11-20, 21-30 and 31-39.

As mentioned previously, increasing the number of iterations of smoothing results in smoother warps. As discussed in [67], the choice of the tradeoff between the similarity measure and regularization is important for segmentation accuracy. Estimating the optimal registration regularization tradeoff is an active area of research [1], [18], [48], [65], [67], [68] that we do not deal with in this paper. Here, we simply repeat the above experiments using {6, 8, 10, 12, 14} iterations of smoothing. For brevity, we will focus the discussion on using 10 iterations of smoothing and comment on results obtained with the other levels of smoothing.

We repeat the above procedure of performing co-registration and cross-validation with the FreeSurfer registration algorithm [27] using the default FreeSurfer settings. Once again, we use the same features and parcellation algorithm [20], [28]. As before, the atlas consists of the mean and variance of cortical geometry.

To compare the cortical parcellation results, we compute the average Dice measure, defined as the ratio of cortical surface area with correct labels to the total surface area averaged over the test set. Because the average Dice can be misleading by suppressing small structures, we also compute the Dice measure for each structure.

On the left hemisphere, FreeSurfer achieves an average Dice of 88.9, while Spherical Demons achieves an average Dice of 89.6 with 10 iterations of smoothing. While the improvement is not big, the difference is statistically significant for a onesided t-test with the Dice measure of each subject treated as an independent sample (*p* = 2 × 10^{−6}). Furthermore, the overall Dice is statistically significantly better than FreeSurfer for all levels of smoothing we considered, with the best overal dice achieved with 12 iterations of smoothing.

On the right hemisphere, FreeSurfer obtains a Dice of 88.8 and Spherical Demons achieves 89.1 with 10 iterations of smoothing. Here, the improvement is smaller, but still statistically significant (*p* = 0.01). Furthermore, the overall dice is statistically significantly better than FreeSurfer for all levels of smoothing we considered, except when 6 iterations of smoothing is used (*p* = 0.06). All results we report in the remainder of this section use 10 iterations of smoothing.

We analyze the segmentation accuracy separately for each structure. To compare Spherical Demons with FreeSurfer, we perform a one-sided paired-sampled t-test treating each subject as an independent sample and correct for multiple comparisons using a False Discovery Rate (FDR) of 0.05 [10]. On the left (right) hemisphere, the segmentations of 16 (8) structures are statistically significantly improved by Spherical Demons with respect to FreeSurfer, while no structure is significantly worse.

Fig. 3 shows the percentage improvement of individual structures over FreeSurfer. Fig. 4 displays the average Dice per structure for FreeSurfer and Spherical Demons (10 iterations of smoothing) for the left and right hemispheres. Standard errors of the mean are displayed as red bars. The numbering of the structures correspond to Table II. Parcellation improvements suggest that our registration is at least as accurate as FreeSurfer.

Percentage Improvement over FreeSurfer. Yellow regions indicate structures scoring better than FreeSurfer. Blue regions correspond to decrease in accuracy. Note that none of these blue regions are statistically significant. The boundaries between parcellation **...**

(a) Dice measure for each structure in the left hemisphere. (b) Dice measure for each structure in the right hemisphere. Black columns correspond to FreeSurfer. White columns correspond to Spherical Demons. * indicates structures where Spherical Demons **...**

The structures with the worst Dice are the frontal pole and entorhinal cortex. These structures are small and relatively poorly defined by the underlying cortical geometry. For example, the entorhinal cortex is partially defined by the rhinal sulcus, a tiny sulcus that is only visible on the pial surface. The frontal pole is defined by the surrounding structures, rather than by the underlying cortical geometry.

Brodmann areas are cyto-architectonically defined parcellations of the cerebral cortex [13]. They can be observed through histology and more recently, through ex-vivo high resolution MRI [6]. Unfortunately, much of the cytoarchitectonics cannot be observed with current in-vivo imaging. Nevertheless, most studies today report their functional findings with respect to Brodmann areas, usually estimated by visual comparison of cortical folds with Brodmann's original drawings without quantitative analysis of local accuracy. By combining histology and MRI, recent methods for creating probabilistic Brodmann area maps in the Talairach and Colin27 normalized space promise a more principled approach [2], [24], [54], [55], [71].

In this experiment, we consider a data set that contains Brodmann labels mapped to the corresponding MRI volume. Specifically, we work with postmortem histological images of ten brains created using the techniques described in [54], [71]. The histological sections were aligned to postmortem MR with nonlinear warps to build a 3D histological volume. These volumes were segmented to separate white matter from other tissue classes, and the segmentation was used to generate topologically correct and geometrically accurate surface representations of the cerebral cortex using FreeSurfer [19]. The eight manually labeled Brodmann area maps (areas 2, 4a, 4p, 6, 44, 45, 17 and 18) were sampled onto the surface representations of each hemisphere, and errors in this sampling were manually corrected (e.g., when a label was erroneously assigned to both banks of a sulcus). A morphological close was then performed on each label to remove small holes. We note that Brodmann areas 4a, 4p and 6 were mapped in only eight of the ten subjects. Fig. 5 shows these eight Brodmann areas on the resulting cortical representations for two subjects. Finally, we map the folding patterns and the Brodmann area labels onto a spherical coordinate system [27].

Brodmann areas 17 (V1), 18 (V2), 2, 4a, 4p, 6, 44 and 45 shown on inflated cortical surfaces of two subjects. Notice the variability of BA44 and BA45 with respect to the underlying folding pattern.

It has been shown that nonlinear surface registration of cortical folds can significantly improve Brodmann area overlap across different subjects [25], [68] compared with volumetric registration. Registering the ex-vivo surfaces is more difficult than in-vivo surfaces because the reconstructed volumes are extremely noisy due to the distortions introduced by the histology, resulting in noisy geometric features, as shown in Fig. 6.

Left: example in-vivo surface used in the parcellation study. Right: example ex-vivo surface used in the Brodmann area study.

We consider two strategies for aligning Brodmann areas. For both strategies, we will use 10 iterations of smoothing for Spherical Demons as it proved reasonable in the previous set of experiments. The first strategy involves co-registering the 10 ex-vivo surfaces using cortical geometry by repeatedly building an atlas and registering the surfaces to the atlas, similar to the previous experiment on cortical parcellation. We use either Spherical Demons or FreeSurfer for registration. We refer to the co-registration using Spherical Demons and FreeSurfer as SD10 and FS10 respectively (10 refers to the number of subjects in the study, not the number of smoothing iterations).

The second strategy involves registering the 10 ex-vivo surfaces to the in-vivo “Buckner40” atlas, constructed from 40 in-vivo subjects, that is distributed with the FreeSurfer software. Once again, we use either Spherical Demons or FreeSurfer for the registration. We refer to the co-registration using Spherical Demons and FreeSurfer as SD40 and FS40 respectively.

To measure the quality of alignment of the Brodmann areas after cortical registration, we use an adaptation of the modified Hausdorff distance [21]. For each pair of registered subjects, we project each Brodmann area from the first subject onto the second subject and compute the mean distance between the boundaries, measured on the original cortical surface of the second subject. We obtain a second measurement by projecting each Brodmann area from the second subject onto the first subject. Since we have 10 surfaces, we get 90 ordered pairs and 90 alignment errors for each Brodmann area.

Table III reports the mean alignment errors for each Brodmann area and for each method. The lowest errors for each Brodmann area are shown in **bold**. We see that for almost all Brodmann areas, the best alignment come from SD10 or SD40. Similarly, Fig. 7 shows the median alignment error for each Brodmann area. The error bars indicate the lower and upper quartile alignment errors.

Median alignment errors of Brodmann areas in *mm* for the four registration methods. The error bars indicate the upper and lower quartile alignment errors. “#” indicates that the median errors of SD10 are statistically lower than those of **...**

Mean alignment errors of Brodmann areas in *mm* for the four registration methods. Lowest errors are shown in **bold**. SD refers to Spherical Demons; FS refers to FreeSurfer.

We use permutation testing to evaluate statistical significance of the results. We cannot use the t-test because the 90 alignment errors are correlated - since the subjects are co-registered together, good alignment between subjects 1 and 2 and between subjects 2 and 3 necessarily implies a higher likelihood of good alignment between subjects 1 and 3.

The tests show that SD10 is better than FS10 and SD40 is slightly better than FS40. SD10 and SD40 are comparable. Compared with FS10, SD10 improves the median alignment errors of 5 (4) Brodmann areas on the right (left) hemisphere (FDR = 0.05) and no structure gets worse. Compared with FS40, SD40 statistically improves the alignment of 2 (1) Brodmann areas on the right (left) hemisphere (FDR = 0.05) with no structure getting worse. Permutation tests on the mean alignment errors show similar results, except that FS40 performs better than SD40 for BA4p on the left hemisphere when using the mean statistic. These results suggest that the Spherical Demons algorithm is at least as accurate as FreeSurfer in aligning cortical folds and Brodmann areas.

The Demons algorithms [57], [66] discussed in Section II and the Spherical Demons algorithm proposed in this paper use a regularization term that modulates the final deformation. Motivated by [12], [14], the Diffeomorphic Demons algorithm [66] admits a fluid prior on the velocity fields. This involves smoothing the velocity field updates *before* computing the exponential map to obtain the displacement field updates to be composed with the current transformation. The resulting algorithm is very similar to the fast implementation [12] of Christensen's well-known fluid registration algorithm [16], except that Christensen's algorithm does not employ a higher-order update method like Gauss-Newon. The Spherical Demons algorithm can similarly incorporate a fluid prior by smoothing the velocity field ^{(t)} in Eq. (28) before computing the exponential map to obtain the displacement updates exp(^{(t)}).

An alternative interpretation of the smoothing implementation of Christensen's algorithm comes from choosing a different metric for computing the gradient from the derivatives [9]. The choice of the metric also arises in our problem when computing the gradient as discussed in Appendix A-C. This suggests that the Spherical Demons algorithm can incorporate a fluid prior by modifying the Gauss-Newton update step Eq. (28). Unfortunately, this process introduces coupling among the vertices resulting in the loss of the speed-up advantage of Spherical Demons (see for example the derivations of [34]). The exploration of the performance of the different fluid prior implementations is outside the scope of this paper.

Because the tools of differential geometry are general, the Spherical Demons algorithm can be in principle extended to arbitrary manifolds, besides the sphere. One challenge is that the definition of coordinate charts for an arbitrary manifold is more difficult than that for the sphere. Approaches of utilizing the embedding space [15] or the intrinsic properties of manifolds [40] are promising avenues of future work.

In this paper, we presented the fast Spherical Demons algorithm for registering spherical images. We showed that the two-step optimization of the Demons algorithm can also be applied on the sphere. By utilizing the one parameter subgroups of diffeomorphims, the resulting deformation is invertible. We tested the algorithm extensively in two different applications and showed that the accuracy of the algorithm compares favorably with the widely used FreeSurfer registration algorithm [27] while offering more than one order of magnitude speedup. Both matlab and ITK versions of the Spherical Demons algorithm are publicly available^{2}.

A clear future challenge is to take into account the original metric properties of the cortical surface in the registration process, as demonstrated in previously proposed registration methods [27], [59].

We note that while fast algorithms are useful for deploying the developed tool on large datasets, they can further allow for complex applications that were previously computationally intractable. For example, we have incorporated the ideas behind Spherical Demons into a meta-registration framework that learns registration cost functions which are optimal for specific applications [68].

We thank Hartmut Mohlberg, Katrin Amunts and Karl Zilles for providing the histological dataset. We also like to thank Xavier Pennec, Wanmei Ou and Serdar Balci for helpful discussions and Tamar Riklin Raviv for reading the conference draft of this paper. We also thank Luis Ibanez and Michel Audette for reading the journal draft of this paper and for their help in implementing the ITK version of the Spherical Demons algorithm. Finally, we thank the reviewers for their many helpful suggestions on improving and clarifying the paper.

Support for this research is provided in part by the NAMIC (NIH NIBIB NAMIC U54-EB005149), the NAC (NIH NCRR NAC P41-RR13218), the mBIRN (NIH NCRR mBIRN U24-RR021382), the NIH NINDS R01-NS051826 grant, the NSF CAREER 0642971 grant, the National Institute on Aging (AG02238), NCRR (P41-RR14075, R01 RR16594-01A1), the NIBIB (R01 EB001550, R01EB006758), the NINDS (R01 NS052585-01) and the MIND Institute. Additional support was provided by The Autism & Dyslexia Project funded by the Ellison Medical Foundation. B.T. Thomas Yeo is funded by the A*STAR, Singapore. Short visits of B.T. Thomas Yeo and Nicholas Ayache between MIT and INRIA were partially funded by the CompuTumor associated teams funding.

In this appendix, we provide details on the computation of the spatial gradient of the warped moving image *M* ^{(t)} and the Jacobian of the deformation ^{(t)}. We also compute the gradients of the demons cost function using the derivatives computed in Eq. (20) and Eq. (21), assuming the *l*_{2} inner product space for vector fields and the canonical metric.

In this appendix, we discuss the computation of
${\mathit{\text{m}}}_{n}^{T}$, the spatial gradient of the warped moving image *M* ^{(t)} at the point *x _{n}*. We can think of

$${M}_{s}(x)=I(p(x)),$$

(33)

where *p*(*x*) is the intersection point and *I*(·) is the barycentric interpolation. Let *p*_{1}, *p*_{2}, *p*_{3} denote the vertices of the triangle containing *p*(*x*) and *$\stackrel{\u20d7}{n}$* denote the 3 × 1 normal vector to the triangle. Since *p*(*x*) = *αx* for some *α* and *p*(*x*) − *p*_{1}, *$\stackrel{\u20d7}{n}$* = 0, we can write

$$p(x)=\frac{{p}_{1},\overrightarrow{n}x,\overrightarrow{n}x}{}$$

(34)

and

$$I(p)=\frac{{A}_{1}(p)\phantom{\rule{0.1em}{0ex}}{M}_{s}({p}_{1})+{A}_{2}(p)\phantom{\rule{0.1em}{0ex}}{M}_{s}({p}_{2})+{A}_{3}(p)\phantom{\rule{0.1em}{0ex}}{M}_{s}({p}_{3})}{A},$$

(35)

where *A*_{1}(*p*), *A*_{2}(*p*) and *A*_{3}(*p*) are the areas of the triangles *pp*_{2}*p*_{3}, *pp*_{1}*p*_{3} and *pp*_{1}*p*_{2} respectively. Note that *A* = *A*_{1}(*p*) + *A*_{2}(*p*) + *A*_{3}(*p*). *M _{s}*(

Computing the derivative of the image value at *x* follows easily from the chain rule:

$$\frac{p(x)x=\frac{{p}_{1},\overrightarrow{n}x,\overrightarrow{n}{I}_{3\times 3}-\frac{{p}_{1},\overrightarrow{n}{x,\overrightarrow{n}2}^{}x{\overrightarrow{n}}^{T}}{}}{}}{}$$

(36)

$$\begin{array}{l}\frac{I(p)p=\frac{{p}_{{A}_{1}}}{}}{}\end{array}$$

(37)

where * _{p}A_{i}*(

A complication arises when *x* corresponds to one of the mesh vertices, since the spatial gradient is not defined in this case. The same problem arises in Euclidean space with linear interpolation and the spatial gradient is typically defined via finite central difference. It is unclear what the equivalent definition on a mesh is. Here, for a mesh vertex *x*, we compute the spatial gradient for each of the surrounding triangles and linearly combine the spatial gradients using weights corresponding to the areas of the triangles.

In this appendix, we discuss the computation of
${\mathit{\text{S}}}_{n}^{T}$, the Jacobian of the deformation ^{(t)} at *x _{n}*. We can think of

$${(t)}^{(}$$

(38)

where *p*(*x*) is the same as in the previous section and

$$I(p)=\frac{{A}_{1}(p){(t)}^{(}}{}$$

(39)

The Jacobian is computed via chain rule, just like in the previous section.

In this appendix, we seek to compute the gradients of *f _{n}*(

$${{\overrightarrow{z}}^{1},{\overrightarrow{z}}^{2}{l}_{2}={\{{E}_{k}\phantom{\rule{0.1em}{0ex}}{\overrightarrow{z}}_{k}^{1}\},\{{E}_{k}\phantom{\rule{0.1em}{0ex}}{\overrightarrow{z}}_{k}^{2}\}{l}_{2}}_{}}_{}$$

(40)

$$=\sum _{k=1}^{N}{{E}_{k}\phantom{\rule{0.1em}{0ex}}{\overrightarrow{z}}_{k}^{1},{E}_{k}\phantom{\rule{0.1em}{0ex}}{\overrightarrow{z}}_{k}^{2}R}_{}$$

(41)

$$=\sum _{k=1}^{N}{{\overrightarrow{z}}_{k}^{1},{\overrightarrow{z}}_{k}^{2}R}_{,}$$

(42)

where

- Eq. (40) follows from the equivalence of the tangent bundles
*T*^{2}and*TS*^{2}induced by the coordinate charts{Ψ}._{n} - Eq. (41) is the result of the
*l*_{2}assumption that the inner product of vector fields is given by the sum of the inner products of individual vectors. - Because we assume the canonical metric, each term in the inner product in Eq. (41) is simply the usual inner product between 3 × 1 vectors ${E}_{k}\phantom{\rule{0.1em}{0ex}}{\overrightarrow{z}}_{k}^{1}$ and ${E}_{k}\phantom{\rule{0.1em}{0ex}}{\overrightarrow{z}}_{k}^{2}$. Since the columns of
*E*are orthonormal with respect to the usual inner product and using linearity of the inner product, Eq. (41) implies Eq. (42), i.e., the inner product_{k}*$\stackrel{\u20d7}{z}$*^{1},*$\stackrel{\u20d7}{z}$*^{2}_{l2}can be computed by the sum of the usual inner product between 2 × 1 tangent vectors ${\overrightarrow{z}}_{k}^{1}$ and ${\overrightarrow{z}}_{k}^{2}$.

Let *df _{n}*(

$${\mathit{\text{df}}}_{n}(\overrightarrow{z})=-{\overrightarrow{m}}_{n}^{T}{E}_{n}{\overrightarrow{z}}_{n}.$$

(43)

Recall that the gradient _{l2}*f _{n}* of

$$-{\overrightarrow{m}}_{n}^{T}{E}_{n}{\overrightarrow{z}}_{n}={\mathit{\text{df}}}_{n}(\overrightarrow{z})$$

(44)

$$={{{l}_{2}}_{\phantom{\rule{0.1em}{0ex}}}{l}_{2}}_{}$$

(45)

$$=\sum _{k=1}^{N}{{{l}_{2}}_{\phantom{\rule{0.1em}{0ex}}}R.}_{}$$

(46)

Therefore, the gradient _{l2}*f _{n}* can be written as a 2

Similarly, we denote the gradient of *g _{n}*(

In this appendix, we demonstrate empirically that iterative smoothing provides a good approximation of spherical vector spline interpolation for a relatively uniform distribution of points corresponding to those of the subdivided icosahedron meshes used in this work. Once again, we work with spheres that are normalized to be of radius 100.

Recall that we seek {ϒ⃗* _{n}*} = {

$$\stackrel{=}{^}$$

(47)

where *K* is a 2*N* × 2*N* matrix consisting of *N* × *N* blocks of 2 × 2 matrices: the (*i, j*) block corresponds to *k*(*x _{i}, x_{j}*)

In constrast, the iterative smoothing approximation we propose can be formalized as follows:

$$\stackrel{=}{^}$$

(48)

where *m* is a positive integer and *K′* is a 2*N* × 2*N* matrix consisting of *N* × *N* blocks of 2 × 2 matrices: the (*i, j*) block corresponds to *λ*(*x _{i}, x_{j}*)

We now demonstrate empirically that for a range of values of *γ*, iterations *m* and the relatively uniform distribution of points corresponding to those of the subdivided icosahedron mesh, there exist kernels *k*(*x _{i}, x_{j}*) that are well approximated by iterative smoothing. Technically, the resulting

More specifically, given a configuration of mesh points, iterations *m* and value of *γ*, we seek (*x _{i}, x_{j}*), such that
$K\phantom{\rule{0.2em}{0ex}}{\left(\frac{{\sigma}_{x}^{T}}{{\sigma}_{T}^{2}}{I}_{2N\times 2N}+K\right)}^{-1}$ is “close” to (

- In the first stage, we seek
*k**(*x*) that is_{i}, x_{j}*not*constrained to be a function of the distance between*x*and_{i}*x*, such that_{j}$$K\phantom{\rule{0.2em}{0ex}}{\left(\frac{{\sigma}_{x}^{T}}{{\sigma}_{T}^{2}}{I}_{2N\times 2N}+K\right)}^{-1}-{({K}^{\prime})}^{m}\approx 0$$(49)Rearranging the terms, we get$${\left({I}_{2N\times 2N}-{\left({K}^{\prime}\right)}^{m}\right)}^{-1}{\left({K}^{\prime}\right)}^{m}\frac{{\sigma}_{x}^{2}}{{\sigma}_{T}^{2}}\approx K$$(50)To make the “≈” concrete, we optimize for$${k}^{}$$(51)where ‖ · ‖is the Frobenius norm._{F}The cost function Eq. (51) can be optimized componentwise, i.e., we can solve for*k**(*x*) for each pair_{i}, x_{j}*x*. For_{i}, x_{j}*γ*= 1,*m*= 10 and a subdivided icosahedron mesh with 642 vertices, we plot the resulting*k**(*x*) as a function of the geodesic distance between_{i}, x_{j}*x*and_{i}*x*in Fig. 8._{j} - In the second stage, we perform a least-squares fit of a b-spline function to the estimated
*k**(*x*) to obtain the final estimate of (_{i}, x_{j}*x*). Fig. 8 illustrates an example kernel (_{i}, x_{j}*x*) we obtain (c.f., the kernel illustrated in [31]). We note that an alternative to b-spline interpolation is to fit the coefficients of the general kernel function suggested in Appendix A of [31]. This will guarantee that the estimated kernel corresponds to an energetic norm. We leave exploring this direction to future work._{i}, x_{j}

We now investigate the quality of the estimate (*x _{i}, x_{j}*) by computing:

$${\Vert K\phantom{\rule{0.2em}{0ex}}{\left(\frac{{\sigma}_{x}^{T}}{{\sigma}_{T}^{2}}{I}_{2N\times 2N}+K\right)}^{-1}-{({K}^{\prime})}^{m}\Vert}_{2}^{2}$$

(52)

where ‖ · ‖_{2} is the *l*_{2} matrix operator norm. The difference metric Eq. (52) measures the *maximum l*_{2} difference between smoothed vector fields obtained from iterative smoothing and spherical vector spline interpolation for *any* possible input vector field {* _{n}*} of unit

Fig. 9 displays the difference metric we obtained with different values of *γ* and iterations *m* for meshes ic2, ic3, ic4 and ic5. Each of the meshes is obtained from recursively subdividing a lower resolution mesh: ic2 indicates that the mesh was obtained from subdividing an icosahedron mesh twice. The number of vertices quadruples with each subdivision, so that ic5 corresponds to 10,242 vertices.

We conclude from the figure that the differences between the two smoothing methods are relatively small and increase with mesh resolution. As discussed in the next section, we run Spherical Demons on different mesh resolutions, including ic7. Unfortunately, because of the large non-sparse matrices we are dealing with, we were only able to compute the differences up to ic5. Computing the difference metric for ic5 took an entire week on a machine with 128GB of RAM. However, the plots in Fig. 9 indicate that the differences appear to have converged by ic5.

To better understand the incurred differences, Fig. 10 illustrates the outputs and differences of the two smoothing methods for different inputs on ic4. The first row illustrates an input vector field which is zero everywhere except for a single tangent vector of unit norm. The results of spline interpolation and iterative smoothing correspond to our intuition that smoothing a single tangent vector propagates tangent vectors of smaller magnitudes to the surronding areas. The two methods also produce almost identical results as shown by the clean difference image in the fourth column.

Comparison of spline interpolation and iterative smoothing (*m* = 10, *γ* = 1). (a) Input vector field (b) Spline Interpolation Output (c) Iterative Smoothing Output (d) Difference between the second and third columns (e) *l*_{2} norm of the difference. **...**

The second row of Fig. 10 demonstrates the worst unit norm input vector field as measured by the difference metric Eq. (52). This worst unit norm input vector field corresponds to the largest eigenvector in Eq. (52). The pattern of large differences correspond to the original 12 vertices of the uniform icosahedron mesh. These original 12 vertices are the only vertices in the subdivided icosahedron meshes with five, instead of six neighbors, as shown by the pentagon pattern. The fact that these 12 vertices are local maxima of differences suggest that these vertices are treated differently by the two smoothing techniques.

The last row of Fig. 10 demonstrates an input vector field that represents the deformation of an actual registration performed in Section IV. The norm of the input vector field is 700 times that in the first two rows, but the discrepancies between spline interpolation and iterative smoothing are less than expected. The differences of 90% of the vectors are less than 0.2mm, with larger differences in the neighborhoods of the 12 vertices identified previously. Since we conclude previously that the difference metric appears to have converged after ic4, the discrepancies are likely to be acceptable at ic7, whose mesh resultion is 1mm.

We should emphasize that the discrepancies between spline interpolation and iterative smoothing do not necessarily imply registration errors. The differences only indicate the deviations of the deformations from true local optima of the Demons registration cost function Eq. (3) assuming the estimated kernel. Approximating smoothing kernels by iterative smoothing is an active area of research in medical imaging [17], [32]. Future work would involve understanding the interaction between the number of smoothing iterations *m* and the choice of the weights
$\text{exp}(-\frac{1}{2\mathit{\gamma}})$ on the quality of the spherical registration.

In this section, we demonstrate how an atlas consisting of a mean image and a standard deviation image can be incorporated into the Spherical Demons algorithm. The standard deviation image replaces Σ in Eq. (3). We first discuss a probabilistic interpretation of the Demons objective function and its relationship to atlases. We then discuss the optimization of the resulting probabilistic objective function.

The Demons objective function reviewed in Section II is defined for the pairwise registration of images. To incorporate a probabilistic atlas, we now reformulate the objective function. Consider the following Maximum-A-Posteriori objective function:

$$\begin{array}{c}({,{\Gamma}^{})}^{}=\underset{,\Gamma}{\text{argmax}}\hfill \hfill \end{array}$$

(53)

$$=\underset{,\Gamma}{\text{argmax}}$$

(54)

$$=\underset{,\Gamma}{\text{argmax}}$$

(55)

Assuming a Gaussian noise model, we define

$$\begin{array}{c}p(F,M\Gamma )\hfill & \phantom{\rule{1em}{0ex}}=\underset{\frac{1}{\sqrt{2\pi}(\sigma (\Gamma ,{x}_{n}))}}{\overset{\text{exp}}{n=1N}}\hfill \end{array}$$

(56)

$$\text{log}\phantom{\rule{0.1em}{0ex}}p(\Gamma )=\frac{1}{Z({\sigma}_{x}^{2})}\text{exp}\phantom{\rule{0.2em}{0ex}}(-\frac{1}{{\sigma}_{x}^{2}}\sum _{n=1}^{N}{\Vert {\stackrel{n}{\to}}_{-}\Vert 2}^{})\phantom{\rule{0.5em}{0ex}},$$

(57)

$$p()=\frac{1}{Z\left({\sigma}_{T}^{2}\right)}\text{exp}\phantom{\rule{0.2em}{0ex}}(-\frac{1}{{\sigma}_{T}^{2}}\text{Reg}())\phantom{\rule{0.5em}{0ex}},$$

(58)

where Reg() is defined via the energetic norm as discussed in Section III-C and for reasons that will soon be clear, we are being purposefully agnostic about the form of *σ*(Γ, *x _{n}*). The objective function in Eq. (55) becomes

$$\begin{array}{c}({,{\Gamma}^{})}^{=}\hfill \end{array}$$

(59)

which is the instantiation of the Demons objective function Eq. (3), except for the extra term
${\sum}_{n=1}^{N}\phantom{\rule{0.1em}{0ex}}\text{log}\phantom{\rule{0.1em}{0ex}}\sigma \phantom{\rule{0.1em}{0ex}}(\Gamma ,{x}_{n})$. Note that we have omitted the partition functions
$Z\left({\sigma}_{x}^{2}\right)$ and
$Z\left({\sigma}_{T}^{2}\right)$ because *σ _{x}* and

As before, *σ*(Γ, *x _{n}*) is the standard deviation of the intensity at vertex

However, recent work [1], [53] suggests that treating the atlas as a moving image might be more correct theoretically. This involves setting the moving image *M* to be the mean image. In this case, *σ*(Γ, *x _{n}*) =

We now discuss the optimization in Eq. (59). Note that the introduction of the new term
${\sum}_{n=1}^{N}\phantom{\rule{0.1em}{0ex}}\text{log}\phantom{\rule{0.1em}{0ex}}\sigma \phantom{\rule{0.1em}{0ex}}(\Gamma ,{x}_{n})$ only affects Step 1 of the Spherical Demons algorithm. By parameterizing Γ^{(t)} = ()^{(t)} exp({E_{n}*$\stackrel{\u20d7}{z}$ _{n}*}), we get

$$\begin{array}{c}\left\{{\overrightarrow{z}}_{n}^{(t)}\right\}\hfill \\ =\underset{{\overrightarrow{z}}_{n}}{\text{argmin}}\sum _{n=1}^{N}\frac{{(F\left({x}_{n}\right)-M\left\{{(t)}^{\text{exp}\left(\left\{{E}_{n}{\overrightarrow{z}}_{n}\right\}\right)}\left({x}_{n}\right)\right)2{(\sqrt{2}\sigma \left\{{(t)}^{\text{exp}\left(\left\{{E}_{n}{\overrightarrow{z}}_{n}\right\}\right)}\left({x}_{n}\right)\right)2}^{}}^{\phantom{\rule{3em}{0ex}}+\frac{1}{{\sigma}_{x}^{2}}\sum _{n=1}^{N}{\Vert {\stackrel{n}{\to}}_{+}^{{G}_{n}^{2}}2}^{}\phantom{\rule{3em}{0ex}}+\sum _{n=1}^{N}\phantom{\rule{0.1em}{0ex}}\text{log}\phantom{\rule{0.1em}{0ex}}\sigma \{{(t)}^{\text{exp}\left(\left\{{E}_{n}{\overrightarrow{z}}_{n}\right\}\right)}\left({x}_{n}\right)\hfill \hfill}}{}\hfill \end{array}$$

(60)

$$\underset{{\overrightarrow{z}}_{n}}{\text{argmin}}\sum _{n=1}^{N}{f}_{n1}^{2}(\overrightarrow{z})+{f}_{n2}^{2}(\overrightarrow{z})+\phantom{\rule{0.1em}{0ex}}\text{log}\phantom{\rule{0.2em}{0ex}}{f}_{n3}(\overrightarrow{z}).$$

(61)

The second term is the same as before, while the first term has become more complicated. Using the product rule and the techniques described in Appendix A, we can find the first derivatives of the first and second terms and estimate their second derivatives using the Gauss-Newton method. The difficulty lies in the third term, which is not quadratic and is even strictly concave, so we have to make further approximations.

Consider the problem of optimizing a one-dimensional function *f*(*x*). Let the current estimate of *x* be *x*_{0}. Newton's optimization [49] involves the following update:

$$\Delta x=-{({f}^{\u2033})}^{-1}({x}_{0})\phantom{\rule{0.2em}{0ex}}{f}^{\prime}({x}_{0}),$$

(62)

where *f*′(*x*_{0}) and *f*″(*x*_{0}) are the gradient and the Hessian of *f* evaluated at *x*_{0} respectively. When *f*″ is negative (positive), the update *x* increases (decreases) the objective function, regardless of whether one is attempting to increase or decrease the objective function! The Gauss-Newton approximation of the Hessian for minimizing non-linear quadratic functions actually stabilizes the Newton's method by ensuring the estimated Hessian is positive.

To optimize Eq. (61) with Newton's method, we need to compute the gradient and the Hessian. Because we are using the *l*_{2} inner product and the canonical metric (see Appendix A-C), the gradient and the Hessian correspond to the first and second derivatives. The first derivative or gradient corresponds to

$$\frac{f{\overrightarrow{z}}_{k}=2{f}_{n1}(\overrightarrow{z})\frac{{f}_{n1}{\overrightarrow{z}}_{k}+2{f}_{n2}(\overrightarrow{z})\frac{{f}_{n2}{\overrightarrow{z}}_{k}+\frac{1}{{f}_{n3}}\frac{{f}_{n3}{\overrightarrow{z}}_{k}}{}}{}}{}}{}$$

(63)

and the second derivative corresponds to

$$\begin{array}{c}\frac{{2}^{f}{{\overrightarrow{\upsilon}}^{\prime}}_{k}^{2}=2{(\frac{{f}_{n1}{\overrightarrow{z}}_{k}}{)}2+2{f}_{n1}(\overrightarrow{z})\frac{{2}^{{f}_{n1}}{{\overrightarrow{\upsilon}}^{\prime}}_{k}^{2}+2{(\frac{{f}_{n2}{\overrightarrow{z}}_{k}}{)}2+}^{}\phantom{\rule{4.5em}{0ex}}+2{f}_{n2}(\overrightarrow{z})\frac{{2}^{{f}_{n2}}{{\overrightarrow{\upsilon}}^{\prime}}_{k}^{2}-{(\frac{{f}_{n3}{\overrightarrow{z}}_{k}}{)}2+\frac{1}{{f}_{n3}}\frac{{2}^{{f}_{n3}}{{\overrightarrow{\upsilon}}^{\prime}}_{k}^{2}}{}}^{}}{}\hfill}{}}^{}}{}\hfill \end{array}$$

(64)

$$\approx 2\phantom{\rule{0.2em}{0ex}}{(\frac{{f}_{n1}{\overrightarrow{z}}_{k}}{)}2+2\phantom{\rule{0.2em}{0ex}}{(\frac{{f}_{n2}{\overrightarrow{z}}_{k}}{)}2-{(\frac{{f}_{n3}{\overrightarrow{z}}_{k}}{)}2\phantom{\rule{0.2em}{0ex}}.}^{}}^{}}^{}$$

(65)

where the last approximation was made using the Gauss-Newton method. Not surprisingly, the third term corresponding to log is negative, which can introduce instability in the Gauss-Newton update. Consequently, we drop the last term, resulting in:

$$\frac{{2}^{f}{{\overrightarrow{\upsilon}}^{\prime}}_{k}^{2}\approx 2\phantom{\rule{0.2em}{0ex}}{(\frac{{f}_{n1}{\overrightarrow{z}}_{k}}{)}2+2\phantom{\rule{0.2em}{0ex}}{(\frac{{f}_{n2}{\overrightarrow{z}}_{k}}{)}2\phantom{\rule{0.2em}{0ex}}.}^{}}^{}}{}$$

(66)

Note that the resulting update Eq. (62) is always in the direction of descent since the estimated second derivative is always positive. Theoretically, it is necessary to do a line search along the Gauss-Newton update direction to ensure convergence. In practice, we find that the objective function decreases reliably for each full Newton's step.

While *υ* and Φ* _{υ}*(

Since the velocity field *υ* is stationary in the case of the one parameter subgroup of diffeomorphisms, *υ* is clearly continuous (and in fact *C ^{∞}*) in

To compute the final deformation of an image, we have to estimate exp(*υ*) at least at the set of image grid points. We can compute exp(*υ*) by numerically integrating the smoothly interpolated velocity field *υ* with Euler integration. In this case, the estimate becomes arbitrarily close to the true exp(*υ*) as the number of integration time steps increases. With a sufficiently large number of integration steps, we expect the estimate to be invertible and the resulting transformation to be diffeomorphic.

The parameterization of diffeomorphisms by a stationary velocity field is popular due to the “scaling and squaring” approach [3] for computing exp(*υ*). Instead of Euler integration, the “scaling and squaring” approach iteratively composes displacement fields. Because we are working on the sphere *S*^{2}, the “scaling and squaring” procedure discussed in [3] has to be slightly modified:

$${\Phi}_{\frac{1}{{2}^{K}}\upsilon}(x)={\left\{{\Psi}_{n}\phantom{\rule{0.2em}{0ex}}\left({E}_{n}\frac{1}{{2}^{K}}\upsilon ({x}_{n})\right)\right\}}_{n=1,,N}$$

(67)

$$\begin{array}{c}{\Phi}_{\frac{1}{{2}^{K-1}}\upsilon}(x)={\Phi}_{\frac{1}{{2}^{K}}\upsilon}\left({\Phi}_{\frac{1}{{2}^{K}}\upsilon}(x)\right)\\ & {\Phi}_{\upsilon}(x)={\Phi}_{\frac{1}{2}\upsilon}\left({\Phi}_{\frac{1}{2}\upsilon}(x)\right)\phantom{\rule{0.2em}{0ex}},\end{array}$$

(68)

where Ψ* _{n}* is the local coordinate chart defined in Eq. (10), such that Ψ

While “scaling and squaring” converges to the true answer as *K* approaches ∞ in the continuous case, in the discrete case, composition of the displacement fields requires interpolation of displacement fields, introducing errors in the process. In particular, suppose Φ* _{t0υ}* (

As discussed in Appendix A-B, we employ barycentric interpolation, followed by normalization to ensure the warp stays on the unit sphere. In practice, we find that the resulting transformation is indeed diffeomorphic. Technically speaking, since we use linear interpolation for the displacement field, the transformation is only homeomorphic rather than diffeomorphic. This is because the transformation is continuous, but not differentiable across mesh edges. However, we follow the convention of [3], [4], [66] who call their transformation diffeomorphic even though they are homeomorphic.

^{1}There are two versions of the code (matlab and ITK) available at http://sites.google.com/site/yeoyeo02/software/sphericaldemonsrelease. The matlab code is used in the experiments of this paper. The preliminary ITK code [35], [36], [37] can also be found at http://hdl.handle.net/10380/3117.

^{2}The matlab code was used for this paper. The ITK code is still preliminary. Please check website http://sites.google.com/site/yeoyeo02/software/sphericaldemonsrelease for updates.

B.T. Thomas Yeo, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, USA.

Mert R. Sabuncu, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, USA.

Tom Vercauteren, Mauna Kea Technologies, Paris, France.

Nicholas Ayache, Asclepios Group, INRIA, Sophia Antipolis, France.

Bruce Fischl, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, USA; Department of Radiology, Harvard Medical School, Charlestown, USA and the Divison of Health Sciences and Technology, Massachusetts Institute of Technology, Cambridge, USA.

Polina Golland, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, USA.

1. Allassonniere S, Amit Y, Trouvé A. Toward a Coherent Statistical Framework for Dense Deformable Template Estimation. Journal of the Royal Statistical Society, Series B. 2007;69(1):3–29.

2. Amunts K, Malikovic A, Mohlberg H, Schormann T, Zilles K. Brodmann's Areas 17 and 18 Brought into Stereotaxic Space - Where and How Variable? NeuroImage. 2000;11:66–84. [PubMed]

3. Arsigny V, Commowick O, Pennec X, Ayache N. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 4190. LNCS; 2006. A Log-Euclidean Framework for Statistics on Diffeomorphisms; pp. 924–931. [PubMed]

4. Ashburner J. A Fast Diffeomorphic Image Registration Algorithm. NeuroImage. 2007;38:95–113. [PubMed]

5. Ashburner J, Andersson J, Friston K. High-Dimensional Image Registration using Symmetric Priors. NeuroImage. 1999;9:619–628. [PubMed]

6. Augustinack J, Van der Kouwe A, Blackwell M, Salat D, Wiggins C, Frosch M, Wiggins G, Potthast A, Wald L, Fischl B. Detection of Entorhinal Layer ii Using 7 Tesla Magnetic Resonance Imaging. Annals of Neurology. 2005;57(4):489–494. [PMC free article] [PubMed]

7. Avants B, Gee J. Geodesic Estimation for Large Deformation Anatomical Shape Averaging and Interpolation. NeuroImage. 2004;23:139–150. [PubMed]

8. Bakircioglu M, Joshi S, Miller M. Landmark Matching on Brain Surfaces via Large Deformation Diffeomorphisms on the Sphere. Proc SPIE Medical Imaging. 1999;3661:710–715.

9. Beg M, Miller M, Trouvé A, Younes L. Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. International Journal on Computer Vision. 2005;61(2):139–157.

10. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society B. 1995;57(1):256–300.

11. Birkhoff G, Rota G. Ordinary Differential Equations. John Wiley and Sons Inc; 1978.

12. Bro-Nielsen M, Gramkow C. Fast Fluid Registration of Medical Images. Proceedigs Visualization in Biomedical Computing. 1996;1131:267–276.

13. Brodmann K. Vergleichende Lokalisationslehre der Grohirnrinde in Ihren Prinzipien Dargestellt auf Grund des Zellenbaues. 1909.

14. Cachier P, Bardinet E, Dormont D, Pennec X, Ayache N. Iconic Feature Based Non-Rigid Registration: The PASHA algorithm. Compute Vision and Image Understanding. 2003;89(2-3):272–298.

15. Chefd'hotel C. Scale Space and Variational Methods in Computer Vision. Vol. 4485. LNCS; 2007. A Method for the Transport and Registration of Images on Implicit Surfaces; pp. 860–870.

16. Christensen G, Rabbit R, Miller M. Deformable Templates Using Large Deformation Kinematics. IEEE Transactions on Image Processing. 1996;5(10):1435–1447. [PubMed]

17. Chung M, Robbins S, Dalton K, Davidson R, Alexander A, Evans A. Cortical Thickness Analysis in Autism with Heat Kernel Smoothing. NeuroImage. 2005;25(4):1256–1265. [PubMed]

18. Commowick O, Stefanescu R, Fillard P, Arsigny V, Ayache N, Pennec X, Malandain G. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 3750. LNCS; 2005. Incorporating Statistical Measures of Anatomical Variability in Atlas-To-Subject Registration for Conformal Brain Radiotherapy; pp. 927–934. [PubMed]

19. Dale A, Fischl B, Sereno M. Cortical Surface-Based Analysis I: Segmentation and Surface Reconstruction. NeuroImage. 1999;9(2):179–194. [PubMed]

20. Desikan R, Segonne F, Fischl B, Quinn B, Dickerson B, Blacker D, Buckner R, Dale A, Maguire R, Hyman B, Albert M, Killiany R. An Automated Labeling System for Subdividing the Human Cerebral Cortex on MRI Scans into Gyral Based Regions of Interest. NeuroImage. 2006;31(3):968–980. [PubMed]

21. Dubuisson M, Jain A. A Modified Hausdorff Distance for Object Matching. Proceedings of the 12th IAPR International Conference on Pattern Recognition. 1994;1:566–568.

22. Durrleman S, Pennec X, Trouvé A, P, Ayache N. Inferring brain variability from diffeomorphic deformations of currents: an integrative approach. Medical Image Analysis. 2008;12(5):626–637. PMID: 18658005. [PMC free article] [PubMed]

23. Eckstein I, Joshi A, Jay Kuo C, Leahy R, Desbrun M. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 4791. LNCS; 2007. Generalized Surface Flows for Deformable Registration and Cortical Matching; pp. 692–700. [PMC free article] [PubMed]

24. Eickhoff S, et al. A new SPM Toolbox for Combining Probabilistic Cytoarchitectonic Maps and Functional Imaging Data. NeuroImage. 2005;25(4):1325–1335. [PubMed]

25. Fischl B, Rajendran N, Busa E, Augustinack J, Hinds O, Yeo BTT, Mohlberg H, Amunts K, Zilles K. Cortical Folding Patterns and Predicting Cytoarchictecture. Cerebral Cortex. 2008;18(8):1973–1980. [PubMed]

26. Fischl B, Sereno M, Dale A. Surface-Based Analysis II: Inflation, Flattening, and a Surface-Based Coordinate System. NeuroImage. 1999;9(2):195–207. [PubMed]

27. Fischl B, Sereno M, Tootell R, Dale A. High-resolution Inter-subject Averaging and a Coordinate System for the Cortical Surface. Human Brain Mapping. 1999;8(4):272–284. [PubMed]

28. Fischl B, van der Kouwe A, Destrieux C, Halgren E, Segonne F, Salat D, Busa E, Seidman L, Goldstein J, Kennedy D, Caviness V, Makris N, Rosen B, Dale A. Automatically Parcellating the Human Cerebral Cortex. Cerebral Cortex. 2004;14:11–22. [PubMed]

30. Geng X, Kumar D, Christensen G. Proceedings of the International Conference on Information Processing in Medical Imaging. Vol. 3564. LNCS; 2005. Transitive Inverse-Consistent Manifold Registration; pp. 468–479. [PubMed]

31. Glaunès J, Vaillant M, Miller M. Landmark Matching via Large Deformation Diffeomorphisms on the Sphere. Journal of Mathematical Imaging and Vision. 2004;20:179–200.

32. Hagler D, Saygin A, Sereno M. Smoothing and Cluster Thresholding for Cortical Surface-Based Group Analysis for fMRI Data. NeuroImage. 2006;33(4):1093–1103. [PMC free article] [PubMed]

33. Hernandez M, Bossa M, Olmos S. Registration of Anatomical Images Using Geodesic Paths of Diffeomorphisms Parameterized with Stationary Velocity Fields. Proceedings of the Workshop on Mathematical Methods in Biomedical Image Analysis, International Conference on Computer Vision. 2007:1–9.

34. Hernandez M, Olmos S. Gauss-Newton Optimization in Diffeomorphic Registration. Proceedings of the International Symposium on Biomedical Imaging: From Nano to Macro. 2008:1083–1086.

35. Ibanez L, Audette M, Yeo BTT, Golland P. Rotational Registration of Spherical Surfaces Represented as QuadEdge Meshes. Insight Journal. 2009. Jan-Jun. http://hdl.handle.net/10380/3063.

36. Ibanez L, Audette M, Yeo BTT, Golland P. Spherical Demons Registration of Spherical Surfaces. Insight Journal. 2009. Jul-Dec. http://hdl.handle.net/10380/3117.

37. Ibanez L, Yeo BTT, Golland P. Iterative Smoothing of Field Data in Spherical Meshes. Insight Journal. 2009. Jul-Dec. http://hdl.handle.net/10380/3091.

38. Klein A, Hirsch J. Mindboggle: A Scatterbrained Approach to Automate Brain Labeling. NeuroImage. 2005;24:261–280. [PubMed]

39. Kühnel W. Differential Geometry: Curves-Surfaces-Manifolds. second. American Mathematical Society; 2006.

40. Lefèvre J, Baillet S. Optical Flow and Advection on 2-Riemannian Manifolds: A Common Framework. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2008;30(6):1081–1092. [PubMed]

41. Lohmann G, Von Cramon D. Automatic Labelling of the Human Cortical Surface using Sulcal Basins. Medical Image Analysis. 2000;4:179–188. [PubMed]

42. Mangin JF, Frouin V, Bloch I, Régis J, Lopez-Krahe J. Automatic Construction of an Attributed Relational Graph Representing the Cortex Topography using Homotopic Transformations. Proceedings of SPIE. 1994;2299:110–121.

43. Marsland S, McLachlan R. Proceedings of the International Conference on Information Processing in Medical Imaging. Vol. 4548. LNCS; 2007. A Hamiltonian Particle Method for Diffeomorphic Image Registration; pp. 396–407. [PubMed]

44. Munkres J. Analysis on Manifold. Westview Press; 1991.

45. Nielsen M, Florack L, Deriche R. Regularization, Scale-Space, and Edge Detection Filters. Journal of Mathematical Imaging and Vision. 1997;7:291–307.

46. Nielsen M, Johansen P, Jackson AD, Lautrup B. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 2489. LNCS; 2002. Brownian Warps: A Least Committed Prior for Non-rigid Registration; pp. 557–264.

47. Olver P. Applications of Lie Groups to Differential Equations. 2nd. Springer-Verlag; 1993.

48. Pennec X, Stefanescu R, Arsigny V, Fillard P, Ayache N. Riemannian Elasticity: A Statistical Regularization Framework for Nonlinear Registration [PubMed]

49. Press W, Teukolsky S, Vetterling W, Flannery B. Numerical Recipes in C: The Art of Scientific Computing. second. Cambridge University Press; 1992.

50. Qiu A, Miller M. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 4791. LNCS; 2007. Cortical Hemisphere Registration via Large Deformation Diffeomorphic Metric Curve Mapping; pp. 186–193. [PMC free article] [PubMed]

51. Rettmann M, Han X, Xu C, Prince J. Automated Sulcal Segmentation Using Watersheds on the Cortical Surface. NeuroImage. 2002;15:329–344. [PubMed]

52. Riviere D, Mangin JF, Papadopoulos-Orfanos D, Martinez J, Frouin V, Regis J. Automatic Recognition of Cortical Sulci of the Human Brain Using a Congregation of Neural Networks. Medical Image Analysis. 2002;6:77–92. [PubMed]

53. Sabuncu M, Yeo BTT, Vercauteren T, Van Leemput K, Golland P. Assymetric Image-Template Registration. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) 2009;5761:565–573.

54. Schleicher A, Amunts K, Geyer S, Morosan P, Zilles K. Observer-Independent Method for Microstructural Parcellation of Cerebral Cortex: A Quantitative Approach to Cytoarchitectonics. NeuroImage. 1999;9:165–177. [PubMed]

55. Schormann T, Zilles K. Three-Dimensional Linear and Non-linear Transformations: An Integration of Light Microscopical and MRI Data. Human Brain Mapping. 1998;6:339–347. [PubMed]

56. Tao X, Prince J, Davatzikos C. Using a Statistical Shape Model to Extract Sulcal Curves on the Outer Cortex of the Human Brain. IEEE Transactions on Medical Imaging. 2002;21:513–524. [PubMed]

57. Thirion J. Image Matching as a Diffusion Process: an Analogy with Maxwell's Demons. Medical Image Analysis. 1998;2(3):243–260. [PubMed]

58. Thompson P, Toga A. A Surface-Based Technique for Warping 3-Dimensional Images of the Brain. IEEE Transactions on Medical Imaging. 1996;15(4):1–16. [PubMed]

59. Thompson P, Woods R, Mega M, Toga A. Mathematical/Computational Challenges in Creating Deformable and Probabilistic Atlases of the Human Brain. Human Brain Mapping. 2000;9(2):81–92. [PubMed]

60. Tosun D, Prince J. Proceedings of the International Conference on Information Processing in Medical Imaging. Vol. 3565. LNCS; 2005. Cortical Surface Alignment Using Geometry Driven Multispectral Optical Flow; pp. 480–492. [PubMed]

61. Trouvé A. Diffeomorphisms, Groups and Pattern Matching in Image Analysis. International Journal on Computer Vision. 1998;28(3):213–221.

62. Tu Z, Zheng S, Yuille A, Reiss A, Dutton R, Lee A, Galaburda A, Dinov I, Thompson P, Toga A. Automated Extraction of the Cortical Sulci Based on a Supervised Learning Approach. IEEE Transactions on Medical Imaging. 2007;26(4):541–552. [PubMed]

63. Twining C, Davies R, Taylor C. Proceedings of the International Conference on Information Processing in Medical Imaging. Vol. 4584. LNCS; 2007. Non-Parametric Surface-Based Regularisation for Building Statistical Shape Models; pp. 738–750. [PubMed]

64. Van Essen D, Drury H, Joshi S, Miller M. Functional and Structural Mapping of Human Cerebral Cortex: Solutions are in the Surfaces. Proceedings of the National Academy of Sciences. 1996;95(3):788–795. [PubMed]

65. Van Leemput K. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 4190. LNCS; 2006. Probabilistic Brain Atlas Encoding Using Bayesian Inference; pp. 704–711. [PubMed]

66. Vercauteren T, Pennec X, Perchant A, Ayache N. Diffeomorphic demons: Efficient non-parametric image registration. NeuroImage. 2009 Mar;45(1, Supp.1):S61–S72. PMID: 19041946. [PubMed]

67. Yeo BTT, Sabuncu M, Desikan R, Fischl B, Golland P. Effects of registration regularization and atlas sharpness on segmentation accuracy. Medical Image Analysis. 2008;12(5):603–615. [PMC free article] [PubMed]

68. Yeo BTT, Sabuncu M, Golland P, Fischl B. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 5761. LNCS; 2009. Task-Optimal Registration Cost Functions; pp. 598–606. [PMC free article] [PubMed]

69. Yeo BTT, Sabuncu M, Vercauteren T, Ayache N, Fischl B, Golland P. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 5241. LNCS; 2008. Spherical Demons: Fast Surface Registration; pp. 745–753. [PMC free article] [PubMed]

70. Yeo BTT, Vercauteren T, Fillard P, Pennec X, Golland P, Ayache N, Clatz O. DT-REFinD: Diffusion Tensor Registration with Exact Finite-Strain Differential. IEEE Transactions on Medical Imaging. 2009 In Press. [PubMed]

71. Zilles K, Schleicher A, Palomero-Gallagher N, Amunts K. Quantitative Analysis of Cyto- and Receptor Architecture of the Human Brain. Elsevier; 2002.

72. Zou G, Hua J, Muzik O. Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) Vol. 4791. LNCS; 2007. Non-rigid Surface Registration using Spherical Thin-Plate Splines; pp. 367–374. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |