A unified framework for focal intensity change detection and deformable image registration. Application to the monitoring of multiple sclerosis lesions in longitudinal 3D brain MRI

1. Introduction

Multiple sclerosis (MS) is an auto-immune neurodegenerative disease characterized by the inflammation of the myelin coating that surrounds the nerves. As a consequence, the transmission of nervous impulses is impaired, causing motor, cognitive and sensorial disabilities. The evolution of MS is characterized by the apparition of focal lesions in the brain and in the spinal cord, and by a progressive atrophy of brain tissues. Both phenomena can be monitored thanks to Magnetic Resonance Imaging (MRI) (Kaunzner and Gauthier, 2017). In the clinical routine, the evolution of the lesion load and of the brain atrophy is generally assessed qualitatively. However, the precise quantification of lesion changes over time may be of great interest to finely characterize the course of the pathology and to evaluate at the early stage the effect of a therapeutic strategy (McNamara et al., 2017). Since the manual delineation of lesion changes in MRI is a tedious and time consuming task, which is prone to both intra- and inter-observer variability, there is a great need for efficient and reliable automated tools (Altay et al., 2013).

Most change detection methods dedicated to lesion monitoring rely on a sequential scheme that first consists in removing all changes that are not of interest in order to detect in a second step only the evolution of lesions (Radke et al., 2005). To correct for global intensity changes induced by the difference of MRI acquisition setups, algorithms for bias field inhomogeneity correction (Song et al., 2017) and histogram-based intensity normalization procedure (Shinohara et al., 2014) are generally considered. Then, geometrical discrepancies due to variation in patient positioning, acquisition-related geometrical distortion and brain atrophy are corrected thanks to registration algorithms involving either rigid, affine or deformable transformations. Finally, the remaining changes corresponding to the evolution of lesions are detected. This final step generally consists in thresholding an intensity-based (Bosc et al., 2003) or deformation-based (Rey et al., 2002) feature map. The threshold can be chosen according to some statistical modeling in order to control the expected number of false positive detections (Rousseau et al., 2007).

The main flaw of the sequential procedure is that it implicitly assumes that each correction step can be performed while not being influenced by the changes that remain to be corrected in the next steps. As a consequence, the sequence order of the correction procedures should be carefully chosen. Moreover, for each correction step, a trade-off should be found between its performance (i.e., the ability of the method to accurately and specifically correct a given kind of change) and its robustness (i.e., the ability not to be biased by another kind of remaining changes). This observation advocates for a unified formulation of the change detection problem allowing to estimate all the different kinds of changes jointly.

In this paper, we address more specifically the interplay between deformable image registration and focal intensity change detection. When deformable registration is performed in the presence of appearing lesions, the estimated transformation tends to make these new lesions disappear in order to minimize the dissimilarity between the two images. This is the reason why the most common practice is to consider only rigid or affine registration in order not to alter lesion shape. However, such linear transforms can only compensate for difference in patient positioning but are not able to capture the complex deformations induced by brain atrophy, which typically occurs in MS. These remaining deformations may yield to spurious detections in atrophied areas, especially in the cortex and around the ventricles.

We propose to account for the intertwining of deformable registration and focal intensity change detection by estimating them jointly. To this end, we show that these two separate tasks can be formulated as a single optimization problem involving a unique energy that models their coupling. Basically, areas corresponding to detected changes are ignored in the registration similarity criterion, which prevents the lesion elimination effect described above. Solving this issue allows us to use of deformable registration, which in turn prevents from detecting spurious changes in atrophied areas. We propose an efficient alternating optimization scheme to solve this unified optimization problem. We focus on demonstrating the benefits of this joint formulation in the particular case of a standard intensity-based data similarity criterion. Nevertheless, the proposed approach is versatile and could easily be extended to more elaborated modeling choices. Experimental analysis is performed on the BrainWeb synthetic dataset and on two annotated real datasets. We first demonstrate in each case the benefits of considering a deformable registration as compared to an affine registration only in order to reduce the number of false detections. Then, we highlight the benefit of the proposed joint formulation as compared to the standard sequential scheme in terms of change detection accuracy. A preliminary version of this work has been published as a conference paper in Dufresne et al. (2020). In this paper, we provide a more extensive experimental analysis, which helps to better characterize the behavior of the proposed method and better understand why it outperforms the sequential approach.

The paper is organized as follows. In Section 2, the sequential approach, which will be considered as the reference baseline method of this work, is described and its limitations are discussed. In Section 3, we introduce the proposed joint model and the optmization strategies that has been set up estimate both change detection map and deformation field. In Section 4, we give implementation details. Finally, we present and discuss experimental results in Section 5.2.

2. The sequential approach

The conventional sequential approach consists of three main steps. First, the images are corrected for global intensity variations, then they are spatially registered, and finally the focal intensity changes due to the evolution of lesions are detected. In this section, we give a brief overview of the common practices in the registration and change detection steps. Our goal is not to cover a comprehensive scope of the field but to formulate the general principles underlying existing methods. In the remainder of this article, we will denote I₁, I₂: Ω → ℝ the baseline and follow-up MRI images, where Ω⊂ℝ³ is the image domain.

2.1. Registration

The registration problem can be formulated as:

where w: Ω → ℝ³ represents the transformation, ρ(·) is a data similarity term, and Ψ(·) is a regularizer weighted by a scalar λ₁ > 0. An overview of deformable registration methods in medical imaging can be found in Sotiras et al. (2013).

Several transformation models can be considered relying either on a parametric representation (e.g., rigid, affine, polynomial, or Bspline-based) or on a non-parametric deformable mapping (i.e., a displacement vector is estimated for each voxel).

The role of the data term is to penalize dissimilarity between I₁ and I₂ warped with the estimated transformation. For monomodal registration, it is common to use intensity-based measures such as the sum of squared intensity differences or the cross correlation. In the multimodal case, mutual information is one of the most widely used similarity metric (Kaunzner and Gauthier, 2017).

In the context of deformable registration, considering the data term only can lead to an ill-posed problem. To overcome this issue, the data term has to be balanced with an additional regularization term Ψ(w) that enforces some constraints on the deformation field. For instance, penalizing the ℓ₂ or ℓ₁ norm of the gradient of w helps to promote smooth solutions.

2.2. Change detection

The change detection step generally consists in thresholding a map of feature differences between registered baseline and follow-up images (see the survey Lladó et al., 2012). These maps can be calculated directly from either intensity-based (Sweeney et al., 2013; Ganiler et al., 2014; Cabezas et al., 2016), or deformation-based (Rey et al., 2002; Cabezas et al., 2016; Salem et al., 2018) features, and sometimes integrate other kind of information (Elliott et al., 2013; Sweeney et al., 2013).

In the perspective of integrating both registration and change detection in a single joint optimization problem, we advocate that they should both rely on the same data similarity. Consequently, we model the binary change map c : Ω → {0, 1} defined at each voxel x as follow:

${\lambda }_2\in {ℝ}^{+*}$ being the detection threshold. The thresholding scheme (Equation 2) can be reformulated as the following optimization problem:

Since simple thresholding can yield noisy results, most MS lesion change detection methods also integrate a denoising step in post-processing to obtain the final change map. The denoising can be realized jointly with the change detection by integrating a regularization term Φ(·) in Equation (3):

Where λ₃ weights the regularization term.

2.3. Limitation of the sequential approach

The main limitation of the sequential approach is illustrated in Figure 1 involving the baseline (Figure 1) and follow-up (Figure 1) MRI acquisitions of a patient suffering from MS. One can observe in the follow-up scan the apparition of a new lesion and a slight enlargement of the ventricle reflecting the brain atrophy process. In the case of affine registration (Figures 1), we can see on the subtraction image that the lesion is well detected, but that spurious detection occur around the ventricles and in the cortical regions due to brain tissue atrophy. Using a deformable registration (Figures 1) helps to remove these spurious detection by compensating the ventricles enlargement and the cortical atrophy. However, it also tends to make the new lesion disappear (Figure 1), thus altering the shape of the corresponding detection (Figure 1). The goal of the proposed joint approach (Figures 1) is to perform an accurate atrophy correction while preserving the shape of appearing and evolving lesions, even when the lesion-to-tissue contrast is quite low.

3. Joint approach

3.1. General formulation

To overcome the limitations of the sequential approach, we advocate a joint modeling of registration and change detection. The two steps are fundamentally intertwined, since registration aims at finding correspondences between images, while change detection determines regions that does not admit correspondences. Therefore, both tasks should be defined with the same objective function to work in synergy. We formulate the following joint minimization problem that achieves this goal by unifying the principles described previously:

With this model, the data term is cancelled in change regions (where c(x) = 1), so that the estimation of the transformation is only driven by the regularization term, thus producing smoothed deformation field in these areas. By this way, it prevents the lesion disappearing effect observed in Figure 1.

3.2. Modeling choices

The formulation (Equation 5) is versatile and could be instantiated with a variety of data and regularization terms. In this paper, the goal is to demonstrate the superiority of the joint formulation over the sequential approach under standard modeling choices, which are detailed in the sequel.

First, we assume that, thanks to the intensity normalization step done as a preprocessing, intensities of both images are comparable, thus allowing us to consider a data term based on intensity difference. Consequently, we consider the following standard similarity measure:

Where σ is a normalization constant defined by the median absolute deviation of the intensity differences between I₁ and I₂. This data term has been used for motion estimation (Bruhn et al., 2005) and is representative of intensity-based features commonly used in change detection methods (Sweeney et al., 2013; Ganiler et al., 2014).

Secondly, we assume that the deformations induced by brain tissue atrophy are complex but still locally smooth. This is why we consider a non-parametric representation of the transformation field w while introducing a first order Tikhonov regularization term:

Where ∇· is the gradient operator.

Finally, we assume that the detected changes should be spatially coherent. Consequently, the change map c is regularized with a standard binary Potts model:

Where δ is the Kronecker function equal to 1 if its argument is true and ${N}(\mathrm{\text{x}})$ is the 6-neighborhood of x.

3.3. Optimization

To solve the optimization problem (Equation 5), we rely on an alternating minimization strategy: at each iteration, we successively minimize with respect to (w.r.t.) each variable, while keeping the other fixed. We detail in this section the optimization strategies dedicated to each of the two subproblems.

Notice that, since the problem is nonconvex, convergence toward the global minimum cannot unfortunately be guaranteed. However, the results presented in Section 5.2 on several datasets suggest that the optimization process converges in practice toward a satisfying solution.

Minimization w.r.t. w To make the problem tractable, we consider the linearized version of the data term (Equation 6) obtained by replacing I₂(y − w(x)) with its Taylor expansion around y:

Where I_t(y) = I₂(y)−I₁(y) is the temporal derivative. Since the Taylor expansion is valid only for small deformations, we embed the estimation in a coarse-to-fine scheme, which is a common practice in registration and motion estimation (Hill et al., 2001).

After Taylor development of I₂(y − w(x)), the new data term ρ_l in Equation (9) is the composition of a quadratic function with a linear function of w, which yields a convex term. Since the regularization term (Equation 7) is also convex, the whole optimization problem is convex. By substituting ρ by ρ_l in Equation (5), the optimisation problem can be addressed with a variety of efficient optimization methods.

We chose to consider the alternated direction method of multipliers (ADMM) framework (Boyd et al., 2011). To this end, we introduce a splitting variable z that decouples the two terms of Equation (5) that depend on w, and we formulate the problem in the constrained form:

The ADMM algorithm is based on the minimization of the augmented Lagrangian associated with Equation (10) w.r.t. w and z, and a gradient ascent on the dual variable (Boyd et al., 2011). It leads to the following iterative updates of w and z (see Fortun et al., 2018 for a similar derivation with different data and regularization terms):

Where ${\mathrm{\text{prox}}}_f(\mathrm{\text{x}})=\underset{\mathrm{\text{y}}}{\mathrm{\text{argmin}}}\frac{1}{2}||\mathrm{\text{x}}-\mathrm{\text{y}}{||}_2^2+f(\mathrm{\text{y}})$ denotes the proximity operator of a function f. The subproblem (Equation 11) is voxel-wise and quadratic, and it admits a simple closed form solution. The subproblem (Equation 12) is equivalent to a denoising operation with the regularizer Ψ(·), and it also has a closed form linear solution that can be computed efficiently in the Fourier domain. μ is the parameter associated with the quadratic penalty in the Augmented Lagrangian associated with Equation (10). The ADMM algorithm is derived from this Augmented Lagrangian and the update (Equations 11, 12) are its minimization w.r.t. w and z. Intuitively, μ controls how fast the constraint w = z is imposed through the optimization process. Thus, even if it is not strictly speaking a step size, it has a similar impact on the convergence speed.

Notice that the ADMM framework is flexible enough to cope with different data and regularization terms with low computational cost. The requirement is to be able to design a splitting of the cost function such that the proximity operators of each iteration have computationally efficient solutions. Examples of admissible models comprise data terms based on the ℓ₁ penalty function or cross-correlation (Vogel et al., 2013), and regularizations by total variation or Nuclear norm of the Jacobian (Bostan et al., 2014).

Minimization w.r.t. c When w is fixed, the estimation of c amounts to a binary segmentation problem with Potts regularization:

We solve it with a graph-cut method (Boykov et al., 2001), which is able to find an exact solution with very low computational cost.

4. Implementation details

4.1. Pre-processing

Before applying the change detection framework, the input images require to be pre-processed as follows. First, images are corrected for bias field inhomogeneity using the N4 algorithm (Tustison et al., 2010). Then, a global scaling of the intensity is performed in order to enforce the median value of the intensities inside the brainmask to be equal to 100. Images are then resampled to 1 mm isotropic resolution. These two steps aim at harmonizing all input images in terms of intensity range and spatial resolution. The follow-up scan is then rigidly registered on the baseline image using ANTs library (Avants et al., 2011)¹ with default parameters and mutual information metric. Differential bias field inhomogeneity is corrected thanks to the method described in Lewis and Fox (2004), while considering a 21 × 21 × 21 median filter size.

4.2. Post-processing

The detection maps are post-processed by discarding the connected components smaller than 3mm³ and the detections outside brain parenchyma (i.e., the union of gray and white matter). Notice that the binary change detection map can also be computed to reflect either positive or negative intensity changes only.

Brain parenchyma masks are computed from T1-weighted images using the FAST method provided in the FSL library (Zhang et al., 2001) or alternatively from FLAIR images using SAMSEG² brain parcellation tool (Cerri et al., 2021), in the case where no T1-weighted image is available.

4.3. Hyperparameter setting

Three hyperparameters have to be set in the proposed joint formulation (Equation 5): λ₁, controlling the spatial regularization of the deformation field, λ₂ acting as a threshold for the map of intensity differences, and λ₃, controlling the spatial regularization of the change map. Here, we suggest strategies to find out relevant parameters settings.

The value of λ₁ should be chosen to optimally estimate longitudinal brain atrophy, since it is the main source of brain deformation in MS. Thus, we consider a subset of 21 images from the dataset OASIS-3 (LaMontagne et al., 2019) that contains longitudinal Alzheimer and normal aging MRI data that exhibit various pattern of longitudinal brain atrophy. We determine the optimal value of λ₁ by selecting the one that leads to the best registration performance on this subset of OASIS-3. To this end, we derive a registration quality metric from the provided segmentation maps of brain structures obtained with Freesurfer. Concretely, the Dice score is computed for each structure between the segmentation maps of the baseline image and of the registered follow-up image. The global registration quality metric is then computed as the sum over all the regions of the median Dice score observed for each region. This procedure leads us to find λ₁ = 70 as an optimal value.

The values of λ₂ and λ₃ have to be set to find the best compromise regarding: (i) The expected intensity difference, (ii) the noise level that corrupts the images, and (iii) the spatial extent of the changes. Here, we suggest an approach to find out optimized setting for each of the two considered databases (see Sections 5.1.2, 3.2). In practice, λ₂ and λ₃ have been fixed to maximize the overall performance of the affine sequential approach (see Section 5.1.5) in terms of local Dice Similarity Coefficient (local DSC, see Section 5.1.4) for each dataset. Considering the local DSC ensures to focus on the ability of the detection scheme (Equation 14) to recover the detected changes while not being influenced by false positive detections that can occur in other parts of the brain. Considering the affine transformation model ensures that the registration step does not to alter the geometry of evolving regions. This procedure leads us to find λ₂ = 16 and λ₃ = 5 as optimal setting for LesjakDB dataset and λ₂ = 25 and λ₃ = 3 for MSSEG-2 dataset.

4.4. Convergence and stopping criteria

The iterations of the alternated minimization of Equation (5) and of the ADMM algorithm (Equations 11–13) are stopped when a stopping criterion is verified or when a maximum number of iterations is reached. The stopping criterion is a threshold on the norm of the relative changes between two consecutive iterations, and is set to 10⁻³ for the alternated minimization and 2.10⁻³ for ADMM. The maximum number of iterations is set to 5 for the alternated minimization and 300 for ADMM.

5. Experimental evaluation

5.1. Evaluation framework

In this section, we report results obtained on one synthetic dataset and two publicly available real patients datasets. The synthetic dataset offer the advantage to have an unambiguously defined ground truth change detection map, while controlling the amount of noise, bias field inhomogeneity and brain atrophy that corrupt the images. The real datasets are used to evaluate the proposed approach in conditions that are closer to the clinical routine, while considering different acquisition conditions and various pathological evolution. The first real patients dataset, denoted in the sequel as LesjakDB (Lesjak et al., 2016), is dedicated to assess the ability of methods to detect every kinds of MS lesion evolutions (shrinkage, growth, new and disappearing), whereas the second dataset, denoted in the sequel as MSSEG-2 (COM, 2021), only focus on the ability to detect new appearing lesions.

5.1.1. Synthetic dataset

We evaluate the proposed method on T2-weighted synthetic volumes generated with the Brainweb simulator (Cocosco et al., 1997), while considering the normal anatomical model (i.e., without lesion) and two multiple sclerosis anatomical models with moderate and severe lesion load. The images are simulated at a 1 mm³ isotropic resolution (image size: 181 x 217 x 181) with bias field inhomogeneity (20%). To simulate realistic brain atrophy for the follow-up image, we applied a deformation field that has been estimated using a deformable registration (Avants et al., 2008) from two T1-weighted MRI scans acquired 4 years apart of a patient suffering from MS that exhibits a significant brain atrophy evolution (in-house dataset). It should be noted that these real data images were first affinely registered onto the brainweb image to ensure the estimated deformable registration to be consistent with the underlying anatomy. The visual inspection confirms that the simulated image exhibits a realistic atrophy pattern. Gaussian additive noise was added with a standard deviation fixed at 5% of the mean intensity in the brightest tissue (cerebrospinal fluid in the T2-weighted simulation). We consider several scenarios of simulated longitudinal acquisition that are summarized in Table 1.

5.1.2. Real dataset LesjakDB: All kinds of lesion evolution

LesjakDB dataset (Lesjak et al., 2016) is composed of 20 longitudinal MRI acquisitions of MS patients with two timepoints. The median time between the baseline and follow-up studies was 311 days, ranging from 81 to 723 days. Each MRI acquisition consists in a 2D T1-weighted, a 2D T2-weighted and 2D-FLAIR sequences. Change detection was conducted on the FLAIR images only. The FLAIR image size is 256 × 256 × 49 with an anisotropic spatial resolution of 0.9 × 0.9 × 3 mm. Ground truth change detection maps are also provided, which were obtained from manual annotations done by two expert raters. We adjusted some of the ground truth annotations that did not match the real lesion changes. The annotated changes include appearing, growing, shrinking and disappearing lesions. Ground-truth detection maps are compared to binary detection maps that include both positive and negative intensity changes.

5.1.3. Real dataset MSSEG-2: Only appearing lesions

MSSEG-2 dataset (COM, 2021) is composed of 100 pairs of FLAIR MRI scans from MS patients acquired on various MR scanners. The provided ground-truth is limited to new appearing lesions, and was build from the consensus of manual annotations delineated by four experts. The dataset is separated into training (40 patients) and testing (60 patients) sets. Since the proposed approach does not require any training step, we consider the whole dataset for testing. However, we distinguish two subgroups of data, namely MSSEG-2-Change corresponding the 61 subjects that exhibit at least one new appearing lesion and MSSEG-2-NoChange corresponding the 39 subjects that do not exhibit any new appearing lesion. Since the provided ground-truth is limited to new appearing lesions, they are compared only to the positive binary change detection maps obtained with the different methods. Notice, that the proposed method framework does not discriminate appearing from evolving lesion. Consequently, lesion evolutions, which are not labeled in the ground truth detection maps, are erroneously considered as false positive detection, thus introducing a bias in some of the evaluation metrics.

5.1.4. Metrics

We report four metrics to evaluate the performance of the methods to detect changes, namely the Dice Similarity Coefficient (DSC), the Positive Predictive Value (PPV), the True Positive Ratio (TPR) and the local DSC. Let TP, TN, FP, and FN be the number of voxels from estimated change detection map that correspond to True Positive, True Negative, False Positive and False Negative, respectively.

The DSC is defined as:

and reflects the overall good overlap between the detection map and the ground truth.

The PPV is defined as:

and reflects the proportion of relevant detections among all the detected changes.

The TPR is defined as:

and reflects the proportion of the ground-truth changes that have been detected.

The local DSC correspond the DSC computed on a restricted area defined as the dilation with a 4-voxel radius spherical structuring element to the ground truth. This metric enables us to focus the evaluation on the local spatial accuracy of the detection method

In addition to the voxel-wise metrics, we also report lesion-wise metrics, namely the Lesion True Positive Ratio (L-TPR) and the Lesion Positive Predictive Value (L-PPV). These metrics have been evaluated thanks to the animaSegPerfAnalyzer validation tool while considering the same hyperparameters as in Commowick et al. (2018).

Since all these metrics are not relevant for data that do not exhibit any changes, we consider in that specific case the number of detected connected components as well as the volume of detected changes to characterize the false positive detections.

5.1.5. Variants used for comparison

To demonstrate the benefits of the proposed joint modeling, we consider three variants of the change detection framework:

• joint: the proposed joint change detection and registration method described in Section 3.

• sequential: The sequential counterpart of the proposed method, which successively performs deformable registration and change detection. For the two steps, we use the same model and optimization algorithms as in substeps of the joint approach described in Section 3.3.

• affine: The sequential approach where the deformable registration has been replaced by affine registration, which corresponds to the most common case. The affine registration was estimated using ANTs library (Avants et al., 2011) ³ with default parameters and mutual information metric. Then, the thresholding and smoothing of the change map routine follows model (Equation 14).

5.2. Results

5.2.1. Synthetic dataset

First, a qualitative visual comparison of the three methods is provided in Figure 2 for the lesion appearance with atrophy scenario. The affine method succeeds to detect almost all the lesion areas, but it suffers from false positive detection around the ventricles due to brain atrophy. Both the sequential and the joint methods compensated for brain atrophy deformation since none of them exhibit false detections around the ventricles. However, the sequential method failed to detect the whole lesion areas due to the over-compensation of lesion changes. This limitation is overcome by the joint approach that succeeds to detect the entire lesion areas.

A quantitative comparison of the three methods under four scenarios is provided in Table 2. First, we consider the appearance of lesion without atrophy. Unsurprisingly, this scenario is the most favorable for the affine method since there is no geometric difference to compensate. The sequential approach yields to significantly lower values of DSC and local DSC. This is due to the lesion over-compensation effect, as supported by the observed low TPR value (i.e., lack of sensitivity) and high PPV value (i.e., high specificity). Finally, the joint approach overcomes the shortcoming of the sequential approach and have performances similar to the affine method, with a slight tendency to underestimate the detected area. Similar observations can be made about the second scenario involving lesion growth without atrophy.

The conclusions are drastically different for the two scenarios involving simulated atrophy. The performance of the affine method is significantly hampered by the numerous false detections due to the atrophy. This is illustrated by the significant decrease of the DSC and PPV values compared to the cases without atrophy, while the TPR and local DSC values are less modified. The sequential approach succeeds to compensate for the simulated brain atrophy, as highlighted by the high PPV value, but still underestimates the changes to detect, as indicated by the low TPR value. The joint approach clearly outperforms the two previous approaches in terms of detection accuracy, as objectified by the significantly higher DSC value.

The behavior of the joint approach can be monitored through the iterations of the alternating optimization scheme (see Figure 3). We can see that the DSC increases across the iterations, and the convergence is reached in a few iterations. Concerning the computational burden of the joint approach, it is about 24min on one single core (Intel(R) Xeon(R) Gold 6130 CPU @ 2.10GHz) for an experiment on the synthetic dataset (image size: 181 x 217 x 181).

5.2.2. LesjakDB

First, a qualitative visual comparison of the three methods is provided in Figure 4. We can draw similar conclusions as for the synthetic dataset (see Figure 2). The affine method demonstrates a high sensitivity (i.e., the lesion evolution is well detected) but a lack of specificity (i.e., numerous false positive detections are detected around the ventricles and in the posterior part of the cortex). Conversely, the sequential method has high specificity but lacks sensibility. The joint approach provides the best visual results, thus illustrating its ability to achieve both high sensitivity and high specificity. Figures 4 show the Jacobian of the deformation fields obtained by the sequential and joint methods, respectively. The specific pattern characterized by the alternance of both high and low values of the jacobian (see areas highlighted by the red squares in Figure 4) reflects the high local contraction and dilation induced by the deformation field to make the lesion disappear, thus explaining the lack of sensitivity of the detection results.

The quantitative evaluation shown in Table 3 (first row) and in Figure 5 confirms the conclusions of the visual analysis. The high sensitivity of the affine method is objectified at the voxel level by a statistically significantly higher TPR than the two other methods. At the lesion level, all the three methods exhibit similar L-TPR values, thus emphasizing their ability to detect the same amount of changing areas. Both the sequential and joint methods yield significantly higher PPV and L-PPV compared to the affine method, which illustrates their ability to reduce the number of false detections induced by brain atrophy at both voxel and lesion levels. This result highlights the benefit of using deformable registration in the context of MS lesion monitoring. The significantly lower TPR achieved by the sequential method compared to the joint method is the consequence of the lesion over-compensation effect. Finally, the joint approach significantly outperforms the two other approaches in term of voxel-wise global accuracy (see DSC).

Figure 6 highlights the variability of the performance of the methods across the subjects. It is interesting to notice that, although the performance of the detection methods greatly varies from one subject to the other, the ranking among the three methods appears to be highly consistent across the subjects. When investigating for the factors that may explain the observed variability, it appears that the volume of the ground-truth seems to play a prominent role: the larger is the volume to detect, the better is the performance of the change detection algorithm (see Figure 7).

5.2.3. MSSEG-2

Similarly as for the synthetic and LesjakDB datasets, the qualitative visual comparison of the two approaches based on deformable registration highlights the lack of sensitivity of the sequential method due the lesion over-compensation effect (see Figure 8). The resulting change detection map (purple) is too small compared to the ground-truth (underlying transparent red) due to the deformable registration that significantly shrinks the lesion. With the joint approach, the shape of the lesion is almost preserved in the warped follow-up image and the change detection map (green) matches almost perfectly the ground-truth.

The quantitative evaluation on the subset MSSEG-2-Change is reported in the second row of Table 3 and in the upper part of Figure 9.

The fact that both the sequential and joint approaches lead to significantly higher PPV values as compared to the affine approach advocates the use of deformable registration to reduce the number of false detections. The benefit of considering the joint over the sequential approach to overcome the lesion overcompensation effect is clearly demonstrated by the significantly higher TPR and local DSC values obtained with joint method.

It is also interesting to notice that the lesion over-compensation effect does not affect the special case of disappearing lesion. Indeed, when registering an image without lesion on a image with a lesion, the dissimilarity in the area of the disappearing lesion cannot be corrected by the transport of intensity of the registration (this is in fact only the case for non symmetric image registration method, see Noblet et al. (2004) for further explanations). To illustrate this phenomenon, we consider the MSSEG-2-Change inverse experiment (see the third row of both Table 3 and the bottom part of Figure 9) that consist in swapping the baseline and the follow-up image, so that the ground-truth now correspond to disappearing lesions. The same conclusion can be drawn from the DSC, PPV, and TPR as compared to the MSSEG-2-Change experiment. The most interesting point concern the local DSC that focuses the evaluation on the disappearing lesion. In that case, there is no significant difference any more between sequential and joint approaches contrary to the MSSEG-2-Change experiment, showing the absence of lesion over-compensation effect in the specific scenario of detecting disappearing lesions.

Note that all the results presented above in this section are evaluated on the 61 subjects of MSSEG-2-Change (i.e., subject presenting at least one new appearing lesion). Indeed, the presented metrics cannot be computed anymore for the 39 subjects of MSSEG-2-NoChange since the ground-truth change detection map is empty. This is why we only report the volume of detected changes for this subset of MSSEG-2 (see Figure 10). We can notice that both sequential and joint approaches lead to significantly lower volume of detected changes as compared to the affine, which appears in line with previous findings that support the use of deformable registration to reduce the number of false detections. Also note that the joint method yields consistently to slightly higher volume of detected changes as compared to the sequential method. This is also the consequence of the lesion over-compensation effect that affects the sequential appproach.

6. Conclusion and perspectives

We have presented a method that unifies registration and change detection for the analysis of longitudinal brain MRI. It is based on the joint modeling of these two tasks as the minimization of a single objective function, for which we have developed an efficient alternating optimization method. The proposed approach has been evaluated in the context of the follow-up of multiple sclerosis lesion, which requires deformable registration to capture characteristic brain atrophy, but also with the potential caveat to shrink appearing lesions. In this context, the conventional sequential detection pipeline leads to large detection inaccuracies around new appearing lesions. We have demonstrated on simulated and real data that the proposed joint approach is able to combine the ability of deformable registration to correct brain atrophy, and a good preservation of the lesions shape to ensure accurate change detection.

The implementation presented in this paper of the proposed joint model relies in fact on quite simple modeling assumptions. The versatility of the optimization approach opens the way for more sophisticated models that could be handled in the same framework. Alternative data fidelity terms such as cross correlation or mutual information, and regularizers such as total variation could be considered to potentially improve the performance of the method. Another perspective is to improve the convergence of the alternating optimization strategy to ensure a better robustness to local minima. To this end we could consider a fuzzy change detection map to turn its estimation into a continuous optimization problem. This would allow us to use more robust optimization approaches such as Proximal Alternating Linearized Minimization (Bolte et al., 2014). Finally, while we have addressed in this paper the registration problem, the unification principle could be extended to other steps of the change detection pipeline. In particular, the intensity normalization and the bias field inhomogeneity correction of the MRI acquisitions are crucial pre-processing tasks that are impacted by the presence of evolving lesions. Therefore, the integration of these two tasks in a single unified model would be a natural extension of the proposed framework.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found at: https://brainweb.bic.mni.mcgill.ca/selection_ms.html; http://lit.fe.uni-lj.si/tools.php?lang=eng; https://portal.fli-iam.irisa.fr/msseg-2/.

Author contributions

The proposed method was coded by ED. ED, DF, and VN equally contributed to the writing of the article and the interpretation of the results. SK provided a clinical expertise for the analysis of the results and the real annotated databases. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by the Région Grand Est and Philips Healthcare.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

References