Home  About  Journals  Submit  Contact Us  Français 
Image registration is typically formulated as an optimization problem with multiple tunable, manually set parameters. We present a principled framework for learning thousands of parameters of registration cost functions, such as a spatiallyvarying tradeoff between the image dissimilarity and regularization terms. Our approach belongs to the classic machine learning framework of model selection by optimization of crossvalidation error. This second layer of optimization of crossvalidation error over and above registration selects parameters in the registration cost function that result in good registration as measured by the performance of the specific application in a training data set. Much research effort has been devoted to developing generic registration algorithms, which are then specialized to particular imaging modalities, particular imaging targets and particular postregistration analyses. Our framework allows for a systematic adaptation of generic registration cost functions to specific applications by learning the “free” parameters in the cost functions. Here, we consider the application of localizing underlying cytoarchitecture and functional regions in the cerebral cortex by alignment of cortical folding. Most previous work assumes that perfectly registering the macroanatomy also perfectly aligns the underlying cortical function even though macroanatomy does not completely predict brain function. In contrast, we learn 1) optimal weights on different cortical folds or 2) optimal cortical folding template in the generic weighted sum of squared differences dissimilarity measure for the localization task. We demonstrate stateoftheart localization results in both histological and functional magnetic resonance imaging data sets.
IN medical image analysis, registration is necessary to establish spatial correspondence across two or more images. Traditionally, registration is considered a preprocessing step [Fig. 1(a)]. Images are registered and are then used for other image analysis applications, such as voxelbased morphometry and shape analysis. Here, we argue that the quality of image registration should be evaluated in the context of the application. In particular, we propose a framework for learning the parameters of registration cost functions that are optimal for a specific application. Our framework is therefore equivalent to classic machine learning approaches of model selection by optimization of crossvalidation error [33], [43], [58].
Image registration is typically formulated as an optimization problem with a cost function that comprises an image dissimilarity term and a regularization term [Fig. 1(a)]. The parameters of the cost function are frequently determined manually by inspecting the quality of the image alignment to account for the characteristics (e.g., resolution, modality, signaltonoise ratio) of the image data. During this process, the final task is rarely considered in a principled fashion. Furthermore, the variability of the results due to these tunable parameters is rarely reported in the literature. Yet, recent work has shown that taking into account the tradeoff between the regularization and similarity measure in registration can significantly improve population analysis [40] and segmentation quality [10], [79].
In addition to improving the performance of applications downstream, taking into account the endgoal of registration could help resolve ambiguities and the illposed nature of image registration.
In this paper, we propose a taskoptimal registration framework that optimizes parameters of any smooth family of registration cost functions on a training set, with the aim of improving the performance of a particular task for a new image [Fig. 1(b)]. The key idea is to introduce a second layer of optimization over and above the usual registration. This second layer of optimization assumes the existence of a smooth cost function or crossvalidation error metric [g in Fig. 1(b)] that evaluates the performance of a particular task given the output of the registration step for a training data set. The training data provides additional information not present in a test image, allowing the taskspecific cost function to be evaluated during training. For example, if the task is segmentation, we assume the existence of a training data set with ground truth segmentation and a smooth cost function (e.g., Dice overlap measure) that evaluates segmentation accuracy. If the registration cost function employs a single parameter, then the optimal parameter value can be found by exhaustive search [79]. With multiple parameters, exhaustive search is not possible. Here, we establish conditions for which the space of local minima is locally smooth and demonstrate the optimization of thousands of parameters by gradient descent on the space of local minima, selecting registration parameters that result in good registration local minima as measured by the taskspecific cost function in the training data set.
We validate our framework on two datasets. The first dataset consists of 10 ex vivo brains with the BAs of each subject obtained via histology [4], [84] and mapped onto the cortical surface representation of each subject obtained via MRI [18]. The second dataset consists of 42 in vivo brains with functional region MT+ (V5) defined using functional magnetic resonance imaging (fMRI). Here, our task is defined to be the localization of BAs and MT+ in the cortical surface representation via the registration of the cortical folding pattern. While it is known that certain cytoarchitectonically or functionallydefined areas, such as V1 or BA28, are spatially consistent with respect to local cortical geometry, other areas, such as BA44, present a challenge for existing localization methods [18], [20].We learn the weights of the weighted sum of squared differences (wSSD) family of registration cost functions and/or estimate an optimal macroanatomical template for localizing the cytoarchitectural and functional regions using only the cortical folding pattern. We demonstrate improvement over existing methods [18].
An alternative approach to overcome the imperfect correlation between anatomy and function is to directly use the functional data for establishing acrosssubject functional correspondence [54], [56]. However, these approaches require extra data acquisition (such as fMRI scans) of all future test subjects. In contrast, our method aims to learn the relationship between macroanatomy and function (or cytoarchitectonics) in a training data set containing information about both macroanatomy and function (or cytoarchitectonics). We use this information to localize function (or cytoarchitectonics) in future subjects, for which only macroanatomical information is available.
Our approach belongs to the class of “wrapper methods” for model or feature selection in the machine learning literature [27], [34]. In particular, our model selection criterion of applicationspecific performance is equivalent to the use of crossvalidation error in various model selection algorithms [33], [43], [58]. Unlike feature selection methods that operate in a discrete parameter space, we work in a continuous parameter space. Consequently, standard algorithms in the “wrapper methods” literature do not apply to this problem.
Instead, our resulting optimization procedure borrows heavily from the mathematical field of continuation methods [2]. Continuation methods have been recently introduced to the machine learning community for computing the entire path of solutions of learning problems (e.g., SVM or Lasso) as a function of a single regularization parameter [16], [28], [46]. For example, the cost function in Lasso [67] consists of the tradeoff between a leastsquares term and a L_{1} regularization term. Leastangles regression (LARS) allows one to compute the entire set of solutions of Lasso as a function of the tradeoff parameter [16]. Because we deal with multiple (thousands of) parameters, it is impossible for us to compute the entire solution manifold. Instead, we trace a path within the solution manifold that improves the taskspecific cost function. Furthermore, registration is not convex (unlike SVM and Lasso), resulting in several theoretical and practical conundrums that we have to overcome, some of which we leave for future work.
The wSSD similarity measure implicitly assumes an independent Gaussian distribution on the image intensities, where the weights correspond to the precision (reciprocal of the variance) and the template corresponds to the mean of the Gaussian distribution. The weights can be set to a constant value [6], [31] or a spatiallyvarying variance can be estimated from the intensities of registered images [19]. However, depending on the wSSD regularization tradeoff, the choice of the scale of the variance is still arbitrary [79]. With weaker regularization, the training images will be better aligned, resulting in lower variance estimates.
Recent work in probabilistic template construction resolves this problem by either marginalizing the tradeoff under a Bayesian framework [1] or by estimating the tradeoff with the minimum description length principle [71]. While these methods are optimal for “explaining the images” under the assumed generative models, it is unclear whether the estimated parameters are optimal for applicationspecific tasks. After all, the parameters for optimal image segmentation might be different from those for optimal group analysis. In contrast, Van Leemput [74] proposes a generative model for image segmentation. The estimated parameters are therefore Bayesianoptimal for segmentation. When considering one global tradeoff parameter, a more direct approach is to employ crossvalidation of segmentation accuracy and to perform an exhaustive search over the values of the tradeoff parameter [79]. This is infeasible for multiple parameters.
By learning the weights of the wSSD, we implicitly optimize the tradeoff betweeen the dissimilarity measure and regularization. Furthermore, the tradeoff we learn is spatially varying. Previous work on learning a spatially varying regularization prior suffers from the lack of ground truth (nonlinear) deformations. For example, [10], [25], [35] assume that the deformations obtained from registering a set of training images can be used to estimate a registration regularization to register new images. However, a change in the parameters of the registration cost function used by these methods to register the training images would lead to a different set of training deformations and thus a different prior for registering new images. Furthermore, the methods are inconsistent in the sense that the learned prior applied on the training images will not result in the same training deformations obtained previously.
While there has been efforts in obtaining ground truth humanannotated deformation fields [37], the images considered typically have welldefined correspondences, rather than for example, the brain images of two different subjects. As suggested in the previously presented examples (Fig. 2), the concept of “ground truth deformations” may not always be welldefined, since the optimal registration may be a function of the application at hand. In contrast, image segmentation is generally better defined in the sense that ground truth segmentation is usually known. Our problem therefore differs from recent work on learning segmentation cost functions [42], [70], [83]. In this paper, we avoid the need for ground truth deformations by focusing on the application of registrationbased segmentation, where ground truth segmentations are better defined and available. However, our framework is general and can be applied whenever a postregistration application can be well quantified by a smooth applicationspecific performance cost function.
This paper is organized as follows. In the next section, we introduce the taskoptimal registration framework. We specialize the framework to align hidden labels in Section III. We present localization experiments in Section IV and discuss outstanding issues in Section V. This paper extends a previously presented conference article [80] and contains detailed derivations, discussions and experiments that were omitted in the conference version.
In this section, we present the taskoptimal registration framework for learning the parameters of a registration cost function. Given an image I, let f (w, Γ) denote a smooth registration cost function, with parameters w and a spatial transformation Γ. For example
where T is the template image, λ is the tradeoff between the image dissimilarity measure and the regularization on the transformation Γ, I Γ denotes the deformed and resampled image I. f is therefore also a function of the image I, which we suppress for conciseness. The optimal transformation Γ* minimizes the cost function for a given set of parameters w
We emphasize that Γ* is a function of w since a different set of parameters w will result in a different solution to (2.2) and thus will effectively define a different image coordinate system.
The resulting deformation Γ* is used to warp the input image or is itself used for further tasks, such as image segmentation or voxelbased morphometry. We assume that the task performance can be measured by a smooth cost function (or crossvalidation error metric) g, so that a smaller value of g (Γ*(w)) corresponds to better task performance. g is typically a function of additional input data associated with a subject (e.g., manual segmentation labels if the task is automatic segmentation), although we suppress this dependency in the notation for conciseness. This auxiliary data is only available in the training set; g cannot be evaluated for the new image.
Given a set of N training subjects, let ${\mathrm{\Gamma}}_{n}^{*}(w)$ denote the solution of (2.2) for training subject n for a fixed set of parameters w and g_{n}(${\mathrm{\Gamma}}_{n}^{*}(w)$) denote the task performance for training subject n using the deformation ${\mathrm{\Gamma}}_{n}^{*}(w)$ and other information available for the nth training subject. A different set of parameters w would lead to different task performance g_{n}(${\mathrm{\Gamma}}_{n}^{*}(w)$). We seek the parameters w* that generalize well to a new subject: registration of a new subject with w* yields the transformation Γ*(w*) with a small taskspecific cost g(Γ*(w*)). One approach to solve this functional approximation problem [17] is regularized risk minimization. Let Reg(w) denote regularization on w and define
Regularization risk minimization seeks
The optimization is difficult because while we assume g_{n} to be smooth, the input to g_{n} (·) is itself the local minimum of another nonlinear cost function f. Furthermore, evaluating the cost function G for only one particular set of parameters w requires performing N different registrations!
In this section, we provide theoretical characterizations of the optimization problem in (2.4). If Γ*(w) is defined strictly to be a global registration optimum, then Γ*(w) is clearly not a smooth function of w, since a small change in w can result in a big change in the global registration optimum. This definition is also impractical, since the global optimum of a nonlinear optimization problem cannot be generally found in practice. Instead, we relax the definition of Γ*(w) to be a local minimum of the registration cost function for fixed values of w. Here, we derive conditions in which Γ*(w) is locally a smooth function of w, so we can employ gradient descent to optimize (2.4).
Let Γ*(w_{0}) denote a local minimum of the registration cost function for a fixed w = w_{0}. Suppose we perturb w by an infinitestimally small δw, so that Γ*(w_{0}) is no longer the registration local minimum for w = w_{0} + δw. We consider two representations of this change in local minimum.
Additive deformation models arise when the space of deformations is a vector space, such as the space of displacement fields or positions of Bspline control points. At each iteration of the registration algorithm, deformation updates are added to the current deformation estimates. The additive model is general and applies to many nonconvex, smooth optimization problems outside of registration. Most registration algorithms can in fact be modeled with the additive framework.
In some registration algorithms, including that used in this paper, it is more natural to represent deformation changes through composition rather than additions [7], [61], [75]. For example, in the diffeomorphic variants of the demons algorithm [75], [81], [82], the diffeomorphic transformation Γ is represented as a dense displacement field. At each iteration, the transformation update is restricted to be a one parameter subgroup of diffeomorphism parameterized by a stationary velocity field. The diffeomorphic transformation update is then composed with, rather than added to, the current estimate of the transformation, thus ensuring that the resulting transformation is diffeomorphic.
Let Γ*(w_{0} + δw) = Γ*(w_{0}) + δΓ*(w_{0}, δw) denote the new locally optimal deformation for the updated set of parameters w_{0} + δw. The following proposition characterizes the existence and uniqueness of δΓ*(w_{0}, δw) as δw is varied. In particular, we show that under some mild conditions, δΓ*(w_{0}, δw) is a welldefined smooth function in the neighborhood of (w_{0}, Γ*(w_{0})). In the remainder, we use ${\partial}_{x},{\partial}_{x}^{2}$ and ${\partial}_{x,y}^{2}$ to denote the corresponding partial derivatives.
Proposition 1: If the Hessian^{1} ${\partial}_{\mathrm{\Gamma}}^{2}f({w}_{0},\mathrm{\Gamma})$ is positive definite at Γ = Γ*(w_{0}), then there exists an ε > 0, such that for all δw, ‖δw‖ < ε, a unique continuous function δΓ*(w_{0}, δw) exists with δΓ*(w_{0}, 0) = 0. Furthermore, δΓ* has the same order of smoothness as _{Γ} f.
Proof: We define the vectorvalued function h(w, Γ) _{Γ} f (w, Γ). Since Γ*(w_{0}) is a local minimum of f(w_{0},Γ), we have
At (w_{0}, Γ*(w_{0})), the Hessian matrix ${\partial}_{\mathrm{\Gamma}}^{2}f({w}_{0},\mathrm{\Gamma})={\partial}_{\mathrm{\Gamma}}^{2}h(w,\mathrm{\Gamma})$ is positive definite by the assumption of the proposition and is therefore invertible. By the Implicit Function Theorem [51], there exists an ε > 0, such that for all δw, ‖δw‖ < ε, there is a unique continuous function δΓ*(w_{0}, δw) such that h(w_{0} + δw, Γ*(w_{0}) + δΓ*(w_{0}, δw)) = 0 and δΓ*(w_{0}, 0) = 0. Furthermore, δΓ*(w_{0}, δw) has the same order of smoothness as h.
Because the Hessian of f is smooth and the eigenvalues of a matrix depend continuously on the matrix [72], there exists a small neighborhood around (w_{0}, Γ*(w_{0})) in which the eigenvalues of ${\partial}_{\mathrm{\Gamma}}^{2}f(w,\mathrm{\Gamma})$ are all greater than 0. Since both sufficient conditions for a local minimum are satisfied (zero gradient and positive definite Hessian), Γ*(w_{0}) + δΓ*(w_{0}, dw) is indeed a new local minimum close to Γ*(w_{0}).
Observe that the conditions in Proposition 1 are stronger than those of typical nonlinear optimization problems. In particular, we do not just require the cost functions f and g to be smooth, but also that the Hessian ${\partial}_{\mathrm{\Gamma}}^{2}f({w}_{0},\mathrm{\Gamma})$ be positive definite at the local minimum. At (w_{0}, Γ*(w_{0})), by definition, the Hessian ${\partial}_{\mathrm{\Gamma}}^{2}f({w}_{0},\mathrm{\Gamma})$ is positive semidefinite, so the positive definite condition in Proposition 1 should not be too restrictive. Unfortunately, degeneracies may arise for local minima with a singular Hessian. For example, let Γ be the 1 × 2 vector [a b] and f (Γ, w) = (ab − w)^{2}. Then for any value of w, there is an infinite number of local minima Γ*(w) corresponding to ab = w. Furthermore, the Hessian at any local minimum is singular. In this case, even though f is infinitely differentiable, there is an infinite number of local minima near the current local minimum Γ*(w_{0}), i.e., δΓ*(w_{0}, δw) is not a welldefined function and the gradient is not defined. Consequently, the parameters w of local registration minima whose Hessians are singular are also local minima of the taskoptimal optimization (2.4). The proof of Proposition 1 follows the ideas of the Continuation Methods literature [2]. We include the proof here to motivate the more complex composition of deformations model.
Let Γ*(w_{0}) be the registration local minimum at w_{0} and δΓ (υ) denote an update transformation parameterized by υ, so that δΓ (0) corresponds to the identity transform. For example, υ could be a stationary [75], [81], [82], nonstationary [8] velocity field parameterization of diffeomorphism, positions of spline control points [52] or simply displacement fields [59]. In the composition model, Γ*(w_{0}) is a local minimum if and only if there exists an ε > 0, such that f (w_{0}, Γ*(w_{0})) < f (w_{0}, Γ*(w_{0}) δΓ (υ)) for all values of ‖υ‖ < ε.
Let Γ*(w_{0}) δΓ(υ*(w_{0}, δw)) denote the new locally optimal deformation for the new parameters w_{0} + δw. In general, there might not exist a single update transformation δΓ(υ*(w_{0}, δw)) that leads to a new local minimum under a perturbation of the parameters w, so that there is no correponding version of Proposition 1 for the general composition model. However, in the special case of the composition of diffeomorphisms model [75], [81], [82] employed in this paper, the following proposition characterizes the existence and uniqueness of υ*(w_{0}, δw) as δw is varied.
Proposition 2: If the Hessian ${\partial}_{\upsilon}^{2}f({w}_{0},{\mathrm{\Gamma}}^{*}({w}_{0})\circ \delta \mathrm{\Gamma}(\upsilon ))$ is positive definite at υ = 0, then there exists an ε > 0, such that for all δw, ‖δw‖ < ε, a unique continuous function υ*(w_{0}, δw) exists, such that υ*(w_{0}, δw) is the new local minimum for parameters w_{0} + δw and υ*(w_{0}, 0) = 0. Furthermore, υ*(w_{0}, δw) has the same order of smoothness as f.
Proof: The proof is a more complicated version of Proposition 1 and so we leave the details to Appendix A.
Just like in the case of the additive deformation model, the parameters w of local registration minima that do not satisfy the conditions of Proposition 2 are also local minima of the taskoptimal optimization (2.4). In the next section, we derive exact and approximate gradients of G.
We now discuss the optimization of the regularized task performance G.
In the previous section, we showed that at (w_{0}, Γ*(w_{0})) with a positive definite Hessian, δΓ*(w_{0}, δw) is a smooth welldefined function such that Γ*(w_{0}) + δΓ*(w_{0}, δw) is the new local minimum at w_{0} + δw. Therefore, we can compute the derivatives of Γ*(w) with respect to w at w_{0}, allowing us to traverse a curve of local optima, finding values of w that improve the taskspecific cost function for the training images. We first perform a Taylor expansion of _{Γ} f (w, Γ) at (w_{0}, Γ*(w_{0}))
where we dropped the term _{Γ} f (w, Γ)_{w0, Γ*(w0)} = 0. For δΓ = δΓ*(w_{0},δw), the lefthand side is equal to 0 and we can write
Therefore, by taking the limit δw 0, we get
Equation (2.8) tells us the direction of change of the local minimum at (w_{0}, Γ*(w_{0})). In practice, the matrix inversion in (2.8) is computationally prohibitive for highdimensional warps Γ. Here, we consider a simplification of (2.8) by setting the Hessian to be the identity
Since −_{Γ} f is the direction of gradient descent of the cost function (2.2), we can interpret (2.9) as approximating the new local minimum to be in the same direction as the change in the direction of gradient descent as w is perturbed.
Differentiating the cost function in (2.4), using the chain rule, we obtain
We note the subscript n on f indicates the dependency of the registration cost function on the nth training image.
In the previous section, we have shown that at (w_{0}, Γ*(w_{0})), assuming the conditions of Proposition 2 are true, υ*(w_{0}, δw) is a smooth welldefined function such that Γ*(w_{0}) δΓ (υ*(w_{0}, δw)) is the new local minimum. Therefore, we can compute the derivatives of υ* with respect to w. As before, by performing a Taylor expansion, we obtain
Appendix B provides the detailed derivations. Differentiating the cost function in (2.4), using the chain rule, we get
Once again, the subscript n on f indicates the dependency of the registration cost function on the nth training image.
Algorithm 1 summarizes the method for learning the taskoptimal registration parameters. Each line search involves evaluating the cost function G multiple times, which in turn requires registering the training subjects, resulting in a computationally intensive process. However, since we are initializing from a local optimum, for a small change in w, each registration converges quickly.
Data: A set of training images {I_{n}}  
Result: Parameters w that minimize the regularized task performance G [see (2.4)]  
Initialize w^{0}.  
repeat  
Step  1.  Given current values of w, estimate ${\mathrm{\Gamma}}_{n}^{*}(w)={\text{arg min}}_{{\mathrm{\Gamma}}_{n}}$, f_{n} (w, Γ_{n}), i.e., perform registration of each training subject n. 
Step  2.  Given current estimates (w, {Γ_{n} (w)}), compute the _{w}G gradient using either

Step  3.  Perform line search in the direction opposite to _{w}G [47]. 
Since nonlinear registration is dependent on initialization, the current estimates (w, Γ*(w)), which were initialized from previous estimates, might not be achievable when initializing the registration with the identity transform. The corresponding parameters w might therefore not generalize well to a new subject, which are typically initialized with the identity transform. To put this more concretely, suppose our current estimates of w and the registration local minima are (w = 5, Γ* (5) = 2). Next, we perform the gradient decent step and update w accordingly. For argument’s sake, let our new estimates of w and the registration local minima be (w = 5.1, Γ* (5.1) = 1.9). Note that this particular value of Γ* (5.1) = 1.9 is achieved by initializing the registration with Γ = 2. Had we initialized the registration with the identity transform (such as for a new subject), then Γ* (5.1) might instead be equal to 2.1, with possibly poorer application performance than (w = 5, Γ* (5) = 2). To avoid this form of overfitting, after every few iterations, we reregister the training images by initializing with the identity transform, and verify that the value of G is better than the current best value of G computed with initialization from the identity transform.
The astute reader will observe that the preceding discussion on “Addition Model” makes no assumptions specific to the taskoptimal registration problem. The framework can therefore also be applied to learn the cost functions in other applications that are formulated as nonlinear optimization problems solved by gradient descent.
We now instantiate the taskoptimal registration framework for localizing hidden labels in images. We demonstrate schemes for either 1) learning the weights of the wSSD family of registration cost functions or 2) estimating an optimal template image for localizing these hidden labels. We emphasize that the optimal template is not necessarily the average of the training images, since the goal is not to align image intensities across subjects, but to localize the hidden labels.
Suppose we have a set of training images {I_{n}} with some underlying ground truth structure manually labeled or obtained from another imaging modality (e.g., Brodmann areas from histology mapped onto cortical surface representations). We define our task as localizing the hidden structure in a test image. In the traditional pairwise registration approach [Fig. 3(a)], a single training subject is chosen as the template. After pairwise registration between the template and test images, the ground truth label of the template subject is used to predict that of the test subject. The goal of predicting the hidden structure in the test subject is typically not considered when choosing the training subject or registration algorithm. For hidden labels that are poorly predicted by local image intensity (e.g., BA44 discussed in Section IA), blind alignment of image intensities lead to poor localization.
In contrast, we pick one training subject as the initial template and use the remaining training images and labels [Fig. 3(b)] to learn a registration cost function that is optimal for aligning the labels of the training and template subjects—perfect alignment of the labels lead to perfect prediction of the labels in the training subjects by the template labels. After pairwise registration between the template and test subject using the optimal registration cost function, the ground truth label of the template subject is used to predict that of the test subject.
We limit ourselves to spherical images (i.e., images defined on a unit sphere), although it should be clear that the discussion readily extends to volumetric images. Our motivation for using spherical images comes from the representation of the human cerebral cortex as a closed 2D mesh in 3D. There has been much effort focused on registering cortical surfaces in 3D [14], [15], [24], [30], [65]. Since cortical areas—both structure and function—are arranged in a mosaic across the cortical surface, an alternative approach is to warp the underlying spherical coordinate system [19], [48], [60], [66], [69], [73], [79], [81]. Warping the spherical coordinate system establishes correspondences across the surfaces without actually deforming the surfaces in 3D. We assume that the meshes have already been spherically parameterized and represented as spherical images: a geometric attribute is associated with each mesh vertex, describing local cortical geometry.
To register a given image I_{n} to the template image T, we define the following cost function:
where transformation Γ_{n} maps a point x_{i} on the sphere S^{2} to another point Γ_{n} (x_{i}) ε S^{2}. The first term corresponds to the wSSD image similiarity. The second term is a percentage metric distortion regularization on the transformation Γ_{n} where _{i} is a predefined neighborhood around vertex i and d_{ij} is the original distance between the neighbors d_{ij} = ‖x_{i} − x_{j}‖ [79]. The weights {λ_{i}}’s are generalizations of the tradeoff parameter λ, allowing for a spatiallyvarying tradeoff between the image dissimilarity term and regularization: a higher weight ${\lambda}_{i}^{2}$ corresponds to placing more emphasis on matching the template image at spatial location x_{i} relative to the regularization. The parameterization of the weights as ${\lambda}_{i}^{2}$ ensures nonnegative weights.
In this work, we consider either learning the weights ${\lambda}_{i}^{2}$ or the template T for localizing BA labels or functional labels by aligning cortical folding pattern. Since the weights of the wSSD correspond to the precision of the Gaussian model, by learning the weights of wSSD, we are learning the precision of the Gaussian model and hence the uncertainty of the sulcal geometry. Optimizing ${\lambda}_{i}^{2}$ leads to placing nonuniform importance on matching different cortical folds with the aim of aligning the underlying cytoarchitectonics or function. For example, suppose there is a sulcus with functional regions that appear on either side of the sulcus depending on the subject. The algorithm may decide to place low weight on the “poorly predictive” sulcus. On the other hand, optimizing T corresponds to learning a cortical folding template that is optimal for localizing the underlying cytoarchitectonics or functional labels of the training subjects. In the case of the previously mentioned “unpredictive sulcus,” the algorithm might learn that the optimal cortical folding template should not contain this sulcus.
We choose to represent the transformation Γ_{n} as a composition of diffeomorphic warps {Φ_{k}} parameterized by a stationary velocity field, so that Γ_{n} = Φ_{1} Φ_{K} [75], [81], [82]. We note that our choice of regularization is different from the implicit hierarchical regularization used in Spherical Demons [81] since the Demons regularization is not compatible with our derivations from the previous section. Instead of the efficient 2Step Spherical Demons algorithm, we will use steepest descent. The resulting registration algorithm is still relatively fast, requiring about 15 min for registering fullresolution meshes with more than 100k vertices, compared with 5 min of computation for Spherical Demons on a Xeon 2.8GHz single processor machine.
In general, a smooth stationary velocity field υ parameterizes a diffeomorphism Φ via a stationary ODE: _{t}Φ(x,t) = (Φ(x, t)) with an initial condition Φ(x,0) = x. The solution at t = 1 is denoted as Φ(x,1) = Φ(x) = exp()(x), where we have dropped the time index. A solution can be computed efficiently using scaling and squaring [5]. This particular choice of representing deformations provides a computationally efficient method of achieving invertible transformations, which is a desirable property in many medical imaging applications. In our case, the velocity field is a tangent vector field on the sphere S^{2} and not an arbitrary 3D vector field.
To register subject n to the template image T for a fixed set of parameters w, let ${\mathrm{\Gamma}}_{n}^{0}$ be the current estimate of ${\mathrm{\Gamma}}_{n}^{*}$. We seek an update transformation exp () parameterized by a stationary velocity field
Let _{i} be the velocity vector tangent to vertex x_{i}, and = {_{i}} be the entire velocity field. We adopt the techniques in the Spherical Demons algorithm [81] to differentiate (3.1) with respect to , evaluated at = 0. Using the fact that the differential of exp ()at = 0 is the identity [44], i.e., [D exp(0)] = , we conclude that a change in velocity _{i} at vertex x_{i} does not affect exp()(x_{n}) for n ≠ i up to the first order derivatives. Defining $\nabla {I}_{n}({\mathrm{\Gamma}}_{n}^{0}({x}_{i}))$ to be the 1 × 3 spatial gradient of the warped image I_{n} ($({\mathrm{\Gamma}}_{n}^{0}(\xb7))$) at x_{i} and $\nabla {\mathrm{\Gamma}}_{n}^{0}({x}_{i})$ to be the 3 × 3 Jacobian matrix of ${\mathrm{\Gamma}}_{n}^{0}$ at x_{i}, we get the 1 × 3 derivative
We can perform gradient descent of the registration cost function f_{n} using (3.2) to obtain ${\mathrm{\Gamma}}_{n}^{*}$, which can be used to evaluate the regularized task performance G to be described in the next section. We also note that (3.2) instantiates _{} f_{n} within the mixed derivatives term in the taskoptimal gradient (2.16) for this application.
We represent the hidden labels in the training subjects as signed distance transforms on the sphere {L_{n}} [36]. We consider a pairwise approach, where we assume that the template image T has a corresponding labels with distance transform L_{T} and set the taskspecific cost function to be
A low value of g_{n} indicates good alignment of the hidden label maps between the template and subject n, suggesting good prediction of the hidden label.
We experimented with a prior that encourages spatially constant weights and template, but did not find that the regularization lead to improvements in the localization results. In particular, we considered the following smoothness regularization on the registration parameters depending on whether we are optimizing for the weights λ_{i} or the template T:
A possible reason for this lack of improvement is that the reregistration after every few line searches already helps to regularize against bad parameter values. Another possible reason is that the above regularization assumes a smooth variation in the relationship between structure and function, which may not be true in reality. Unfortunately, the relationship between macroanatomical structure and function is poorly understood, making it difficult to design a more useful regularization that could potentially improve the results. In the experiments that follow, we will discard the regularization term of the registration parameters (i.e., set Reg(w) = 0). We also note that Reg(w) is typically set to 0 in machine learning approaches of model selection by optimization of crossvalidation error [33], [43], [58].
To optimize the task performance G over the set of parameters w, we have to instantiate the taskoptimal gradient specified in (2.16). We first compute the derivative of the taskspecific cost function with respect to the optimal update *. Once again, we represent * as the collection $\left\{{\stackrel{\u20d7}{\upsilon}}_{i}^{*}\right\}$, where ${\stackrel{\u20d7}{\upsilon}}_{i}^{*}$ is a velocity vector at vertex x_{i}. Defining $\nabla {L}_{n}{({\mathrm{\Gamma}}_{n}^{*}({x}_{i}))}^{T}$ to be the 1 × 3 spatial gradient of the warped distance transform of the nth subject L_{n} $({\mathrm{\Gamma}}_{n}^{*}(\xb7))$ at x_{i}, we get the 1 × 3 derivative
Weight Update: To update the weights {λ_{j}} of the wSSD, we compute the derivative of the registration local minimum update * with respect to the weights. Using the approximation in (2.14), we obtain the 3 × 1 derivative of the velocity update with respect to the weights of the wSSD cost function
Here δ (k, i) if k = i and is zero otherwise. Since (3.10) is in the same direction as the first term of the gradient descent direction of registration [negative of (3.2)], increasing ${\lambda}_{k}^{2}$ will improve the intensity matching of vertex x_{k} of the template. Substituting (3.10) and (3.6) into (2.16) provides the gradient for updating the weights of the wSSD cost function.
Template Update: To update the template image T used for registration, we compute the 3 × 1 derivative of the velocity update with respect to the template T
Since (3.14) is in the same direction as the first term of the gradient descent direction of registration [negative of (3.2)], when T (x_{k}) is larger than I_{n} $({\mathrm{\Gamma}}_{n}^{*}({x}_{k}))$, increasing the value T (x_{k}) of will warp vertex x_{k} of the template further along the direction of increasing intensity in the subject image. Conversely, if T (x_{k}) is smaller than I_{n} $({\mathrm{\Gamma}}_{n}^{*}({x}_{k}))$, decreasing the value of T (x_{k}) will warp vertex x_{k} of the template further along the direction of decreasing intensity in the subject image. Substituting (3.14) and (3.6) into (2.16) provides the gradient for updating the template used for registration. Note that the template subject’s hidden labels are considered fixed in template space and are not modified during training.
We can in principle optimize both the weights {λ_{i}} and the template T. However, in practice, we find that this does not lead to better localization, possibly because of too many degreesoffreedom, suggesting the need to design better regularization of the parameters. A second reason might come from the fact that we are only using an approximate gradient rather than the true gradient for gradient descent. Previous work [82] has shown that while using an approximate gradient can lead to reasonable solutions, using the exact gradient can lead to substantially better local minima. Computing the exact gradient is a challenge in our framework. We leave exploration of efficient means of computing better approximations of the gradient to future work.
We now present experiments on localizing BAs and fMRIdefined MT+ (V5) using macroanatomical cortical folding in two different data sets. For both experiments, we compare the framework with using uniform weights [31] and FreeSurfer [19].
We consider the problem of localizing BAs in the surface representations of the cortex using only cortical folding patterns. In this study, ten human brains were analyzed histologically postmortem using the techniques described in [57] and [84]. 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 a freely available suite of tools [21]. Six manually labeled BA maps (V1, V2, BA2, BA44, BA45, MT) 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. Finally, the left and right hemispheres of each subject were mapped onto a spherical coordinate system [19]. The BAs on the resulting cortical representations for two subjects are shown in Fig. 2(b). We do not consider BA4a, BA4p, and BA6 in this paper because they were not histologically mapped by the experts in two of the ten subjects in this particular data set (even though they exist in all human brains).
As illustrated in Fig. 2(c) and discussed in multiple studies [3], [4], [18], we note that V1, V2, and BA2 are wellpredicted by local cortical geometry, while BA44, BA45, and MT are not. For all the BAs however, a spherical morph of cortical folding was shown to improve their localization compared with only Talairach or nonlinear spatial normalization in the Euclidean 3D space [18]. Even though each subject has multiple BAs, we focus on each structure independently. This allows for an easier interpretation of the estimated parameters, such as the optimal template example we provide in Section IVA3. A clear future direction is to learn a registration cost function that is jointly optimal for localizing multiple cytoarchitectural or functional areas.
We compare the following algorithms.
We note that both the taskoptimal and uniformweights methods use a pairwise registration framework, while FreeSurfer uses an atlasbased registration framework. Under the atlasbased framework, all the ex vivo subjects are registered to an atlas (Fig. 4). To use the BA of a training subject to predict a test subject, we have to compose the deformations of the training subject to the atlas with the inverse deformation of the test subject to the atlas. Despite this additional source of error from composing two warps, it has been shown that with carefully constructed atlases, using the atlasbased strategy leads to better registration because of the removal of template bias in the pairwise registration framework [6], [23], [26], [31], [32], [39], [79].
We run the taskoptimal and uniformweights methods on a lowresolution subdivided icosahedron mesh containing 2562 vertices, whereas FreeSurfer results were computed on highresolution meshes of more than 100k vertices. In our implementation, training on eight subjects takes on average 4 h on a standard PC (AMD Opteron, 2GHz, 4GB RAM). Despite the use of the lowresolution mesh, we achieve stateoftheart localization accuracy. We also emphasize that while training is computationally intensive, registration of a new subject only requires one minute of processing time since we are working with lowresolution meshes.
Fig. 5 displays the mean and standard errors from the 90 trials of leavetwoout. On average, taskoptimal template performs the best, followed by taskoptimal weights. Permutation tests show that taskoptimal template outperforms FreeSurfer in five of the six areas, while taskoptimal weights outperforms FreeSurfer in four of the six areas after corrections for multiple comparisons (see Fig. 5 for more details). For the Broca’s areas (BA44 and BA45) and MT, this is not surprising. Since local geometry poorly predicts these regions, by taking into account the final goal of aligning BAs instead of blindly aligning the cortical folds, our method achieves better BA localization. FreeSurfer and the uniformweights method have similar performance because a better alignment of the cortical folds on a finer resolution mesh does not necessary improve the alignment of these areas.
Since local cortical geometry is predictive of V1, V2, and BA2, we expect the advantages of our framework to vanish. Surprisingly, as shown in Fig. 6, taskoptimal template again achieve significant improvement in BAs alignment over the uniformweights method and FreeSurfer. Taskoptimal weights is also significantly better than the uniformweights method, but only slightly better than FreeSurfer. Permutation tests show that taskoptimal template outperforms FreeSurfer in five of the six areas, while taskoptimal weights is outperforms FreeSurfer in three of the six areas after corrections for multiple comparisons (see Fig. 6 for more details). This suggests that even when local geometry is predictive of the hidden labels and anatomybased registration achieves reasonable localization of the labels, tuning the registration cost function can further improve the task performance. We also note that in this case, FreeSurfer performs better than the uniformweights method on average. Since local cortical folds are predictive of these areas, aligning cortical folds on a higher resolution mesh yields more precise alignment of the cortical geometry and of the BAs.
We note that the FreeSurfer Buckner40 atlas utilizes 40 in vivo subjects consisting of 21 males and 19 females of a widerange of age. Of these, 30 are healthy subjects whose ages range from 19 to 87. 10 of the subjects are Alzheimer’s patients with age ranging from 71 to 86. The average age of the group is 56 (see [12] for more details). The T1weighted scans were acquired on a 1.5T Vision system (Siemens, Erlangen Germany), with the following protocol: two sagittal acquisitions, FOV = 224, matrix = 256 × 256, resolution = 1 × 1 × 1.25 mm, TR = 9.7 ms, TE = 4 ms, Flip angle = 10°, TI = 20 ms and TD = 200 ms. Two acquisitions were averaged together to increase the contrasttonoise ratio. The histological data set includes five male and five female subjects, with age ranging from 37 to 85 years old. The subjects had no previous history of neurologic or psychiatric diseases (see [4] for more details). The T1weighted scans of the subjects were obtained on a 1.5T system (Siemens, Erlangen, Germany) with the following protocol: flip angle 40°, TR = 4 ms, TE = 5 ms and resolution = 1 × 1 × 1.17 mm. While there are demographic and scanning differences between the in vivo and ex vivo data sets, the performance differences between FreeSurfer and the taskoptimal framework cannot be solely attributed to this difference. In particular, we have shown in previous work that FreeSurfer’s results are worse when we use an ex vivo atlas for registering ex vivo subjects ([81,Table III]). Furthermore, FreeSurfer’s results are comparable with that of the uniformweights baseline algorithm, as well as previously published results [18], where we have checked for gross anatomical misregistration. We emphasize that since the goal is to optimize Brodmann area localization, the learning algorithm might take into account the idiosyncrasies of the registration algorithm in addition to the relationship between macroanatomy and cytoarchitecture. Consequently, it is possible that the performance differences are partly a result of our algorithm learning a registration cost function with better local minima, thus avoiding possible misregistration of anatomy.
Fig. 7 illustrates representative localization of the BAs for FreeSurfer and taskoptimal template. We note that the taskoptimal boundaries (red) tend to be in better visual agreement with the ground truth (yellow) boundaries, such as the right hemisphere BA44 and BA45.
Fig. 8 illustrates an example of learning a taskoptimal template for localizing BA2. Fig. 8(a) shows the cortical geometry of a test subject together with its BA2. In this subject, the central sulcus is more prominent than the postcentral sulcus. Fig. 8(b) shows the initial cortical geometry of a template subject with its corresponding BA2 in black outline. In this particular subject, the postcentral sulcus is more prominent than the central sulcus. Consequently, in the uniformweights method, the central sulcus of the test subject is incorrectly mapped to the postcentral sulcus of the template, so that BA2 is misregistered. Fig. 8(b) also shows the BA2 of the test subject (green) overlaid on the cortical geometry of the template subject after registration to the initial template geometry. During taskoptimal training, our method interrupts the geometry of the postcentral sulcus in the template because the uninterrupted postcentral sulcus in the template is inconsistent with localizing BA2 in the training subjects. The final template is shown in Fig. 8(c). We see that the BA2 of the subject (green) and the taskoptimal template (black) are wellaligned, although there still exists localization error in the superior end of BA2.
In the next section, we turn our attention to a fMRI data set. Since the taskoptimal template performed better than the taskoptimal weights, we will focus on the comparison between the taskoptimal template and FreeSurfer.
We now consider the application of localizing fMRIdefined functional areas in the cortex using only cortical folding patterns. Here, we focus on the socalled MT+ area localized in 42 in vivo subjects using fMRI. The MT+ area defined functionally is thought to include primarily the cytoarchitectonicallydefined MT and a small part of the medial superior temporal (MST) area (hence the name MT+). The imaging paradigm involved subjects viewing an alternating 16 s blocks of moving and stationary concentric circles. The structural scans were processed using the FreeSurfer pipeline [21], resulting in spherically parameterized cortical surfaces [11], [19]. The functional data were analyzed using the general linear model [22]. The resulting activation maps were thresholded by drawing the activation boundary centered around the vertex with maximal activation. The threshold was varied across subjects in order to maintain a relatively fixed ROI area of about 120 mm^{2} (±5%) as suggested in [68]. The subjects consist of 10 females and 32 males, with age ranging from 21 to 58 years old. 23 of the 42 subjects are clinically diagnosed with schizophrenia, while the other 19 subjects are healthy controls. Imaging took place on a 3T MR scanner (Siemens Trio) with echoplanar (EP) imaging capability. Subjects underwent two conventional highresolution 3D structural scans, constituting a spoiled GRASS (SPGR) sequence (128 sagittal slices, 1.33 mm thickness, TR = 2530 ms, TE = 3.39 ms, Flip angle = 7°, voxel size = 1.3 × 1 × 1.3 mm). Each functional run lasted 224 s during which T2*weighted echoplanar (EP) images were acquired (33 × 3mmthick slices, 3 × 3 × 3 mm voxel size) using a gradient echo (GR) sequence (TR = 2000 ms; TE = 30 ms; Flip angle = 90°). To maximize training data, no distinction is made between the healthy controls and schizophrenia patients.
In this experiment, we use each of the 10 ex vivo subjects as a template and the remaining nine subjects for training a taskoptimal template for localizing MT. We then register each taskoptimal template to each of the 42 in vivo subjects and use the template subject’s MT to predict that of the test subjects’ MT+. The results are 420 Hausdorff distances for each hemisphere. For FreeSurfer, we align the 42 in vivo subjects to the Buckner40 atlas. Once registered in this space, we can use MT of the ex vivo subjects to predict MT+ of the in vivo subjects.
Fig. 9 reports the mean and standard errors of the Hausdorff distances for both methods on both hemispheres. Once again, we find that the taskoptimal template significantly outperforms the FreeSurfer template (p < 10^{−5} for both hemispheres). We note that the errors in the in vivo subjects (Fig. 9) are significantly worse than those in the ex vivo subjects (Fig. 5). This is not surprising since functionally defined MT+ is slightly different from cytoarchitectonically defined MT. Furthermore, the ex vivo surfaces tend to be noisier and less smooth than those acquired from in vivo subjects [81]. Since our framework attempts to leverage domain specific knowledge about MT from the ex vivo data, one would expect these mismatches between the data sets to be highly deterimental to our framework. Instead, FreeSurfer appears to suffer more than our framework.
To understand the effects of the training set size on localization accuracy, we perform crossvalidation within the fMRI data set. For each randomly selected template subject, we consider 9, 19, or 29 training subjects. The resulting taskoptimal template is used to register and localize MT+ in the remaining 32, 22, or 12 test subjects, respectively. The crossvalidation trials were repeated 100, 200, and 300 times, respectively, resulting in a total of 3200, 4400, and 3600 Hausdorff distances. This constitutes thousands of hours of computation time. For FreeSurfer, we perform a pairwise prediction of MT+ among the in vivo subjects after registration to the Buckner40 atlas, resulting in 1722 Hausdorff distances per hemisphere.
Fig. 10 reports the mean and standard errors of the Hausdorff distances for FreeSurfer and taskoptimal template on both hemispheres. We see that the FreeSurfer alignment errors are now commensurate with the ex vivo results (Fig. 5). However, the taskoptimal template still outperforms FreeSurfer (p < 10^{−5} for all cases).We also note that the accuracy of MT+ localization improves with the size of the training set. The resulting localization error with a training set of 29 subjects is less than 7 mm for both hemispheres. For all training set sizes, the localization errors are also better than the ex vivo MT experiment (Fig. 5).
The experiments in the previous section demonstrate the feasibility of learning registration cost functions with thousands of degreesoffreedom from training data. We find that the learned registration cost functions generalize well to unseen test subjects of the same (Sections IVA and IVB2), as well as different imaging modality (Section IVB1). The almost linear improvement with increasing training subjects in the fMRIdefined MT+ experiment (Fig. 10) suggests that further improvements can be achieved (in particular in the histological data set) with a larger training set. Unfortunately, histological data over a whole human hemisphere is difficult to obtain, while fMRI localization experiments tend to focus on single functional areas. Therefore, a future direction of research is to combine histological and functional information obtained from different subjects and imaging modalities during training.
Since our measure of localization accuracy uses the mean Hausdorff distance, ideally we should incorporate it into our taskspecific objective function instead of the SSD on the distance transform representing the BA. Unfortunately, the resulting derivative is difficult to compute. Furthermore, the gradient will be zero everywhere except at the BA boundaries, causing the optimization to proceed slowly. On the other hand, it is unclear how aligning the distance transform values far from the boundary helps to align the boundary. Since distance transform values far away from the boundary are larger, they can dominate the taskspecific objective function g. Consequently, we utilize the distance transform over the entire surface to compute the gradient, but only consider the distance transform within the boundary of the template BA when evaluating the task performance criterion g.
The idea of using multiple atlases for segmentation has gained recent popularity [29], [49], [50], [53], [55], [76]. While we have focused on building a single optimal template, our method can complement the multiatlas approach. For example, one could simply fuse the results of multiple individuallyoptimal templates for image segmentation. A more ambitious task would be to optimize for multiple jointlyoptimal templates for segmentation.
In this work, we select one of the training subjects as the template subject and use the remaining subjects for training. The taskspecific cost function g evaluates the localization of the hidden labels via the template subject. During training (either for learning the weights or template in the registration cost function), the Brodmann areas of the template subject are held constant. Because the fixed Brodmann areas are specific to the template subject, the geometry of the template subject should in fact be the best and most natural initialization. It does not make sense to use the geometry of another subject (or average geometry of the training subjects) as initialization for the template subject’s Brodmann areas, especially since the geometry of this other subject (or average geometry) is not registered to the geometry of the template subject. However, the use of a single subject’s Brodmann (or functional) area can bias the learning process. An alternative groupwise approach modifies the taskspecific cost function g to minimize the variance of the distance transforms across training subjects after registration. In this case, both the template geometry and Brodmann (functional) area are estimated from all the training subjects and dynamically updated at each iteration of the algorithm. The average geometry of the training subjects provided a reasonable template initialization. However, our initial experiments in the ex vivo data set do not suggest an improvement in task performance over the pairwise formulation in this paper.
While this paper focuses mostly on localization of hidden labels, different instantiations of the taskspecific cost function can lead to other applications. For example, in group analysis, the taskspecific cost function could maximize differences between diseased and control groups, while minimizing intragroup differences, similar to a recent idea proposed for discriminative Procrustes alignment [38].
In this paper, we present a framework for optimizing the parameters of any smooth family of registration cost functions, such as the image dissimilarityregularization tradeoff, with respect to a specific task. The only requirement is that the task performance can be evaluated by a smooth cost function on an available training data set. We demonstrate stateoftheart localization of Brodmann areas and fMRIdefined functional regions by optimizing the weights of the wSSD imagesimilarity measure and estimating an optimal cortical folding template. We believe this work presents an important step towards the automatic selection of parameters in image registration. The generality of the framework also suggests potential applications to other problems in science and engineering formulated as optimization problems.
The authors would like to thank P. Parillo for discussion on the optimization aspects of this paper. The authors would also like to thank C. Brun, S. Durrleman, T. Fletcher, and W. Mio for helpful feedback on this work.
This work was supported in part by the NAMIC (NIH NIBIB NAMIC U54EB005149), in part by the NAC (NIH NCRR NAC P41RR13218), in part by the mBIRN (NIH NCRR mBIRN U24RR021382), in part by the NIH NINDS R01NS051826 Grant, in part by the NSF CAREER 0642971 Grant, in part by the National Institute on Aging (AG02238), in part by the NCRR (P41RR14075, R01 RR1659401A1), in part by the NIBIB (R01 EB001550, R01EB006758), in part by the NINDS (R01 NS05258501), and in part by the MIND Institute. Additional support was provided by The Autism & Dyslexia Project funded by the Ellison Medical Foundation. The work of B. T. Thomas Yeo was supported by the A*STAR, Singapore.
In this appendix, we prove Proposition 2: If the Hessian ${\partial}_{\upsilon}^{2}f({w}_{0},{\mathrm{\Gamma}}^{*}({w}_{0})\circ \delta \mathrm{\Gamma}(\upsilon ))$ is positive definite at υ = 0, then there exists an ε > 0, such that for all δw, ‖δw‖ < ε, a unique continuous function υ*(w_{0}, δw) exists, such that υ*(w_{0}, δw) is the new local minimum for parameters w_{0} + δw and υ*(w_{0}, 0) = 0. Furthermore, υ*(w_{0}, δw) has the same order of smoothness as f.
In the next section, we first prove that the Hessian ${\partial}_{{\upsilon}_{1}}^{2}f({w}_{0},{\mathrm{\Gamma}}^{*}({w}_{0})\circ \delta \mathrm{\Gamma}({\upsilon}_{1})){}_{{\upsilon}_{1}}=0$ is equal to the mixderivatives matrix ${\partial}_{{\upsilon}_{1},{\upsilon}_{2}}^{2}f({w}_{0},{\mathrm{\Gamma}}^{*}({w}_{0})\circ \delta \mathrm{\Gamma}({\upsilon}_{1})\circ \delta \mathrm{\Gamma}({\upsilon}_{2})){}_{{\upsilon}_{1}={\upsilon}_{2}=0}$ under the composition of diffeomorphisms model [75], [81], [82]. We then complete the proof of Proposition 2.
We will only provide the proof for when the image is defined in ^{3} so as not to obscur the main ideas behind the proof. To extend the proof to a manifold (e.g., S^{2}), one simply need to extend the notations and bookkeeping by the local parameterizing the velocity fields υ_{1} and υ_{2} using coordinate charts. The same proof follows.
Let us define some notations. Suppose the image and there are M voxels. Let $\stackrel{\u20d7}{x}$ be the ^{3M} rasterized coordinates of the M voxels. For conciseness, we define for the fixed parameters w_{0}
Therefore, p is a function from ^{3M} to . Under the composition of diffeomorphisms model, δΓ (υ) is the diffeomorphism parameterized by the stationary velocity field υ defined on the M voxels, so that δΓ (υ)(·) is a function from ^{3M} to . To make the dependence of δΓ (υ) on υ explicit, we define
and so is a function from ^{3M} × ^{3M} to ^{3M}. In other words, we can rewrite
and
Now that we have gotten the notations out of the way, we will now show that
Hessian: We first compute the 1× 3M Jacobian via the chain rule
From the above equation, we can equivalently write down the jth component of the 1 × 3M Jacobian
where ^{n} and ${\upsilon}_{1}^{j}$ denote the nth and j th components of and υ_{1}, respectively. Now, we compute the (i, j) th component of the 3M × 3M Hessian using the product rule
Because _{υ1} _{υ1=0} is the identity matrix and the 1 × 3M Jacobian _{υ1} p((υ_{1}, $\stackrel{\u20d7}{x}$))_{υ1=0} = (_{}p)(_{υ1} )_{υ1=0} = 0 (because derivative is zero at local minimum), we get _{}p _{υ1=0} = 0, and so the second term in (A.10) is zero.
To simplify the first term of (A.10), we once again use the fact that _{υ1} _{υ1=0} is the identity matrix, and so the summand is zero unless k = i and n = j. Consequently, (A.10) simplifies to
or equivalently
MixDerivatives Matrix: We first compute the 1 × 3M Jacobian via the chain rule
From the above equation, we can equivalently write down the jth component of the 1 × 3M Jacobian
Now, we compute the (i, j)th component of the 3M × 3M mixderivatives matrix using the product rule
Like before, we have _{}p_{υ1=υ2=0} = 0, and so the second term is zero. Because _{υ1} _{υ1=0} is the identity, ${\partial}_{{\upsilon}_{1}^{i}}{\Upsilon}_{1}^{n}$ is zero unless k = i. Since ${\Upsilon}_{1}^{n}({\upsilon}_{1}=0,\stackrel{\u20d7}{x})=\stackrel{\u20d7}{x},{\partial}_{{\stackrel{\u20d7}{x}}^{j}}{\Upsilon}_{1}^{n}$, is also equal to zero unless n = j. Therefore, we get
or equivalently
We now complete the proof of Proposition2. Let h(w, υ_{1}) _{υ2} f (w, Γ*(w_{0}) δΓ (υ_{1}) δΓ (υ_{2}))_{υ2=0}. Since δΓ(0) = Id, we have
where the last equality comes from the definition of Γ*(w_{0}) being a local minimum for the composition model.
Since the mixderivatives matrix _{υ1} h(w, υ_{1})_{υ1=0} is invertible by the positivedefinite assumption of this proposition, by the Implicit Function Theorem, there exists an ε > 0, such that for all δw, ‖δw‖ < ε, there is a unique continuous function υ* (w_{0}, δw), such that h(w_{0} + δw, υ* (w_{0}, δw)) = 0 and υ* (w_{0}, 0) = 0. Furthermore, υ* (w_{0}, δw) has the same order of smoothness as f.
Let $k(w,{\upsilon}_{1})={\partial}_{{\upsilon}_{2}}^{2}f(w,{\mathrm{\Gamma}}^{*}({w}_{0})\circ \delta \mathrm{\Gamma}({\upsilon}_{1})\circ \delta \mathrm{\Gamma}({\upsilon}_{2})){}_{{\upsilon}_{2}=0}$. Then k (w_{0}, 0) is positive definite at υ_{1} = 0 by the assumption of the proposition. By the smoothness of derivatives and continuity of eigenvalues, there exists a small neighborhood around (w_{0}, υ_{1} = 0) in which the eigenvalues of k (w, υ_{1}) are all greater than zero. Therefore, Γ*(w_{0}) δΓ (υ*(w_{0}, δ_{w})) does indeed define a new local minimum close to Γ*(w_{0}).
To compute _{w} υ*, we perform a Taylor expansion
and rearranging the terms for υ_{1} = υ*, we get
^{1}Here, we assume that the transformation Γ is finite dimensional, such as the parameters of affine transformations, positions of spline control points or dense displacement fields defined on the voxels or vertices of the image domain.
B. T. Thomas Yeo, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA (Email: ude.tim.liasc@samohty).
Mert R. Sabuncu, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA. Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA 02129 USA (Email: ude.tim.liasc@ucnubasm)
Tom Vercauteren, Mauna Kea Technologies, 75010 Paris, France (Email: moc.hcetaekanuam@neretuacrev.mot)
Daphne J. Holt, Massachusetts General Hospital Psychiatry Department, Harvard Medical School, Charlestown, MA 02139 USA (Email: gro.srentrap@tlohd)
Katrin Amunts, Department of Psychiatry and Psychotherapy, RWTH Aachen University and the Institute of Neuroscience and Medicine, Research Center Jülich, 52425 Jülich, Germany (Email: ed.nehcaaku@stnumak).
Karl Zilles, Institute of Neuroscience and Medicine, Research Center Jülich and the C.&O. VogtInstitute for Brain Research, University of Düsseldorf, 52425 Jülich, Germany (Email: ed.hcileujzf@selliz.k)
Polina Golland, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA (Email: ude.tim.liasc@anilop)
Bruce Fischl, Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA. Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA 02129 USA. Department of Radiology, Harvard Medical School and the Divison of Health Sciences and Technology, Massachusetts Institute of Technology, Cambridge, MA 02139 USA (Email: ude.dravrah.hgm.rmn@lhcsif)
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. 