Search tips
Search criteria 


J Med Imaging (Bellingham). 2016 January; 3(1): 014005.
Published online 2016 March 16. doi:  10.1117/1.JMI.3.1.014005
PMCID: PMC4793085

Deformation-based atrophy computation by surface propagation and its application to Alzheimer’s disease


Obtaining regional volume changes from a deformation field is more precise when using simplex counting (SC) compared with Jacobian integration (JI) due to the numerics involved in the latter. Although SC has been proposed before, numerical properties underpinning the method and a thorough evaluation of the method against JI is missing in the literature. The contributions of this paper are: (a) we propose surface propagation (SP)—a simplification to SC that significantly reduces its computational complexity; (b) we will derive the orders of approximation of SP which can also be extended to SC. In the experiments, we will begin by empirically showing that SP is indeed nearly identical to SC, and that both methods are more stable than JI in presence of moderate to large deformation noise. Since SC and SP are identical, we consider SP as a representative of both the methods for a practical evaluation against JI. In a real application on Alzheimer’s disease neuroimaging initiative data, we show the following: (a) SP produces whole brain and medial temporal lobe atrophy numbers that are significantly better than JI at separating between normal controls and Alzheimer’s disease patients; (b) SP produces disease group atrophy differences comparable to or better than those obtained using FreeSurfer, demonstrating the validity of the obtained clinical results. Finally, in a reproducibility study, we show that the voxel-wise application of SP yields significantly lower variance when compared to JI.

Keywords: registration, atrophy computation, Jacobian integration, cube propagation

1. Introduction

Alzheimer’s disease is known to manifest itself through the degeneration of certain nerve cells.1 This degeneration may be captured through structural magnetic resonance imaging (MRI) as “atrophy.”25 Using deformation-based computational anatomy tools, two kinds of atrophy estimates can be derived from structural MRI scans: (a) regional atrophy—atrophy in disease-relevant brain regions and (b) voxel-wise atrophy—per voxel volume change. The latter has a wider scope since not only can it be used to compute (a), but also to perform voxel-wise statistics on populations to observe general atrophy trends in the disease. Additionally, deformation-based atrophy estimation methods are appealing due to its less reliance on accurate delineations of regions of interest which allows for computing atrophy in regions that are difficult to segment both manually and automatically, for instance amygdala.

Estimating deformation-based atrophy involves using image registration to derive deformation fields that describe spatial transformation needed to match brain structures. Since the deformation field encodes relative positioning of each voxel in the brain, their gradients encode information such as the volume change via the Jacobian matrix. The determinant of this matrix yields local volume change information that can be later integrated to obtain regional atrophy measures. This method is popularly termed as Jacobian integration (JI).69 The numerical precision of JI relies on the smoothness of the deformation with respect to the discretization grid (See Sec. 5.1). As an alternate to JI, atrophy can also be computed using volumetric meshing (VM).1014 Here, each voxel is subdivided into simplices (with closed form expressions for their volume), and the volume of each of these simplices are summed after propagation of their vertices through the deformation field. Due to the completeness in the representation of a deformed voxel, this method is more precise than JI. We will refer to this method as simplex counting (SC). Figure 1 refers to the classes of deformation-based atrophy estimation methods considered in this paper. While SC has been proposed previously in the literature,10,1114 the mathematical properties that underpins its improvement over JI and a thorough evaluation against JI is missing in the literature. In lieu, this paper primarily deals with assessing both JI and geometric methods in the context of atrophy estimation. To this end, we particularly focus on the following:

Fig. 1

Classification of atrophy estimation methods considered in this study. VM—volumetric meshing; SP—surface propagation; SC—simplex counting; and JI—Jacobian integration.

  • • We propose surface propagation (SP)—an improvement in the manner of SC meshing that significantly reduces its computational complexity.
  • • We derive the orders of approximation of JI (using midpoint, trapezoidal, and Simpsons rule) and SP/SC.
  • • We inspect the influence of using SP/SC for atrophy estimation in Alzheimer’s disease (AD).

SP is a time-discrete, surface flux-based method that estimates the volume of an object (a cubic voxel or a segmentation) from the surface of the object only. The essence of SP is the same as SC, i.e., dividing the volume into simplicies. However, SP is constructed in such a way that the interior volumes cancel out and the contribution to the volume is only from the surface. This is because each voxel is meshed only on the exterior in the first place unlike SC where subdividing the interior of a voxel is paramount. One may derive similarities to the Gauss divergence theorem where SP relates to the surface integral and SC to the volume integral. Except for the empirical experiment, SP is considered as the representative of the general volumetric meshing methods, i.e., observations derived from the head-to-head comparison with JI can be extended to SC too.

SP is compared head-to-head with JI in three experiments: (1) robustness to noise in an artificial dataset; (2) performance in separating diagnostic groups in a standardize dataset from the Alzheimer’s disease neuroimaging initiative (ADNI); (3) scan-rescan performance in a small, in-house dataset (voxel-wise SP is used in this case). In addition, we compare the obtained results in (2) to those from FreeSurfer, in order to check the clinical validity of the combination of the registration methodology and SP. Preliminary applications of voxel-wise SP can be found in Refs. 15 and 16.

2. Related Work

Deformation fields used to extract voxel-wise atrophy information can either be obtained from a displacement-based17,18 or a flow-based1922 registration scheme. Given a registration scheme, the Jacobian matrix relevant for computing voxel-wise atrophy can be analytically computed if the transformation model is parametric, for instance B-spline-based freeform deformation schemes. For methods that encode transformation in each voxel (like Demons), the Jacobian needs to be approximated using an appropriate finite differencing scheme. The geometric atrophy estimation method (SP/SC) is attractive because it only relies on the final position of voxels after the deformation. Hence, atrophy can be estimated in a consistent fashion regardless of the underlying deformation model.

Regional atrophy can be obtained by integrating the Jacobian determinant values over an ROI. However, there are specialized methods available to compute region specific atrophy. Although such methods are not a focus of this paper, we will briefly discuss them.

One of the most popular automated way of computing atrophy is using the boundary shift integral (BSI).23 In this method, atrophy is captured as a difference in intensities (residing in an intensity window) resulting from a shift in the boundaries of the structure of interest. There is an implicit assumption of robust intensity normalization to this method. The available implementation of BSI computes whole brain and hippocampal atrophy. BSI is a global measure as opposed to a voxel-wise measure from deformation fields. In addition, BSI depends on hard segmentation whereas soft segmentations can be used to obtain regional atrophy from deformation fields.

Regional atrophy may also be computed by computing the difference in the number of voxels residing in the segmentation of a region in successive MRI scans. Some of the popular methods in this category are FreeSurfer24,25 and SIENA(X).26,27 Recent studies28,29 have shown that voxel-wise JI produced atrophy estimates with less variance and increased statistical power when compared to the popular segmentation methods. A more thorough overview of existing atrophy estimation methods can be found in Ref. 30.

Another class of registration-based atrophy estimation method is flux based.21 This method has been published in conjunction with stationary velocity field-based registration methods. Figure 2 shows an overview of obtaining volume change information from different classes of transformation models. In case the transformation model is flow based, i.e., velocity fields (v) integrated to obtain deformation (ϕ), then one can either estimate the flux first (through surface integration, SI) and then exponentiate the same or exponentiate the deformation and then use SP to obtain volume change. However, if one has access to only the deformation then in order to compute the volume change via the flux, one needs to first compute the logarithm of the deformation to get the velocity field and then compute the flux and volume change. Computing the Log is an additional step that may lead to numerical imprecision. We do not explore flux estimation from the velocity field in this study because in case of SP, regardless of the way a deformation is obtained, the volume change information can be obtained in a single step as shown in Fig. 2.

Fig. 2

Volume change computation from registration, v is the velocity fields, ϕ is the deformation, Exp is the Lie group exponential, SI is the surface integration, and Ω is the region of interest.

3. Outline

We will begin by introducing the registration scheme used in the paper in Sec. 4. Following this, in Secs. 5.1 and 5.1.1, we will discuss JI and the orders of approximation of JI. In Sec. 5.2, we will introduce the contribution of this paper—surface propagation and in Sec. 5.2.1, we will derive the order of approximation of SP, which applies to SC too. Further, in Sec. 6, we will present the experiments conducted to illustrate the use of SP in atrophy estimation and compare it with JI. Finally, in Secs. 7 and 8, we will discuss the observations from the experiments and provide concluding remarks.

4. Registration

In this section, we present a brief overview of the registration methodology used to estimate the deformation field. The registration is formulated as a multiscale (image smoothing) rigid transformation followed by a multiscale (image smoothing), multiresolution (control point resolution) nonrigid transformation.16,31,32

In the registration pipeline, the cost functions for both the rigid and nonrigid registrations are on the form


where the similarity measure M is normalized mutual information (NMI),31 I is the moving image, R is the reference or baseline image, λ0 is a user defined regularization constant, S is the regularization term and ϕ(x,y,z):R3R3 is a deformation field that maps a location (x,y,z) in the reference image to the moving image.

In the rigid registration, λ=0 and ϕri is a linear transformation with six degrees of freedom. In the non-rigid registration, λ>0 and ϕnr is a free-form deformation based on cubic b-splines17


where l=i+x/δ, m=j+y/δ, n=k+z/δ, u=(x/δ)x/δ, v=(y/δ)y/δ, w=(z/δ)z/δ, δ is the spacing between the control points, pl,m,n are b-spline control point coefficients at l, m, n, and Bi, Bj, Bk are the one-dimensional cubic b-spline tensors.

The nonrigid registration is driven by the cost function (1) with a regularization based on the discrete Laplacian in the spline’s control points (pl,m,n), given by


The reason for using this regularization is twofold: (1) It penalizes oscillations in displacements leading to a smoother deformation field and (2) using a second-order term makes the regularization unbiased in first order transformations (affine image deformations) and thereby globally unbiased w.r.t. atrophy. This setting ensures the smoothness of the Jacobian and in practice often nonfolding. A composition of multiple transformations is used, ϕ(x,y,z)=(ϕnϕn1ϕ1)(x,y,z), because this enables modeling larger deformations without self intersections. A limited memory Broyden–Fletcher–Goldfarb–Shanno33 scheme is used for optimization. The optimization was performed in MATLAB using the min-Func package.34

To reduce memory consumption, only a subset of the image voxels are considered in evaluating the similarity measure M. At the finest grid of the multiscale, multiresolution deformation, we chose a denser sampling of the evaluation points in a band region around the cerebral boundaries and the ventricular boundaries. The band is generated by the following mathematical morphology operations: (IsE)(IsE), where E is the structuring element [three-dimensional (3-D) sphere of radius 6 mm], is dilation, is erosion, and Is is the binary image consisting of the brain segmentations (minus ventricles, cerebrospinal fluid, and brain stem). In addition, dense evaluation points in a box (five voxels beyond the boundary) containing the hippocampus is included in the finest grid. Note that the registration method is not an out-of-box freeform deformation17 implementation. The difference is in how multiscale resolution of control points is handled. While Ref. 17 uses knot-splitting, we use a composition of warps approach.

5. Atrophy Computation

Given a deformation field ϕ, the volume change can be computed in order to quantify the amount of atrophy. Here, we present two methods for atrophy computation from the deformation field, the widely used JI (for different integration schemes), and SP, the alternative we propose.

5.1. Jacobian Integration

Given an ROI Ω and its volume vol(Ω), the volume of the transformed region ϕnr(Ω) is given as


where J(x,y,z)=ϕnr(x,y,z) is the Jacobian of the transformation field, and |·| denotes the determinant. In this work, we use a warp of a warp approach. Each intermediate warp (n) may be defined as follows


where fn and pn are the B-Spline and its parameter of the n’th warp of ϕnr. The integral of the determinant of Jacobian of such a transformation may be written as


Here, n represents also the intermediate warps (lower n have coarser resolution of control points and higher n have finer resolution). Note that each warp is parameterized separately. Hence, the partial derivatives can be analytically evaluated at the sampling points. It is possible to analytically evaluate the integral for a single warp,35 however, the analytical integral for a composition of warps is nontrivial: If ϕ is a third-order spline, then its Jacobian is a matrix of second-order polynomials, its determinant is a sixth-order polynomial, and the product of n determinants is a 6n polynomial. Hence, JI is usually approximated as a Riemannian sum of the left-hand-side of Eq. (6):


where x, y, z are coordinates that are followed through the composition of warps, h is the sampling distance and i, j, kZ3.

5.1.1. Errors of approximation for Jacobian integration

The integral of Jacobian determinants in a composition of warps setting is not trivial to compute, hence the Jacobian determinant values at the voxel center are summed to obtain an approximation. Alternatively, it can be computed using the Lagrange’s linear and quadratic spline interpolants in the form of trapezoidal or Simpson’s rule for a better approximation. The truncation errors of these integrations schemes are O(h2), O(h3), O(h3), O(h4) for left point, midpoint, trapezoidal and Simpson’s rule, where h is the sampling distance. For detailed discussion on truncation errors, we direct the readers to Ref. 36.

5.2. Surface Propagation

SP is a modified SC approach, where each face of a voxel is divided into two triangles, such that neighboring voxels share triangles. Unlike SC, the interior of the voxels is not subdivided. Since the interior is not subdivided, the contributions from adjacent voxel faces cancel out and the contribution to the volume is from the triangles on the surface of the voxel. We will illustrate the theory behind surface propagation by considering a single voxel as an object. If one desires voxel-wise change maps similar to JI and SC, such a representation is useful.

The key element in SP is the Gauss’ divergence theorem, which relates surface integrals to the volume of an object. Consider a smooth (C1) surface S=Ω as the boundary of a ROI, Ω, and a smooth deformation ϕ:R3R3, then Gauss’ divergence theorem relates the closed surface integral to the volume integral


where dv=dxdydz is the unit volume element, da is the unit surface area element, n=(nx,ny,nz) is the (outward) surface unit normal, and ·ϕ is the divergence of ϕ.

Now consider an axis-aligned flow such that ϕ^=(x,0,0), ·ϕ^=1 and ϕ^·n=xnx. Now Eq. (8) can be rewritten as


Consider a subdivision of Ω into a set of nonoverlapping volumes ω such that


where ω has the shape of deformed voxel, and each side is triangulated such that triangle i has area Ai. Considering one subvolume, ω˜i, if we define the x coordinate of the barycenter of each surface triangle as xi=Aixda/|A| and the unit normal along the x-axis as nxi then the discrete volume form over ω, is


where ω˜ is the surface of the sub-volume, N is the number of triangles and Axi=nxiAi is the area of the I’th triangle projected down to the y-z plane. Since the normals of triangles on faces of adjacent voxels cancel each other out, the only triangles that contribute to the volume are outermost faces of the outermost voxels. Therefore, Eq. (11) will boil down to


where Ns is the number of triangles on the outermost faces of the the outermost voxels.

Given all the surface triangles are planar, substitution of the integral by a sum in Eq. (12) is exact. However, in a nonrigid framework this is not the case, as the planar surface triangles are nonlinearly deformed with nonplanar shapes as shown in Fig. 3. Hence, the triangular model after deformation is only an approximation, the order of which will be derived in Sec. 5.2.1.

Fig. 3

A linear approximation of a deformed surface by a triangle. The volume between the planar triangle and the surface is the error of approximation in both SP and JI; h and k are the sampling distances of the discrete grid, ϕ˜ is the transformation ...

An illustration of the method for one voxel as an object can be seen in Fig. 4. In the figure, the surfaces of a voxel are represented by N triangles that undergo a transformation. The total volume change of the transformed voxel is the sum of the volume under each of the N transformed triangles. The volume under a transformed triangle is computed by multiplying the projected area, Axi (on y-z plane in the illustration), of the I’th triangle by the x-component of the transformed centroid of the triangle ϕnr(xi). For instance, if (x1,y1,z1),(x2,y2,z2),(x3,y3,z3) are three points at one face of the p’th voxel forming a triangle, ϕnr(x1,y1,z1),ϕnr(x2,y2,z2),ϕnr(x3,y3,z3) is their new position. Once the new position of these three points is known, the x coordinate of the centroid of triangle formed by these points is multiplied with the projected area of this triangle onto the y-z plane. For regional measures, SP is applied where the vertices of the triangles on the surface are propagated through the deformation and the volume under each of the deformed triangles are summed (the interior triangles do not contribute to the volume).

Fig. 4

An illustration of how the volume change of a single voxel is computed in surface propagation. Each surface of the voxel is represented by two planar triangles (marked in red) that undergoes a transformation. The volume under each transformed triangle ...

5.2.1. Error of approximation for cube propagation

In SP and SC, the error of approximations are nearly the same and only different by a constant factor. For all practical purposes, this difference is negligible. In a voxel deformed by a transformation function ϕnr and triangulated as described in Sec. 5.2, the error of approximation will be limited to the surface, i.e., the volume between the actual surface and the linear approximation of it in the form of a set of triangles. The surface error on an interior volume element is canceled by the same surface element of the neighboring volume element, thus leaving only errors on the surface of an object. We will now describe the order of approximation of SP.

Consider a surface triangle with vertices (v1,v2,v3) with barycentric coordinates (b1,b2,b3), where ||v1v3||=h and ||v2v3||=k are small, and consider a warp ϕ˜ such that takes the original vertices v1=(1,0,0), v2=(0,1,0), and v3=(0,0,1) to ϕ˜(v1)=(h,0,0), ϕ˜(v2)=(0,k,0), and ϕ˜(v3)=(0,0,0). Although the deformation is nonlinear, the Jacobian of the transformation can be easily computed as the sine of the angle θ between the sides v1v3 and v2v3. The illustration of the deformation can be seen in Fig. 3. The local Cartesian coordinates (x,y) of an arbitrary point v inside the triangle may be expressed in terms of barycentric coordinates as follows,


where A=(1/2)hksin(θ) is the area of the triangle with θ as the angle between the sides of lengths h, k and db1=dx/h, db2=dy/k. This integral allows us to compute the volume constrained by the sides of the triangle, where the height at each point, v, is given by the function f(v). The interesting part is that we can express the integral in terms of the vertices of the triangle and barycentric coordinates only. Accordingly, the integral in ϕ˜’s coordinate system will be


The function f(x,y) can be expressed using a Taylor series as follows


By construction f(0,0)=f(h,0)=f(0,k)=0, since they lie on the surface f, and therefore we have fx(0,0)=(h/2)fxx(0,0)+O(h2) and fy(0,0)=(k/2)fyy(0,0)+O(k2). Thus, the integral on the triangle given in Eq. (16) can be rewritten as


Thus, assuming O(h)=O(k), then the error of approximation of a surface by triangulation is O(h4).

6. Experiments and Results

6.1. Synthetic Noise Experiment

In order to inspect the robustness of SP and JI, we test the precision of each of the methods in estimating the true volume change in the presence of noise. We consider a deformation of a grid in 3-D space using a field proportional to the gradient of the Gaussian plus random noise. The deformation is given as


where ϕi is the displacements, η>0 is a parameter used to control the magnitude of the noise component, σ is the Gaussian width, f is a frequency picked from a uniform discrete distribution 1 to 10 Hz, and χ is a random phase also picked from a uniform distribution [π,π]. σ=0.9  mm is the width of the Gaussian. The parameter η was varied from 0 to 0.01. This deformation and noise model is simplistic and only serves the purpose of investigating the deterioration of the various integration schemes when the signal to noise ratio is varied. The difference in the true volume and the mean approximated volume is measured as a function of η. Jacobian determinant values were computed at every deformed grid point, and the integral was approximated using left point, midpoint, trapezoidal, and Simpson’s rules. For each η, the experiment was repeated 100 times, and the mean of the approximated volume was recorded. The test involved keeping the grid size constant to 9  mm3. The actual deformation variation recorded on real scan-rescan data (the data used in the experiment in Sec. 6.3) in various substructures spatially scaled to the size of this experiment achieve a noise level of 2×103 for ventricles to 9×103 for the hippocampus.

From Fig. 5, we observe the following: (a) as we begin with no noise in the deformation field, i.e., the beginning of the horizontal axis, JI using Simpson’s rule expectedly has the highest order of precision among the integration schemes considered, (b) as the noise increases, i.e., as we move to the right on the horizontal axis, the accuracy in Simpson’s, trapezoidal and midpoint decreases. As expected, SP and SC are identical, and they are more stable in the presence of noise. Hence, the choice of the integration scheme depends on the signal to noise ratio in the deformation field. Due the similarity, SP will be considered as a representative of VM.

Fig. 5

Precision of atrophy estimation schemes in the presence of noise. The vertical axis denotes the mean absolute difference between the estimated volume change and the analytically derived, noise-free change. The numerically integrated change was performed ...

6.2. Application on a Standardized Alzheimer’s Disease Neuroimaging Initiative Dataset

In this section, we will test the various integrations schemes on a standardized ADNI subset for the diagnostic power of the atrophy computation. We apply SP to data from ADNI in order to inspect estimated atrophy scores in normal controls (NC) and AD patients, and to evaluate the diagnostic performance compared to JI. The clinical validity of the obtained atrophy scores is also checked by comparing the scores qualitatively to those from FreeSurfer.

6.2.1. Data

Data used in the preparation of this article were obtained from the ADNI database.37 The ADNI was launched in 2003 by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering (NIBIB), the Food and Drug Administration, private pharmaceutical companies and nonprofit organizations, as a $60 million, 5-year public–private partnership. The primary goal of ADNI has been to test whether serial magnetic resonance imaging, positron emission tomography, other biological markers, together with clinical data and neuropsychological assessments can be combined to measure the progression of mild cognitive impairment and early Alzheimers disease. Determination of sensitive and specific markers of very early AD progression is intended to aid researchers and clinicians to develop new treatments and monitor their effectiveness, as well as lessen the time and cost of clinical trials, see Ref. 38. The description of the ADNI standardized dataset used here can be found in Ref. 39.

We analyzed baseline and 12-month follow-up 1.5 T1 weighted MRI volumes (169 AD and 101 NC). The demographics of the subjects can be found in Table 1. The raw DICOM images were preprocessed using FreeSurfer’s intensity correction method following which, the images were resampled as isotropic voxels of dimension 2563 with each voxel of size 1  mm3. Segmentations used to define the ROIs (whole brain, ventricles, medial temporal lobe, and hippocampus) for quantification of atrophy were also obtained using cross-sectional FreeSurfer.40

Table 1

Baseline demographics for the standardized ADNI dataset. MMSE: mini mental state examination.

6.2.2. Experimental setup

For each subject in the ADNI dataset, the baseline and month 12 follow-up scans were registered using the parameters in Table 2. The number of transformation levels, and thereby number of compositions, was three and five in the rigid and nonrigid registration, respectively. The volume change was computed from the resulting deformation field using SP and JI. JI was computed using both the left point rule, the midpoint rule, the trapezoidal rule, and Simpson’s rule. Note that the ROIs are defined at baseline and the same exact ROIs were used for both the methods.

Table 2

Various parameters used during the registration.

6.2.3. Metrics for evaluation

To characterize the performance of the various integration schemes, we recorded the mean and standard deviation of atrophy in the diagnostic groups as well as computed the diagnostic capability in terms of effect size measured as the Cohen’s D and the area under the receiver operator characteristic curve (AUC). Cohen’s D—disease specific effect size were computed as


where m1, s1 and m2, s2 are the means and standard deviations of NC and AD subjects respectively. To compute the p-value for the pairwise method comparison, we carried out a two-tailed t-test for the null hypothesis of equal measures N1N2=0, where N1 and N2 are independent random measures. We computed a probability distribution for the difference between the Cohen’s D for the two measures and computed p as p(N1>N2)=1cdfN1N2(0) and p(N2>N1)=cdfN1N2(0).41 The p-values for comparing the AUCs were computed using the DeLong test.42 Contrary to some previous studies, we do not exclude any potential outliers when computing statistics.

6.2.4. Detailed comparison of surface propagation and Jacobian integration (left point)

A detailed comparison between SP and JI computed using the left point rule is performed because this is the most widely used version of JI. In addition, we observe that none of the other integration schemes yield significantly better results, when compared to the left point rule, see Sec. 6.2.5. Only diagnostic group separation capabilities according to Cohen’s D of the other schemes are tested against SP.

As seen in Table 3, JI generally exhibits higher mean and standard deviation than SP in both AD and NC. Both AUC and Cohen’s D are generally higher for SP over JI. In a head-to-head pairwise comparison, SP shows significant improvement over JI (p<0.01) in separating NC and AD based on whole brain and hippocampal atrophy, both according to Cohen’s D and AUC. The numbers are similar or higher, however not significantly different, in ventricles and medial temporal lobe. After Bonferroni correcting the p-value threshold across ROIs (p<0.05/4), SP still remained significant over JI in both hippocampus and whole brain. Further, after Bonferroni correcting the p-value threshold across ROIs and statistical measures (p<0.05/8), SP remained significant over JI in whole brain.

Table 3

Statistics based on atrophy estimated in different brain structures using SP, JI, and FreeSurfer. Mean and standard deviation are atrophy in %. * indicates that SP has significantly different Cohen’s D or AUC compared to JI after Bonferroni correction ...

6.2.5. Comparison of surface propagation and Jacobian integration using higher order integration schemes

Using higher order integration schemes to compute JI do not yield any improvement, which is reflected in the Cohen’s D values in Table 4. As was the case when using the left point rule to compute JI, SP yields a significantly better separation between the diagnostic groups when atrophy is computed in the whole brain and in the hippocampus. The numbers a similar or slightly higher in medial temporal lobe and ventricles, yet not significantly different.

Table 4

Cohen’s D for separating AD and NC groups using different JI schemes and SP. SP is significantly better than these schemes in hippocampus and whole brain. WB: whole brain, Hip: hippocampus, Vent: ventricles, MTL: medial temporal lobe.

Note, in order to perform a fair comparison between SP and JI using Simpson’s rule, we subdivided each voxel into 33=27 smaller voxels effectively adapting Simpson’s rule to a single voxel. Since the analytical expression of the transformation is available, the deformations are simply read out at subvoxel corners and centers.

6.3. Scan-Rescan Test Using In-House Data

To asses the variability of JI using the left point rule and of SP in measuring volume changes, we apply the same registration framework on an in-house dataset comprising nine pairs of scans from a coffee-table experiment and compute atrophy in the ROIs. We use the voxel-wise version of SP in assessing variance of the estimated volume changes. The motivation of this experiment is that we expect no structural changes to be seen when scanning a subject within a very small interval using the same scanning protocol.

Both JI and SP gave nonsignificant atrophies in the considered ROIs. The uncertainties in atrophy measurements were <1% in all the structures. If we zoom in and compare the variance in atrophy over 1000 random regions of each image, each sized 203 voxels, they differed. SP had a lower atrophy standard deviation of 0.79% than JI, which had a 0.97% standard deviation. Further, in a two-sample f-test to test the equality in the variances, SP was significantly different from JI (p<0.01). The mean atrophy over these random regions were negligible and nonsignificant (atrophy=0).

6.4. Application on Deformation Fields from Localized Cross Correlation-Demons

In order to test how SP and JI compare when using deformation fields from a widely known and freely available registration method, we also applied SP and JI (left point) to deformation fields computed using localized cross correlation (LCC)-Demons43 on the ADNI data used above (Table 3). Note, since the deformation field is not available in analytical form, JI is computed using finite differencing and the deformation field is linearly interpolated to the corners of the cube in SP. We conducted experiments using the recommended intrasubject parameters43 (σlaplacian=1.0 and regularization/data term ratio 0.15), and with a less regularized setting (σlaplacian=0.5 and regularization/data term ratio 0.05), were the latter was optimised for class-separation in the hippocampus. Here, we report the highest class-separation in terms of AUC (SP / JI): whole Brain (WB) 0.642 / 0.635, hip 0.59 / 0.58 Vent 0.78 / 0.78, medial temporal lobe (MTL) 0.65 / 0.66. For WB, Vent, and MTL these were achieved used the recommended regularization while the less regularised run gave highest separation for the hippocampus. The only significant (Delong and Delong) difference was achieved for WB, where SP is slightly more separating than JI. Notice that all AUCs are below those reported using deformation fields from the registration method in Sec. 4 (see Table 3).

6.5. Clinical Validity by Comparison to FreeSurfer

In addition to all the above experiments, we also present atrophy numbers from another popular segmentation software, Freesurfer, in order to inspect the validity of the obtained results (see Table 3). The mean atrophies of SP (and of JI) over different brain regions are close to those obtained by both cross-sectional40 and longitudinal FreeSurfer,25 albiet slightly on the lower side. The separation between AD and NC appears to be better than FreeSurfer. We do not perform any pairwise comparison to FreeSurfer, since FreeSurfer does not involve estimating atrophy from a deformation field which is the purpose of this paper.

7. Discussion

We propose SP as a method to estimate atrophy from a deformation field. In SP, the volume of a deformed object is obtained by summing the x coordinate of the surface triangle centroids weighted by the projected areas. We tested both SP and JI using deformations arising from a freeform deformation scheme, which uses NMI as a similarity measure. With SP, the surface-only simplification is not available with SC schemes, since the interior of the voxels need to be subdivided. The numerical properties of SC and SP are identical as validated by the synthetic experiment. For all practical applications, SP was used as a representative of VM methods for comparisons with JI. The main findings of this paper are that (1) synthetic experiments demonstrated that SP/SC is more robust than JI in the presence of moderate to severe noise; (2) SP produced whole brain and medial temporal lobe atrophy numbers that were significantly better than JI at separating between NC and Alzheimer’s disease patients in a standardized dataset from the ADNI; (3) SP displayed significantly less variance at a local patch level when compared to JI using a small, in-house scan-rescan dataset; (4) SP produced disease group atrophy differences comparable to or larger than those obtained using FreeSurfer in the standardized ADNI dataset, demonstrating the validity of the obtained clinical results.

JI is usually obtained by summing the transformations’ Jacobian determinant evaluated at voxel centers (left point rule). Higher order integration schemes like the midpoint, trapezoidal, and Simpson’s rule can also be used to approximate the integral. Since a parametric registration scheme is used, the analytical expression of the transformation is available, analytical derivatives can be directly obtained at the sampling points. However, in the context of nonparametric registration schemes, the partial derivatives in JI need to be approximated with schemes like finite differencing in addition to the approximation in the integral. In case of SP and SC, the volume change information can be directly computed, since only the displacements are required. The synthetic experiments demonstrated that SP and SC are more stable and precise in the presence of noise as compared to JI using various integration schemes. We have derived the order of approximation of SP to be O(h3), where h is the maximum edge length of a surface triangle. This results is applicable for all methods using mesh-propagation for volume change computation, and is the same as the order of approximation of JI using the Trapezoidal integration scheme. Using Simpson’s rule, the order of approximation for JI may become higher, but the synthetic example demonstrates that only in regimes with very low noise are these more accurate than mesh-propagation.

One of the key features of SP is the manner of meshing. Only the surface of a cube is triangulated to estimate the volume. Since the normal of triangles on the interface between two voxels are in opposite directions, they cancel out leaving only the triangles on the outermost faces of voxels on the boundary contributing to the volume of the object. This property is not available with SC because the interior of the voxel needs to be subdivided too. Utilizing only the surface to estimate volume has applications beyond just gain in computational time. For instance, when regularization terms such as penalizing the determinant of Jacobian is used to preserve volume, rigidity is enforced on each voxel of the image. This may prevent larger deformations from being achieved. If the surface version of CP is used, regularity is enforced only on the global volume of a region allowing points inside the region to move more freely.

In the experiment on ADNI data, baseline and 12-month follow-up images were registered using the registration scheme described in Sec. 4. Regional volume change was then computed from the resulting deformation field using SP and JI formulated using the various integration schemes. Pairwise head-to-head comparison with SP was performed with JI and its variations. SP-based measures present a significant improvement over JI-based measures in both hippocampus and whole brain. Figure 6 shows two outliers from Jacobian measures. After we discount for these two outliers, there is a marginal increase in Cohen’s D to 0.64. However, SP was still significantly better than JI. Among the outlier’s seen in Fig. 6, we inspected them to find the reason for an exaggerated estimation of atrophy. We observed a differential bias (i.e., residual bias between the N3 bias corrected baseline and followup) in the images which lead to such an unrealistic atrophy estimation. Despite this, SP appears to measure atrophy closer to the NC mean, whereas the Jacobian measure is implausibly high.

Fig. 6

Scatter plot of JI-based atrophy scores versus SP-based atrophy scores for whole brain (top-left), hippocampus (top-right), ventricles (bottom-left), and medial temporal lobe (bottom-right). The green line is the identity line.

In the brain regions with the most spatially varying deformation field (i.e., the hippocampus and whole brain), JI-based methods suffered the most. Whereas in region with the lowest variance, (i.e., the ventricles), the statistics were nearly identical. This could be explained as follows: (a) Since the registration is a composition of warps, the determinant of Jacobian is a very high-degree polynomial. Involvement of higher order terms tend to amplify any noise in the data with a relatively high variance thus truncating any gain from the use of the higher order terms. This may be one of the reasons why using higher order integration schemes did not yield any statistically significant improvements in the measurements over the traditional left point JI. (b) Homogeneity in the deformation field diminish the contributions from higher order terms making the integration schemes, regardless of their order, perform similar. It can be hypothesized that high regularity in the method and subsequently in the deformation field results in a similar performance of both JI and SP. This is reflected in the near identical results based on ventricular measures, where the deformation field varied the least. Similarly, the two methods give very similar results when based on the Gaussian smoothed deformation fields obtained with LCC-Demons, even though SP is slightly more (and significantly so) separating than JI in the whole brain measurements.

We choose dense sampling of evaluation points, since we intend to have high precision in clinically relevant regions such as hippocampus and cerbral cortex. The other relatively uniform regions were sparsely sampled. This type of sampling can introduce skewed distribution estimates for both the images potentially influencing NMI. It would be possible to compensate for such an effect by applying appropriate weights in the histogram and its gradients, and this is a topic for future work.

8. Conclusion

Much work has gone into improving various aspects of image registration such as image normalization, regularization, similarity measures. However, there has been very little effort in assessing the precision of integration methods that estimate volume changes from the resulting deformation field. Therefore, this study was mainly aimed at comparing JI-based methods to VM-based methods. In the process, we presented surface propagation—a modification of simplex counting with similar numerical properties but less computational complexity. For noisy deformation fields, we saw that JI-based schemes lose precision due to an incomplete representation of volume. In the scan rescan experiments, voxel-wise SP displayed significantly lower variability in measuring no-change than JI in randomly sampled regions. The improvements presented by SP over JI in separating diagnostic groups may also be extrapolated to SC, since SP and SC have the same numerical properties. In addition, the registration-based method together with SP presented in this study, enables a precise quantification of whole brain and regional atrophy based on successive structural MRI scans. The generalization of deformation-based morphometry to include other measures than volume change by the Jacobian determinant is interesting for several applications, e.g., Ref. 44. This is a natural next step for cube-propagation.


We gratefully acknowledge the Laboratory for Computational Neuroimaging at the Athinoula A. Martinos Center for Biomedical Imaging for providing the software (FreeSurfer) used for the segmentations in this paper. Data collection and sharing for this project was funded by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904). ADNI is funded by the National Institute on Aging, the NIBIB, and through generous contributions from the following: Alzheimer’s Association, Alzheimer’s Drug Discovery Foundation, BioClinica, Inc., Biogen Idec Inc., Bristol-Myers Squibb Company, Eisai Inc., Elan Pharmaceuticals, Inc., Eli Lilly and Company, F. Hoffmann-La Roche Ltd. and its affiliated company Genentech, Inc., GE Healthcare, Innogenetics, N.V., IXICO Ltd., Janssen Alzheimer Immunotherapy Research and Development, LLC, Johnson and Johnson Pharmaceutical Research and Development LLC, Medpace, Inc., Merck and Co., Inc., Meso Scale Diagnostics, LLC, NeuroRx Research, Novartis Pharmaceuticals Corporation, Pfizer Inc., Piramal Imaging, Servier, Synarc Inc., and Takeda Pharmaceutical Company. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health ( The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Disease Cooperative Study at the University of California. This research was also supported by NIH grants P30 AG010129, K01 AG030514, and the Dana Foundation.

This study was funded by the Danish Research Foundation (Den Danske Forskningsfond), the Danish National Advanced Technology Foundation (Højteknologifonden) and Biomediq A/S. Data used in preparation of this article were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database ( As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found at:



Akshay Pai is a postdoctoral researcher at the University of Copenhagen. He recieved his PhD in medical image analysis from the University of Copenhagen in 2015. His research interests include automatic and quantitative modeling of neurodegenerative diseases.


Biographies for the other authors are not available.


1. McKhann G., et al. , “Clinical diagnosis of Alzheimer’s disease: report of the NINCDS-ADRDA work group under the auspices of Department of Health and Human Services task force on Alzheimer’s disease,” Neurology 34, 939–944 (1984). [PubMed]
2. Scahill R., et al. , “A longitudinal study of brain volume changes in normal aging using serial registered magnetic resonance imaging,” Arch. Neurol. 60, 989–994 (2003). [PubMed]
3. Hampel H., et al. , “Biomarkers for Alzheimer’s disease: academic, industry and regulatory perspectives,” Nat. Rev. Drug Discovery 9, 560–574 (2010). [PubMed]
4. Fox N., et al. , “Effects of Abeta immunization (AN1792) on MRI measures of cerebral volume in Alzheimer disease,” Neurology 64, 1563–1572 (2005). [PubMed]
5. Frisoni G. B., et al. , “The clinical use of structural MRI in Alzheimer disease,” Nat. Rev. Neurol. 6, 67–77 (2010). [PMC free article] [PubMed]
6. Boyes R., et al. , “Cerebral atrophy measurements using Jacobian integration: comparison with the boundary shift integral,” NeuroImage 32, 159–169 (2006). [PubMed]
7. Hua X., et al. , “Tensor-based morphometry as a neuroimaging biomarker for Alzheimer’s disease: an MRI study of 676 AD, MCI, and Normal subjects,” NeuroImage 43, 458–469 (2008). [PMC free article] [PubMed]
8. Hua X., et al. , “Unbiased tensor-based morphometry: improved robustness and sample size estimates for Alzheimer’s disease clinical trials,” NeuroImage 66, 648–661 (2013). [PMC free article] [PubMed]
9. Thompson P., et al. , “Growth patterns in the developing brain detected by using continuum mechanical tensor maps,” Nature 404(6774), 190–193 (2000). [PubMed]
10. Holland D., Dale A. M., and the Alzheimer’s Disease Neuroimaging Initiative, “Nonlinear registration of longitudinal images and measurement of change in regions of interest,” Med. Image Anal. 15, 489–497 (2011). [PMC free article] [PubMed]
11. Grandy J., “Efficient computation of volume of hexahedral cells,” Technical Report UCRL-ID-128886, Lawrence Livermore National Lab, California: (1997).
12. Burger M., Ruthotto L., Modersitzki J., “A hyperelastic regularization energy for image registration,” SIAM J. Sci. Comput. 35(1), B132–B148 (2011).
13. Carr H., Moller T., Snoeyink J., “Artifacts caused by simplicial subdivision,” IEEE Trans. Visual Comput. Graphics 12(2), 231–242 (2006). [PubMed]
14. Yushkevich P. A., et al. , “Bias in estimation of hippocampal atrophy using deformation-based morphometry arises from asymmetric global normalization: an illustration in ADNI 3T MRI data,” NeuroImage 50, 434–445 (2010). [PMC free article] [PubMed]
15. Pai A., et al. , “Evaluation of bias in brain atrophy estimation,” in Novel Imaging Biomarkers in Alzheimer’s Disease, MICCAI workshop, (2012).
16. Pai A., et al. , “Cube propagation for focal brain atrophy estimation,” in IEEE Symp. on Biomedical Imaging, (2013).
17. Rueckert D., et al. , “Non-rigid registration using free-form deformations: application to breast MR images,” IEEE Trans. Med. Imaging 18(8), 712–721 (1999). [PubMed]
18. Loeckx D., et al. , “Nonrigid image registration using conditional mutual information,” IEEE Trans. Med. Imaging 29(1), 19–29 (2010). [PubMed]
19. Vercauteren T., et al. , “Non-parametric diffeomorphic image registration with the demons algorithm,” in Medical Image Computing and Computer Assisted Intervention (MICCAI), Vol. 4792, pp. 319–326 (2007). [PubMed]
20. Beg M., et al. , “Computing large deformation metric mappings via geodesic flows of diffeomorphisms,” Int. J. Comput. Vision 61, 139–157 (2005).
21. Lorenzi M., et al. , “4D registration of serial brain MRI images: a robust measure of changes applied to Alzheimer’s disease,” in Miccai Workshop on Spatio-Temporal Image Analysis for Longitudinal and Time-Series Image Data, Beijing, China (2010).
22. Christensen G., Rabitt R., Miller M., “Deformable templates using large deformation kinematics,” IEEE Trans. Med. Imaging 5(10), 1435–1447 (1996). [PubMed]
23. Freeborough P., Fox N., “The boundary shift integral: an accurate and robust measure of cerebral volume changes from registered repeat MRI,” IEEE Trans. Med. Imaging 16, 623–629 (1997). [PubMed]
24. Dale A., Fischl B., Sereno M., “Cortical surface-based analysis. i. Segmentation and surface reconstruction,” NeuroImage 9, 179–194 (1999). [PubMed]
25. Reuter M., et al. , “Within-subject template estimation for unbiased longitudinal image analysis,” NeuroImage 61(4), 1402–1418 (2012). [PMC free article] [PubMed]
26. Smith S. M., et al. , “Normalized accurate measurement of longitudinal brain change,” J. Comput. Assisted Tomogr. 25(3), 466 (2001). [PubMed]
27. Smith S. M., “Fast robust automated brain extraction,” Hum. Brain Mapp. 17(3), 143–155 (2002). [PubMed]
28. Anderson V. M., et al. , “Gray matter atrophy rate as a marker of disease progression in {AD},” Neurobiol. Aging 33(7), 1194–1202 (2012). [PMC free article] [PubMed]
29. Nakamura K., et al. , “Jacobian integration method increases the statistical power to measure gray matter atrophy in multiple sclerosis,” NeuroImage: Clin. 4, 10–17 (2014). [PMC free article] [PubMed]
30. Cash D. M., et al. , “Assessing atrophy measurement techniques in dementia: results from the miriad atrophy challenge,” NeuroImage 123, 149–164 (2015). [PMC free article] [PubMed]
31. Darkner S., Sporring J., “Locally orderless registration,” IEEE Trans. Pattern Anal. Mach. Intell. 35(6), 1437–1450 (2012)(PrePrint). [PubMed]
32. Lillholm M., et al. , “Evaluation of WBAA with registration-based cube propagation for brain atrophy quantification,” in Proc. of the 99th Annual Meeting of the Radiological Society of North America (RSNA ‘13), (2013).
33. Liu D. C., et al. , “On the limited memory BFGS method for large scale optimization,” Math. Program. 45, 503–528 (1989).
35. Shackleford J. A., et al. , “Analytic regularization of uniform cubic B-spline deformation fields,” Lect. Notes Comput. Sci. 7511, 122–129 (2012). [PubMed]
36. Sporring J., “Truncation error for simplex propagation,” Technical Report, ISSN: 0107-8283, University of Copenhagen, Copenhagen: (2013).
39. Wyman B., et al. , “Standardization of analysis sets for reporting results from ADNI MRI data,” Alzheimer’s and Dementia 9(3), 332–337 (2012). [PMC free article] [PubMed]
40. Fischl B., et al. , “Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain,” Neuron 33, 341–355 (2002). [PubMed]
41. Holland D., et al. , “Unbiased comparison of sample size estimates from longitudinal structural measures in ADNI,” Hum. Brain Mapp. 33(11), 2586–2602 (2012). [PMC free article] [PubMed]
42. DeLong E. R., DeLong D. M., Clarke-Pearson D. L., “Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach,” Biometrics 44, 837–845 (1988). [PubMed]
43. Lorenzi M., et al. , “LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm,” NeuroImage 81, 470–483 (2013). [PubMed]
44. Lepore N., et al. , “Generalized tensor-based morphometry of HIV/AIDS using multivariate statistics on deformation tensors,” IEEE Trans. Med. Imaging 27(1), 129–141 (2008). [PMC free article] [PubMed]

Articles from Journal of Medical Imaging are provided here courtesy of Society of Photo-Optical Instrumentation Engineers