Fast and Robust Unsupervised Identification of MS Lesion Change Using the Statistical Detection of Changes Algorithm

SDC Lesion-Detection Algorithm

Given 2 MR images I₁ and I₂ of the same brain acquired at 2 time points, the voxel-subtraction signal d = I₂ − I₁ (Fig 1) is assumed to follow a Gaussian distribution N(μ, σ²) with mean μ and SD σ, in which σ can be estimated from the set of nonlesion WM voxels on the subtraction image. Most of these voxels belong to the intersection of the 2 WM masks obtained by brain segmentation tools (such as FSL; http://www.fmrib.ox.ac.uk/fsl⁵) from T1-weighted structural images. These masks typically exclude large lesions that are hypointense on T1WI and therefore consist mainly of nonlesion voxels (On-line Fig 1).

• Download figure
• Open in new tab
• Download powerpoint

Fig 1.

Schematic of the proposed SDC lesion-change detection algorithm on the T2-weighted FLAIR subtraction image. One unchanged lesion (yellow arrow) and 1 new lesion (red arrow) are correctly identified by the SDC. The algorithm automatically generates red ROIs, which encompass the detected areas of change on the second FLAIR image to help the human reader quickly identify lesion changes.

The MR imaging signal-to-noise property can be used to formulate an optimal SDC of lesions as a composite statistical test between 2 hypotheses of the following likelihood functions: In this work, the SDC test statistic was computed over a 3-voxel connected neighborhood based on the currently accepted minimum MS lesion size requirement of 3 mm (3 voxels in 1-mm³ isotropic images)⁶ and on the assumption that the subtraction signals within this small neighborhood are similar. Denoting the subtraction signals at the i-th voxel and its neighbor voxels as d_i1, …, d_i3 and assuming μ > 0 (positive change), the test statistic t_i can be computed from the log-likelihood ratio test⁷ and compared with a threshold γ to make a decision: Here γ was chosen to control the false-positive rate P_FP = P(t_i > γ | H₀). According to the Neyman-Pearson lemma, this test provides the best detection power for a given P_FP regardless of the unknown mean μ (uniformly most powerful detector).⁷

The test statistic is maximized over all possible neighborhoods surrounding the voxel, to increase the sensitivity of lesion detection: where V_i denotes a 3-voxel connected neighborhood system of the i-th voxel (On-line Fig 2). Intuitively, this test statistic encodes in probabilistic terms the expectation that a bright voxel on the subtraction image is more likely to be identified as “changed” if at least 2 of its neighboring voxels also have high signals.

MR Imaging Experiment

This was a retrospective study of 30 patients with MS with 2 consecutive brain MRIs (mean scan interval, 267 ± 104 days; range, 15–410 days) performed on 3T scanners (Magnetom Skyra, VE11A software; Siemens, Erlangen, Germany). The imaging protocol consisted of an MPRAGE T1WI sequence for brain structure (TR/TE/TI = 2300/2.3/900 ms, 1 mm³ isotropic) and a T2WI FLAIR sequence for lesion detection (TR/TE/TI = 7600/446/2450 ms, 1 mm³ isotropic). After skull removal and bias field correction, FLAIR images were coregistered into the half-way space using the FMRIB Linear Image Registration Tool algorithm (FLIRT; http://www.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT)⁵ to ensure that the degree of blurring introduced by coregistration was similar between images because this similarity improves subtraction. To account for changes in image contrast or dynamic range (eg, due to different receiver gain settings or slight changes in imaging parameters), we performed image-intensity normalization before subtraction. The robust intensity range (second and 98th percentiles, denoted as m and M, respectively) was computed for each image. The image intensity of the second image I₂ was then scaled linearly to match that of the first image I₁ as follows: I_2,scaled = αI₂ + β, where α = (M₁ − m₁) / (M₂ − m₂) and β = [(M₁ − αM₂) + (m₁ − αm₂)] / 2. In addition, brain GM, WM, and CSF masks were obtained from the T1WI using the FMRIB Automated Segmentation Tool (FAST; https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FAST) algorithm.⁵ The SDC test statistic (Equations 1 and 2) was then computed and thresholded to generate a change mask (Fig 1). The false-positive rate P_FP was set to 0.0001, which means that, on average, 50 of approximately 500,000 WM voxels may be incorrectly labeled as “changed.” To reduce the number of false-positives, we imposed additional constraints on lesion size (≥3 voxels), location (lesions located within 2 voxels of the CSF border had to be part of a larger lesion that extended outside this border), and intensity on the second FLAIR image (>2 SDs above the mean normal-appearing WM intensity, ie, WM voxels that do not appear bright on FLAIR were excluded).

For comparison, the lesion prediction algorithm (LPA), part of the Lesion Segmentation Tool toolbox (LST; http://www.applied-statistics.de/lst.html),⁸ was used to compute the lesion masks from FLAIR images. This algorithm consists of a binary classifier in the form of a logistic regression model trained on the data of 53 patients with MS.⁸ The lesion masks were then subtracted to obtain the change mask without any human revision. Like the SDC, lesion changes of <3 voxels were excluded.

Statistical Analysis

A neuroradiologist with 6 years of experience reviewed the 2 FLAIR and the subtraction images with the help of computer-generated color ROIs that encompassed the detected lesion changes (Fig 2). These were labeled as “true-positive” or “false-positive.” The reader also reviewed the images outside these ROIs to count the number of missed (false-negative) and unchanged (true-negative) lesions. Lesion changes detected by the SDC and LPA were presented in randomized order (both by subject and by detection algorithm) to the reader, who was blinded to the algorithm. The image review time was recorded for each subject and algorithm. A 2-tailed paired-sample t test was used to compare the mean review time per subject of SDC and LPA. The sensitivity and specificity of each method were calculated using the generalized estimating equation logistic regression, which accounts for the correlation among the measurements within the same subject.⁹

• Download figure
• Open in new tab
• Download powerpoint

Fig 2.

Comparison of lesion changes identified by the LPA and SDC (indicated by the red ROIs generated by the algorithms and superimposed on the T2 source FLAIR images). Each column shows images acquired at a different time points. The LPA identifies more false-positives (yellow arrows) yet misses a new small lesion (red arrow). The SDC correctly classifies positive lesion changes in concordance with the human expert.