|Home | About | Journals | Submit | Contact Us | Français|
One of the most challenging problems in modern neuroimaging is detailed characterization of neurodegeneration. Quantifying spatial and longitudinal atrophy patterns is an important component of this process. These spatiotemporal signals will aid in discriminating between related diseases, such as frontotemporal dementia (FTD) and Alzheimer’s disease (AD), which manifest themselves in the same at-risk population. Here, we develop a novel symmetric image normalization method (SyN) for maximizing the cross-correlation within the space of diffeomorphic maps and provide the Euler-Lagrange equations necessary for this optimization. We then turn to a careful evaluation of our method. Our evaluation uses gold standard, human cortical segmentation to contrast SyN’s performance with a related elastic method and with the standard ITK implementation of Thirion’s Demons algorithm. The new method compares favorably with both approaches, in particular when the distance between the template brain and the target brain is large. We then report the correlation of volumes gained by algorithmic cortical labelings of FTD and control subjects with those gained by the manual rater. This comparison shows that, of the three methods tested, SyN’s volume measurements are the most strongly correlated with volume measurements gained by expert labeling. This study indicates that SyN, with cross-correlation, is a reliable method for normalizing and making anatomical measurements in volumetric MRI of patients and at-risk elderly individuals.
Frontotemporal dementia (FTD) prevalence may be higher than previously thought and may rival Alzheimer’s disease (AD) in individuals younger than 65 years . Because FTD can be challenging to detect clinically, it is important to identify an objective method to support a clinical diagnosis. MRI studies of individual patients are difficult to interpret because of the wide range of acceptable, age-related atrophy in an older cohort susceptible to dementia. This has prompted MRI studies that look at both the rate and the anatomic distribution of change [2–7].
Manual, expert delineation of image structures enables in vivo quantification of focal disease effects and serves as the basis for important studies of neurodegeneration . Expert structural measurements from images also provide the gold-standard of anatomical evaluation. The manual approach remains, however, severely limited by the complexity of labeling 2563 or more voxels. Such labor is both time consuming and expensive to support, while the number of individual experts available for such tasks is limited. A third significant difficulty is the problem of inter-rater variability which limits the reliability of manual labeling . While rarely available for large-scale data processing, an expert eye remains valuable for limited labeling tasks that give a basis for algorithmic evaluation.
Deformable image registration algorithms are capable of functioning effectively in time-sensitive clinical applications  and high-throughput environments and are used successfully for automated labeling and measurement research tasks. One challenge is reliable performance on non-standard data, as in studies of potentially severe neurodegenerative disorders. These types of images violate the basic assumptions of small deformations and/or simple intensity relationships used in many existing image registration methods.
Diffeomorphic image registration algorithms hold the promise of being able to deal successfully with both small [10–16] and large deformation problems [17–26]. State of the art methods also give full space-time optimizations, are symmetric with respect to image inputs and allow probabilistic similarity measures . Inverse consistent image registration (ICIR) is an additional common alternative to diffeomorphic mapping. Inverse consistency was first introduced by Thirion as an extension to his Demons algorithm  but was popularized by Johnson and Christensen  and others . Symmetric methods are distinct from ICIR in that symmetric algorithms, first, guarantee that results are identical regardless of the order of input data and, second, use exact inverse transformations guaranteed by diffeomorphisms. Inverse consistency approximates symmetry by including variational penalties in the normalization optimization algorithm. Depending on the weights of the various data, regularization and inverse consistency terms, consistency may be satisfied (or not) at the expense of the other matching criterion. Furthermore, inverse consistent algorithms use approximate inverse transformations . Because the inverse transformations themselves are approximate, the consistency term, as well, is compromised.
Here, we develop a novel symmetric diffeomorphic optimizer for maximizing the cross-correlation in the space of topology preserving maps. The cross-correlation measure has been used in medical image registration before [10,11,31] and more extensively in computer vision. However, this measure has not been investigated for the diffeomorphic case. Furthermore, it has not been used in symmetric normalization or “inverse consistent” studies. Applying our novel normalization formulation to cross-correlation provides the advantage (or option) of symmetrizing the cross-correlation Euler-Lagrange equations. We show that these symmetric Euler-Lagrange equations can be computed with only minor additional computational costs. We then give a careful evaluation of the performance of our symmetric diffeomorphic algorithm for high dimensional normalization of elderly and neurodegenerative cortical anatomy. We compare the method to an elastic cross-correlation optimizer as well as the Demons algorithm which was shown to outperform other methods in a careful evaluation of inter-subject brain registration .
Thirion’s Demons algorithm  is known to perform well in inter-subject deformable image registration. The method uses an approximate elastic regularizer to solve an optical flow problem, where the “moving” image’s level sets are brought into correspondence with those of a reference or “fixed” template image. In practice, the algorithm computes an optical flow term which is added to the total displacement (initially zero). The total displacement is then smoothed with a Gaussian filter. The process repeats for a set number of iterations for each resolution in a multi-resolution optimization scheme. The method is freely available in the Insight ToolKit and has been optimized by the ITK community (www.itk.org).
Dawant et al. used the Demons algorithm for segmenting the caudate nucleus, the brain and the cerebellum for a morphometric comparison of normal and chronic alcoholic individuals . Their evaluation of the algorithm found reasonable agreement between automated and manual labeling. They also showed results on the automated labeling of hippocampus but did not evaluate performance. Their comparison used the Dice statistic (overlap ratio),
which measures both difference in size and location between two segmentations, R1 and R2. The #(R) operator counts the number of pixels in the region, R. This sensitive measure varies in the range [0, 1] where values greater than 0.6 for smaller structures and 0.8 for larger structures are considered good by some authors [34,8].
The range of acceptable values for the Dice statistic are, of course, highly dependent upon the application. Both the amount of certainty that one has in the “gold standard” dataset and, also, the specific use of the segmentations determine a reasonable operating range. Our goal, in this paper, is to use manually segmented structures as a foundation for comparing automated normalization methods. In this respect, it is relative performance (measured with respect to the Dice statistic) that is of critical value.
A diffeomorphism is a differentiable map with a differentiable inverse [35,36]. We typically restrict our solutions to the diffeomorphic space Diff0 with homogeneous boundary conditions. That is, we assume that rigid and scaling transformations have been factored out and the image border maps to itself1. Shortest paths between elements in this space are termed geodesic. Diffeomorphic methods were introduced into medical computer vision  for the purpose of providing a group theoretical, large deformation space-time image registration framework. Current developments in large deformation computational anatomy by Miller, Trouve and Younes extended the methods to include photometric variation and to use Euler-Lagrange equations . However, the standard version of these methods, Beg’s Large Deformation Diffeomorphic Metric Matching (LDDMM) , do not formulate the transformation optimization symmetrically. They are only symmetric in theory and their implementation requires parallel computation. Through personal communication, however, we understand that a symmetrization operator, based on the transformation Jacobian, is being included in current developments by Younes and Trouv’e [37,38]. However, these methods do not guarantee symmetry for similarity metrics other than the intensity difference.
Our current work extends the Lagrangian diffeomorphic registration technique described in . This new formulation has symmetry properties required for a geodesic connecting two images, I and J, in the space of diffeomorphic transformations and guarantees symmetry regardless of the chosen similarity measure. This formulation accounts for the natural symmetry in the problem: both images move along the shape (diffeomorphism) manifold. Symmetric diffeomorphisms guarantee two properties that are intrinsic to the notion of a geodesic path: the path from I to J is the same as it is when computed from J to I, regardless of similarity metric or optimization parameters. Symmetry is required for distance estimates and makes results independent of arbitrary decisions about which image is “fixed” or “moving.”
Our method is also unique in that it guarantees sub-pixel accurate, invertible transformations in the discrete domain by directly including invertibility constraints in the optimization. While diffeomorphisms are theoretically guaranteed to be invertible, interpolation errors can lead to invertibility errors that increase linearly with the number of interpolation steps. Our solution, on the other hand, directly minimizes this error by exploiting the invertibility guaranteed by diffeomorphisms. Finally, the method is efficient enough to use on single-processor machines and in processing large datasets.
We define a diffeomorphism of domain Ω, generally, for transforming image I into a new coordinate system by I = I ○ (x, t = 1). The parameters of these transformations are time, t, a spatial coordinate, x, and a velocity field, υ(x, t) on Ω, which is a square-integrable, continuous vector field . The correspondence maps, , are gained by integrating the velocity fields in time, ; the distance is then , where L defines the linear operator regularizing the velocity. The functional norm, || · ||L, induces regularity on the velocity field via a linear differential operator such as L = a2 + bId (a, b constants). Such a diffeomorphism gives a dense map in both space and time and is illustrated in figure 1. A basic fact of diffeomorphisms allows them to be decomposed into two parts, 1 and 2. We exploit this fact to define a variational energy that explicitly divides the image registration diffeomorphisms into two halves such that I and J contribute equally to the path and deformation is divided between them. Assume that x indicates the identity position of some anatomy in image I and z indexes the identity position of the same anatomy in image J. We assume, also, that the diffeomorphism maps homologous anatomy in these images. The prior knowledge that this diffeomorphic map should apply evenly to both images can be captured by including the constraint D(Id, 1(x, 0.5)) = D(Id, 2(z, 0.5)) directly in the formulation of the problem. The result is a method that finds correspondences with equal consideration of both images. Note that below we will derive the equations assuming intensity difference as a similarity measure, for simplicity. However, in actuality, we have a variety of statistical image similarity measures (robust intensity difference, cross-correlation, mutual information) at our disposal, as in , or employ user landmarks as in . After this introductory section, we will develop our method for the cross-correlation.
Define the image registration optimization time, t [0, 1] where t indexes both 1 and 2, though in opposite directions. The similarity seeks 1 such that 1(x, 1)I = J. Recalling the basic definition of diffeomorphisms allows us to write any geodesic through composing two parts, each of which is a geodesic (any sub-part of a geodesic is a geodesic). We apply this fact in the second step of our derivation below,
which converts the similarity term from |1(x, 1)I − J| to |1(x, t)I − 2(z, 1 − t)J|2. A visualization of these components of is in Figure 1. The forward and backward optimization problem is then, solving to time t = 0.5,
Minimization with respect to 1 and 2 provides the symmetric normalization (SyN) solution and also solves a 2-mean problem. Landmarks may also be included, as in the Lagrangian Push Forward method , by dividing the similarity term, as done with the image match terms above. The constraint D(1(x, 0.5)) = D(2(x, 0.5)) is built into the fact that we integrate the solution from 0 to 0.5. Because the velocities are of (approximately) constant arc length and are, at each iteration, of exactly the same length, the length of 1, integrated over [0, 0.5], is equivalent to the length of 2 integrated over [0, 0.5].
As noted in the introduction, this method is quite distinct from inverse consistent image registration (ICIR) . ICIR uses vector fields, h and g, to define correspondence from I to J and J to I, respectively. Therefore, in total, ICIR uses four vector fields in its approach to normalization, h, g, h−1 and g−1. The inverses of these fields are not guaranteed to exist (as the optimization is not performed in the diffeomorphic space) and no exact method is used to compute the inverses. Instead, judging from the brief discussion in , the inverse is only guaranteed to be exact at a few points in the domain and the inversion algorithm itself is not well-specified or exact. This inexactness means that the variational term measuring the difference between two vector fields, ||h(x)−g−1(x)||2, only coarsely estimates consistency. Furthermore, this means that ICIR must include two terms to compute consistency, that is, both ||h(x)−g−1(x)||2 and ||h−1(x)−g(x)||2. If inverses were exact (or very close to exact) only one of the above terms would be enough as minimizing one would imply minimizing the other.
SyN, alternatively, provides an inverse that is guaranteed to be everywhere sub-pixel accurate. Furthermore, SyN avoids a basic computational redundancy present in ICIR - the solutions h(x) and g(x) overlap in time. SyN’s pair of solutions, 1 and 2, on the other hand, do not overlap in time. They are only two parts of a longer path. This distinction is shown in figure 2. However, SyN is still able to compute the measure of inverse consistency by simply composing our time 0.5 maps together and comparing 1(x, 1) and in the I domain (a similar computation may be done in the J domain). Recall that, by definition and diffeomorphic computations, we are able to guarantee and . We therefore have an inverse consistency error of zero, up to the inaccuracy induced by the single interpolation required in computing the time 1 maps, e.g. . A key to the symmetry of our method is the ability to compute sub-pixel accurate inverses such that, for each i = 1, 2, we have (to a user-selected numerical precision) and . Typical precision is 1.0 × 10−6 L1 norm and 0.2 L∞ norm, measured in terms of voxel/pixel coordinates. Further details on the numerical methods employed in optimizing this energy may be found in  and will be summarized below.
We now propose to symmetrically solve the following image matching problem: Find a spatiotemporal mapping, Diff0, such that the cross correlation between the image pair is maximized. This formulation of the problem represents a change of philosophy when compared to Bajcsy  and, later, Gee’s  elastic matching, which also used cross-correlation. Elastic image registration methods seek to balance a regularization term and a similarity term allowing one to find a constrained deformable solution. The approach recommended here departs from this strategy by allowing an unconstrained optimization of the similarity term within the space of diffeomorphisms. This strategy is used when finding optimal mappings in lower-dimensional (e.g. affine) transformation groups. The main advantage of an unconstrained search within the space of diffeomorphisms is its simplicity: one allows the method to maximize the similarity until a local maximum or limit on computation time is reached. 2 The disadvantage is one requires a similarity metric that, when optimized in Diff, provides a useful solution. Therefore, design and choice of the similarity metric requires great care.
The mutual information (MI) similarity metric garnered significant interest in recent years [42,43,14,44]. Intuitively speaking, MI estimates the optimal matching between images by inferring how much global information is shared in the image pair, as estimated from the pair’s joint histogram. At the same time, MI prevents over-fitting by penalizing clustering of the marginal image probabilities . The globality of this approach makes it extremely useful for robust rigid registration but may limit performance in deformable registration, in particular in cases where non-stationary noise patterns or intensity inhomogeneity requires a locally adaptive similarity. This problem has been addressed by Studholme et al  and in our previous work . One difficulty with a locally varying estimate to the MI is that joint probabilities need a large number of samples for reliable statistics. Therefore, as locality in the MI estimate increases, its statistical reliability decreases.
Cross-correlation (CC), on the other hand, adapts naturally to situations where locally varying intensities occur and is suitable for some multi-modality problems. The CC depends only on estimates of the local image average and variance which may be accurately/exactly measured with relatively few samples. Furthermore, the cross-correlation has shown historically to perform well in many real-world computer vision applications where one requires robustness to unpredictable illumination, reflectance, etc. An example of the robustness of our method to strong MRI inhomogeneity (bias field) is shown in figure 3. For these reasons, we revisit the classical cross-correlation as a similarity metric for use in our emerging diffeomorphic image registration.
We now provide a symmetric formulation of the diffeomorphic image registration problem as driven by CC, along with the Euler-Lagrange equations for this problem. First, define I1 = I(1(x, 0.5)) and J2 = J (2(x, 0.5)). We also define a variable to represent each image with its local mean subtracted as Ī(x) = I1(x) − μI1 (x) and (x) = J2 − μJ2 (x). We compute μ over a local nD window centered at each position x, where D is the image dimension. We usually choose n = 5. The cross-correlation is then
where the inner product is also taken over a nD window. Note that our CC is implicitly a function of 1 and 2 through its dependence on I1 and J2. We now are able to define the variational optimization problem in similar fashion to equation 3,
The problem, here, is the same as before but with the cross-correlation as the driving similarity term. We now take the variation of this function with respect to 1 at time 0.5 and 2 at time 0.5. Following Beg’s derivation, adapted for our symmetric normalization and particular similarity metric, we find two Euler-Lagrange equations,
The gradients given in equations 6 and 7 include that fact that the velocities exist in the space of smooth vector fields given by the linear operator L as well as the determinant of each i transformation Jacobian, |Di|. This equation is similar to both that derived for LDDMM  and the derivative for the cross-correlation given by Hermosillo . This equation differs from Beg’s Euler-Lagrange equation in that we have an E-L equation for both 1 and 2 instead of just 1. In addition, Beg used the intensity difference metric given in equation 3. The derivation of this equation differs from that in Hermosillo  by the presence of the derivative of the diffeomorphism Ī/i, which leads to the Jacobian term, instead of a vector field (small deformation) derivative. Furthermore, we have represented the CC in a different arrangement of terms that suits our own derivation of the CC variation and computational implementation. This arrangement of terms allows one to precompute, for each iteration, the locally varying values of A, B, C and store them for use in the local computation of the derivative.
Our novel symmetric formulation, therefore, does not add significant additional cost to the normalization. The main additional cost is that of smoothing the derivative estimate to gain the υi, not in the actual estimate of the similarity term. Both the 1 and 2 derivatives require the terms A, B and C. Therefore, in estimating equation 6 and equation 7, we use the procedure detailed in algorithm 1 below.
Algorithm 1 : Computing Cross-Correlation Derivatives
(1) Deform I by 1(0.5) and J by 2(0.5).
(2) Calculate Ī and from the result of step (1).
(3) Calculate and store images representing A, B and C.
These steps enable us to loop over the image domain to rapidly compute equation 6 and equation 7 at the same time from these precomputed variables. Note, however, that it could be possible to modify Lewis’s fast normalized correlation measure  to speed up this process even further.
Equations 6 and 7 gives the update to the velocities at time 0.5 and consequently the diffeomorphisms 1 and 2. However, we also require the i to satisfy our o.d.e. and exact invertibility constraints. We satisfy these constraints by, given a new velocity estimate, first updating the i through the o.d.e. and then integrating backward in time to find i −1 and verify invertibility. This approach is a standard one and detailed, algorithmically, in the LPF method . For completeness, we include an abbreviated explanation here.
First, assume arbitrary and υ, related by the standard o.d.e. d(x, t)/dt = υ((x, t), t) used to generate diffeomorphisms. The update method for our diffeomorphism comes from discretization of this o.d.e., such that,
This discretization is used to update both i from time 0 to 0.5. Second, algorithm 2 gives a method for generating inverses when an arbitrary diffeomorphism is updated by a small time-step through a velocity field, as in the previous equation. The algorithm typically converges within one to a few iterations, particularly if the time step is small. The existence of a solution is guaranteed by the integrability condition established for diffeomorphic image registration , while uniqueness comes from the uniqueness theorem of o.d.e.s .
Algorithm 2 : Inversion Method
The algorithm uses a fixed point method to push the inverse of forward by a small amount performing a gradient descent on −1() = Id, enforcing −1() = Id to a sub-pixel level. The same approach was used in the Lagrangian Push Forward algorithm . We use a temporary variable, ψ, to represent the input di3eomorphism that, on output, will be the inverse of . Input (x) = y, ψ−1(ỹ) = x where ỹ ≠ y and output ψ−1(y) = x = −1(y). At convergence, ỹ = y and ||ψ−1() − Id||∞ < ε2r where r is the image resolution and ε2 a small constant, typically 0.1.
1: while ||ψ−1((x)) − x||∞ > ε2r do
2: Compute ν−1(x) = ψ−1((x)) − x.
3: Find scalar γ such that ||ν−1||∞ = 0.5r.
4: Integrate ψ−1 s.t. ψ−1(ỹ, t)+ = γ ν−1(ψ−1(ỹ, t)).
5: end while
We now summarize the SyN method in algorithm 3. Our implementation essentially updates each i with the current estimate to the velocity, follows the update to each i by calling algorithm 2 to generate inverse maps, and then re-estimates the velocity from the new estimate to the i. This formulation and diffeomorphic representation guarantee SyN’s sub-pixel invertibility and algorithmic independence to input permutations. These methods allow us to symmetrically match images to the degree that discrete diffeomorphisms are invertible. Furthermore, while extended here to cross-correlation, similar techniques may be used to efficiently symmetrize almost any other similarity measure. We now leverage ITK for a fair implementation and evaluation of SyN.
Algorithm 3 : Overview of Symmetric Normalization Method
(1) Initialize and .
(2) Repeat the following steps until convergence:
(3) Compute the cross correlation as described in algorithm 1.
(4) Compute each νi by smoothing the result of step (3) in this table. One may also use the modified midpoint method for each velocity, as in the LPF algorithm , to give smoothness in time.
(5) Update each i by νi through the o.d.e. as described in equation 8. This step automatically adjusts the time step-size such that the maximum length of the updates to the i is sub-pixel and approximately constant over iterations. We explicitly guarantee ||ν1(·, t) || = ||ν2(·,t) ||. We also update the estimate to the geodesic distance by trapezoidal rule, as in the LPF method.
(6) Use algorithm 2 to get the inverses of the i.
(7) Generate the time 1 solutions from and .
The Demons algorithm is freely available in the standard ITK distribution and has been quantitatively evaluated by the ITK community. We have implemented SyN within our extended version of the ITK deformable image registration framework, described in . Therefore, we are in a position to measure performance gains by varying first the similarity metric and then the the transformation model, keeping an identical underlying code base.
For this part of the study, we replace the itkDemonsRegistrationFunction, in the itkDemonsRegistrationFilter, with our own itkCrossCorrelationRegistrationFunction. The interface and operation of this modified deformable registration algorithm are identical to the Demons algorithm, except that the image forces come from the small deformation version of equation 6, as may be found in . In addition to switching the itkDemonsRegistrationFunction for our cross-correlation image forces, we also modified the time step used during this elastic registration process. The time-step was increased from 1, for Demons, to 4 for cross-correlation. The Demons force, based on optical flow, is extremely aggressive compared to our analytical cross-correlation derivative. We found, for our dataset, that time-steps beyond 4 led to non-positive Jacobians in a significant fraction of our normalizations.
As mentioned previously, we have implemented SyN in our extended Insight ToolKit. For this study, SyN will use the itkPDEDeformableRegistrationFilter as implemented in ITK for its Gaussian smoothing operator 3. This smoothing operator will be applied only to the incremental output of our itkCrossCorrelationRegistrationFunction. The Demons and elastic cross-correlation algorithms both apply Gaussian smoothing to the total displacement field, in accordance with small deformation transformation models. The presence of a large deformation assumption, along with the use of a space-time parameterization of the transformation, are the main (not the only) differences between diffeomorphic and elastic registration algorithms.
The key contrast between the first pair of methods (Demons, elastic cross-correlation) is the similarity metric. The only difference between the second two (elastic cross-correlation, SyN cross-correlation) is in the transformation model. We are therefore investigating, first, if cross-correlation provides better normalization than Demons optical flow, for our neurodegenerative and elderly T1 MRI dataset. The second comparison evaluates if changing from the elastic to the symmetric diffeomorphic transformation model yields additional improvement. Our hypothesis was that each change would yield benefits. This hypothesis was born out in our experimental data, as suggested by the preliminary study data shown in figure 4.
Our database consists of 20 T1 MRI images (0.85 × 0.85 × 1.5 mm, GE Horizon Echospeed 1.5 T scanner) from 10 normal elderly and 10 frontotemporal dementia patients. The 10 frontotemporal dementia individuals are a different set than used in our previous study . Each of the 20 images, along with the elderly template, was manually labeled with the protocol described in . This protocol was shown to be highly reproducible for both small and large structures via six-month intra-rater reliability and inter-rater reliability measurements. Left hippocampus labeling, for example, showed a 0.92 intra-rater overlap ratio (equation 1) and 0.83 average for inter-rater overlap. As the hippocampus is relatively small, these values are reasonable. Finally, as the labeling was intended to be only on the cortex (not cerebrospinal fluid), we masked the label images by the inverse of the cerebrospinal fluid segmentations. This also makes our evaluation more sensitive to differences in the accuracy of sulcal and gyral alignment.
We now use these two comparative pairings of three algorithms to study performance for characterizing the volumetric differences between elderly and frontotemporal dementia cortex, hippocampus, amygdala and cerebellum. An example comparison of the SyN and Demons methods’ normalization abilities is in figure 4. This study reveals the ability of these methods, SyN, elastic cross-correlation and Demons, to reproduce results gained from an expert user’s labeling of our 20 image dataset. This test is performed by using each method to map the template labels (an elderly individual also labeled by our expert user) to each individual. Twenty rigid registrations and sixty non-rigid registrations were required (twenty subjects were each deformably registered by three algorithms). The atlas labelings are then warped by the composed rigid and non-rigid transformation into the space of the patient image. We then compute Dice overlap ratios between the manual and automatic structural segmentations for each structure.
Each method was applied in an identical four-level multi-resolution scheme and ran until convergence or a fixed (maximum) number of iterations was reached. We allowed up to 100 iterations at the first level, 100 iterations at the second level, 100 iterations at the third level and 20 iterations at the full resolution. The relative running times for each algorithm varied depending upon the particular cases being run. However, the relative average running times (in terms of the average Demons running time) were 1 (Demons), 4.2 (elastic cross-correlation) and 5.5 (SyN). The mean runtime of the Demons method on a 2.0 GHz Intel Mobile Pentium processor was 20.4 minutes. The relative similarity of the elastic time cost to the cost of SyN also indicates that our symmetric approach does not create significant additional computational cost.
A visualization of six individuals from our dataset and the normalization of the template to these individuals is shown in figure 5. Associated intensity error images are shown in figure 6. Label error images are shown in figure 7. The eight labeled individual structures are shown, along with the summary statistics for our results, in figure 8. Both cross-correlation algorithms produced segmentation results above the minimum threshold of 0.6 for all structures, as shown in figure 8. We also compared the minimum Jacobian of the elastic cross-correlation and SyN cross-correlation methods and did not find significant differences (T= −1.67725, P < 0.101703). This indicates that SyN’s results are not significantly less constrained at a local level. The smallest difference in performance between the Demons and elastic cross-correlation methods was on the cerebellum. The highest difference was on the frontal lobe. The largest performance difference in SyN and elastic cross correlation was on the temporal and frontal lobe. Performance gains expectedly focus on frontal and temporal lobe due to the known effect of FTD on these two areas. The shape changes caused by this complex disease are likely difficult to capture with a constrained method such as elastic normalization. The smallest differences between the elastic and diffeomorphic methods were found in the hippocampus and amygdala. This is likely to indicate that the similarity metric does not provide a rich enough feature space over which to optimize the correspondence of these structures. Therefore, a better optimizer operating with that similarity metric will not be likely to provide an advantage.
A likely improvement in performance would come from using a method such as the STAPLE algorithm  to bootstrap results. Similarly, an optimal template would augment results . However, both of these enhancements would increase performance in a consistent manner across all our tested registration methods and would be unlikely to change the relative performance of our algorithms.
The second test we performed was to evaluate whether the volume measurements obtained by the normalization are strongly correlated with the measurements of the expert user. We performed this analysis only on regions (frontal, temporal, parietal lobes) where one may reasonably expect FTD to induce a difference from the normal Elderly population and where labeling performance was good. The purpose of this test is to determine if conclusions made by analyzing the output of our automated normalization methods are at all consistent with the expert user’s analysis. As we are using the expert labeling as our gold standard, the “better” method should produce volume measures that correlate more strongly with the volume measures gained by the expert. The volume, in cubic centimeters, of each structure was calculated by summing the voxel volumes that were given the appropriate label for that structure. As shown in table 1, SyN outperformed elastic cross-correlation (and the Demons method) according to the degree to which the automated volumes correlate with the manual volumes for each of three structures. SyN also had the least overall discrepancy in the measured volume, that is, Σi |VAlg − VMan|, where VMan is manually measured volume and VAlg refers to algorithmically measured volume. These results are shown in table 2.
We also plot the estimated volume versus the manually computed volume for the temporal lobe in figure 9, for frontal lobe in figure 10 and for parietal lobe in figure 11. We use linear regression to fit the estimated volume with a line in order to assess the closeness of the method’s estimate to the manual estimate. The Pearson correlation of SyN volumes with the manually measured volumes was 0.86 for temporal lobe, 0.89 for frontal lobe and 0.71 for parietal lobe. The correlation of the elastic method volumes with the manually measured volumes was 0.69 for temporal lobe, 0.67 for frontal lobe and 0.42 for parietal lobe. The correlation of Demons volumes with the manually measured volumes was 0.79 for temporal lobe, 0.67 for frontal lobe and 0.66 for parietal lobe. One interesting finding, here, is that the Demons method, compared to the elastic method, produces stronger correlations with manual labelings (in terms of volume measurements) but produces smaller overlap ratios.
Our final analysis tests the significance of the difference in volumes between the FTD and elderly members of our image population. The P-values of these results, obtained by non-parametric permutation testing, are shown in table 3. None of the structures showed a significant volumetric difference between groups. However, SyN appears to reflect the significance obtained by the manual rater more closely than the other methods. At the same time, we do not find these results to be strong enough to warrant definitive conclusions. Comparing a larger dataset (or possibly a different set of control and subject images) would likely show stronger differences than this dataset, as we found in our previous, preliminary study .
Thus, we can see how an apparently small, but consistent difference in performance (as measured by overlap ratio) can have an impact on the validity of the study outcome. That is, the performance gains in overlap ratio translate to more accuracy in making clinically meaningful measurements, such as volume. Although the quality of our results are reasonable, by some standards, these algorithms, as applied here, cannot claim to accurately reproduce manual labelings. The reasons for this are well-known. The primary reason is that expert knowledge is not directly encoded in these methods. Secondarily, the uncertainty inherent to the problem of neuroanatomical labeling limits the accuracy and reproducibility of both manual and automated segmentation. Finally, it is not yet known the extent to which brains of different individuals, when represented as magnetic resonance images, are diffeomorphic to each other. This problem is even less well understood with elderly and patient brains. Note also that, when measuring the trends in the difference of elderly and FTD structure volumes, the Demons, elastic correlation and SyN significance tend to estimate larger volumes than the manual results. This is likely caused by two things: segmentation bias towards the template and the fact that registration-based segmentations are smoothed (elastic and Demons more than SyN), while the manual segmentations are not. Both of these types of errors may be observed in figure 7. Figure 7 also shows that the labeler may have used different “styles” of labeling or changed the decision-making process over the dataset. All of these confounds impact the overlap ratios and correlations that we find here. We accentuate that our goal, here, is not specifically to reproduce the expert labels, but to use these labeled structures as a neuroanatomically-based evaluation of relative algorithmic performance.
We first described the symmetric normalization formulation. We then extended this formulation to use the cross-correlation similarity function providing, in addition, Euler-Lagrange equations for the variational problem in the symmetric diffeomorphic context. We contrasted our method with the popular inverse consistent image registration technique, which was outperformed by the Demons method in an unbiased comparative evaluation of brain segmentation and alignment . We then provided algorithmic details of the SyN methodology. Finally, we leveraged our ITK implementation to compare the SyN cross-correlation optimizer with the Demons and an elastic cross-correlation optimizer. This enabled us to explore the effect of similarity metric and optimizer separately and purely, as we use an identical linear operator for regularization and an identical code base, the Insight ToolKit. The relative similarity of the elastic computation time to the computational cost of SyN also indicates that our symmetric approach does not require significant additional computational expense. We found, in short, that the cross-correlation behaves well on elderly and neurodegenerative data. In particular, it outperforms the Demons optical flow on our dataset. We also found that the symmetric diffeomorphic optimizer outperforms the elastic optimizer on the same dataset, while guaranteeing that the topology of the images is preserved and, importantly, that the algorithm’s performance does not depend upon the order in which one inputs the images.
This careful comparison establishes the distinct advantage of SyN for segmenting elderly and neurodegenerative cerebrum, cerebellum and, in particular, the temporal and frontal lobe. Note that, in addition to better performance, SyN provides a dense space-time map and transformation inverses. The differences in performance are consistent, statistically significant and have a major impact on study outcome. One can extrapolate even larger differences between SyN and algorithms with lower dimensionality than either Demons or SyN. For this reason, along with the theoretical advantages that translate into practical benefits, we promote diffeomorphic algorithms and the cross-correlation similarity metric in neuroimaging research, in particular when studying nonstandard datasets, such as FTD and AD.
We thank the reviewers for greatly improving the contents of this paper. Much of this work was supported by NIH grant R01-EB006266.
1Note that extending the background space of an image allows almost any diffeomorphism of the image to be captured in Diff0.
2The contribution of the regularization term to the total energy is small compared to the similarity. However, the regularization term in the diffeomorphic matching problem remains important. It guarantees that one finds a path of minimal length.
3The Gaussian kernel, Gσ, is an estimate to the Green’s kernel for the linear operator L = 2 +Id  and is adequate for providing the velocity field regularity necessary for integrating o.d.e.s.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.