Perivascular Spaces Segmentation in Brain MRI Using Optimal 3D Filtering

Perivascular Spaces (PVS) are a feature of Small Vessel Disease (SVD), and are an important part of the brain’s circulation and glymphatic drainage system. Quantitative analysis of PVS on Magnetic Resonance Images (MRI) is important for understanding their relationship with neurological diseases. In this work, we propose a segmentation technique based on the 3D Frangi filtering for extraction of PVS from MRI. We used ordered logit models and visual rating scales as alternative ground truth for Frangi filter parameter optimization and evaluation. We optimized and validated our proposed models on two independent cohorts, a dementia sample (N = 20) and patients who previously had mild to moderate stroke (N = 48). Results demonstrate the robustness and generalisability of our segmentation method. Segmentation-based PVS burden estimates correlated well with neuroradiological assessments (Spearman’s ρ = 0.74, p < 0.001), supporting the potential of our proposed method.

Abstract-Perivascular Spaces (PVS) are a recently recognised feature of Small Vessel Disease (SVD), also indicating neuroinflammation, and are an important part of the brain's circulation and glymphatic drainage system. Quantitative analysis of PVS on Magnetic Resonance Images (MRI) is important for understanding their relationship with neurological diseases. In this work, we propose a segmentation technique based on the 3D Frangi filtering for extraction of PVS from MRI. Based on prior knowledge from neuroradiological ratings of PVS, we used ordered logit models to optimise Frangi filter parameters in response to the variability in the scanner's parameters and study protocols. We optimized and validated our proposed models on two independent cohorts, a dementia sample (N=20) and patients who previously had mild to moderate stroke (N=48). Results demonstrate the robustness and generalisability of our segmentation method. Segmentation-based PVS burden estimates correlated with neuroradiological assessments (Spearman's ρ = 0.74, p < 0.001), suggesting the great potential of our proposed method.

I. INTRODUCTION
P ERIVASCULAR SPACES (PVS), also known as Virchow-Robin spaces, are fluid-filled spaces that follow the typical course of cerebral penetrating vessels. PVS have the same Magnetic Resonance Imaging (MRI) contrast characteristics of Cerebrospinal Fluid (CSF), that is they appear hypointense (dark) on T1-weighted (T1) and hyperintense (bright) on T2-weighted images (T2) [1], [2]. They appear as small 3D tubular structures that, depending on the viewing plane, are linear or round, with a diameter generally smaller than 3mm (see Fig. 1) [3].
Enlargement of perivascular spaces is associated with other morphological features of Small Vessel Disease (SVD) such as white matter hyperintensities and lacunes [4]; cognitive impairment [5] and inflammation [6]. Most studies use visual rating scales to assess PVS burden [7], [8], but are prone to observer variation, particularly in the Centrum Semiovale (CS), due to the coexistence of PVS with other neuroradiological features of SVD that confound their identification in this region.
Efforts have been made to computationally assess PVS [3], [9]. Recent semi-automatic methods are based on thresholding and require user intervention either for the choice of parameters or for manual editing of the resulting masks [10], [11]. One of the most promising approaches proposed for PVS automatic segmentation uses the Frangi filter [12] parameterised through a Random Forest scheme [13] that learns discriminative PVS characteristics from manually segmented ground truth on MR images acquired at 7T [14], [15]. However, MRI in clinical research and practice is mostly performed in scanners with field strengths at 1.5T or 3T, and the reference standards available are visual ratings performed by neuroradiologists, which makes the learning-based approach proposed by Park et al. [15] limited for practical use. Moreover, it is difficult to assess enlarged PVS burden at high field 7T MRI since normal PVS and deep medullary veins with similar intensities to PVS confound visualization. Our current goal is to present a segmentation approach for enlarged PVS that arXiv:1704.07699v1 [cs.CV] 25 Apr 2017 can be used widely in current clinical research studies, to further elucidate their pathological significance and assess their potential role in neurological disorders.
We propose a novel application of ordered logit models, usually used in statistics as a regression model for ordinal dependent variables, as this model provides a good estimate for capturing the sources of influence that explain the ordinal dependent variables (i.e. in this case the PVS visual rating scores) considering the uncertainty (i.e. subjectivity, interobserver variability) in the measurement of such data [16]. We use this model to estimate the parameters of the Frangi filter [12] to obtain the maximum likelihood of a vessel-like structure to be a PVS in the CS, by also estimating the count of PVS that most likely falls in the class corresponding to the category given by the neuroradiologist in this brain region.
We calibrated different ordered logit model, according to the rating scale available for every dataset. We optimized the parameters of the Frangi filter to deal with T1 and T2 modalities, and combined the resulting filtered images. Validation was carried out on different cohorts, using images acquired in different sites.
The remainder of the paper is organised as follows: Section II introduces the data used in this paper. In Sec. III we describe the basis of our optimization approach. Section IV discusses some implementation details and the model validation. Following are the validation results which are presented in Section V. Finally the discussion and conclusions are shown in Sec. VI.

II. MATERIALS
Two datasets were used for developing, testing and validating the method:  [17]. The characteristics of the sequences relevant for PVS assessment are summarized in Table I.  [12], largely used for enhancing blood vessels, for instance in retinal images [18]. Given the absence of an accurate computational "ground truth" (i.e. manual labels of each PVS by experts), we propose a modelling technique to use the available information (i.e. PVS burden assessed using visual rating scales) to optimize the filter parameters. For this scope, an ordered logit model [16] has been used to simulate the relationship between the number of PVS and the rating categories, taking into account the uncertainty in the measurements. The flowchart of the proposed optimization process is illustrated in Fig. 2.
In the remainer of this section we summarize the two visual rating scales used in our work and review the fundamentals of Frangi filtering and ordered logit modelling.

A. Visual Ratings of PVS
Two established visual rating scales for PVS severity were used in the present work.
On the basis of the measured contrast-to-noise ratios (CNRs), the PVS scores proposed by Patankar et al. [19] were based principally on the appearances seen on inversion recovery images. PVS should be scored in the centrum semiovale as 0 (none), 1 (less than five per side), 2 (more than five on one or both sides).
Two slightly modified versions of these rating methods, as described in [10], were also used in this work. Coregistered MRIs were used for assessment, with T2 for primary identification, T1 for confirmation, and Proton Density (PD) for rejection as required. To reduce double-counting, a slice increment of 3 was implemented as a standardized rating protocol. To reduce ceiling effects and account for a greater range of PVS, the Patankar scale was standardized across regions: 0 (none), 1 (one to five), 2 (six to ten), 3 (eleven to fifteen), 4 (sixteen or more). CS was defined as the White Matter (WM) projections superior to the ventricles, and while both scales normally include sub-insular regions lateral to the lentiform nucleus, we only included the thalamus, lentiform nucleus, caudate nucleus, and internal capsule for evaluation.

B. Frangi filter
Frangi [12] analyses the second order derivatives of an image I, defined in the Hessian matrix H s (v) as: where λ 1 , λ 2 and λ 3 are the ordered eigenvalues (|λ 1 | ≤ |λ 2 | ≤ |λ 3 |) of the Hessian matrix, , and α, β, c are thresholds which control the sensitivity of the filter to the measures R A , R B and S.
For a bright tubular structure in a 3D image we expect: For a dark structure λ 2 , λ 3 ≥ 0 and the conditions in Eq. 2 should be reversed.
Given a set of scales s ∈ [s min , s max ], the responses are combined as: where s min and s max are the minimum and maximum scales at which relevant structures are expected to be found [12].

C. Ordered Logit Models
An ordered logit model defines the relationship between an ordinal variable (y) which can vary between 0 and m(m ∈ N + ), and the vector of independent variables (x) by using a latent continuous variable (y ) defined in an one-dimensional space characterized by threshold points (µ 0 , . . . , µ m−1 ) as described in equation: where β and µ i are parameters to be estimated, is the error component which has a logistic random distribution with expected value equal to 0 and variance equal to π/ √ 3, that accounts for the measurement error. This modelling approach provides a relevant methodology for capturing the sources of influence (independent variables) that explain an ordinal variable (dependent variable) taking into account the measurement uncertainty of such data [16].
Since y is not a deterministic quantity, it is only possible to define the probability to belong to each class: where L is the logistic cumulative distribution function.
In our work, the ordinal variable (y) is the rating class (from 0 to 4) and the independent variables (x) is the number of PVS.

A. Model Calibration
The ordered logit model has been calibrated by maximizing a likelihood function based on a synthetic dataset generated in 3 steps. In the first step 1000 numbers of PVS Count (P C i , i = 1, · · · , 1000) have been generated using a lognormal distribution (see Fig. 3), that reflects the observed PVS distribution in known datasets [11]. In the second step, the uncertainty has been simulated for each P C i casting a New value of PVS Count (N P C i ) using a normal distribution with mean equal to P C i and standard deviation equal to one. Therefore, the probability that N P C i is included between P C i − 3 and P C i + 3 is 0.997. In the third step, a Rating Class (RC ij ) has been assigned to each generated N P C i . Assuming m classes, the log-likelihood function can be written as: where RC ij is equal to one if the i th generated number belong to the j th rating class and it is equal to zero otherwise.  For the modified Patankar scale [10] a rating class from 0 to 4, being 0(none), 1(1−5), 2(6−10), 3(11−15), 4(> 15) PVS, has been assigned to each generated number. The estimated parameters for the modified Patankar rating scale are β = 1.906, µ 0 = 2.269, µ 1 = 9.569, µ 2 = 18.995, µ 3 = 28.639, and the model is illustrated in Fig. 5.

B. Image Preprocessing
Images were preprocessed to generate the Region-of-Interest (ROI) masks. A fuzzy C-means clustering algorithm was applied to T1 images [20]. This is an unsupervised iterative clustering technique that effectively assigns each voxel to one of 4 membership classes: background, Cerebrospinal Fluid (CSF), Gray Matter (GM), and White Matter (WM). After a series of morphological and thresholding operations, the CSF and GM re-labelled voxels were combined to generate the final CSFGM mask which was used for false positive minimization. The CS was identified as the region of WM, superior to the lateral ventricles, present in each of the cerebral hemispheres under the cerebral cortex. In this paper we focused on the CS rather then the basal ganglia (BG), due to the availability of the ROI masks.
PVS masks were also available. These were obtained using Lesion Explorer [21], which implements 2 false positive minimization strategies: i) in order to reduce errors from minor imaging artifacts and improve differentiation from lacunar infarcts, candidate PVS are required to satisfy acceptance criteria from both T1 and T2, and rejection criteria from PD, and ii) to address potential registration errors and partial volume effects, the cortical GM segmentation from the CSFGM mask was dilated by 1 voxel. This resulted in a relatively conservative estimate of the overall PVS burden and thus, limited its utility as a ground truth for segmentation optimization.

C. Parameter Optimization
In order to apply the 3D Frangi filtering, the MRI volumes were first resliced to make 1mm isotropic voxels using linear interpolation. Then volumes have been filtered according to Eq. (2) and (3) and voxels having F(v) larger than a threshold t were kept. PVS were identified using 3D connected component analysis as the tubular structures with lengths between 3 and 50mm [3], [11]. This provided the initial PVS binary masks.
For each slice we calculated the PVS density as the area of the PVS mask divided by the area of the CS mask (obtained as described in Sec. IV-B). We automatically selected the slice in the CS with highest density of PVS. This slice corresponded to the representative slice having the highest number of PVS selected by the radiologist for assessing the Wardlaw visual ratings [7]. The count of PVS in this slice was derived automatically with 2D connected component labelling. Similarly, the total number of PVS in the entire CS was obtained with 3D connected component labelling. This count of PVS corresponded to the count performed by the radiologist for the Patankar ratings [19].
A log-likelihood function has been defined to optimize the segmentation parameters: Frangi filter scales s min , s max and threshold t. In this work, we used the default configuration for the other Frangi filter parameters (α = 0.5, β = 0.5, c = 500), as in our previous work, [22] we noted that optimizing these parameters produced essentially similar results, at the cost of a much higher computational time.
Based on the count of PVS (x i (s min , s max , t)) for each case i we obtained the probabilities of each case i to belong to the five rating classes (P (y = j|x i ), j = 0, . . . , 4) using the ordered logit model. The PVS visual rating category provided by an expert radiologist was then used to select a probability for each i case (P i ). The sum of the logarithms of these selected probabilities is the log-likelihood function to maximize: where N is the number of cases.

D. Model Validation
Segmentation procedures are commonly evaluated by assessing the voxel-wise spatial agreement between two binary masks, one obtained by the automatic method and a manual one. In our case, the manual segmentation of PVS was not available, as it would have been a very tedious and time consuming task to manually annotate these tiny structures in a reasonable size dataset. Quantitative comparison with other methods [9], [11] was unfeasible as they have been applied to MR images having different resolution, acquired using different protocols in different cohorts.
The performance of the models was therefore validated using Spearman's ρ (statistical analysis were performed using MATLAB Robust correlation toolbox [23]). Correspondence of PVS total count and volume vs. visual ratings was also assessed to test generalizability.

A. Experiments
For developing and optimizing the segmentation approach, the imaging datasets of 20 subjects were selected from a sample of the Sunnybrook Dementia Study (SDS) [10].
The optimization procedure has been applied to T1weighted and T2-weighted MRI sequences of the SDS dataset.
Frangi filter scales (s min and s max ) and two thresholds (t 1 and t 2 , one for each modality) have been optimized. The range of the parameters that undergo the optimization process has been defined as in Table II. PVS were assessed by three raters: two experienced neuroradiologists using the two modified visual rating scales as described in [10], and a third rater strictly following the guideline of the original Wardlaw [7] and Patankar [19] rating methods. The two ratings (modified Wardlaw and Patankar) of the first raters were close to the conservative estimate of PVS burden obtained as described in sec. IV-B. The third rater counted all visible PVS in T1 and T2, including the very small ones discarded by the first raters. All raters were blind to each other. Visual ratings are summarized in Table III. Three sets of experiments were performed as indicated in Table III. The symbol indicates the scale/rater used for optimization, while the symbol specifies those used for validation.
For illustration Fig. 6a and Fig. 6b show the surface plots for the parameter optimization using the modified Wardlaw rating scale. Fig. 6c and Fig. 6d show the surface plots for the parameter optimization using the modified Patankar rating scale. The optimal parameters obtained with the 2 models are very similar (s min = 1.4, s max = 3.2, t 1 = 0.96, t 2 = 0.35 for the first model, s min = 1.4, s max = 3.2, t 1 = 0.95, t 2 = 0.35 for the second one).
From the plots we can observe that the most significant parameter of the Frangi filter is the minimum scale (s min ). From the figures it is also clear that the Frangi filter was needed. Indeed for any combination of s min and s max the threshold values play a smaller role.
The optimal parameters obtained with the Wardlaw model using PVS assessed by the third rater where slightly different from the previous ones (s min = 0.2, s max = 2, t 1 = 0.96, t 2 = 0.1). The plots of the log-likelihood function are shown in Fig. 8.
The trend of plots confirms the validity of the model demonstrating that the model was able to adapt to the rater, and finds the best parameters to segment the PVS accounted by that rater.

B. Qualitative Evaluation
Magnified views of PVS segmentation using a thresholdbased method and the proposed method are shown in Fig. 9. It is clear that the proposed method detected most of the PVS, including the tiny ones, thanks to the enhancement of tubular structure performed by the Frangi filtering using the appropriate scale. The threshold based method missed them, as it was forced to be conservative in order to distinguish PVS from confounding tissue boundaries.
Examples of segmented PVS for two representative SDS cases having few and many PVS are shown in Fig. 10 and Fig. 11. For each case, we show T1, T2 and the PVS overlay in red.
Volume rendering of the segmented PVS for two cases having few and many PVS are shown in Fig. 12 for visual qualitative evaluation.

C. Quantitative Evaluation
When comparing segmentation results with the modified Wardlaw and Patankar visual ratings of the second rater, a high correlation was found for both methods (Spearman's ρ = 0.58, p = 0.006 and ρ = 0.71, p = 0.0004 respectively). However, low and no significant correlation was found with total PVS number in or volume CS, suggesting low generalizability. This replicates our previous analysis [10]. For the segmentation results obtained with the optimal parameters of the model optimized with the original Wardlaw scale a stronger correlation was found (Spearman's ρ = 0.74, p = 0.0002). In addition PVS total count and volume correlates with visual rating scores (Spearman's ρ = 0.67, p = 0.001 and ρ = 0.53, p = 0.015, respectively).

D. Application to alternative acquisitions
To validate the new PVS segmentation method we applied it to MRI of cases of the Mild Stroke Study (MSS). Visual ratings using the Wardlaw rating [7] were available for all the cases [4]. Automatic brain, cerebrospinal fluid (CSF) and normalappearing white matter extraction were performed on T1weighted MRI using optiBET [24] and FSL-FAST [25] respectively. All subcortical structures were segmented, also automatically, using other tools from the FMRIB Software Library (FSL) and an age-relevant template as per the pipeline described elsewhere [17]. After identifying the lateral ventricles as the CSF-filled structures with boundaries with the subcortical structures, the CS was identified as the region of  normal-appearing white matter, superior to the lateral ventricles, present in each of the cerebral hemispheres under the cerebral cortex. T1-weighted sequence and CS region were linearly registered to the T2-cube images [26], The optimization procedure has been applied to T2-cube MRI sequences of 20 patients, and tested on 48 patients of the same study.
The results of this experiment suggest high generalizability of the method. CS is more difficult to rate visually than BG, so future application of this method to BG may be more straightforward.

VI. CONCLUSIONS
We presented an automatic method for 3D segmentation of PVS in conventional MRI. The 3D Frangi filter enhances and captures the 3D geometrical shape of PVS, thus this method shows promise for identifying and quantifying PVS that run both longitudinally and transversally in the CS, avoiding the  double-counting limitations of slice-based methods. The novelty of this work is the fact that the ordered logit model allows use of the visual ratings for segmentation optimization in absence of alternative computational ground truth. The ordered logit model could deal with the measurement uncertainty and the unequal class intervals of the rating scores. One limitation of this method is that it relies on the image preprocessing step for the ROI masks. If the masks provided by this step are not accurate, the method can detect as PVS boundary of grey matter and giri. Another limitation of this method is that it requires high resolution and quasi isotropic structural MRI. Noisy images have been excluded for this study, otherwise any noise spot of tubular shape can be wrongly segmented as PVS. This can be overcome by a learning method. However, learning methods require GT, and not just visual ratings assessment.
The automatically segmented PVS count and volume agree with visual ratings. The method is fully automatic and therefore free from inter-and intra-rater variability. However, much more testing is required in a wider range of subjects including those with high burden of other ageing and neuroinflammation features. The quantitative assessment of PVS volume and count is more suitable for longitudinal studies than visual ratings, that tend to be susceptible to ceiling/flooring effects. The accurate segmentation of PVS will allow the analysis of their spatial distribution, orientation and density. Moreover, it will enable the study of the spatial and volumetric relationships of PVS with other markers of SVD, e.g. acute lacunar infarcts, white matter hyperintensities, lacunes, and microbleeds. Additionally, this method shows promise for use in longitudinal studies where PVS burden can be assessed in relation to mea-sures of cerebral blood brain barrier permeability, perfusion and cerebrovascular reactivity. Quantitative measurements will better characterize the severity of PVS in ageing people and their associations with dementia, stroke and vascular diseases. This is the first work to propose a multicentre study of PVS segmentation. It shows excellent multi-centre reproducibility.