Keywords

1 Introduction

Recent randomized clinical trials revolutionized stroke care by confirming the benefit of modern endovascular therapy (EVT) in patients with acute ischemic stroke (AIS) due to proximal anterior circulation occlusion [1]. It is important for the EVT therapy being tested within the clinical trials to show the capability of saving brain tissue by measuring post-treatment infarct volume in addition to showing benefit across clinical scales, such as the modified Rankin Scale (mRS) and the National Institute of Health Stroke Scale (NIHSS). A recent study shows that the beneficial effect of EVT on functional outcome could be explained by a reduction in post-treatment infarct volume measured manually [2]. Post-treatment infarct of AIS patients typically includes ischemic infarct only. But in around 10% of AIS patients [3] hemorrhage occurs after endovascular interventions, including intracerebral hemorrhage (ICH) (Fig. 1) and subarachnoid hemorrhage (SAH). These are feared complications. In those patients, post-treatment ischemic infarct is accompanied by hemorrhage, leading to two components within and around the final infarct (Fig. 1(a)). Non-contrast CT scanning is widely available and is the most commonly used diagnostic tool for AIS patients. Manually contouring ischemic infarct and hemorrhage from the follow up NCCT images is tedious, time consuming, and observer dependent. Therefore, an auto- or semi-automated infarct segmentation from clinically used follow-up NCCT images would be desired.

However, it is challenging to segment ischemic infarct and hemorrhage from NCCT images, as shown in Fig. 1(a), due to: (1) low signal to noise ratio of clinically used NCCT images compared to MR images, (2) low contrast of brain soft tissue making anatomical structures less differentiable, (3) anisotropy of NCCT images with big thickness (5 mm) while having high in-plane resolution (\(0.625>0.625\) mm\(^2\)), (4) interference from atrophy, leukoaraiosis, and partial volume effect at the area around cerebrospinal fluid (CSF) that has similar intensity distributions to ischemic infarct, and intra-ventricular calcification and artefact showing the same appearance as hemorrhage. While previous studies mostly focused on segmenting either infarct [4] or hemorrhage [5] individually, these techniques cannot be directly applied on brain scans of the patients with these infarcts and hemorrhage occurring at the same time. Zimmerman et al. [6] manually measured intracerebral hemorrhage and infarct. Bardera et al. [7] and Gillebert et al. [4] developed two different models to segment hemorrhage and infarct sequentially.

Fig. 1.
figure 1

A post-treatment NCCT image of an AIS patient with manual segmentations of hemorrhage and ischemic infarct super-imposed in (a); (b): initialization to the algorithm, red: infarct, green: hemorrhage, blue: ventricle, yellow: background; (c) scheme of the proposed contour evolution for a single contour.

Contributions: In this work, segmenting hemorrhage and ischemic infarct is integrated into the same framework, for which, a novel joint segmentation approach is proposed for extracting these two different lesions simultaneously from the follow up NCCT images. The proposed approach makes use of multi-region time-implicit contour evolution combined with random forest (RF) learned semantic information. We show that the evolution of multiple contours can be globally optimized through convex optimization technique, leading to an efficient and parallelized algorithm.

2 Method

We propose a segmentation approach to simultaneously partition the image into 4 disjoint regions \(\mathcal {C}_i\), \(i=1 \ldots 4\), denoting ischemic infarct, hemorrhage, CSF, and other tissues, using multi-region time-implicit level-set evolution. It should be noted that the CSF was segmented as an additional region in order to remove the false positives caused by the partial volume effect of CSF.

2.1 Preprocessing

The input non-contrast CT images are first skull-stripped using an in-house software. Then, some voxels are sampled manually for each region, as shown in Fig. 1(b), to initialize the subsequent contour evolution. The user-sampled voxels are not only used to estimate prior intensity probability functions (PDFs) \(f_i(x)\) for each region, but also hard coded in the algorithm serving as constraints supervising the optimization scheme.

2.2 Multi-region Time-Implicit Level-Set Evolution

A fully time-implicit level-set scheme using global optimization [8, 9] has demonstrated great advantages on both implementation and computation compared to classical level-set methods [10]. In this work, we employ the evolution of 4 time-implicit level-sets \(\mathcal {C}_i\), \(i=1 \ldots n\) to simultaneously partition the image into 4 disjoint regions, under the constraint:

$$\begin{aligned} \varOmega \, = \, \cup _{i=1}^n \mathcal {C}_i \, , \quad \mathcal {C}_k \, \cap \, \mathcal {C}_l \, = \, \emptyset \, , \quad \forall k \, \ne \,l. \end{aligned}$$
(1)

where \(\varOmega \) means image domain, and \(n=4\) and \(i=1 \ldots 4\). Let \(\mathcal {C}_i^t\), \(i=1 \ldots n\), be the i-th region at the current time frame t, which moves to position \(\mathcal {C}_i^{t+1}\) at the next time frame \(t+1\). For each region \(\mathcal {C}_i^t\), \(i=1 \ldots n\), we define \(\mathcal {C}_i^+\) and \(\mathcal {C}_i^-\) indicating expansion and shrinkage of \(\mathcal {C}_i^t\) at time t with respect two costs of \(c_i^+(x)\) and \(c_i^-(x)\) (Fig. 1(c)).

With these definitions, the evolution of each region is intended to minimize the discrepancy between shrinkage and expansion over the discrete time frame while minimizing contour length:

$$\begin{aligned} \min _{\mathcal {C}_i} \; \sum _{i=1}^n\, \Big \{\int _{\mathcal {C}_i^-} c_i^-(x) \, dx \, + \, \int _{\mathcal {C}_i^+} c_i^+(x) \, dx \Big \}\, + \, \sum _{i=1}^n\,\int _{\partial \mathcal {C}_i} g(s) \, ds \end{aligned}$$
(2)

subject to (1), where g(s) is the weighting function along the contour boundaries. In general, the level-set \(\mathcal {C}\) evolution is driven by distance functions and image features for image segmentation, for which, the cost functions \(c_i^-(x)\) and \(c_i^+(x)\), w.r.t. region changes are typically defined as:

$$\begin{aligned} c^-_i(x) \, = \, c_i^+(x) \, = \, \mathrm {dist}(x, \partial \mathcal {C}_i^t)/h \, + F_i(x) \end{aligned}$$
(3)

where \(h > 0\) is constant, \(\mathrm {dist}(.)\) is the distance defining distances from x to the contour \(\partial \mathcal {C}_i^t\), and the image feature \(F_i(x)\) is derived according to the specified intensity distribution [11]. In this study, we use geodesic distance function \(gdist(x, \partial \mathcal {C})\) rather than traditional Euclidean distance function dist(.), as geodesic distance defines the distances from x to the contour \(\partial \mathcal {C}^t\) using the shortest path along the image intensities, which could provide more adaptive spatial contexts [12].

In particular, in addition to image features derived from images, an additional cost \(r_i(x)\) for each voxel is incorporated into the contour evolution, which includes semantic information derived from a random forest trained model [4]. The semantic information used for the random forest training and testing includes intensity information, and the mean and standard deviation of intensity in local regions, as well as the difference of those features between the current voxel and quasi-symmetry one from its contralateral side. An infarct occurrence spatial probability map obtained from training images served as an additional feature.

The cost functions are accordingly defined as follows:

$$\begin{aligned} \begin{aligned} c_i^+(x) \, = c_i^-(x) \, =&\, -\omega _1 \,\log f_i(x) \, - \, \omega _2 \, \log r_i(x) \,\\&+ \omega _3 \, \frac{1}{h} \,\mathrm {gdist}(x, \partial \mathcal {C}_i^t) \, \quad \forall x \, \in \, \mathcal {C}_i^t \end{aligned} \end{aligned}$$
(4)

where the weighting parameters \(\omega _1,\omega _2,\omega _3 > 0\), \(\omega _1+\omega _2+\omega _3= 1\) weight the contributions from the intensity PDFs \(f_i(x)\), random forest learned probability \(r_i(x)\), and geodesic distance for each voxel, respectively. The parameters \(\omega _i\) were empirically defined and kept as constants through evaluation stage. Specifically, we have implemented two independent random forest models trained with datasets from another study for classifying hemorrhage and ischemic infarct voxels, sequentially. There were no random forest models specially derived for CSF and other tissues that means only intensity and distance priors were used for segmenting these two regions.

2.3 Spatially Continuous Potts Model

The introduced variational problem (2) can be equally reformulated as a Potts problem [13]. For this purpose, we define two cost functions \(D_i^s(x)\) and \(D_i^q(x)\) w.r.t. the current contour \(\mathcal {C}_i^t\), \(i = 1. . .n\), at time t:

$$\begin{aligned} D_i^s \, := \, \left\{ \begin{array}{ll} c_i^-(x) \, , &{} \text {where} \,\, x \in \mathcal {C}_i^t \\ 0 \, , &{} \text {otherwise} \end{array} \right. \end{aligned}$$
(5)
$$\begin{aligned} D_i^q \, := \, \left\{ \begin{array}{ll} c_i^+(x) \, , &{} \text {where}\,\, x\notin \mathcal {C}_i^t \\ 0 \, , &{} \text {otherwise} \end{array}. \right. \end{aligned}$$
(6)

If \(u_i(x) \in \{0,1\}\), \(i=1\ldots n\), is defined as the indicator function of the region \(\mathcal {C}_i\), the disjoint constraint in (1) can be represented by

$$\begin{aligned} \sum _{i=1}^n\,u_i(x)=1; \, \quad u_i(x) \in \{ 0, 1\} \, \quad \forall x \in \varOmega \end{aligned}$$
(7)

Via the cost functions (5) and (6). It can be proven that the variational formulation (2) can be represented as the Potts model:

$$\begin{aligned} \min _{u_i(x) \in \{0,1\}} \; \sum _{i=1}^n\left\langle u_i, D_i^q - D_i^s \right\rangle \, + \, \sum _{i=1}^n\int _{\varOmega } g(x) \left|\nabla u_i\right| \, dx \end{aligned}$$
(8)

subject to the contour disjointness constraint (7), where the weighted length term in (2) is encoded as the weighted total-variation function (the second term in (8)).

We introduce convex relaxation technique to optimize the energy function (8) by solving a continuous min-cut/max-flow problem [9, 14], which implies that the position of the contour \(\mathcal {C}_i^{t+1}\) at the next time step \(t+1\) is globally optimal, and the given contour \(\mathcal {C}_i^{t}\) at the current time frame is always moving to its globally optimal position. In addition, the new contour position at each evolution step is computed in a time-implicit manner, which allows a large time-step and parallelization of the algorithm, substantially speeding up calculations.

3 Experiments and Results

Image Acquistion: The size of 3D NCCT images is \(512\times 512\times (18-37)\) with a voxel spacing from \(0.37\times 0.37\times 3 \) mm\(^3\) to \(0.49\times 0.49\times 5\) mm\(^3\). Sixteen stroke patients with hemorrhage and infarct after treatment from a clinical trial [1] were used for this study. Hemorrhage and infarct in each image were manually segmented slice by slice on parallel axial view using ITK-SNAP by a trained observer and verified by an experienced clinician. Manual segmentations were used to evaluate the proposed segmentation approach. In addition, another 30 patients (different form the used 16 patients) were used to train the random forest model and derive the lesion probability r(x).

Evaluation Metrics: Our segmentation method was evaluated by comparing the algorithm segmented results with manual segmentation in terms of the overlap of hemorrhage, infarct, and those two regions as a whole, respectively, using volume-based metrics: dice similarity coefficient (DSC); and distance-based metrics: the mean absolute surface distance (MAD) and maximum absolute surface distance (MAXD) [8, 9, 14, 15]. In order to show the efficacy of the proposed approach, two semi-automated segmentation methods: graph cut (GC) [16] and multi-region segmentation (MLS) [17], and a fullly automated random forest method followed by the optimization of Markov random field (RF+MRF) [5] were compared. For fair comparison, the initializations to our proposed method were also used as initializations to MLS and GC. It should be noted that the major part of the proposed method was implemented in Matlab, while the derivation of lesion probability map using random forest was implemented in Python. The methods used for comparison above were re-implemented based on the papers and publicly available code.

Results: The quantitative results in Table 1 show that the proposed approach yielded mean DSCs of 76.4% for ischemic infarct, 56.5% for hemorrhage, and 85.0% for the two regions with favorable values of MAD and MAXD. A statistical analysis, using one-way ANOVA followed by Tukey-Kramer correction regarding DSC, MAD, and MAXD,shows that the proposed method significantly outperformed the fully automated method RF+MRF and two other semi-automated methods :MLS and GC (\(p<0.001\) for all three comparisons using all three metrics). The Spearman correlation analysis between manually and algorithm segmented infarct volume shows that the Spearman correlation coefficients are 0.954 (\(p<0.0001\)) for ischemic infarct, 0.964 (\(p<0.0001\)) for hemorrhage, and 0.986 (\(p<0.0001\)) for the two regions, respectively (Fig. 2).

Fig. 2.
figure 2

Segmentation accuracy in terms of DSC in different regions using 16 patients.

Table 1. Overall performance results using 16 patient images in terms of DSC (%), MAD (mm), and MAXD (mm).

4 Discussion and Conclusion

A novel semi-automated segmentation approach is proposed, which is capable of extracting hemorrhage and ischemic infarct simultaneously from the follow-up NCCT images of post-treatment AIS patients. The quantitative results using 16 patients demonstrate its efficacy in terms of DSC, MAD, and MAXD. Compared to other state-of-the-art segmentation methods, the proposed approach improves segmentation accuracy greatly for ischemic infarct, hemorrhage, and the two regions as a whole.

There are few methods segmenting post-treatment ischemic infarct accompanied by hemorrhage from follow-up NCCT images. Most of studies segmented these two distinct pathologies individually. Compared to other ischemic infarct segmentation methods, the mean DSC of 76.4% obtained in this study is higher than a median of 74% (42–90%) [18] and 58–62% [4] in previous studies. Compared to other hemorrhage segmentation methods, the DSC of \( 56.5\pm 13.2\)% for the hemorrhage generated in this study is higher than \(55\pm 24\)% reported in [19], but lower than 62–78% [4] and 89.9% [5] in previous studies. However, in addition to the volume biased DSC value, the Spearman correlation coefficient of 0.964 (\(p<0.0001\)) demonstrated a strong correlation between the algorithm and manually segmented hemorrhage volume. Additionally, the proposed approach generated favorable DSC, MAD, and MAXD values provided both infarct and hemorrhage as a whole.

One limitation of this study is that the proposed approach was only validated on a small number of images. Extensive evaluation on a larger dataset is required to make it clinically applicable. In addition, given that the fully automated RF+MRF yielded poorer performance, the proposed method required user sampled voxels in 3–5 slices out of 28–44 slices in a CT image as initialization to improve performance even though it required additional \({<}1\) min. The introduced observer variability should be assessed in the future.

To summarize, we have proposed a segmentation approach to quantifying infarct volume of post-treatment AIS patients who present with ischemic infarct and hemorrhage simultaneously. The quantitative evaluations showed excellent accuracy and strong correlation with the gold standard of manual segmentation, suggesting that the proposed technique has the potential of being used to assess imaging outcomes in AIS patients after treatment.