1 Introduction

Cerebrovascular segmentation is significant to present the three-dimensional (3D) structure of cerebral vessels and is beneficial to clinical diagnosis and treatment process. Traditional segmentation mainly depends on the manual operation of clinicians, which has the disadvantages of poor reproducibility and low efficiency. Recently, Deep-learning methods [1] realized the optimal segmentation on many medical image analysis tasks, while it is very hard to obtain vascular ground-truths (GT) from complex image context. High-efficient statistical modeling benefits not only clinicians but also deep-learning enthusiasts. Thus, it is our main research motivation.

For cerebrovascular segmentation from Time-of-flight Magnetic Resonance Angiography (TOF-MRA) data, Moccia et al. [2] completely reviewed the existing segmentation methodology, including active contour model, threshold-based segmentation, Hessian-matrix based model, and statistical model. The last two methods draw our recent interest in developing a more accurate and robust model.

Hessian-Matrix Based Model:

It is often used to differentiate and enhance the ball-, tube-, and plane-like structures with the analysis of the second-order intensity derivatives or the Hessian-matrix eigenvalues. Hessian-based feature, i.e., eigenvalues and eigenvectors have been designed to construct a few of vessel functions by some researchers [3]. Especially, the recent vessel function proposed by Jerman et al. [4] presents satisfying contrast and the ratio of signal to noise on the resultant filtering response image. Although Hessian-matrix based model highlights the vessel-like object, it cannot directly lead an accurate segmentation.

Statistical Model:

An adaptive algorithm for cerebrovascular segmentation from TOF-MRA data is proposed by Wilson and Noble [5], where the finite mixture models (FMM) is a weighted combination of two Gaussian distributions and a uniform one. Wen et al. [6] split FMM into a Rayleigh distribution and two Gaussian ones. However, they did not consider the neighborhood constraint around a voxel, which resulted in vascular fragments and pseudo-vessels. Hassouna et al. [7] employed Maximum a Posterior & Markov random field (MAP-MRF) to relieve the FMM-processing noise and keep the segmentation continuity. Zhou et al. [8] further enriched the vascular structure by using statistical modeling on MAP–MRF and multi-pattern neighborhood system.

The above statistical methods only make limited improvement on their designated experimental dataset. As far as we know, no one has proposed a self-adapting statistic model to segment cerebral vessel in the TOF-MRA data from MR equipment with different scanner-types, imaging parameters, and individual difference of patients.

Our work contributes to the challenge in four aspects. Firstly, we present an efficient data standardization strategy to make the algorithm focus on the region without a skull, which helps statistic model adapt to segment cerebral vessel in the TOF-MRA data from anywhere and any patients. Secondly, a knowledge-based expectation maximization (EM) algorithm is proposed to estimate the Gaussian mixture models (GMM) parameters, which helps to obtain better parameters of the vascular distribution. Thirdly, the vascular probability feature map is extracted and a novel neighborhood energy function for MRF is obtained to further improve the segmentation performance. Lastly, our method is validated on three different datasets (total of 139 clinical data). Due to data being influenced by different factors, such as scanners, imaging parameters, and patients, our proposed method greatly improves the Dice similarity coefficients by 7.15%–38.20% and improves the sensitivity by 9.45%–48.41% than the existing methods. In terms of vascular network coverage, and environmental complexity, the visual comparison in between our method with the above traditional methods are shown in Fig. 1.

Fig. 1.
figure 1

The cerebrovascular segmentation results using different methods. The maximum intensity projections of the TOF-MRA without skull (a1 ~ a2). Subfigures (b–e) are the segmentation results of the proposed method, Zhou’s [8], Lu’s [9], and Wilson and Nobel [5].

Full size image

2 Method

The proposed method consists of four steps: (1) data standardization; (2) GMM and initialization; (3) GMM parameters estimation using knowledge-based EM algorithm; (4) a novel neighborhood constraint energy function for MAP-MRF.

2.1 Data Standardization

As a fact, most of cerebrovascular voxels distribute in the middle and high intensity, and the proportion of cerebrovascular voxels in intracranial volume is relatively small (about 1% ~ 5%). Another fact comes from numerous experiential validations. Namely, GMM accuracy is mostly dominated by the non-vascular region, whereas little by fitting the vascular regions. The existing methods often modeled the whole data space of original TOF-MRA data, where the non-vascular region inevitably facilitates the fitting errors and computational redundancy. In this respect, skull-stripping step is necessary, for which we use the FSL-BET [10] tool to greatly reduce the non-vascular region and stabilize the histogram of intracranial volume.

To improve the robustness, the histogram specification is employed to reduce the difference of the histograms of various TOF-MRA data. We use the proposed method without histogram specification to process all the TOF-MRA data. Then we select a TOF-MRA data that has the best visual vessel segmentation result and regard its intensity distribution as the specific histogram. After that, we standardize the histogram of the other TOF-MRA data to the specific histogram. Since discrete transformation may have a one-to-many relationship, a 3 × 3 × 3 Gaussian kernel (commonly used with an empirical variance of 0.4) is used to convolve with the TOF-MAR data, which increases the vascular intensity continuity as shown in Fig. 2.

Fig. 2.
figure 2

The process of data standardization illustrates the resultant histograms of (a) the original TOF-MRA data; (b) skull-stripped data; (c) the specific one and (d) the transformed one.

Full size image

2.2 GMM and Initialization

After data standardization, the resultant intensity distribution (shown in Fig. 2-a) facilitates to divide the voxels of skull-stripping data into one vessel class (cerebral vasculature) and two background classes (brain tissues and cerebrospinal fluid). Therefore, we use two Gaussian distributions to model the background classes, and another one for the vessel class, namely

$$ \left\{ {\begin{array}{*{20}l} {M\left( y \right) = \sum\nolimits_{i = 1}^{3} {w_{Gi} f_{Gi} \left( {y |u_{i} ,\sigma_{i} } \right)} } \hfill \\ {\sum\nolimits_{i = 1}^{3} {w_{Gi} = 1} } \hfill \\ {f_{Gi} \left( {y |u_{i} ,\sigma_{i} } \right) = \frac{1}{{\sqrt {2\pi } \sigma_{i} }}exp\left( {\frac{{ - \left( {y - u_{i} } \right)^{2} }}{{2\sigma_{i}^{2} }}} \right)\quad \,i \in \left[ {1,2,3} \right] } \hfill \\ \end{array} } \right. $$
(1)

where the functions \( f_{Gi} \left( {y |u_{i} ,\sigma_{i} } \right) \left( {i = 1,2,3} \right) \) are Gaussian distributions; \( f_{G3} \left( {y |u_{3} ,\sigma_{3} } \right) \) corresponds to the cerebrovascular region; \( f_{G2} \left( {y |u_{2} ,\sigma_{2} } \right) \) corresponds to the gray and white matters; \( f_{G1} \left( {y |u_{1} ,\sigma_{1} } \right) \) corresponds to cerebrospinal fluid and lateral ventricle; \( w_{Gi} \left( {i = 1,2,3} \right) \) are the weights of these categories; \( y \) is the observed intensity of each voxel in TOF-MRA data. K-means algorithm is used to initialize the GMM parameters using three points (i.e., a quarter of the peak value, peak value, and twice peak value) as the start points. The curve line of the GMM initialization is shown in Fig. 3-b.

Fig. 3.
figure 3

After data standardization, the intensity of skull-stripped TOF-MRA data (a) and the initial curve line of GMM (b).

Full size image

2.3 Knowledge-Based EM Algorithm for GMM Parameters Estimation

As we all know, the proportion of non-vascular voxels is much larger than that of cerebrovascular voxels in whole TOF-MRA volume, which will make the EM algorithm unable to estimate the parameters of vascular distribution fairly. Hence, we collect some vascular points using region growing algorithm, which would provide some knowledge for estimating vascular distribution. Then all the voxels \( D \) of skull-stripped volume are divided into un-labeled and labeled voxels (i.e., \( D_{u} \) and \( D_{l} \)). To make full use of labeled data and better fit the vascular distribution, a knowledge-based machine learning algorithm is used to estimate the GMM parameters and is iteratively updated as follows.

$$ \left[ {u_{i} } \right]_{k + 1} = \frac{{\mathop \sum \nolimits_{{y_{j} \in D_{u} }} \left[ {p\left( {G_{i} |y_{j} } \right)} \right]_{k} y_{j} + \mathop \sum \nolimits_{{y_{j} \in D_{li} }} y_{j} }}{{\mathop \sum \nolimits_{{y_{j} \in D_{u} }} \left[ {p\left( {G_{i} |y_{j} } \right)} \right]_{k} + N\left( {D_{li} } \right)}} $$
(2)
$$ \left[ {\sigma_{i}^{2} } \right]_{k + 1} = \frac{{\mathop \sum \nolimits_{{y_{j} \in D_{u} }} \left[ {p\left( {G_{i} |y_{j} } \right)} \right]_{k} \left( {y_{j} - \left[ {u_{i} } \right]_{k} } \right)^{2} + \mathop \sum \nolimits_{{y_{j} \in D_{li} }} \left( {y_{j} - \left[ {u_{i} } \right]_{k} } \right)^{2} }}{{\mathop \sum \nolimits_{{y_{j} \in D_{u} }} \left[ {p\left( {G_{i} |y_{j} } \right)} \right]_{k} + N\left( {D_{li} } \right)}} $$
(3)
$$ [w_{i} ]_{k + 1} = \frac{{\mathop \sum \nolimits_{{y_{j} \in D_{u} }} \left[ {p\left( {G_{i} |y_{j} } \right)} \right]_{k} + N\left( {D_{li} } \right)}}{N\left( D \right)} $$
(4)

where \( D_{li} \) is the labeled dataset of the \( i^{th} \) distribution. \( N\left( \cdot \right) \) is the number of voxels in the given set. According to Bayesian criterion, the posterior probability \( p\left( {G_{i} |x_{j} } \right) \) is calculated by:

$$ p\left( {G_{i} |y_{j} } \right) = \frac{{w_{Gi} f_{Gi} \left( {y_{j} |u_{i} ,\sigma_{i} } \right)}}{{\mathop \sum \nolimits_{i = 1}^{3} w_{Gi} f_{Gi} \left( {y_{j} |u_{i} ,\sigma_{i} } \right)}} $$
(5)

According to the maximum a posteriori (MAP) classification, if a voxel \( y_{j} \) meets \( p\left( {G_{3} |y_{j} } \right) - max\left\{ {p\left( {G_{i} |y_{j} } \right);i = 1,2} \right\} > 0 \), we infer it to be a vascular point initially.

2.4 Neighborhood Constraint with Vessel Prior

Probability Feature Map.

Hessian-based multi-scale filter [4] highlights the tubular structure and suppress background efficiently, but its filtering response cannot represent the probability of vessels. Without loss of generality, we facilitate the filtering response to be probabilistic using estimated GMM parameters. The probability feature map \( V_{f} \) is defined as:

$$ V_{f} \left( v \right) = \frac{1}{{1 + \left( {\frac{\tau }{v}} \right)^{2} }} $$
(6)

where \( v \) is multi-scale filtering response, and the parameter \( \tau \) is obtained by solving the following equations:

$$ \aleph \left( \tau \right) = w_{G3} \int_{V} {1 dv;\,\aleph \left( \tau \right)} = \int_{\tau }^{1} {} \int_{V} {\varepsilon \left( {v,\tau } \right)dvd\tau ;\,\varepsilon \left( {v,\tau } \right)} = \left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 1 & {v \ge \tau } \\ \end{array} } \\ {\begin{array}{*{20}c} 0 & {v < \tau } \\ \end{array} } \\ \end{array} } \right. $$
(7)

where \( w_{G3} \) is the estimated weight of the vascular Gaussian distribution in Eq. (1), and \( \aleph \left( \tau \right) \) presents the number of voxels that meet the condition. Figure 4 shows the obvious difference between the probability feature map and the multi-scale filtering response.

Fig. 4.
figure 4

Illustration of differences in between three kinds of image data with their local information on (a1) TOF-MRA data, (b1) multi-scale filtering response, and (c1) probability feature map. The subfigures in (a2 ~ c2) show the intensity-value boards within the blue boxed region of that in (a1 ~ c1). Note, the intensity-value in b2, c2 is multiplied by 100 for equilibrium display.

Full size image

Neighborhood Constraint.

In this study, we organize the novel six-neighborhood energy function by GMM process result and probability feature map. The new potential energy function \( U\left( {x_{r } ,r} \right) \) comes from two parts, not only comes from the initial label field obtained by GMM corresponding to \( \varphi_{1} \left( {x_{r} , x_{{ r^{*} }} } \right) \), but also comes from probability feature map corresponding to \( \varphi_{2} \left( { V_{{f_{r} }} , V_{{f_{{r^{*} }} }} } \right) \), which is defined as:

$$ U\left( {x_{r} ,r} \right) = \sum\nolimits_{{ r^{*} \in \eta_{r} }} {(\alpha_{1} \varphi_{1} \left( {x_{r} ,x_{{ r^{*} }} } \right) + \alpha_{2} \varphi_{2} \left( { V_{{f_{r} }} , V_{{f_{{r^{*} }} }} } \right))} $$
(8)

where \( r \) is the index of a point and \( r^{*} \in \eta_{r} \) is the index around six-neighborhood. where \( x_{r} \), \( x_{{ r^{*} }} \) takes a value from the label set \( L = \left\{ {L_{v} ,L_{b} } \right\} \) with \( L_{v} \) and \( L_{b} \) being the vessel class and the background class, respectively. \( V_{{f_{r} }} \) is the value of the \( r^{th} \) position in a probability feature map. Given that scale factors meet \( \alpha_{1} + \alpha_{2} = 1 \), and \( \alpha_{1} \) and \( \alpha_{2} \) is set to be 0.5. \( \varphi_{1} \left( {x_{r} ,x_{{ r^{*} }} } \right) \) and \( \varphi_{2} \left( { V_{{f_{r} }} , V_{{f_{{r^{*} }} }} } \right) \) are defined as:

$$ \varphi_{1} \left( {x_{r} ,x_{{ r^{*} }} } \right) = \left\{ {\begin{array}{*{20}c} {0,} & { if x_{r} = x_{{ r^{*} }} } \\ {1,} & { otherwise} \\ \end{array} } \right. $$
(9)
$$ \varphi_{2} \left( { V_{{f_{r} }} , V_{{f_{{r^{*} }} }} } \right) = \left| {V_{{f_{r} }} - V_{{f_{{r^{*} }} }} } \right| $$
(10)

Then, according to the MAP–MRF segmentation approach [8], the posterior probability of the vessel and background class is expressed as:

$$ p\left( {x_{r} = L_{v} |y_{r} } \right) \propto f_{G3} \left( {y_{r} |u_{3} ,\sigma_{3} } \right)exp\left( { - U\left( {x_{r} = L_{v} ,r} \right)} \right) $$
(11)

And the posterior probability of the is expressed by:

$$ p\left( {x_{r} = L_{b} |y_{r} } \right) \propto \frac{{\mathop \sum \nolimits_{i = 1}^{2} w_{Gi} f_{Gi} \left( {y_{r} |u_{i} ,\sigma_{i} } \right)}}{{\mathop \sum \nolimits_{i = 1}^{2} w_{Gi} }}exp\left( { - U\left( {x_{r} = L_{b} ,r} \right)} \right) $$
(12)

According to the MAP criterion, if a voxel meets \( p\left( {x_{r} = L_{v} |y_{r} } \right) > p\left( {x_{r} = L_{b} |y_{r} } \right) \), we infer it to be a vascular point.

3 Experiment and Result

3.1 Materials and Implementation Environment

Our experimental materials consist of three datasets which are acquired from several MR scanners with different imaging parameters for different individuals that have healthy or diseased cerebrovascular structures. (1) MIDAS-109: 109 normal TOF-MRA datasets are collected using SIEMENS ALLEGRA 3.0T MRI scanner (TR = 35.0, TE = 3.56, flip angle = 22) from a public dataset MIDAS [11]. Each volume has the size of 448 × 448 × 128 (voxels), and the resolution is uniformly 0.51 mm × 0.51 mm × 0.80 mm. (2) AVM-20: 20 abnormal TOF-MRA datasets with arteriovenous malformation (AVM) mass, are collected using GE Signa HDx 3.0T MRI scanner (TR = 25.0, TE = 3.5, flip angle = 20) from General Hospital of Southern Theater Command of Chinese People’s Liberation Army. Each volume has the size of 512 × 512 × 332 (voxels), and the resolution is uniformly 0.56 mm × 0.56 mm × 0.55 mm. (3) GZ-10: 10 normal TOF-MRA datasets are collected using GE Signa HDx 3.0T MRI scanner (TR = 30.0, TE = 3.2, flip angle = 15) from General Hospital of Southern Theater Command of PLA of China. Each volume has the size of 512 × 512 × 283 (voxels), and the resolution is uniformly 0.43 mm × 0.43 mm × 0.50 mm.

The algorithms are implemented by using MATLAB2018b, on a PC Intel® CoreTM i7-6850 k CPU @3.60 GHz*12, while the 3D visualization uses a hybrid rendering engine of VTK-8.0 with Visual Studio 2017.

3.2 Quantification Result

To make a quantitative evaluation, ten TOF-MRA data are selected from the above three datasets (Note, five from MIDAS-109, three from AVM-20, two from GZ-10). The cerebrovascular structures are manually segmented with voxel-by-voxel selections from the 10 TOF-MRA data by three neurosurgeons. The manual segmentations are treated as the ground truth to estimate the performance of models. Meanwhile, to derive a quantitative comparison of the segmentation results, three metrics are employed, i.e., dice similarity coefficient \( DSC = \frac{2TP}{2TP + FN + FP} \), sensitivity \( Sen = \frac{TP}{TP + FN} \) and positive predictive value \( PPV = \frac{TP}{TP + FP} \). where \( TP \), \( FP \), \( TN \), and \( FN \) denote the voxel number of true-positive, false-positive, true-negative and false-negative results, respectively.

Our experiments consisted of four methods in three datasets, and the qualitative evaluation results are presented in Table 1. The proposed method gets an average DSC score of 89.12%, an average Sen score of 83.72%, and an average PPV score of 95.66%, which outperforms all other methods. Besides, clinical evaluation proceeded through the observations of neurosurgeons is illustrated in Table 2. Their evaluation reports were divided into three-level cases. i.e., Good, general, and poor cases according to cerebrovascular network coverage, over-segmentation, and the percentage of pseudo-vascular structures. Cerebrovascular segmentation results are shown in Fig. 1, which shows that all these methods show similar performance on the large vessel segmentation. Comparing some segmentation details, it is obvious that our segmentation result with richer vascular network coverage than all the other methods.

Table 1. The qualitative evaluation results of four different methods for three datasets
Full size table
Table 2. Clinical evaluation of the three TOF-MRA data sets
Full size table

4 Conclusion

Our method steps contribute four parts for the MAP-MRF framework, i.e., data standardization, knowledge-based EM estimation of GMM parameters, and Markov high-level model with novel neighborhood constraint energy function, which make our proposed method more effective than the existing ones. Meanwhile, the proposed method is validated on three TOF-MRA datasets that contain 139 data totally and are acquired from different MR equipment, scanning parameters, and individuals. Extensive experimental results demonstrate that our method wins out of other ones in terms of the universality and accuracy for non-homologous TOF-MRA datasets. In addition, our method performs well on TOF-MAR data with AVM, which indicates a particular significance for computer-assisted clinical procedures. In the feature work, we will produce weak label on un-marked cerebrovascular data sets and conduct weak supervised learning, especially, we will further optimize and expand the method to automate vascular annotation and develop the deep-learning on vascular segmentation of more TOF-MRA datasets.