|Home | About | Journals | Submit | Contact Us | Français|
Measures of brain changes can be computed from sequential MRI scans, providing valuable information on disease progression for neuroscientific studies and clinical trials. Tensor-based morphometry (TBM) creates maps of these brain changes, visualizing the 3D profile and rates of tissue growth or atrophy. In this paper, we examine the power of different nonrigid registration models to detect changes in TBM, and their stability when no real changes are present. Specifically, we investigate an asymmetric version of a recently proposed unbiased registration method, using mutual information as the matching criterion. We compare matching functionals (sum of squared differences and mutual information), as well as large deformation registration schemes (viscous fluid and inverse-consistent linear elastic registration methods versus symmetric and asymmetric unbiased registration) for detecting changes in serial MRI scans of 10 elderly normal subjects and 10 patients with Alzheimer's Disease scanned at 2-week and 1-year intervals. We also analyzed registration results when matching images corrupted with artificial noise. We demonstrated that the unbiased methods, both symmetric and asymmetric, have higher reproducibility. The unbiased methods were also less likely to detect changes in the absence of any real physiological change. Moreover, they measured biological deformations more accurately by penalizing bias in the corresponding statistical maps.
In recent years, computational anatomy has become an exciting interdisciplinary field with many applications in functional and anatomic brain mapping, image-guided surgery, and multimodality image fusion (Avants and Gee (2004); Christensen et al. (1996); Chung et al. (2001); Collins et al. (1994); Grenander and Miller (1998); Miller (2004); Shen and Davatzikos (2003); Thompson and Toga (2002)). The goal of image registration is to align, or spatially normalize, one image to another. In multi-subject studies, this reduces subject-specific anatomic differences by deforming individual images onto a population average brain template. When applied to serial scans of human brain, image registration offers tremendous power in detecting the earliest signs of illness, understanding normal brain development or aging, and monitoring disease progression (Camara et al. (2008); Crum et al. (2001); Hogan et al. (2000); Klein et al. (2009); Pieperhoff et al. (2008); Ridha et al. (2006); Wang et al. (2003)). Recently, there has been an expanding literature on various nonrigid registration techniques, with different image matching functionals, regularization schemes, and numerical implementations. Many algorithms were developed that regularize different differential operators, including elastic regularization (Broit (1981)), viscous fluid registration (Christensen et al. (1996)), and large deformation diffeomorphic metric mapping (Beg et al. (2005)), among other works. We have argued previously that most work in tensor-based morphometry (TBM) is interested in relative volume gain or loss, and if that is to be statistically evaluated, it is preferable to be working on distributions that have no bias and minimal skew, in order to obtain correct interpretations.
It is also important to observe that the Unbiased technique, though a novel concept in medical image registration, can be adapted and combined with any non-rigid registration algorithm. As such, we chose to add unbiased registration to work with fluid regularization in order to show its advantages in measuring biological deformations more accurately, by penalizing bias in the corresponding statistical maps and generating more stable and reproducible results than if only fluid regularization were used.
In (Leow et al. (2007); Yanovsky et al. (2007b)), our group systematically examined the statistical properties of Jacobian maps (the determinant of the local Jacobian operator applied to the deformations), and proposed an unbiased large-deformation image registration approach. In this context, unbiased means that we strive to obtain a zero-mean and symmetric log-Jacobian distribution under the null hypothesis, when a pair of images is matched. We argued that this distribution is beneficial when recovering change in regions of homogeneous intensity, and in ensuring symmetrical results when the order of two images being registered is switched. An unbiased algorithm is advantageous as it does not detect gain with more likelihood than loss when signals of equal magnitude log-Jacobian but opposite sign are present. We applied this method to a longitudinal MRI dataset from a single subject, and showed promising results in eliminating spurious signals. We also noticed that different registration techniques, when applied to the same longitudinal dataset, may sometimes yield visually very different Jacobian maps, causing problems in interpreting local structural changes. Given this ambiguity and the increasing use of registration methods to measure brain change, more information is required concerning the baseline stability, reproducibility, and statistical properties of signals generated by different nonrigid registration techniques.
The main contribution of this paper is to provide quality calibrations for different non-rigid registration techniques in tensor-based morphometry (TBM). We systematically investigate and compare the performances of different non-rigid registration techniques including two common matching functionals: L2, or the sum of squared intensity differences, versus mutual information, and four regularization techniques (fluid registration, inverse-consistent linear elastic registration, and the Symmetric and Asymmetric Unbiased registration techniques). Our experiments are designed to determine which registration method is more reproducible and more reliable with less artifactual variability, especially in regions of homogeneous image intensity. Also, for the first time we investigate the Asymmetric version of the Unbiased registration (by contrast with the Symmetric Unbiased model), as well as analyze unbiased models with mutual information based matching functionals (prior work has focused on the case where the summed squared intensity differences is used as the criterion for registration).
Following analyses in the preparatory phase of the Alzheimer's Disease Neuroimaging Initiative (ADNI) (Leow et al. (2006)), the foundation of our calibrations is based on the assumption that, by scanning healthy normal human subjects twice over a 2-week period using the same protocol, serial MRI scan pairs should not show any systematic biological change. Therefore, any regional structural differences detected using TBM over such a short interval may be assumed to be errors. We apply statistical analysis to the profile of these errors, providing information on the reliability, reproducibility and variability of different registration techniques. Moreover, serial images of 10 subjects from the ADNI follow-up phase (images acquired one year apart) were analyzed in a similar fashion and compared to the ADNI baseline data. In images collected one year apart, real anatomical changes are present; neurobiological changes due to aging and dementia include widespread cell shrinkage, regional gray and white matter atrophy and expansion of fluid-filled spaces in the brain. Thus, a good computational technique should be able to differentiate between longitudinal image pairs collected for the ADNI baseline (2-week) and follow-up (1-year) phases. We refer to prior papers for details of the ADNI acquisition protocol, but briefly, all subjects were scanned with a standardized MRI protocol, developed after a major effort evaluating and comparing 3D T1-weighted sequences for morphometric analyses (Jack et al. (2008)).
In the experiments that follow, all scans were collected according to the standard ADNI MRI protocol (http://www.loni.ucla.edu/ADNI/Research/Cores/), which acquires a high-resolution sagittal T1-weighted 3D MP-RAGE sequence for each subject, with a reconstructed voxel size of 0.9375 × 0.9375 × 1.2 mm3. Additional image corrections were also applied, using a processing pipeline at the Mayo Clinic, consisting of: (1) a procedure termed GradWarp for correction of geometric distortion due to gradient non-linearity (Jovicich et al. (2006)), (2) a “B1-correction”, to adjust for image intensity non-uniformity using B1 calibration scans (Jack et al. (2008)), (3) “N3” bias field correction, for reducing intensity inhomogeneity Sled et al. (1998), and (4) geometrical scaling, according to a phantom scan acquired for each subject (Jack et al. (2008)), to adjust for scanner- and session-specific calibration errors. Additional phantom-based geometric corrections were applied to ensure spatial calibration was kept within a specific tolerance level for each scanner involved in the ADNI study (Gunter et al. (2006)).
At this point, we would like to motivate the unbiased approach, which couples the computation of deformations with statistical analyses on the resulting Jacobian maps. As a result, the unbiased approach ensures that deformations have intuitive axiomatic properties by penalizing any bias in the corresponding statistical maps. In the following sections, we describe the mathematical foundations of this approach, define energy functionals for minimization, and perform thorough statistical analyses to demonstrate the advantages of the unbiased registration models.
We first introduce the notation used in this paper. Throughout this paper, we denote the vectors by bold fonts and scalars by regular fonts. Let Ω be an open and bounded domain in n, for arbitrary n. Without loss of generality, assume that the volume of Ω is 1, i.e. |Ω| = 1. Let I1 : Ω → and I2 : Ω → be the two images to be registered. We seek to estimate a transformation g : Ω → Ω such that I2 matches I1 when deformed by g. In this paper, we will restrict this mapping to be differentiable, one-to-one, and onto. We denote the Jacobian matrix of a deformation g to be Dg. The inverse mapping of g is denoted by g−1.
The displacement field u(x) from the position x in the deformed image I2○g(x) back to I2(x) is defined in terms of the deformation g(x) by the expression g(x) = x − u(x) at every point x Ω. Thus, we consider the problems of finding g and u to be equivalent. It is sometimes more convenient to write expressions in terms of either g or u. For instance, we can denote the determinant of the Jacobian matrix of deformation g as either |Dg(x)| or |D(x − u(x))|.
We now describe the construction of the Unbiased Large-Deformation Image Registration. We associate three probability density functions (PDFs) to g, g−1, and the identity mapping id:
By associating deformations with their corresponding global density maps, we can now apply information theory to quantify the magnitude of deformations (Yanovsky et al. (2007a)). In our approach, we choose the Kullback-Leibler (KL) divergence and symmetric Kullback-Leibler (SKL) distance. The KL divergence between two probability density functions, p1(x) and p2(x), is defined as
We define the SKL distance as
The Unbiased method solves for the deformation g (or, equivalently, for the displacement u) minimizing the energy functional E, consisting of the image matching term F and the regularizing term R which is based on KL divergence or SKL distance. The fidelity term F dependents on I2 and I1, as well as the displacement u. The general minimization problem can be written as
Here, λ > 0 is a weighting parameter.
To quantify the magnitude of deformation g, in this paper we introduce a new regularization term RKL, which is an asymmetric measure between Pid and Pg:
This regularization term can be shown to be
Thus, the energy functional in (4) implementing Asymmetric Unbiased registration can be written as
for some distance measure F between I2(x − u) and I1(x). The second term on the right-hand side of the equality in (7) is equivalent to the logarithmic barrier in numerical optimization theory (Nocedal and Wright (1999)) and is well-behaved.
In this section, we describe the regularization functional based on the symmetric KL distance between Pid and Pg:
The energy functional employing Symmetric Unbiased registration can be rewritten as
for some distance measure F. Notice that the symmetric unbiased regularizing functional is pointwise nonnegative, while the asymmetric unbiased regularizer in equation (6) can take either positive or negative values locally. Although in theory, the asymmetric KL regularization may potentially favor voxel expansion over the identity transformation (at least locally), this is not the case globally. Indeed, given a body force of zero everywhere, the deformation with minimum asymmetric KL energy is still the identity transformation. This can be readily appreciated by noticing that the function − log(x), though not strictly non-negative, is still a convex function with respect to its argument, x. Thus, expansion in some regions would induce contraction in others, driving the overall asymmetric KL energy upwards and away from zero. Moreover, although symmetrization was shown to be important in (Christensen and Johnson (2001)), here we show that in practice, the asymmetric unbiased method does not seem to perform much differently than its symmetric version. This further supports our conclusion that the log transformation may be a more fundamental operation than symmetrization in understanding Jacobian maps in the context of medical imaging.
In this paper, the matching functional F takes two forms: the L2 norm (the sum of squared intensity differences) and MI (mutual information). These functionals have each been widely used in the past for nonrigid registration, to measure the intensity agreement between a deforming image and the target image. We briefly describe both distances in this section.
The L2-norm matching functional is suitable when the images have been acquired through similar sensors and thus are expected to present the same intensity range and distribution. The L2 distance between the deformed image I2(x − u) and target image I1(x) is defined as
Computing the first variation of functional FL2 gives the following gradient
Mutual information is a measure of how much information one random variable has about another. The use of mutual information for image registration was first introduced in (Collignon et al. (1995) and Viola and Wells (1995)). One of the main advantages of using mutual information is that it can be used to align images of different modalities, without requiring knowledge of the relationship (joint intensity distribution) of the two registered images. We refer the readers to (D'Agostino et al. (2003); Hermosillo et al. (2001); Wells et al. (1996)) for relevant discussions on mutual information.
To define the mutual information between the deformed image I2(x − u) and the target image I1(x), we follow the notations in (Hermosillo et al. (2001)), where pI1 and are used to denote the intensity distributions estimated from I1(x) and I2(x − u), respectively. An estimate of their joint intensity distribution is denoted as . In this probabilistic framework, the link between two modalities is fully characterized by a joint density.
We let i1 = I1(x), i2 = I2(x − u(x)) denote intensity values at point x Ω. Given the displacement field u, the mutual information computed from I1 and I2 is provided by
We seek to maximize the mutual information between I2(x − u) and I1(x), or equivalently, minimize the negative of :
where |Ω| is the volume of Ω, Qu is defined as
and ψ(ξ1, ξ2) is a two-dimensional Parzen windowing kernel, which is used to estimate the joint intensity distribution from I2(x − u) and I1(x). Here, we use the Gaussian kernel with variance σ2:
In general, we expect minimizers of the energy functional E(u) to exist. Computing the first variation of the functional in (4), we obtain the gradient of E(I1, I2, u), namely uE(I1, I2, u). We define the force field f, which drives I2 into registration with I1, as
Here, R(u) is either RKL(u) or RSKL(u). Explicit expressions for components of uR(u), in both cases, are derived in Appendices A and B for two and three dimensional cases, respectively. Also, the gradient uF(I1, I2, u) depends on the choice of the fidelity term.
Given the force field, the most straightforward way to minimize (4) might seem to involve parameterizing the descent direction by an artificial time τ,
However, in our case, we do not solve Euler-Lagrange equations using the gradient descent method. In order to regularize the flow, we employ the fluid regularization proposed in (Christensen et al. (1996)). Given the velocity field v, the following partial differential equation can be solved to obtain the displacement field u:
The instantaneous velocity as in (D'Agostino et al. (2003)) is obtained by convolving f with Gaussian kernel Gσ of variance σ2:
This equation can be solved efficiently using the Fast Fourier transform (FFT).
To avoid possible confusion, we summarize the methods we will be referring to in our subsequent analyses. In later discussions, minimization of the following energies
using equations (18), (21), (20) will be referred to as L2-Asymmetric Unbiased and L2-Symmetric Unbiased models, respectively. The model above, provided λ = 0, will be referred to as the L2-Fluid model.
Similarly, minimization of
will be referred to as the MI-Asymmetric Unbiased and MI-Symmetric Unbiased models, respectively. Such models, with λ = 0, define the MI-Fluid model.
|Algorithm 1 Unbiased Nonlinear Registration|
Based on the authors' approach in (Leow et al. (2007); Yanovsky et al. (2007b)), we observe that, given that there is no systematic structural change within two weeks, any deviation of the Jacobian map from one should be considered error. Thus, we expect that a better registration technique would yield log |Dg| values closer to 0 (i.e., smaller log Jacobian deviation translates into better methodology). Mathematically speaking, one way to test the performance is to consider the deviation map dev of the logged (i.e., logarithmically transformed) Jacobian away from zero, defined at each voxel as
For two different registration methods A and B, we define the voxel-wise deviation gain of A over B (denoted by SA,B) as
For the ADNI baseline dataset (in which patients are scanned twice with MRI, two weeks apart), two distinct types of t tests are used, a within-subject paired t test and a group paired t test. A within-subject paired t test is conducted for each subject by pooling all voxels inside a region of interest, as defined by the ICBM whole brain mask (the ICBM brain is a standardized population average image, defined by the International Consortium for Brain Mapping (Mazziotta et al. (2001))). This determines whether two methods differ significantly inside the whole brain (for each subject). A group paired t test, on the other hand, is performed across subjects, by computing a voxel-wise t-map of deviation gains. In this case, to statistically compare the performance of two registration methods, we rely on the standard t test on the voxel mean of S. To construct a suitable null hypothesis, we notice that the following relation would hold, assuming B outperforms A
Thus, the null hypothesis in this case would be testing if the mean deviation gain is zero
To determine the ranking of A and B, we have to consider one-sided alternative hypotheses. For example, when testing if B outperforms A, we use the following alternative hypothesis
For n subjects, the voxel-wise T statistic, defined as
thus follows the Student's t distribution under the null hypothesis and may be used to determine the p-value that the null hypothesis is true. If the alternative hypothesis is accepted, we confirm that sequence B outperforms A at point x. Otherwise, we would rank A and B equally if the null hypothesis is not rejected.
For both the ADNI follow-up dataset (in which patients are scanned twice with MRI, one year apart) and ADNI baseline dataset, we create a voxel-wise t map using the local log Jacobian values of the ten subjects, allowing us to test the validity of the zero mean assumption. To simplify the notation, we introduce J to denote J = |Dg|. The following voxel-wise T statistic compared to a two-tailed Student's t distribution may then be used to test the above null hypothesis
We reject the null hypothesis if the p value calculated above exceeds a preset threshold based on a suitable confidence interval. Notice the voxel-wise variance of log J provides us with a way to assess the repeatability of a deformation method, i.e., measuring the voxel-wise spread of the given multiple observations (with higher variance corresponding to lower repeatability).
To determine the overall global effects of different registration methods on the deviation of log Jacobian maps throughout the brain, we performed permutation tests to adjust for multiple comparisons (Bullmore et al. (1999); Nichols and Holmes (2001)). Following the analyses in (Leow et al. (2006)), we resampled the observations by randomly flipping the sign of (i = 1, 2, …, n) under the null hypothesis. For each permutation, voxelwise t tests are computed. We then compute the percentage of voxels inside the chosen ROI (in this case the ICBM mask) with T statistics exceeding a certain threshold. The multiple comparisons corrected p value may be determined by counting the number of permutations whose above-defined percentage exceeds that of the un-permuted observed data. This is comparable to ‘set-level inference’ in the widely-used SPM (Statistical Parametric Mapping) functional image analysis package (Friston et al. (1995)). For example, we say that sequence B outperforms A on the whole brain if this corrected p value is smaller than 0.05 (that is, less than 5% of all permutations have the above-defined percentage greater than that of the original data). All possible (210 = 1024) permutations were considered in determining the final corrected p value.
To visually assess the global significance level of the voxel-wise t tests on deviation gains and log-Jacobian values, we also employed the cumulative distribution function (CDF) plot, as in several prior studies (Brun et al. (2007); Chiang et al. (2007); Lepore et al. (2007); Morra et al. (2009)). In brief, we plot observed cumulative probabilities against the theoretical distribution under the null hypothesis. These CDF plots are commonly created as an intermediate step, when using the false discovery rate (FDR) method to assign overall significance values to statistical maps (Benjamini and Hochberg (1995); Genovese et al. (2002); Storey (2002); Zhu et al. (2007)). As they show the proportion of supra-threshold voxels in a statistical map, for a range of thresholds, these CDF plots (sometimes called Q-Q plots) offer a measure of the effect size in a statistical map. They also may be used to demonstrate which methodological choices influence the effect size in a method that creates statistical maps (Brun et al. (2007); Chiang et al. (2007); Lepore et al. (2007)).
In the case of deviation gains S of a worse technique A over a better technique B in the ADNI baseline data, we expect a CDF curve to lie above the Null line, in the sense that a better technique exhibits less systematic changes. In the case of log-Jacobian values, a better registration technique, on the other hand, should be able to separate the CDF curves between ADNI baseline and follow-up phases (this is what we refer to as the separation of CDF curves in the presence of real physiological changes).
In Sections 7 and 8, we tested the Asymmetric Unbiased and Symmetric Unbiased models and compared the results to those obtained using the Fluid registration model (Christensen et al. (1996); D'Agostino et al. (2003)) and inverse-consistent linear elastic registration (Christensen and Johnson (2001); Leow et al. (2005)). Of note, even though Asymmetric Unbiased and Symmetric Unbiased methods minimize different energy functionals, our experiments showed that they generate very similar maps. For each regularization technique, we employed both L2 and mutual information matching functionals (see equations (22)-(25)).
To obtain a fair comparison, re-gridding was not employed. Re-gridding is a method to relax the energy computed from the linear elasticity prior after a certain number of iterations, which allows large-deformation mappings to be recovered without any absolute penalty on the displacement field (other than via the smoothness constraint on the velocity field which is integrated to give the displacement) (Christensen et al. (1996)). It is essentially a memory-less procedure, as how images are matched after each re-gridding is independent of the final deformation before the re-gridding, rendering the comparison of final Jacobian fields and cost functionals problematic. Moreover, we consider the strategy of re-gridding, through the relaxation of deformation fields over time, to be less rigorous from a theoretical standpoint, as the imposition of a regularizer can be used to secure distributional properties in the resulting statistics (e.g., symmetric log-Jacobian).
For the experiments in this paper, different values of parameter σ (the standard deviation) in equation (21) were chosen. The values we used were σ = 7.0, 9.0, 12.0. Also, different values of parameter λ were employed. For example, in our MI-Unbiased experiments, λ was chosen to be 1.0, 2.0, 5.0, 10.0. As will be shown in the experimental sections, Fluid registration yielded best results with σ = 9.0 and 12.0, and therefore, all fluid registration results were visualized with σ = 9.0, unless otherwise mentioned. As will be shown in the experimental sections, for all values of λ, the unbiased registration outperformed the fluid registration with statistical significance when comparing the value of the data term at convergence (for the 10 subjects used) using a t-test to compare the results of one set of parameters versus the other. Unless otherwise mentioned, MI-Unbiased registration results were generated using values of σ = 9.0 and λ = 5.0. A similar procedure was employed to choose the parameters for L2-based methods. Notice that our procedure for selecting the parameters closely follows the approach employed in (Cachier et al. (2003)), who compared algorithms over a range of regularization constants and picked the best constant for each algorithm. Alternatively, parameters may be selected as in (Yeo et al. (2008)), where the best regularization constant was found using cross-validation.
In this section, nonlinear registration was performed on a dataset that we shall refer to as the “ADNI Baseline” dataset, collected during the preparatory phase of the ADNI project, which includes serial MRI images of ten normal elderly subjects acquired two weeks apart. Each of the ten pairs of scans is represented on a 220 × 220 × 220 grid. Here, the foundation of calibrations is based on the assumption that, by scanning normal control human subjects serially within a two-week period using the same MRI protocol, no systematic structural changes should be recovered.
In our first experiment, we compared methods based on L2 matching (L2-Fluid, L2-Asymmetric Unbiased, and L2-Symmetric Unbiased). Uniform values of λ = 500 and λ = 1000 were used for all deformations using L2-Symmetric Unbiased and L2-Asymmetric Unbiased algorithms, respectively. Since the Asymmetric Unbiased model quantifies only the forward deformation, the weight of the corresponding regularization functional is half the magnitude of that of the Symmetric Unbiased model, and hence, a weighting parameter twice as large should be used.
Figures 1--44 show the results of registering a pair of serial MRI images for one of the subjects (subject 3). The deformation was computed in both directions (time 2 to time 1, and time 1 to time 2) using methods based on L2 matching. In Figure 1, Jacobian maps of deformations are superimposed on brain volumes. Both Asymmetric Unbiased and Symmetric Unbiased methods generate less noisy Jacobian maps with values closer to the identity mapping, which shows the superior stability of the Unbiased approach in the absence of physiological changes. We also visually assessed the inverse consistency of the mappings (Christensen and Johnson (2001)) by concatenating forward and backward Jacobian maps (in an ideal situation, this operation should yield the identity). Again, we observe noticeable visual differences between the results obtained using the unbiased methods and Fluid registration. Figure 2 shows the L2-norm decrease per iteration for all L2-based models. Figure 3 plots the KL divergence and SKL distance measures for each of the L2-based methods. For L2-Fluid method, both KL and SKL measures increase with increasing numbers of iterations. On the other hand, even though the Asymmetric Unbiased method minimizes the KL divergence and the Symmetric Unbiased model minimizes the SKL distance, these two measures stabilize for both unbiased methods. Figure 4 shows the histograms of voxel-wise deviation gains of L2-Fluid over L2-Asymmetric Unbiased as well as L2-Fluid over L2-Symmetric Unbiased. The histograms are skewed to the right, indicating the superiority of both unbiased registration methods over Fluid registration.
Of note, we have also considered a different deviation map, defined as , in place of (26). We performed statistical analyses with this definition of deviation gain, which yielded very similar results. These results are therefore not shown in this paper.
In Figure 5, we compared L2-Fluid and L2-Symmetric Unbiased methods, conducting a within-subject paired t test inside the ICBM mask for each of the ten subjects. In this case, p < 0.0001 for all subjects, indicating that the Symmetric Unbiased registration, when coupled with L2 matching cost functional, produces more reproducible maps with less variability.
Figure 6 shows the mean Jacobian maps obtained using L2-Fluid, L2-Asymmetric Unbiased, and L2- Symmetric Unbiased registration algorithms. Jacobian maps generated using unbiased models have values closer to 1, whereas L2-Fluid model generated noisy mean maps. Figure 7, shows the results when performing 3D voxel-wise paired t tests for the deviation gain of L2-Fluid over L2-Asymmetric Unbiased and L2-Fluid over L2-Symmetric
Unbiased. T maps for the deviation gains are empirically thresholded at 2.28 (p = 0.005 on the voxel level with 9 degrees of freedom) to show statistical significance.
Figure 8 shows results obtained using Multiple Comparison Analysis using permutation testing on deviation gains of L2-Fluid over L2-Symmetric Unbiased. The results indicate that out of 1024 permutations, no permutation yields a larger percentage of voxels with p < 0.05 than the observed data, which indicates that L2-Symmetric Unbiased method outperforms L2-Fluid with p < 0.001.
To emphasize the differences between the distributions of log Jacobian values for Fluid and unbiased (both asymmetric and symmetric) methods, in Figure 9, we plotted the cumulative distribution function of the p-values in deviation gains as defined in (27). In these plots, the interval p [0, 0.005] is the most important. For a null distribution, this cumulative plot falls along the line y = x in xy-plane, as represented by the dashed black line. Larger upward inflections of the CDF curve near the origin are associated with significant deviation gains, indicating that both Asymmetric Unbiased and Symmetric Unbiased methods outperform Fluid method in being less likely to exhibit structural changes in the absence of systematic biological changes.
In our second experiment, we compared the performance of methods based on mutual information matching (MI-Fluid, MI-Asymmetric Unbiased, and MI-Symmetric Unbiased). Uniform values of λ = 5 and λ = 10 were used for all deformations using MI-Symmetric Unbiased and MI-Asymmetric Unbiased algorithms, respectively. Since the Asymmetric Unbiased model quantifies only the forward deformation, the weight of the corresponding regularization functional is half the magnitude of that of the Symmetric Unbiased model, and hence, a weighting parameter twice as large should be used.
Figures 10--1313 show the results of registering a pair of serial MRI images for one of the subjects (subject 3). The deformation was computed in both directions (time 2 to time 1, and time 1 to time 2) using methods based on Mutual Information matching. In Figure 10, Jacobian maps of deformations are superimposed on brain volumes. Both Asymmetric Unbiased and Symmetric Unbiased methods generate less noisy Jacobian maps with values closer to the identity mapping, which shows the superior stability of the Unbiased approach in the absence of physiological changes. We also observed that MI-Asymmetric Unbiased and MI-Symmetric Unbiased methods to produce inverse consistent maps with less variability. Figure 11 shows the mutual information decrease per iteration for all MI-based models. Figure 12 plots the KL divergence and SKL distance measures for each of the MI-based methods. For MI-Fluid method, both KL and SKL measures increase with increasing numbers of iterations. On the other hand, these two measures stabilize for both unbiased methods. Figure 13 shows the histograms of voxel-wise deviation gains of MI-Fluid over MI-Asymmetric Unbiased as well as MI-Fluid over MI-Symmetric Unbiased. The histograms are skewed to the right, indicating the superiority of both unbiased registration methods over Fluid registration.
In Figure 14, we compared MI-Fluid and MI-Symmetric Unbiased methods, conducting a within-subject paired t test for each of the ten subjects. In this case, p < 0.0001 for all subjects, indicating that the Symmetric Unbiased registration, when coupled with mutual information matching cost functional, produces more reproducible maps with less variability.
Figure 15 shows the mean Jacobian maps obtained using MI-Fluid, MI-Asymmetric Unbiased, and MI- Symmetric Unbiased registration algorithms. Jacobian maps generated using unbiased models have values closer to 1, whereas MI-Fluid model generated noisy mean maps. Figure 16, shows the results when performing 3D voxel-wise paired t tests for the deviation gain of MI-Fluid over MI-Asymmetric Unbiased and MI-Fluid over MI-Symmetric Unbiased. T maps for the deviation gains are empirically thresholded at 2.28 (p = 0.005 on the voxel level with 9 degrees of freedom) to show statistical significance.
Figure 17 shows results obtained using Multiple Comparison Analysis using permutation testing on deviation gains of MI-Fluid over MI-Symmetric Unbiased. The results indicate that out of 1024 permutations, no permutation yields a larger percentage of voxels with p < 0.05 than the observed data, which indicates that MI-Symmetric Unbiased method outperforms MI-Fluid with p < 0.001.
In Figure 18, we plotted the cumulative distribution function of the p-values in deviation gains. Larger upward inflections of the CDF curve near the origin are associated with significant deviation gains, indicating that both Asymmetric Unbiased and Symmetric Unbiased methods outperform Fluid method in being less likely to exhibit structural changes in the absence of systematic biological changes.
Lastly, we compared L2 and mutual information cost functionals for both Fluid and Symmetric Unbiased regularization. (Since Asymmetric Unbiased and Symmetric Unbiased regularizations produce similar results, we do not show the results for the asymmetric version). We again conducted within-subject paired t tests (Figure 19) as well as group paired t tests (Figure 20) on the voxel-wise deviation gains for all voxels inside the ICBM brain mask. We showed that MI-Fluid outperforms L2-Fluid with p < 0.0001. However, the result of the comparison of L2-Symmetric Unbiased and MI-Symmetric Unbiased is inconclusive.
In this section, we analyze a dataset we shall call the “ADNI Follow-up” phase dataset, which includes serial MRI images (220 × 220 × 220) of ten subjects acquired one year apart. These data were collected as part of a larger study to track degenerative brain changes in MRI in 800 subjects, ages 55 to 90, including 200 elderly controls, 400 subjects with mild cognitive impairment, and 200 patients with AD. As the images are now one year apart, real anatomical changes are present, which allows methods to be compared in the presence of true biological changes.
In Figures 21 and and22,22, nonlinear registration was performed using Fluid, Asymmetric Unbiased, and Symmetric Unbiased methods. Visually, the Fluid method generates noisy mean Jacobian maps, while maps generated using unbiased methods suggest a volume reduction in gray matter as well as ventricular enlargement. Here, both Asymmetric Unbiased and Symmetric Unbiased methods perform equally well. Figure 23 displays the cumulative distribution of p-values for the voxel-wise log Jacobian t-maps for both ADNI Baseline and ADNI Follow-up datasets. We expect a better method to separate these two CDF curves, indicating that a real biological change has occurred between the two time points. A greater separation is accomplished when Asymmetric Unbiased and Symmetric Unbiased methods are used, while the Fluid method does not differentiate between the two datasets.
In Figures 24 and and25,25, we further compared different registration methods by matching images artificially corrupted with Gaussian noise (zero mean; variance 12.0). Figure 25 displays the cumulative distribution of p-values for the voxel-wise log Jacobian t-maps for both ADNI Baseline and ADNI Follow-up datasets. A greater separation is accomplished when Unbiased methods are used.
We also compared Fluid registration and Unbiased registration methods using different values for the two parameters σ and λ (Figure 26). We used λ = 1.0, 2.0, 5.0, 10.0 (for the Unbiased registration), and σ = 7,9,12 (for both Fluid and Unbiased registration). In general, Fluid registration and Unbiased registration with smaller λ values generated noisier mean maps, while maps generated using Unbiased registration with larger λ values suggest a volume reduction in gray matter as well as ventricular enlargement. As the value of the smoothing parameter σ increased, the resulting Jacobian maps became smoother, making the biological effects, such as reduction in gray matter, harder to detect.
In Figures 27 and and28,28, we compared Unbiased registration methods with the inverse-consistent linear elastic matching of (Christensen and Johnson (2001); Leow et al. (2005)). Figure 28 displays the cumulative distribution of p-values for the voxel-wise log Jacobian t-maps for both ADNI Baseline and Follow-up datasets for Unbiased technique with different values of λ, Fluid registration, and inverse-consistent linear elastic matching. Unbiased methods, especially with larger values of parameter λ, show a bigger separation between the baseline and follow-up curves than Fluid and inverse-consistent linear elasticity methods do, indicating that Unbiased methods, with any λ value, were able to better differentiate between the two datasets.
This is the first paper which aims to provide quality calibrations for different non-rigid registration techniques in TBM. We systematically investigated and compared the performances of different non-rigid registration techniques including common matching functionals (L2-metric and mutual information), and regularization techniques (fluid registration, inverse-consistent linear elastic matching, Asymmetric Unbiased, and Symmetric Unbiased techniques). Experiments were conducted to determine which registration method is more reproducible, more reliable, with less artifactual variability in regions of homogeneous image intensity. We also introduced a novel asymmetric unbiased registration model (the Asymmetric Unbiased model) and for the first time, analyzed unbiased models with mutual information based matching functionals.
We showed that both Asymmetric and Symmetric Unbiased models generate very similar results. Although in theory, the asymmetric KL regularization may potentially favor voxel expansion over the identity transformation, this is not the case globally. Indeed, given a body force of zero everywhere, the deformation with minimum asymmetric KL energy is still the identity transformation. Thus, expansion in some regions would induce contraction in others, driving overall asymmetric KL energy upwards and away from zero. Here we show that in practice, the asymmetric unbiased method does not seem to perform much differently than its symmetric version. This further supports our conclusion that the log transformation may be a more fundamental operation than symmetrization in understanding Jacobian maps in the context of medical imaging.
We compared Fluid registration and Unbiased registration methods with different values of parameters σ and λ. Fluid registration, as well as Unbiased registration with smaller values of λ, generated noisier mean maps, while maps generated using Unbiased registration with larger values of λ suggest a volume reduction in gray matter as well as ventricular enlargement. As the value of smoothing parameter σ increased, the Jacobian map became smoother, making the biological effects, such as reduction in gray matter, harder to detect.
Our analyses showed that both Asymmetric Unbiased and Symmetric Unbiased models perform significantly better than the fluid registration technique and demonstrated that the Unbiased registration, with any choice of parameters, has more power differentiating between regions of change and no change than the inverse-consistent linear elastic matching of (Christensen and Johnson (2001); Leow et al. (2005)).
The fluid registration is indeed a useful nonrigid registration technique for many applications since it generates a close alignment between the images being registered. Even so, as we showed in this paper, the fluid registration has its limitations in tensor-based morphometry. Unbiased regularization, however, ensures that the resulting deformations have intuitive axiomatic properties by penalizing any bias in the corresponding statistical maps. With Unbiased registration, we can generate accurate alignment while ensuring that deformations have intuitive axiomatic properties.
This also leads to a key issue in the field of non-linear registration, as validation studies have been relatively difficult to perform as ground truth is not generally available. Even in cases where test images were transformed using simulated deformations and thus the ground truth is known a priori, one may still wonder whether these simulated deformations correspond to anatomically plausible brain changes. As a result, efforts such as the ADNI project are important, as it provides a platform for researchers to compare nonlinear registration techniques using standardized imaging data. Along this direction, here we reported one of the first studies using the ADNI dataset, and with potentially interesting results for the rest of the registration community.
Importantly, the proposed Unbiased framework, though a novel concept in medical image registration, can be easily adapted and combined with any non-rigid registration algorithm. In this paper, we implemented the unbiased registration by adapting a fluid regularization algorithm in order to show its advantages in measuring biological deformations. The construction of Unbiased nonlinear elastic registration algorithm (Yanovsky et al. (2008a,b)) is another example of how an arbitrary nonrigid registration algorithm can be extended to compute unbiased deformations.
In addition to investigating various regularizers, we compared L2 and Mutual Information matching functionals in detecting changes in tensor-based morphometry. When applied to serial scans obtained using the same protocol, the results were inconclusive when comparing L2-Unbiased and MI-Unbiased (both asymmetric and symmetric) models. However, L2-Fluid performs less favorably than MI-Fluid. In other words, mutual information performs better as a fidelity metric when coupled with Fluid registration, but not in the case of unbiased registration.
To explain this result, we postulate that by not constraining the final deformations (as in Fluid registration), assuming intensity 1-to-1 correspondence (i.e., matching using L2) may lead to local oscillations of the deformation maps, as minimizing L2 forces a local search for the smallest intensity differences. One result of this is a Jacobian map with locally extreme values, translating into spurious signals and, in our case, less reproducibility. On the other hand, the Unbiased method eliminates local oscillations, allowing better realization of true physiological signals even when L2 is used as a data fidelity term. Again, this supports our conclusion that unbiased framework is fundamental in the understanding and realization of physiological signals.
Although various techniques have been extensively applied to detect disease effects and monitor brain changes with TBM, this paper is the first calibration study to systematically compare registration models for tensor-based morphometry. We believe our results are important, as they provide greater insight into the interpretation of TBM results in the future.
This work was supported in part by the National Institutes of Health under Grant U54 RR021813 NIH/NCRR, Grant U01 AG024904, Grant P41 RR13642, and Grant R21 EB001561. The work of Igor Yanovsky was also partially supported by NSF VIGRE Grant DMS-0601395 and CCB-NIH Grant 30886, and was also carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. The work of Paul Thompson was also supported in part by the National Center for Research Resources, the National Institute for Biomedical Imaging and Bioengineering, the National Institute of Aging, and the National Institute for Neurological Disorders and Stroke, Grant R21 RR019771, Grant EB01651, Grant AG016570, Grant NS049194.
In this Appendix, we derive explicit expressions for uR(u) in (18) when Ω 2. Let us denote the components of vector x to be (x1, x2) and the components of vector u be (u1, u2). We also denote jui = ui/xj.
To simplify the notation, we let J = |Dg| = |D(x – u)|. Also, denote L(J) = LKL(J) = − log J, when R = RKL and L(J) = LSKL(J) = (J – 1) log J, when R = RSKL. Note that J : 2×2() → , where M2×2() is the set of 2 × 2 matrices with real elements, and L : → . Jacobian J is a function of jui, for i,j = 1,2, and is given by
We would like to minimize the functional
We find the first Euler-Lagrange equation. For some :
With notation L' = dL/dJ, the first Euler-Lagrange equation becomes:
Thus, minimizing the energy R(u) with respect to u1, for fixed u2, yields the first component of uR(u):
Note that L'KL(J) = −1/J and L'SKL(J) = 1 + log J – 1/J.
Similarly, the Euler-Lagrange equation for the second component of uR(u) can be found to be:
In this Appendix, we derive an explicit expression for uR(u) in (18) when Ω 3. Let us denote the components of vector x to be (x1, x2, x3) and the components of vector u be (u1, u2, u3). Here, we will use the same notation we used in Appendix A.
Jacobian J is a function of jui, for i, j = 1, 2, 3, and is given by
We would like to minimize the functional
For some η, we have
Hence, the first Euler-Lagrange equation becomes:
Thus, minimizing the energy R(u) with respect to u1, for fixed u2 and u3, yields the first component of uR(u):
Similarly, the other two Euler-Lagrange equations can be found to be:
In this Appendix, we derive the gradient uFMI(u) of the mutual information matching functional in (14), adopting the approach of (Chefd'Hotel et al. (2001); Hermosillo et al. (2001)), modeling the joint intensity distribution of deformed image I2(x − u) and image I1(x) as a continuous function using the Parzen windowing method.
We compute the first variation of FMI(u) by perturbing u in the following way
Thus, we have
However, note that
Hence, the second term on the right hand side of the equality in (C.2) reduces to
Equation (C.1) becomes
The joint intensity distribution estimated from I2(x − u) and I1(x) is given by
where |Ω| is a volume of Ω and ψ(ξ1,ξ2) is a two-dimensional Parzen windowing kernel.
The derivative of (C.7) can also be computed:
Let us denote
If we let ε = 0, equation (C.6) gives the first variation of FMI(u):
Here, * denotes a convolution. Thus,
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.