Temporal filtering of longitudinal brain magnetic resonance images for consistent segmentation.

Longitudinal analysis of magnetic resonance images of the human brain provides knowledge of brain changes during both normal aging as well as the progression of many diseases. Previous longitudinal segmentation methods have either ignored temporal information or have incorporated temporal consistency constraints within the algorithm. In this work, we assume that some anatomical brain changes can be explained by temporal transitions in image intensities. Once the images are aligned in the same space, the intensities of each scan at the same voxel constitute a temporal (or 4D) intensity trend at that voxel. Temporal intensity variations due to noise or other artifacts are corrected by a 4D intensity-based filter that smooths the intensity values where appropriate, while preserving real anatomical changes such as atrophy. Here smoothing refers to removal of sudden changes or discontinuities in intensities. Images processed with the 4D filter can be used as a pre-processing step to any segmentation method. We show that such a longitudinal pre-processing step produces robust and consistent longitudinal segmentation results, even when applying 3D segmentation algorithms. We compare with state-of-the-art 4D segmentation algorithms. Specifically, we experimented on three longitudinal datasets containing 4-12 time-points, and showed that the 4D temporal filter is more robust and has more power in distinguishing between healthy subjects and those with dementia, mild cognitive impairment, as well as different phenotypes of multiple sclerosis.

Introduction

Segmentation of brain magnetic resonance (MR) images is an important step in analyzing brain structures. By providing quantitative measures of the shape and size of brain structures, a better understanding of normal aging (Resnick et al., 2003, Thambisetty et al., 2010), as well as diseases such as Alzheimer's (Querbes et al., 2009) and multiple sclerosis (MS) (Shiee et al., 2012) can be gained. Longitudinal segmentation in brain imaging provides additional insights into the dynamics of brain anatomy for monitoring atrophy and other structural changes that may be related to disease progression. For example, longitudinal imaging studies have revealed accelerated brain volume decline in mild cognitive impairment (Driscoll et al., 2009) and accelerated gray matter atrophy in MS progression (Fisher et al., 2008). The goal of this work was to develop an algorithm that improves the accuracy and stability of brain segmentations in longitudinal data.

The primary challenge of a 4D image segmentation method, in contrast to 3D analysis, is ensuring consistency or stability of the results while retaining sensitivity. Longitudinal processing is aimed toward quantifying time-varying changes of a subject. However, the presence of image artifacts and noise can reduce the sensitivity of a 4D segmentation method by overshadowing the time-varying effects. An example is shown in Fig. 1 where four longitudinal T1-weighted SPGR (spoiled gradient recalled) scans of a healthy volunteer are processed with independent 3D Freesurfer (FS) (Dale et al., 1999) and a state-of-the-art 4D segmentation method, longitudinal Freesurfer (or 4D Freesurfer) (Reuter et al., 2012). Each scan is separated by approximately one year. The inconsistency is visually evident in the hard segmentations (orange arrow in Fig. 1 second row), where the cortical gray matter shrinks and grows over time. Although all images are scanned in the same scanner and with the same pulse sequence parameters, small differences in the noise level or intensities across time-points give rise to inconsistency in the 4D segmentations. In comparison, use of our 4D filter produces a more consistent segmentation, where the cortex shrinks gradually, as expected in normal aging (Fig. 1 bottom row). In this paper, we address this issue of segmentation stability in longitudinal imaging studies of aging or diseases where one might expect gradual changes.

Fig. 1

Top row shows SPGR scans of four consecutive time-points of a healthy volunteer (age 68 at T = 1), with each time point being one year apart. Next four rows show magnified views of the cortex on the original images, segmentations from 3D Freesurfer, 4D Freesurfer, and 3D Freesurfer following by our 4D filtering. An arrow shows where the cortex shrinks and grows periodically on both 3D and 4D Freesurfer segmentations, while it is more consistent with the 4D filter.

Several methods have been previously proposed to analyze 4D images in a longitudinally consistent manner. CLASSIC (Xue et al., 2006) is a 4D segmentation algorithm that contains a 4D registration and 4D segmentation loop. First, all 3D time-points are co-registered using a 4D registration (Shen and Davatzikos, 2004, Prastawa et al., 2012, Kim et al., 2013) to account for the anatomical changes and then are segmented by a 4D version of a fuzzy C-means algorithm (Pham, 2001) to account for temporal smoothness in segmentation. Here, smoothness refers to lack of sudden changes or discontinuities in intensities across time. A 4D segmentation pipeline was also proposed in (Dai et al., 2012), where temporal constraints on cortical thickness are enforced to obtain longitudinally smooth cortical thickness measures. In 4D Freesurfer (Reuter et al., 2012), individual time-points are first segmented with 3D Freesurfer (Dale et al., 1999), a mean template is created, and then the time-points are segmented using the mean template as a target in an unbiased fashion. However, 4D FreeSurfer does not take advantage of the fact that much of the brain often remains unchanged over time, which could be used to improve robustness to noise and other artifacts.

In this paper, we propose a patch based 4D temporal filtering algorithm as a pre-processing step to any segmentation method so as to obtain temporally consistent longitudinal brain segmentations. A patch is a small 3D subimage (e.g., 3 × 3 × 3 voxels) centered at a voxel of interest. For a longitudinal dataset where the scans are of similar contrast, time-points are first co-registered using a rigid transformation. Once the images are rigidly registered, we assume that some changes in anatomy can be modeled by smooth temporal changes in image intensities. Patches are used to model the temporal change in order to take advantage of the contextual information within the neighborhood around a voxel (Roy et al., 2013a, Roy et al., 2011, Roy et al., 2014). We then automatically distinguish between two types of patches, (1) patches that show gradual temporal atrophy, but can be corrupted by noise and minor intensity variations, and (2) patches that do not follow any gradual trend. The intensities of only the first type of patches are filtered using an auto-regressive model of first order (AR(1)), under the assumption that the patch intensities of the tth time-point can be predicted from the (t - 1)th time-point. The intensities of the second type of patches are not changed. Once the images are temporally filtered, they can be segmented using any segmentation method. We show that by using temporal filtering as a pre-processing step, stable and robust 4D segmentations are obtained. We have presented preliminary elements of this work in conference form (Roy et al., 2013c, Roy et al., 2013b).

The paper is organized as follows. First, the motivation for intensity based regression of normal anatomical changes is proposed in Section 2.2. Then the auto-regressive model on the intensities is described in Section 2.3. Next, we evaluate the accuracy of the method on simulated longitudinal atrophies and show the improvement in robustness on test–retest scans of real data in Section 3.2. We also show that 4D filtering improves discrimination of dementia or mild cognitive impairment when compared against healthy groups using a variety of segmentation algorithms in 3.3, 3.4. In particular, we show that using 3D FreeSurfer with our proposed approach improves upon 4D Freesurfer especially when the number of time-points is large (Section 3.4). Finally, we explore in Section 3.5 the effect of 4D filtering in finding GM atrophy on subjects with multiple sclerosis (MS).

Materials and methods

Data sets

We experimented on 6 sets of data, listed in Table 1. First, we estimated the parameter f based on the image contrast and scanners, as noted in Section 2.3. We simulated atrophy on a dataset containing two subjects (denoted by Sim-2), one with both SPGR (GE 1.5T, 0.94 × 0.94 × 1.5 mm3, T = 6.9 ms, T = 3.4 ms, flip angle = 8∘) and MPRAGE (GE 1.5T, 0.94 × 0.94 × 1.5 mm3, T = 35 ms, T = 5 ms, flip angle = 45∘) and one with MPRAGE from a Philips 3T (0.82 × 0.82 × 1.1 mm3, T = 10.21 ms, T = 6 ms, flip angle 8∘) scanner. These two subjects were chosen from the BLSA and MS dataset, respectively (explained later in this section).

Table 1

The table describes the details of the datasets used for all the experiments.

	Name	Contrast	#	Scanner	Resolution (mm3)	TR (ms)	TE (ms)	Flip angle
Simulation	Sim-2	SPGR &MPRAGE	2	GE 1.5T & Philips 3T	0.94 × 0.94 × 1.5 & 0.82 × 0.82 × 1.1	6.9 & 35	3.4 & 5.0	8∘ & 45∘
Test–retest	TR-10	MPRAGE	1	Siemens 3T	0.82 × 0.82 × 1.17	10.2	6	8∘
CoRR-29	SPGR	29	GE 3T	1.0 × 1.0 × 1.0	8.06	N/A	8∘
OASIS	OA-49	MPRAGE	49	Siemens 1.5T	1.0 × 1.0 × 1.25	9.7	4	10∘
BLSA	BL-39	SPGR	49	GE 1.5T	0.94 × 0.94 × 1.5	35	5	45∘
MS	MS-59	MPRAGE	59	Philips 3T	0.82 × 0.82 × 1.1	10.21	6	8∘

To estimate the test–retest stability of the method, we use two datasets. For the first set, MPRAGE scans (0.82 × 0.82 × 1.17 mm3, T = 10.2 ms, T = 6 ms, flip angle 8∘) from a 3T Siemens scanner (denoted by TR-10) were acquired every week for 10 consecutive weeks on a healthy male volunteer. Precautions were taken to have identical scanning and patient conditions every time, e.g., scans were done at the same time of the day, and the volunteer was allowed the same amount of sleep and fluid intake on the day of scanning. The second dataset consists of SPGR scans (GE 3T, 1.0 mm3, T = 8.06 ms, T = 450 ms, flip angle 8∘) of 29 healthy subjects, each subject having 10 scans over a month, separated by 3 days. This dataset (denoted by CoRR-29) is obtained from the Consortium for Reliability and Reproducibility (CoRR) database (Zuo et al., 2014), specifically the HNU dataset (Hangzhou Normal University) (http://dx.doi.org/10.15387/fcp_indi.corr.hnu1). Note that although the scans were available for 30 subjects, we discarded one due to an incomplete field of view.

To demonstrate the discriminative power in longitudinal analysis versus cross-sectional analysis, 49 subjects with 3 - 5 time-points (age range 60 - 92) were selected from the OASIS (Open Access Series of Imaging Studies) (Marcus et al., 2007) longitudinal database (dataset denoted by OA-49). Each visit was separated by approximately one year and there was an average of three visits per subject. For every subject and every visit, there were 3–4 repeat MPRAGE scans (1 × 1 × 1.25 mm3, T = 9.7 ms, T = 4 ms) from a Siemens 1.5T scanner. The repeat scans are co-registered and averaged to improve signal-to-noise ratio. Among the 49 subjects, 15 of them were diagnosed with dementia, and the other 34 were characterized as non-demented throughout the span of the study. Note that there were also 14 patients in the database who were initially non-demented, but were diagnosed with dementia at later time-points. We did not include them in our subset.

To evaluate our method when applied to a greater number of time-points, we experimented on 39 subjects from the BLSA (Baltimore Longitudinal Study of Aging) database (Resnick et al., 2000), with 15 of them diagnosed with mild cognitive impairment (MCI) at all time-points (dataset denoted by BL-39). The other 24 subjects are healthy controls. Each subject contains 4–11 visits (average 9 visits per subject), each separated by approximately one year. SPGR images were acquired axially (0.94 × 0.94 × 1.5 mm3, T = 35 ms, T = 5 ms, flip angle = 45∘) for each visit on a GE 1.5T scanner.

The 4D filtering was also applied on a set of MS patients, where we explored the progression of atrophy on different phenotypes of MS. The data set includes 59 patients with MS (dataset denoted by MS-59), with each patient having 3–8 visits (average 4), each visit separated by a year. The average age of the participants was 44 years (range 22–67), with an average disease duration 9 years. Among the 59 patients, 22 were diagnosed with relapsing remitting MS (RRMS), 15 with primary progressive (PPMS), and 22 with secondary progressive MS (SPMS). All subjects had T1-w MPRAGE (0.82 × 0.82 × 1.1 mm3) and FLAIR (0.82 × 0.82 × 2.2 mm3) scans acquired on a 3T Philips scanner (TR/TE = 10.21/6 ms, flip angle 8∘).

To demonstrate the flexibility of our approach, we apply it as a preprocessing step to several different segmentation algorithms. These algorithms each use different underlying methodologies. The Freesurfer (Dale et al., 1999) algorithm employs a Markov random field model in combination with statistical atlases and deformable surfaces. Because it includes both a 3D and 4D implementation, a natural comparison that we focus on is the combination of temporal filtering and 3D Freesurfer against 4D Freesurfer (Reuter et al., 2012). Other algorithms that we tested include Atropos (Avants et al., 2011), which is based on an expectation–maximization approach, FIRST (Patenaude et al., 2011), which is based on active appearance models, TOADS (Bazin and Pham, 2008), which uses fuzzy clustering with topological constraints, and MALP-EM (Ledig et al., 2015), which combines multi-atlas label fusion with expectation–maximization.

Motivation

Most existing 4D segmentation methods involve a nonlinear 4D registration step where either the individual images (Dale et al., 1999) or their segmentations (Xue et al., 2006) are deformably registered to a common template. Thus the temporal change in anatomy is modeled by a geometric “morphing” of the images. However, it has been shown that the same change in anatomy can also be modeled by an intensity change (Miller and Younes, 2001).

An example of how longitudinal anatomical changes can be explained by intensity changes is shown in Fig. 2(a), where three time-points (years 1, 2, and 6) of a healthy volunteer are shown. Each of the time-points is rigidly registered (Jenkinson and Smith, 2001) to a reference space (e.g., an atlas or an average of all time points). Each of the images are scaled so that the mode of the white matter (WM) intensity is at unity. The mode is automatically found by fitting a smooth kernel estimator to the image histogram (Pham and Prince, 1999). Fig. 2(b) shows the temporal intensity profiles of three voxels, one each in deep WM, ventricles, and on the ventricle-WM boundary. As the subject is healthy, the anatomical change associated with normal aging, e.g., enlargement of the ventricles, is manifested in the monotonic decrease in intensity of the voxel on the WM-ventricle boundary (blue line), with perturbations due to image noise. Also, the voxel inside WM (or ventricle) remains WM (or ventricle) for all the time-points, thus its intensity remains in the realm of the corresponding tissue intensities. A similar trend can also be observed for a patch as well. Fig. 2(c) shows the same 5 × 5 patch (in 2D) over the three time-points, demonstrating the shifting boundary of the ventricles. We propose that after rigid registration, the trend in some patches can be attributed to various factors like noise (e.g., green and magenta lines in Fig. 2(b)), actual gradual anatomical changes (e.g., blue line in Fig. 2(b)), or small registration errors. To remove the noise and other artifacts from those patches, we propose to fit the longitudinal intensities to a smooth first order auto-regressive model.

Fig. 2

The top row shows SPGR scans of three time-points of a normal subject having 6 longitudinal scans, separated by a year. The second row shows the intensity profiles of three voxels, one inside ventricle (magenta), one in deep WM (green), and one on the ventricle-WM boundary (blue). The AR(1) fit of the blue line is shown in red. The bottom row shows the intensities of 5 × 5 patches around the blue voxel over the 3 time-points. The leftmost and rightmost voxels in the patches represent WM and ventricle voxels that remain unchanged over the time course, while the voxels in the middle columns of the patches are generally decreasing in intensity, indicating enlargement of the ventricle. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Algorithm

Image patches are defined as p × q × r 3D subimages associated with each voxel of the image. Patches are typically small and centered on the voxel of interest—e.g., p = q = r = 3. It is convenient for writing out the mathematics to describe a patch as the 1D vector of size d × 1, where d = pqr. The voxels within a patch are always ordered in the same way, using a consistent rasterization, to create this 1D representation.

We assume that there are T time-points available for a subject, {s1, … , s}. Each of the time-points {s1, … , s} is rigidly registered to a reference template. Each scan is skull-stripped with the same mask (generated from either the baseline or a mean template) and corrected for any intensity inhomogeneities independently. They are also scaled so that the mode of the WM intensities of each image is unity (Pham and Prince, 1999, Shinohara et al., 2014, Roy et al., 2013d). The images are also assumed to be scanned with the same acquisition protocol. Each image is first decomposed into d × 1 patches, y( , i = 1 , … , N , t = 1 , … , T, where N is the total number of voxels in each image. Since the images are rigidly registered, the trend at the ith voxel is obtained from the collection of patches { y(1), … , y(}.

At every voxel, an observed image patch y( is assumed to consist of a “true patch” x( and additive noise, expressed as.

where the observed patch is a noisy perturbation of the true patch at each voxel. Here ( accounts for both image noise as well as small intensity variations. We assume that there are two types of patches: patches that follow a gradual trend and patches that do not. Our assumption is that aging or neurodegenerative processes typically cause gradual shrinking of brain structures. These constitute the first type of patch, while other patches fall into the second type. We discuss this in greater detail in the remainder of the section.

As shown in Fig. 2, there are some patches for which the atrophy results in smooth intensity change, such as ventricular enlargement. We assume that the gradual anatomical changes can be modeled by smooth changes in intensities; therefore, for those patches, x( are modeled as an AR(1) regression,

Assuming a is the ith true patch for the baseline image (t = 1), then x( can be re-written as,

The d × d matrix M is positive definite and contains the parameters of the AR model. The exponent on M signifies multiple products with itself ((t - 1) power). We note that the intensities are a nonlinear function of the time t; therefore the model has the capability of detecting nonlinear changes in the intensity profile. However, because M itself is not a function of t, x( must follow a monotonic trend. The first order model is based on the assumption that the serial images of a subject are obtained at approximately uniform intervals, which is true for all subjects in our database. However, when a subject is scanned at highly nonuniform intervals, higher order models can be used in the same framework.

There are some patches that do not follow monotonic change in intensities. Such patches might include acute lesions in MS subjects which may appear and disappear in successive time-points, the presence of a new tissue-boundary in developing brains that is not possible to predict from the previous time-point, or tissue variations due to hydration level changes (Nakamura et al., 2015). We distinguish between these two types of tissue patches. The patches that do not follow smooth, monotonic intenity changes are where the fitting error of the AR regression (Eq. (3)) is sufficiently large. Therefore, for such patches, the original intensities must be preserved. Eq. (3) is modified to satisfy this criteria for all patches,

Here, z( is an estimated patch accounting for both gradual and non-gradual atrophies. The scalar w( is a weight that determines if the trend at the ith patch is due to a gradual anatomical change (such as growing ventricles in normal aging, w( → 1) or some atrophy that is not gradual. In the presence of such non-gradual atrophies, the intensities cannot be smoothed with an AR model, thus w( is set to 0.

The weights w( are computed in a data-dependent manner so as to distinguish between the two types of patches. If the model fitting error ||  y( - x(||2 is too large, it indicates the presence of a non-monotonic atrophy that cannot be explained by smooth change in intensities. Following this notion, w( is defined as,

where f is a smoothing parameter that also acts as a soft noise threshold. Note that w( depends on the deviation of the observed patch y( from the true patch x(, when the patch is assumed to follow a gradual atrophy. If this deviation ||(|| = ||  y( - x(|| ≫ f, w( ≈ 0, then the patch is largely unaffected by the smoothing filter, i.e., x( ≈ y(. Estimation of f is described in Section 3.1.

Based on the data y(, we first estimate the weights w( and AR(1) parameters M and a. Then the estimated patches z( are found from Eq. (4) using these estimates. For mathematical simplicity, we make the assumption that M is a positive definite diagonal matrix, M = diag{m, … , m}, and thus M = diag{m, … , m}. Then the filtered intensities x( are found for the ith patch by minimizing the norm of the error with respect to M and a. The total sum of squared error of the estimated patches is given by

Note that although we employ an autoregressive model of the first order, it is evident from Eq. (6) that all time-points are used to estimate the model parameters.

Estimates of M and a are obtained by differentiating ETotal with respect to these variables and setting the gradients to zero. Under the diagonal assumption of M, the update equations are as follows,

where and denote the component of y( and a, respectively, . Initializing and a = y(1), Eqs. (7), (8) are repeated until the values of and a converge. The converged values for M, a, and a chosen value for f (see Section 3.1) are used to provide the estimated patches z( based on Eq. (4). It is noted that the computation in Eqs. (7), (8) is simplified by the diagonal assumption of M, and they are solved for each of the patch dimensions separately. Nevertheless, the patch-based error ||(|| introduces contextual information in the computation of M. The diagonal M can be replaced by an arbitrary d × d positive definite matrix M. However, the number of time-points is usually much less than the patch dimension (i.e. T ≪ d), introducing instability in the computation of an arbitrary positive definite matrix.

Fig. 3 shows four out of ten time-points of a healthy volunteer, each separated by two years, and the corresponding m for the center voxel of each 3 × 3 × 3 patch in the image (d = 27). Voxels that remain WM throughout the time span (deep WM denoted by a white arrow), have a m ≈ 1 indicating their intensity remains stable over the time series. Enlargement of the ventricles (black arrow) is also evident, resulting in m < 1 shown in dark blue, which indicates a decrease in intensity over time, corresponding to the transition from WM voxels to ventricle voxels. Some cortical thinning is also observed with m < 1 near deep sulci (magenta and white arrows), where GM voxels turn into CSF voxels. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

(a) Four time-points of a subject before and after 4D processing, and (b) a map of m using 3 × 3 × 3 patches after temporal filtering (d = 27). The white arrow points to a voxel inside deep WM with m ≈ 1, indicating no gradual change in intensities since the voxel remains as WM over the time-period. The magenta arrow indicates a voxel on the cortex with m > 1, indicating a gradual increase in intensities (GM to WM), since the cortex shrinks over time. The black arrow indicates a voxel near the ventricle-WM boundary, where m < 1. It indicates a gradual increase in ventricle size, where the voxel changes from WM to ventricle over time.

The parameter f acts as a soft threshold for the AR model fitting to the temporal intensities of a patch. By design, it can also be a user input to the algorithm. Larger f implies less error (||(||) tolerance, indicating filtered intensities should follow the autoregressive model more strictly. Ideally, f should vary from patch to patch based on the tissue type and amount of expected atrophy at a patch. For example, with a ventricular CSF patch, the inherent signal-to-noise ratio is lower on T1-weighted images (Sijbers et al., 1998) than a pure WM patch, indicating that the choice of f should be higher for CSF patches than WM. The choice of f should also depend on the acquisition protocol of the image, e.g., SPGR vs MPRAGE and the scanner. However, in this paper, we only address the variation of f with respect to the image contrast and scanner, and choose a global f based on a simulation study on an SPGR and MPRAGE images of different scanners in Section 3.1.

Results

The algorithm was implemented in MATLAB. We used 3 × 3 × 3 patches, i.e., d = 27, for all our experiments. The processing time is about 2 hours on a 12-core 2.92GHz Intel processor for an 11 time-point data set with images of 1 mm3 resolution, requiring about 2-3GB of memory. 3D and 4D Freesurfer require about 24 and 36 h, respectively on the same machine. Before temporal filtering, all images were pre-processed by skull-stripping (Carass et al., 2011), registration to MNI atlas (Mazziotta et al., 1995) (www.mristudio.org), resampling to 1 mm3 isotropic resolution, and inhomogeneity correction using N4 (Tustison et al., 2010).

Simulation study

To estimate the parameter f in Eq. (5), we used an atrophy simulation algorithm (Karacali and Davatzikos, 2006) on the Sim-2 data (Section 2.1) to simulate anatomical changes in the brain. One subject with both SPGR and MPRAGE scans (GE 1.5T) from the BL-39 and one subject with MPRAGE from MS-59 dataset (Philips 3T) were used for this experiment. The simulation algorithm takes an MR image and its 3-class hard segmentation (CSF, GM, and WM) and generates an atrophied image around designated foci with a given radius of atrophy. We used decreasing radii of atrophy (20 to 15 mm) to simulate six time-points where the ventricles grow over time, simulating a common occurrence in late aging. A 3-class fuzzy c-means (Bezdek, 1980) (FCM) based segmentation was used as the input. The six time-points were filtered using varying f ∈ (0, 1) and the resultant images were again individually segmented using FCM. The same experiment was repeated for all of the three acquisitions. The misclassification rates in both cases, averaged over three classes, are obtained for each f. Misclassification rate is defined as the ratio of average number of misclassified voxels (over 6 time-points) and the average number of voxels involved in the simulated atrophy.

As mentioned earlier in Section 2.3, a very small f will not have any effect on the filtering, as smaller f imparts no change in intensities (as f → 0, w( → 0 in Eq. (5)). With increasing f, the temporal intensities are smoothed more and more, thereby becoming more robust to noise, while at the same time smoothing some actual tissue changes. We obtained these optimal values: fSPGR1.5T = 0.31, fMPRAGE1.5T = 0.21, fMPRAGE3T = 0.19, at less than 2% misclassification error. Since the f values for MPRAGE on both scanners are very close, we used fMPRAGE = 0.21 for all the experiments with MPRAGE images described later. A larger value of fSPGR compared to fMPRAGE arises from the observation that the contrast to noise ratio in MPRAGE images is typically higher than in SPGR, indicating a smaller noise threshold is needed for detection of gradual atrophy.

Robustness on test–retest scans

In this section, we evaluate the robustness of the 4D filtering using test–retest scans of healthy volunteers, i.e. TR-10 and CoRR-29 data (Section 2.1). Ideally for a healthy subject, the tissue volumes are not expected to change much within a span of ten weeks or a month. For TR-10 dataset, the volumes of cortical GM, cerebral WM, and ventricles are plotted in Fig. 4 for individual 3D Freesurfer segmentation of the original images (blue triangles), 4D Freesurfer (magenta spheres), and 3D Freesurfer segmentation of temporally filtered images (green stars). Similar to a previously reported study on multi-center test–retest reliability of Freesurfer (Jovicich et al., 2013), 4D FS produces almost half the variability compared to 3D FS on GM and WM, while there is very little difference between the two on ventricles. We also observe that 4D Freesurfer has more variability than the 4D filtering method. The difference is primarily emphasized in ventricle volume, where the coefficient of variation is 1.67% for 4D Freesurfer, while it is 0.53% for our filtering method, shown in Table 2. As the ventricle volume is the most robust statistic for scans spanning over only ten weeks, a significantly smaller coefficient of variation for our temporal filtering indicates greater consistency of the segmentation. A similar decrease is shown for the other tissues. We note that 4D Freesurfer has a consistent bias from individual Freesurfer segmentations with larger ventricles or smaller WM. This is likely due to the difference in the transformation spaces of the methods.

Fig. 4

Ventricle, WM, and cortical GM volumes (in mm3) are plotted w.r.t. the time-points for a healthy subject, scanned weekly for ten consecutive weeks.

Table 2

Coefficient of variation (in percent) of volume changes in ten consecutive weeks for a healthy subject (TR-10 dataset). Bold indicates lowest coefficient compared to other methods.

	Ventricle	Cortical GM	Subcortical GM	Cerebral WM
3D Freesurfer	1.529	1.756	1.523	1.332
4D Freesurfer	1.669	0.757	0.741	0.557
4D Filter + 3D Freesurfer	0.532	0.568	0.589	0.380

A similar analysis was carried out for the CoRR-29 dataset, with 29 healthy volunteers having 10 scans spanned over a month. Table 3 shows the mean and standard deviations of coefficients of variations for the three methods. 4D Filter followed by 3D Freesurfer shows the lowest average coefficient of variation (p < 0.01 using Wilcoxon signed rank test) among the three on ventricles, cortical GM and WM, indicating significantly improved segmentation stability.

Table 3

Mean ± standard deviations of coefficient of variation (in percent) of volume changes in ten consecutive scans are shown for 29 healthy subjects (CoRR-29 dataset). Bold indicates significantly lowest (p < 0.01) coefficient compared to other methods.

	Ventricle	Cortical GM	Subcortical GM	Cerebral WM
3D Freesurfer	3.48 ± 1.13	2.41 ± 1.21	2.07 ± 0.54	1.15 ± 0.33
4D Freesurfer	2.22 ± 0.60	1.40 ± 0.22	0.99 ± 0.33	1.00 ± 0.17
4D Filter ± 3D Freesurfer	1.24 ± 0.49	1.25 ± 0.57	0.98 ± 0.39	0.57 ± 0.22

To estimate the stability of the algorithm with respect to the number of time-points, we also created an augmented set of the TR-10 dataset by randomly sampling T images (T = {4, … , 9}) from the 10 weekly scans. The sampling is done 10 times for each T. Therefore, this augmented dataset contains 60 subsets, each subset containing 4 to 9 scans from the weekly test–retest scans. Then the images were processed with the 4D filtering method. Then both the processed and original images are segmented with FAST (Zhang et al., 2001) to find CSF, GM, and WM segmentations. Since we do not expect the tissues to change much within a span of 10 weeks on a healthy control, Dice coefficients (Dice, 1945) for the three tissue classes were used as a similarity metric and were computed on each subset with respect to the baseline of that subset. The mean Dice coefficients for T = 4 , … , 9 on all three tissues were significantly larger (p < 10-5 with Wilcoxon rank-sum test) than the unfiltered images, indicating significant improvement in segmentation stability after 4D filtering. Also the Dice coefficients between T and (T + 1) time-points are not significantly different (p > 0.01) for any T on any tissue, indicating that the 4D filtering is robust to the variation in number of time-points. Median WM Dice for T = 4 was 0.97, while it was 0.98 for T = 9, compared to 0.91 with unfiltered images.

We also tested our approach with a 3-class Expectation–Maximization (EM) based segmentation method Atropos, which also includes 3-D and 4-D implementations (Avants et al., 2011). Scans of the same healthy subject over 10 weeks are segmented using 4D Atropos with a temporal smoothness weight of 0.3. CSF, GM, and WM tissue volumes are plotted in Fig. 5 for 3D Atropos, 4D Atropos, and 3D Atropos on 4D filtered images. Use of filtering shows an improvement in segmentation stability of our temporal filter over a Markov random field temporal smoothness constraint on the segmentations. Clearly, 3D Atropos shows noisy volume trends, and 4D Atropos shows a slight decreasing trend in WM volume, indicating that the WM volume changes by 5%, while there is a 7% increase in CSF. Our 4D filtering with 3D segmentation produces the most stable segmentation without the spurious volume changes.

Fig. 5

MPRAGE scans of a healthy volunteer, scanned for 10 consecutive weeks, are segmented using both 3D and 4D EM based segmentation Atropos (Avants et al., 2011). We compared the 4D segmentation volumes in mm3 (magenta lines) with individual 3D Atropos segmentations of original data (blue lines) and 4D filtered images (green lines). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Atrophy detection on OASIS

As no longitudinal dataset with manual labels are freely available to test sensitivity of the proposed filtering method, we instead demonstrate that the 4D filter improves discrimination power between groups of healthy controls and patients with neurodegeneration on the OA-49 data (Section 2.1). Since the filtering can be applied to any 3D segmentation method, we filtered and segmented both demented and non-demented subjects with TOADS (Bazin and Pham, 2008), FIRST (Patenaude et al., 2011), and 3D Freesurfer (Dale et al., 1999). 4D Freesurfer (Reuter et al., 2012) has already been shown to distinguish between the healthy and dementia groups on this dataset. Here, we replicated that result with 4D Freesurfer and showed that the combination of 4D filter and 3D Freesurfer performs similarly, if not better, in distinguishing the two groups than 4D Freesurfer. Similar improvement in distinguishing power after filtering was shown for TOADS and FIRST as well. TOADS provides a tissue segmentation of the whole brain into multiple labels, cerebellar and cerebral WM, GM, sulcal CSF, thalamus, caudate, putamen, and ventricles. FIRST generates subcortical GM segmentation into multiple labels, such as left and right thalamus, caudate, putamen, globus pallidus, hippocampus, and amygdala. Similar labels are also obtained from Freesurfer as well. To measure longitudinal atrophy in segmentation volumes, we propose absolute percent volume change per year (APV), defined as

where is a robust baseline volume for every subject, obtained from a linear fit of the longitudinal volumes (Reuter et al., 2012). V denotes the volume at tth time-point, and Δt denotes the age difference between time-points (t + 1) and t. Note that in addition to the absolute volume changes, APV also normalizes with respect to the actual age difference between two scans. For subjects with more than two time-points, multiple APV values were obtained. However since the rates of atrophy may increase or decrease with age, we do not compute an average APV for a subject. Instead we use the APV values of all time-points of all subjects to compare the atrophy rates between demented and non-demented groups using non-parametric statistical testing.

Fig. 6 shows barplots of median APV when both original and filtered images are segmented with TOADS and FIRST. Subcortical GM volumes have been shown to be associated with dementia and early onset of Alzheimer's disease (Driscoll et al., 2009, Lehmann et al., 2010, Jovicich et al., 2009, den Heijer et al., 2010). Visually, the variation in APV decreases after 4D filtering (shown by lower interquartile range), indicating more confidence in discriminating demented vs non-demented groups. We hypothesized that APV for the demented group is larger than the non-demented group. A Mann–Whitney U test indicates that ventricle and putamen volumes for TOADS segmentations are indeed significantly (p < 0.01) larger in the demented group after 4D filtering. The caudate and cerebral WM volumes do not show any significant difference after filtering, although the interquartile ranges decrease, indicating decrease in variances. Similarly, both the left and right hippocampus and amygdala volume changes obtained from FIRST are significantly larger (p < 0.001) after 4D filtering. Left and right globus pallidus volume changes are significantly larger (p < 0.05) as well. Although we did not find any significant difference between demented and healthy subjects for thalamus volumes, the p-value decreases (p = 0.067 to 0.054) after filtering. Note that the statistic obtained from the Mann–Whitney U test is closely related to the area under an ROC (receiver operating characteristics) curve for a binary classification. Therefore, a lower p-value indicates a higher U statistic, which in turn indicates higher area under the curve, indicating improved sensitivity and specificity to distinguish between dementia and non-dementia based on atrophy rates.

Fig. 6

Barplots show average percent volume change (APV) per year on both demented and non-demented groups from OASIS database, when segmented by TOADS and FIRST. The range corresponds to 25th and 75th percentiles. Asterisk indicates significantly (p < 0.01) higher APV in demented group compared to the non-demented group. On TOADS segmentations, both ventricles and putamen show significant atrophy in dementia on 4D filtered images compared to un-filtered images. On FIRST segmentations, left and right hippocampus and amygdala show significant atrophy with 4D filtering, which was absent in the original images.

Next we tested our approach when used with Freesurfer. Fig. 7 shows barplots of absolute percent volume change per year for 3D Freesurfer, 4D Freesurfer, and 4D filtering followed by 3D Freesurfer for both demented and non-demented groups. As shown earlier (Reuter et al., 2012), magnitude of volume changes, as well as their variances, are sometimes higher in 3D Freesurfer than 4D Freesurfer, e.g., left/right thalamus, putamen, and amygdala. Significant increase in atrophy was observed in the demented group on the left/right hippocampus (p < 0.05) for all three methods, which is consistent with earlier results (Reuter et al., 2012). We observed significant increases in rates of volume changes in the demented group in the left/right ventricles and thalamus (p < 0.05) on both 4D Freesurfer and 4D filtering followed by 3D Freesurfer, but not on 3D Freesurfer. We also observed significant increases in left putamen (p < 0.05) volumes, though there is no significant increase on either 4D or 3D Freesurfer (p > 0.10). Atrophy in the putamen has been shown to be associated with dementia (Moller et al., 2015). Therefore, our 4D filter followed by 3D Freesurfer has similar performance, if not better, as 4D Freesurfer in distinguishing neurodegeneration, with significant reduction in computation time (about 20–24 h of 3D Freesurfer processing).

Fig. 7

Barplots show median of absolute percent volume change per year (APV) on left and right subcortical GM structures, when 49 subjects from OASIS dataset are segmented with 3D and 4D Freesurfer, as well as 4D filtering followed by 3D Freesurfer. The range corresponds to 25th and 75th percentiles. Plus indicates significantly (p < 0.05) higher APV on left/right thalamus and ventricles from demented group observed in both 4D filter and 4D Freesurfer, but not on 3D Freesurfer. Asterisk indicates significant (p < 0.05) atrophy on left putamen in demented group observed on 4D filtering, but not 3D or 4D Freesurfer.

Improvement in atrophy detection on BLSA

While the OASIS dataset has an average of three time-points per subject, we also experimented on the BL-39 dataset (Section 2.1), which has an average of nine time-points per subject separated by one year, to show the improvement in sensitivity and specificity in detecting atrophy when average number of time-points is larger. As before, we segmented the images with 3D and 4D Freesurfer, and 4D filtering followed by 3D Freesurfer. Fig. 8 shows barplots of median values and interquartile ranges of APV for these three segmentations. Visually, magnitude and variations of atrophy rates of 3D Freesurfer are generally larger than 4D segmentations (e.g., left/right pallidum and amygdala), indicating longitudinally unstable segmentations. We observed significantly higher atrophy rates in the MCI group than the healthy volunteers for the left/right hippocampus and left/right ventricles (p < 0.05) on all three methods. However, when the atrophy is more subtle and relatively small, both 3D and 4D Freesurfer failed to detect any significant change. Atrophies in thalamus and caudate are known to be associated with Alzheimer's Disease (Ryan et al., 2013). Significant atrophy in the MCI group was observed on the left/right thalamus (p < 0.01), left caudate, right putamen, and right amygdala (p < 0.05), which was not observed on both 4D and 3D Freesurfer segmentations (shown as asterisk in Fig. 8).

Fig. 8

Barplots show median of absolute percent volume change (APV) on left and right subcortical GM structures, when 39 subjects from BLSA dataset are segmented with 3D and 4D Freesurfer, as well as 3D Freesurfer after our 4D filtering. The range corresponds to 25th and 75th percentiles. Asterisk indicates significantly (p < 0.05) higher APV on MCI group than controls observed in the combination of 4D filtering and 3D Freesurfer, when the APV is not significant in both 4D and 3D Freesurfer.

We also computed absolute percent change in cortical thickness per year. Accelerated cortical thinning has been observed in normal aging (Fjell et al., 2014), Alzheimer's Disease, as well as other neurodegenerative diseases (Mak et al., 2015). In the BLSA dataset, we observed in Fig. 9 that the MCI group showed a higher rate of cortical thinning than the controls on all three methods. However, the rate is not significant in 3D and 4D Freesurfer (p > 0.10 in both cases), but it is significant (p = 0.015) when the images are processed by our approach followed by 3D Freesurfer. Also, the variations in the rates of change of thicknesses are lower (smaller interquartile range) in the 4D filter, indicating a more stable segmentation. Therefore, our method is more sensitive to distinguishing diseased brains from normals than 4D Freesurfer when the number of time-points is large. The difference in performance of 4D Freesurfer can be explained by the absence of any explicit 4D smoothness model in segmentations. Therefore the effect of random temporal noise becomes more evident when the number of time-points becomes larger.

Fig. 9

Barplots show median and interquartile ranges of absolute percent change in cortical thickness, when 39 subjects from BLSA dataset are segmented with 3D and 4D Freesurfer, as well as 4D filtering followed by 3D Freesurfer. Asterisk indicates significantly (p = 0.015) higher thinning in cortex in MCI group (n = 15) compared to controls (n = 24) observed in the combination of 4D filtering and 3D Freesurfer, while it is not significant in both 4D and 3D Freesurfer.

Atrophy in MS dataset

We have shown on two datasets, OA-49 and BL-39 in 3.3, 3.4, that our 4D filtering has superior sensitivity in distinguishing healthy subjects from patients with dementia and MCI, especially when the number of time-points is large. In this section, we explore the volumetric changes of GM in different phenotypes of MS on the MS-59 data (Section 2.1). In the previous experiments, we showed the advantage of the 4D filter with intensity based segmentation methods (FIRST, TOADS, and Freesurfer). Here we use a registration based label fusion algorithm MALP-EM (Ledig et al., 2015) after with and without the filtering as a pre-processing step.

Since the 4D intensity model does not account for lesions, the WM lesions in MPRAGE images are first segmented (Shiee et al., 2009), and then inpainted with WM intensities (Battaglini et al., 2012). The inpainted images are filtered with the 4D filtering, followed by segmentations with a recent multi-atlas label fusion method MALP-EM (multi-atlas label propagation using expectation maximization). MALP-EM produces whole brain labeling by first registering multiple atlases to a subject, then transferring the corresponding atlas labels to subject space, and combining the labelmaps in subject space via expectation maximization (EM). In this section, we show that the 4D filter in conjunction with the label fusion algorithm can detect changes in atrophy between different phenotypes of MS, compared to unfiltered images. The results are validated by corroborating with previous findings.

Increased decline in total brain and GM volume has been observed in MS patients (Battaglini et al., 2009) compared to controls, while atrophies in putamen, thalamus (Eshaghi et al., 2014), cerebellum, and ventricles (Ramasamy et al., 2009) are shown to be associated with the progression of the disease in cross-sectional studies. In different types of MS, rates of atrophy in different tissues are dominant in nature. E.g., ventricles, putamen, hippcampus show significant atrophy in RRMS, while caudate nucleus, along with regions of frontal, patietal, temporal and occipital cortex, show significant atrophies in SPMS and PPMS (Pagani et al., 2005). In a cross-sectional study, lower GM volume was observed in RRMS compared to SPMS (Roosendaal et al., 2011). It was also observed that total brain volume change is higher in SPMS compared to both RRMS and PPMS (Stefano et al., 2010). Specifically, brain volume change was shown to be a powerful statistic that can predict (p < 0.001) relapses in RRMS (Horakova et al., 2009). Similarly, significantly higher volume loss was observed in the thalamus, left superior temporal gyrus, and left superior frontal sulcus on SPMS (Riccitelli et al., 2010) compared to PPMS.

Fig. 10(a)–(b) show barplots of median (and interquartile ranges) absolute percent volume changes per year (APV) on the 59 subjects with MS, where the segmentations are obtained from MALP-EM after 4D filtering. Significant atrophy is observed in left/right superior frontal (p < 0.001), left superior occipital, left superior parietal, and right superior temporal (p < 0.05) gyri in both SPMS and PPMS compared to RRMS (cf. Pagani et al., 2005) after filtering (Fig. 10(b)). Significant atrophy in SPMS was also observed in left/right thalamus (p = 0.01) and left superior temporal gyrus (p = 0.05) compared to PPMS (cf. (Riccitelli et al., 2010)). Higher rate of change was observed for total brain volume in SPMS compared to RRMS (p = 0.068) (cf. Stefano et al., 2010), although it is similar to PPMS (p = 0.45). Only right superior occipital gyrus and left thalamus shows significant atrophy in RRMS compared to SPMS (p = 0.009 and 0.013 respectively) with raw images (Fig. 10(a)), indicating that 4D filter has more power in distinguishing between different types of MS. Note that although right superior occipital gyrus is significant in unfiltered images with 4D FS (p = 0.009 RRMS vs SPMS), it becomes non-significant in filtered images (p = 0.285).

Fig. 10

Barplots show median (and interquartile range) of absolute percent volume change per year (APV) on 59 subjects from MS dataset, with label-fusion segmentations on (a) un-filtered images, and (b) 4D filtered images. Blue asterisk indicate structures where significant differences (p < 0.05) are observed between at least one of the MS phenotypes. Brain indicates total brain volume. PPMS, RRMS, and SPMS stands for primary progressive, relapsing remitting, and secondary progressive MS. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We also segmented raw and filtered images with 4D Freesurfer and 3D Freesurfer, respectively. Significant atrophy was observed in the combination of 4D filter and 3D Freesurfer segmentations in left and right thalamus (p < 0.05) in RRMS compared to both SPMS and PPMS, corroborating with the fact that the thalamus has been shown to have significant neurodegeneration in RRMS (Bergsland et al., 2012). Additionally, left ventricles in RRMS (p = 0.03) have higher atrophy rates than PPMS (cf. Pagani et al., 2005), while left hippocampus has higher atrophy rates (p = 0.02) in PPMS than RRMS. Total brain volumes from 3D Freesurfer does not show any significant difference between any groups (p > 0.10). The only significant difference on 4D Freesurfer segmentations of raw images were observed on right thalamus, which was more atrophied in PPMS compared to both RRMS and SPMS.

Discussion and conclusions

Our method employs an auto-regressive (AR(1)) filter on the intensities of images collected longitudinally. Higher order autoregressive models can certainly be included in the framework, where intensity at the tth time-point can depend on all the previous time-points, or in a simpler case, (t - 1)th through (t - N)th time-points. Eq. (6) can easily be modified to account for these higher order models. However, we use a simple model to demonstrate that the 4D filtering can be used as a pre-processing step to any 3D segmentation algorithm, while being similar or superior to alternative approaches. We have shown that our 4D filter can detect tissue atrophy in dementia or MCI, as well as different types of MS.

The autoregressive model promotes the temporal intensity trend of a patch to be smooth. Other approaches have imposed a temporal smoothness penalty on segmentations via 4D Markov random field priors. In contrast, we enforce smoothness on the intensity values, so that the algorithm is not tied to any particular segmentation method. This can be a valuable advantage in many circumstances where the properties of a particular segmentation algorithm are desired. For example, when combined with the TOADS algorithm as shown in Section 3.3, our filtering can lead to stable topology-preserved longitudinal segmentations. However, like other 4D segmentations, a limitation of our method is that if a new scan is observed, the filtered images at all time-points must be re-computed.

The temporal auto-regressive filter in Eq. (4) has two components, one to enforce smooth temporal intensities (Eq. (3)), and a weight to penalize smoothness if the variation is greater than a noise threshold (w( in Eq. (5)). Other smoothing filters, e.g., a moving average filter, can be used in place of the AR model, but a key component of our work is the weight w( which ensures that the patches are smoothed only when the variation in the temporal dimension is small enough. The filtering depends on a suitable choice of the threshold f. In theory, f can vary spatially and also depend on the expected rate of change in the underlying anatomy. In future work, we will explore the effect of a temporally and spatially varying f.

In addition to the temporal stability of the brain segmentation, bias in the longitudinal analysis is a common problem, that can arise from the asymmetric interpolations when the baseline (first time point) is used as a target for 4D nonlinear registrations (Yushkevich et al., 2010) in deformation based morphometry. To account for the bias in registrations, two (or more) time-points of a subject can be transformed to a “halfway” space to remove the directional bias, as done in SIENA (Smith et al., 2001). 4D Freesurfer creates an unbiased “median” template of all the time-points and registers all the time-points to the template. In this work, we simply used the MNI atlas as the target for rigid registrations of all the time-points before applying the 4D filtering (as described in Section 2.3), since the focus of this paper is on the robustness of the segmentations. The rigid registration to the MNI atlas before filtering can be replaced by a rigid registration to an unbiased median template, generated by ANTS (Avants et al., 2008) or 4D Freesurfer.

4D FS has no explicit noise model; it only transforms all the time-points into a common unbiased space. Therefore, we would expect a very similar atrophy detection performance of 4D filter followed by 4D FS as the combination of 4D filter and 3D FS, since the noise has already been reduced by the filter. For the test–retest study (TR-10 and CoRR-29), we expect that combination of 4D FS and 4D filter to increase the robustness, i.e. decrease the coefficients of variations, similar to the comparison between 3D vs 4D FS.

In summary, we have described a temporal filtering algorithm to obtain stable longitudinal segmentations that can be used as a pre-processing step to any segmentation method. We used both intensity based (FIRST, TOADS, Freesurfer), as well as registration based (MALP-EM) segmentation algorithms to show that our method is significantly more stable than other approaches while remaining sensitive to actual longitudinal changes.