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. Tensorbased 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 inverseconsistent 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 2week and 1year 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, imageguided 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 multisubject studies, this reduces subjectspecific 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 tensorbased 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 nonrigid 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 largedeformation image registration approach. In this context, unbiased means that we strive to obtain a zeromean and symmetric logJacobian 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 logJacobian 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 nonrigid registration techniques in tensorbased morphometry (TBM). We systematically investigate and compare the performances of different nonrigid registration techniques including two common matching functionals: L^{2}, or the sum of squared intensity differences, versus mutual information, and four regularization techniques (fluid registration, inverseconsistent 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 2week 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 followup 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 fluidfilled spaces in the brain. Thus, a good computational technique should be able to differentiate between longitudinal image pairs collected for the ADNI baseline (2week) and followup (1year) 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 T1weighted 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 highresolution sagittal T1weighted 3D MPRAGE sequence for each subject, with a reconstructed voxel size of 0.9375 × 0.9375 × 1.2 mm^{3}. 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 nonlinearity (Jovicich et al. (2006)), (2) a “B1correction”, to adjust for image intensity nonuniformity 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 sessionspecific calibration errors. Additional phantombased 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 I_{1} : Ω → and I_{2} : Ω → be the two images to be registered. We seek to estimate a transformation g : Ω → Ω such that I_{2} matches I_{1} when deformed by g. In this paper, we will restrict this mapping to be differentiable, onetoone, 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 I_{2}○g(x) back to I_{2}(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 LargeDeformation 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 KullbackLeibler (KL) divergence and symmetric KullbackLeibler (SKL) distance. The KL divergence between two probability density functions, p_{1}(x) and p_{2}(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 I_{2} and I_{1}, 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 R_{KL}, which is an asymmetric measure between P_{id} and P_{g}:
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 I_{2}(x − u) and I_{1}(x). The second term on the righthand side of the equality in (7) is equivalent to the logarithmic barrier in numerical optimization theory (Nocedal and Wright (1999)) and is wellbehaved.
In this section, we describe the regularization functional based on the symmetric KL distance between P_{id} and P_{g}:
As shown in (Leow et al. (2007); Yanovsky et al. (2007b)), the regularization term is linked to statistics on Jacobian maps as follows
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 nonnegative, 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 L^{2} 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 L^{2}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 L^{2} distance between the deformed image I_{2}(x − u) and target image I_{1}(x) is defined as
Computing the first variation of functional F_{L}_{2} 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 I_{2}(x − u) and the target image I_{1}(x), we follow the notations in (Hermosillo et al. (2001)), where p^{I1} and ${p}_{\text{u}}^{{I}_{2}}$ are used to denote the intensity distributions estimated from I_{1}(x) and I_{2}(x − u), respectively. An estimate of their joint intensity distribution is denoted as ${p}_{\text{u}}^{{I}_{1},{I}_{2}}$. In this probabilistic framework, the link between two modalities is fully characterized by a joint density.
We let i_{1} = I_{1}(x), i_{2} = I_{2}(x − u(x)) denote intensity values at point x Ω. Given the displacement field u, the mutual information computed from I_{1} and I_{2} is provided by
We seek to maximize the mutual information between I_{2}(x − u) and I_{1}(x), or equivalently, minimize the negative of $M{I}_{\text{u}}^{{I}_{1},{I}_{2}}$:
The gradient of (14) is derived in Appendix C and is given by
where Ω is the volume of Ω, Q_{u} is defined as
and ψ(ξ_{1}, ξ_{2}) is a twodimensional Parzen windowing kernel, which is used to estimate the joint intensity distribution from I_{2}(x − u) and I_{1}(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(I_{1}, I_{2}, u), namely _{u}E(I_{1}, I_{2}, u). We define the force field f, which drives I_{2} into registration with I_{1}, as
Here, R(u) is either R_{KL}(u) or R_{SKL}(u). Explicit expressions for components of _{u}R(u), in both cases, are derived in Appendices A and B for two and three dimensional cases, respectively. Also, the gradient _{u}F(I_{1}, I_{2}, 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 EulerLagrange 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
and
using equations (18), (21), (20) will be referred to as L^{2}Asymmetric Unbiased and L^{2}Symmetric Unbiased models, respectively. The model above, provided λ = 0, will be referred to as the L^{2}Fluid model.
Similarly, minimization of
and
will be referred to as the MIAsymmetric Unbiased and MISymmetric Unbiased models, respectively. Such models, with λ = 0, define the MIFluid model.
We are now ready to give an algorithm for minimizing either one of energy functionals (22) through (25) above.
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 voxelwise deviation gain of A over B (denoted by S^{A,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 withinsubject paired t test and a group paired t test. A withinsubject 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 voxelwise tmap 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 onesided alternative hypotheses. For example, when testing if B outperforms A, we use the following alternative hypothesis
For n subjects, the voxelwise T statistic, defined as
where
and
thus follows the Student's t distribution under the null hypothesis and may be used to determine the pvalue 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 followup dataset (in which patients are scanned twice with MRI, one year apart) and ADNI baseline dataset, we create a voxelwise 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 voxelwise T statistic compared to a twotailed Student's t distribution may then be used to test the above null hypothesis
where
and
We reject the null hypothesis if the p value calculated above exceeds a preset threshold based on a suitable confidence interval. Notice the voxelwise variance of log J provides us with a way to assess the repeatability of a deformation method, i.e., measuring the voxelwise 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 ${S}_{i}^{A,B}$ (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 abovedefined percentage exceeds that of the unpermuted observed data. This is comparable to ‘setlevel inference’ in the widelyused 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 abovedefined percentage greater than that of the original data). All possible (2^{10} = 1024) permutations were considered in determining the final corrected p value.
To visually assess the global significance level of the voxelwise t tests on deviation gains and logJacobian 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 suprathreshold voxels in a statistical map, for a range of thresholds, these CDF plots (sometimes called QQ 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 logJacobian values, a better registration technique, on the other hand, should be able to separate the CDF curves between ADNI baseline and followup 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 inverseconsistent 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 L^{2} and mutual information matching functionals (see equations (22)(25)).
To obtain a fair comparison, regridding was not employed. Regridding is a method to relax the energy computed from the linear elasticity prior after a certain number of iterations, which allows largedeformation 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 memoryless procedure, as how images are matched after each regridding is independent of the final deformation before the regridding, rendering the comparison of final Jacobian fields and cost functionals problematic. Moreover, we consider the strategy of regridding, 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 logJacobian).
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 MIUnbiased 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 ttest to compare the results of one set of parameters versus the other. Unless otherwise mentioned, MIUnbiased registration results were generated using values of σ = 9.0 and λ = 5.0. A similar procedure was employed to choose the parameters for L^{2}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 crossvalidation.
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 twoweek period using the same MRI protocol, no systematic structural changes should be recovered.
In our first experiment, we compared methods based on L^{2} matching (L^{2}Fluid, L^{2}Asymmetric Unbiased, and L^{2}Symmetric Unbiased). Uniform values of λ = 500 and λ = 1000 were used for all deformations using L^{2}Symmetric Unbiased and L^{2}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 144 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 L^{2} 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 L^{2}norm decrease per iteration for all L^{2}based models. Figure 3 plots the KL divergence and SKL distance measures for each of the L^{2}based methods. For L^{2}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 voxelwise deviation gains of L^{2}Fluid over L^{2}Asymmetric Unbiased as well as L^{2}Fluid over L^{2}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 $de{v}_{2}(\text{x})=\leftD\text{g}(\text{x})1\right$, 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 L^{2}Fluid and L^{2}Symmetric Unbiased methods, conducting a withinsubject 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 L^{2} matching cost functional, produces more reproducible maps with less variability.
Figure 6 shows the mean Jacobian maps obtained using L^{2}Fluid, L^{2}Asymmetric Unbiased, and L^{2} Symmetric Unbiased registration algorithms. Jacobian maps generated using unbiased models have values closer to 1, whereas L^{2}Fluid model generated noisy mean maps. Figure 7, shows the results when performing 3D voxelwise paired t tests for the deviation gain of L^{2}Fluid over L^{2}Asymmetric Unbiased and L^{2}Fluid over L^{2}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 L^{2}Fluid over L^{2}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 L^{2}Symmetric Unbiased method outperforms L^{2}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 pvalues 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 xyplane, 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 (MIFluid, MIAsymmetric Unbiased, and MISymmetric Unbiased). Uniform values of λ = 5 and λ = 10 were used for all deformations using MISymmetric Unbiased and MIAsymmetric 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 101313 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 MIAsymmetric Unbiased and MISymmetric Unbiased methods to produce inverse consistent maps with less variability. Figure 11 shows the mutual information decrease per iteration for all MIbased models. Figure 12 plots the KL divergence and SKL distance measures for each of the MIbased methods. For MIFluid 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 voxelwise deviation gains of MIFluid over MIAsymmetric Unbiased as well as MIFluid over MISymmetric Unbiased. The histograms are skewed to the right, indicating the superiority of both unbiased registration methods over Fluid registration.
In Figure 14, we compared MIFluid and MISymmetric Unbiased methods, conducting a withinsubject 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 MIFluid, MIAsymmetric Unbiased, and MI Symmetric Unbiased registration algorithms. Jacobian maps generated using unbiased models have values closer to 1, whereas MIFluid model generated noisy mean maps. Figure 16, shows the results when performing 3D voxelwise paired t tests for the deviation gain of MIFluid over MIAsymmetric Unbiased and MIFluid over MISymmetric 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 MIFluid over MISymmetric 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 MISymmetric Unbiased method outperforms MIFluid with p < 0.001.
In Figure 18, we plotted the cumulative distribution function of the pvalues 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 L^{2} 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 withinsubject paired t tests (Figure 19) as well as group paired t tests (Figure 20) on the voxelwise deviation gains for all voxels inside the ICBM brain mask. We showed that MIFluid outperforms L^{2}Fluid with p < 0.0001. However, the result of the comparison of L^{2}Symmetric Unbiased and MISymmetric Unbiased is inconclusive.
In this section, we analyze a dataset we shall call the “ADNI Followup” 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 pvalues for the voxelwise log Jacobian tmaps for both ADNI Baseline and ADNI Followup 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 pvalues for the voxelwise log Jacobian tmaps for both ADNI Baseline and ADNI Followup 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 inverseconsistent linear elastic matching of (Christensen and Johnson (2001); Leow et al. (2005)). Figure 28 displays the cumulative distribution of pvalues for the voxelwise log Jacobian tmaps for both ADNI Baseline and Followup datasets for Unbiased technique with different values of λ, Fluid registration, and inverseconsistent linear elastic matching. Unbiased methods, especially with larger values of parameter λ, show a bigger separation between the baseline and followup curves than Fluid and inverseconsistent 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 nonrigid registration techniques in TBM. We systematically investigated and compared the performances of different nonrigid registration techniques including common matching functionals (L^{2}metric and mutual information), and regularization techniques (fluid registration, inverseconsistent 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 inverseconsistent 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 tensorbased 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 nonlinear 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 nonrigid 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 L^{2} and Mutual Information matching functionals in detecting changes in tensorbased morphometry. When applied to serial scans obtained using the same protocol, the results were inconclusive when comparing L^{2}Unbiased and MIUnbiased (both asymmetric and symmetric) models. However, L^{2}Fluid performs less favorably than MIFluid. 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 1to1 correspondence (i.e., matching using L^{2}) may lead to local oscillations of the deformation maps, as minimizing L^{2} 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 L^{2} 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 tensorbased 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 DMS0601395 and CCBNIH 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 _{u}R(u) in (18) when Ω ^{2}. Let us denote the components of vector x to be (x_{1}, x_{2}) and the components of vector u be (u_{1}, u_{2}). We also denote _{j}u_{i} = u_{i}/x_{j}.
To simplify the notation, we let J = Dg = D(x – u). Also, denote L(J) = L_{KL}(J) = − log J, when R = R_{KL} and L(J) = L_{SKL}(J) = (J – 1) log J, when R = R_{SKL}. Note that J : _{2×2}() → , where M_{2×2}() is the set of 2 × 2 matrices with real elements, and L : → . Jacobian J is a function of _{j}u_{i}, for i,j = 1,2, and is given by
We would like to minimize the functional
We find the first EulerLagrange equation. For some $\eta \in {C}_{c}^{1}(\Omega )$:
With notation L' = dL/dJ, the first EulerLagrange equation becomes:
Thus, minimizing the energy R(u) with respect to u_{1}, for fixed u_{2}, yields the first component of _{u}R(u):
Note that L'_{KL}(J) = −1/J and L'_{SKL}(J) = 1 + log J – 1/J.
Similarly, the EulerLagrange equation for the second component of _{u}R(u) can be found to be:
In this Appendix, we derive an explicit expression for _{u}R(u) in (18) when Ω ^{3}. Let us denote the components of vector x to be (x_{1}, x_{2}, x_{3}) and the components of vector u be (u_{1}, u_{2}, u_{3}). Here, we will use the same notation we used in Appendix A.
Jacobian J is a function of _{j}u_{i}, for i, j = 1, 2, 3, and is given by
We would like to minimize the functional
For some η, we have
Hence, the first EulerLagrange equation becomes:
Thus, minimizing the energy R(u) with respect to u_{1}, for fixed u_{2} and u_{3}, yields the first component of _{u}R(u):
Similarly, the other two EulerLagrange equations can be found to be:
and
In this Appendix, we derive the gradient _{u}F_{MI}(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 ${p}_{\text{u}+\varepsilon \eta}^{{I}_{1},{I}_{2}}\left({i}_{1},{i}_{2}\right)$ of deformed image I_{2}(x − u) and image I_{1}(x) as a continuous function using the Parzen windowing method.
We compute the first variation of F_{MI}(u) by perturbing u in the following way
Thus, we have
However, note that
and
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 I_{2}(x − u) and I_{1}(x) is given by
where Ω is a volume of Ω and ψ(ξ_{1},ξ_{2}) is a twodimensional 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 F_{MI}(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.
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. 