- Email: [email protected]
Automated brain extraction from head CT and CTA images using convex optimization with shape propagation
Computer Methods and Programs in Biomedicine 176 (2019) 1–8
Contents lists available at ScienceDirect
Computer Methods and Programs in Biomedicine journal homepage: www.elsevier.com/locate/cmpb
Automated brain extraction from head CT and CTA images using convex optimization with shape propagation Mohamed Najm, Hulin Kuang, Alyssa Federico, Uzair Jogiat, Mayank Goyal, Michael D. Hill, Andrew Demchuk, Bijoy K. Menon, Wu Qiu∗ Department of Clinical Neuroscience, University of Calgary, Calgary, Alberta, Canada
a r t i c l e
i n f o
Article history: Received 26 March 2019 Revised 20 April 2019 Accepted 28 April 2019
Keywords: Brain extraction Convex optimization Geodesic level-set Non-contrast computer tomograph Computer tomograph angiogram
a b s t r a c t Background and objective: Non-Contrast Computer Tomography (NCCT) and CT angiography (CTA) are the most used and widely acceptable imaging modalities in clinical practice for the diagnosis and treatment of acute ischemic stroke (AIS) patients. Brain extraction of CT/CTA images plays an essential role in stroke imaging research. There is no robust automated brain extraction method in the literature that is well established for both NCCT and CTA images. Thus, a validated and automated brain extraction tool for CT imaging would be of great value for both research and clinical practice. Methods: The proposed brain extraction method is based on the contour evolution technique, which extracts brain tissues from acquired NCCT and CTA images in a slice-by-slice fashion. Specifically, the proposed approach makes use of a novel propagation framework, which is initialized by a localized slice with the largest brain section in axial views, followed by a geodesic level-set evolution for automatically extracting the brain section in each slice. In particular, the segmented contour propagated from the previous slice is reused to penalize the defined object function for contour evolution to enforce the shape continuity between any two adjacent contours. We show that the defined contour evolution function can be solved iteratively by globally optimal convex optimization. Results: The proposed brain extraction approach is quantitatively evaluated using 40 NCCT and CTA images acquired from 20 AIS patients and drawn from 4 different vendors, compared to manual segmentations using Dice and Jaccard coefficient metrics. The quantitative results show that the proposed segmentation algorithm is consistently accurate for both NCCT and CTA images using Dice metric. The proposed method is further validated on 1736 NCCT and CTA images of 1331 AIS patients acquired from three multi-national multi-centric clinical trials. A visual check performed on these data demonstrates a low failure rate of 0.4% for 1331 NCCT images and a zero-failure rate for 405 CTA images. Conclusions: Both quantitative and qualitative evaluation suggest that the proposed brain extraction approach for NCCT and CTA images can be used for different clinical imaging settings, thus serving to improve current image analysis in the field of neuroimaging. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Automated skull stripping or brain extraction is the process of isolating brain tissue on radiologic images. This technique has become a valuable tool for image processing in neuroimaging research [37]. The literature on brain extraction tools for magnetic resonance imaging (MRI) is well established, and the brain extraction tool (BET) from the FMRIB Software Library (FSL) is a robust and validated tool for brain extraction [37]. However, there is a gap in research pertaining to computer tomography (CT) brain extraction [22]. This is especially poignant given that CT is ∗
Corresponding author. E-mail address: [email protected] (W. Qiu).
https://doi.org/10.1016/j.cmpb.2019.04.030 0169-2607/© 2019 Elsevier B.V. All rights reserved.
more widely used, compared to MR imaging, in clinical practice for the diagnosis and treatment of stroke patients. Post processing of CT images is a necessary step in ischemic stroke lesion and haemorrhage segmentation on non contrast CT (NCCT) [2,3,5,35], collateral grading on CT angiography (CTA) [13], automated Alberta Stroke Program Early CT Score (ASPECTS) scoring [17] etc. Thus, a validated and automated brain extraction tool for CT imaging would be of great value for both research and clinical practice. There have been adaptations to the BET in order to utilize CT images, however these methods have not been validated against the gold standard of manual segmentation [32,38]. Mandell et al. [19] developed a brain extraction method for CT images that was validated; however, this method was not completely automated. Muschelli et al. [22] developed an automated and validated
2
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
method of brain extraction for CT images based on the BET with successful results. However, the algorithm created for this method was found to be parameter dependent (fraction intensity (FI)), and was only validated on NCCT, not on CTA images. Recently, deep learning technique based on convolutional neural network (CNN) has been applied on this application, demonstrating promising results [1]. This method however, was only validated on 22 NCCT images, and no results were reported on CTA images. In addition, like any other deep learning based segmentation methods [36], it requires a large amount of training images labeled appropriately, which is time and labour intensive. In this study, we propose a contour evolution based segmentation approach [18,23] for extracting brain tissues from acquired NCCT and CTA images in a slice by slice fashion. The proposed approach is based on a novel propagation framework, which is initialized by a localized slice with the largest brain section in axial views, followed by a variational time-implicit geodesic levelset evolution for automatically extracting the brain in each slice. The proposed contour evolution approach incorporates intensity probability density functions (PDFs) of brain tissue as well as a shape constraint propagated from the slice previously segmented. Specifically, the segmented contour propagated from the previous slices, reused as a shape constraint, is used to penalize the convex optimization method to assist segmenting the other slices, making use of the shape continuity between any two adjacent slices. We show that the defined contour evolution function can be iteratively solved by global optimal convex optimization. Moreover, the algorithm of convex optimization can be parallelized to obtain a high computational performance. The proposed skull stripping approach fills a gap in the established literature with extensive validation using clinically used NCCT and CTA images. The purpose of this paper is to provide an algorithm that is automated and validated against manual segmentation, and has high accuracy and robustness in brain extraction for both CTA and NCCT images.
Fig. 1. Approach overview.
2. Method 2.1. Algorithm overview Fig. 1 shows the pipeline of the proposed brain extraction method, which consists of the following steps: (1) The input 3D head CT or CTA axial images are first thresholded by a pre-defined thresholding values, followed by a component analysis in each slice along the z-axis. The generated component profile could help to localize an initial 2D slice with approximately the largest brain cross section. (2) 2D morphological operations and component analysis are used to segment a brain section from the localized initial slice. (3) The obtained segmentation result in the initial slice is then propagated to the head apex and base, respectively, for segmenting its two neighbouring slices, in which the segmented result in the initial slice is reused as both initial contour and the propagation shape constraint for the proposed convex optimization segmentation algorithm. This step is repeated until all 2D slices in a 3D image are segmented. (4) Post processing, including hole filling and small island removal, is finally used to refine the segmented brain.
Fig. 2. Component analysis along the z-axis.
The cases with a large oblique angle compared to a NCCT atlas in MNI space [32] were reoriented by registering them onto the atlas using affine transformation in the toolbox of NiftyReg [34]. Note that the affine transformation might change intracranial volume, so the developed technique might not be suitable for brain volume measurement [22]. Each NCCT/CTA image was thresholded using a thresholding Hounsfield unit (HU) value of 100 for NCCT images and 400 for CTA images. The thresholded images were then smoothed using Gaussian filter with a kernel size of σ = 1 mm3 . Intensity normalization typically performed for preprocessing MR images was not used in this application. [44] 2.3. Initial slice localization and segmentation
2.2. Preprocessing Although orientation correction is not necessary for most of cases, it is required for the cases with large gantry tilt. The cases with a large gantry tilt were recovered using dcm2nii tool [33].
Component analysis was performed on each slice after thresholding, therefore leading to the largest component in each slice reserved, which was assumed as the approximated brain cross sections. Fig. 2 shows an example of component analysis, in
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
which, the area of the largest component in each slice as well as the area difference of the largest component between any two adjacent slices are plotted. The slice with the largest area of the brain section was localized as the initial slice for the propagation strategy. After a 2D erosion with a disk-structuring element in the initial section, 2D connected component analysis was applied on the initial slice to extract the brain section, therefore generating two intensity probability density functions (PDFs) for foreground (brain) and background (skin, skull, muscle, etc.). These two PDFs serve as intensity priors to segment the remained 2D slices, which will be used in the contour evolution algorithm described below. 2.4. Propagation based segmentation using geodesic contour evolution with shape constraint 2.4.1. Time-implicit level-set evolution Recent studies [7,28] demonstrated a fully time-implicit levelset scheme using global optimization, which shows advantages over the classical level-set approach [10] in terms of both implementation and computation. Given an outer force f and meancurvature k, the level-set problem driven by mean-curvature [7]:
∂t C = −k + f
min C
C C t
1 d(x, ∂ C t ) dx + ds − f dx , h ∂C C
(2)
for each discrete time frame from t to t + h, where h > 0 is constant, and d(. ) indicates the signed distance of x to ∂ Ct , and C C t is the symmetric difference between the contour C and C t . The function (1) or (2) can be globally minimized by the continuous min-cut algorithm [9,24,42] or graph-cuts [6,8]. This implies that the level-set function can be evolved to a globally optimal position during each discrete time-frame. For image segmentation problem, the evolution of level-set is typically driven by the distance function as well as image feature. Generally, the cost functions c− (x ) and c+ (x ) associated with contour region changes are defined as the combination of image feature cost and distance function, such that
c (x ) = c (x ) = −
+
ω1 d(x, ∂ C )/h + ω2 f (x ) t
min C
c (x ) dx + −
C−
C+
c (x ) dx +
+
∂C
g(s ) ds
(4)
2.4.2. Propagation shape constraint for contour evolution In this study, in order to enforce the consistency between the current evolved brain contour and the propagated brain contour from previous slice, the difference between these two contours is also taken into account during the contour evolution. Specifically, an additional energy term is proposed to penalize the difference between the evolved contour and the propagated contour to account for the consistency of brain contours along the z-axis of 3D NCCT/CTA image. In particular, the distance between the voxel on the current contour ∂ C t and the contour C∗ propagated from the previous slice is minimized, assisting the evolution of the current contour, for which, the optimized energy function can be formulated as follows:
1
min C
C C t
x to the contour ∂ C t . Geodesic distance function measures the shortest path along the image intensities while considering gradient information, therefore giving more adaptive spatial contexts compared to Euclidean distance [11,28,30,39,40]. The cost function (3) can be extended as:
c− (x ) = c+ (x ) = − ω1 log f (x ) +
1 gd(x, ∂ C t ) + gd(x, ∂ C ∗ ) dx + ds h1 h2 ∂C
(5)
Geodesic distance function gd(x, C ) is used in this study, instead of Euclidean distance function d(. ), measuring the distance from
ω2
+ ω3
1 gd(x, ∂ C t ) h1 1 gd(x, ∂ C ∗ ) h2
(6)
∀x ∈ C t
in which the parameters ω1 , ω2 , ω3 > 0 (ω1 + ω2 + ω3 = 1) weight the contributions from three aspects: image intensity, geodesic distance, and shape constraint. f (x ) is denoted as intensity PDFs derived from the localized initial slice in Section 2.3. If L(x ) ∈ {0, 1} is the indicator function of the region C, the minimization of (5) can be equivalently formulated as a spatially continuous min-cut problem:
min
L(x )∈{0,1}
1 − L, Ds + L, Dt +
g(x )|∇ L| dx ,
(7)
where Ds and Dt are cost functions for the object and background pixels, respectively, which are associated with the current contour evolution [25,26,28–30], and defined as:
Ds :=
(3)
where f (x ) denotes the specified image feature, such as intensity distributions. Thus, the corresponding optimization function can be formulated as:
Fig. 3. Contour evolution with propagation constraint. (a) Propagation direction; (b) Contour evolution with shape constraint in single slice.
(1)
can be iteratively solved by minimizing the following energy function:
3
Dt :=
c− ( x ) , 0,
where x ∈ C t otherwise
(8)
c+ ( x ) , 0,
where x ∈ / Ct . otherwise
(9)
An efficient continuous max-flow algorithm is introduced to minimize the nonlinear energy function (7) [28,30]. If the binary constraint Li (x ) ∈ {0, 1} in (7) is relaxed by its convex relaxation Li (x ) ∈ [0, 1], optimizing the energy function (7) is converted into a continuous min-cut problem, which can be globally minimized [41]. The continuous min-cut problem indicates that the given contour can be evolved to its globally optimal position at each time step. In addition, the evolution of the contour at each step is time-implicitly calculated, which is beneficial for algorithm parallelization. We refer to Yuan et al. [41] for the configuration of the continuous max-flow used in this study. (Fig. 3) 3. Experiments 3.1. Image acquisition Forty NCCT and CTA images of 20 AIS patients at admission (median age: 73 years; Interquartile Range (IQR): 67–81 years; 60% male), one NCCT and one CTA image for each patient, were acquired for quantitative evaluations. The National Institutes of Health Stroke Scale (NIHSS) score at admission for the 20 patients is a median value of 17 (IQR, 12–19). Manual segmentations of the
4
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8 Table 1 Image information of 40 NCCT and CTA images acquired from 20 AIS patients used for quantitative evaluation. NCCT
CTA
Vendor
Image size
Image spacing (mm3 )
Image size
Image spacing (mm3 )
Siemens GE Philips Toshiba
512 × 512 × (18–33) 512 × 512 × (28–32) 512 × 512 × (30–49) 512 × 512 × (29–35)
0.41 × 0.41 × 3.1 − 0.47 × 0.47 × 5 0.45 × 0.45 × 5 0.37 × 0.37 × 3.1 − 0.49 × 0.49 × 5 0.43 × 0.43 × 5
512 × 512 × (37–109) 512 × 512 × (177–293) 512 × 512 × (268–378) 512 × 512 × (143–160)
0.36 × 0.36 × 1.5 − 0.40 × 0.40 × 4 0.45 × 0.45 × 0.63 − 0.49 × 0.49 × 0.63 0.41 × 0.41 × 0.5 − 0.49 × 0.49 × 5 0.43 × 0.43 × 1
40 images served as ground truth. The detailed image information is listed in Table 1. The small sample size used for quantitative comparison of the manual and automatic segmentations did not allow for exposure to many different image types. Thus, we opted to test the proposed algorithm on a larger dataset with a greater variety within the images. We then analyzed 1736 CT/CTA head (without neck) images, comprising of 1331 NCCT and 405 CTA baseline images. These images were acquired from 1331 stroke patients with different levels of stroke severity at admission. The patients were drawn from three multi-national multi-center hospital-based cohort studies that focus on the diagnosis and treatment of acute ischemic stroke patients. Those trials were Prove-IT (Precise and Rapid Assessment of Collaterals Using Multi-phase CTA in the Triage of Patients with Acute Ischemic Stroke for IA Therapy) [20], INTERRSECT (Identifying Nonocclusive and occlusive Thrombus characteristics associated with recanalization and ideal parenchymal and arterial imaging findings predictive of Early neurological Recovery after Recanalization using Serial CT angiography) [21], and ESCAPE (Endovascular treatment for Small Core and Anterior circulation Proximal occlusion with Emphasis on minimizing CT to recanalization times) [14]. These acquired images were scanned across more than 5 vendors in 22 centers in 5 countries, representing most of image types currently clinically used. 3.2. Quantitative evaluation metrics Our data sets were manually segmented by consensus readings of two experts. The brain were manually contoured through the axial planes on parallel slices using ITK-SNAP [43]. 2R ∩R We used region-based metrics, Dice = R M+R A [28,30] and JacR ∩R
M
A
card = RM ∪RA [15] coefficients, to evaluate the performance of the M A proposed segmentation approach, where RM and RA define the regions enclosed by the manual surface M and algorithm-segmented surface A. In addition, a visual check done by an imaging scientist with more than 10 years experience was performed for these 1736 images to qualitatively determine if the results generated by the proposed automatic segmentation was a ‘failure’, of ‘intermediate’, or of ‘good’ quality. ‘Good’ quality was defined as > 90% accuracy, ‘intermediate’ was defined as 70 − 90% accuracy, and ‘failure’ was defined as < 70% accuracy. In order to increase the reliability of visual check, a post doc stroke fellow was asked to review brain extraction results of 250 cases randomly selected from 1736 images. The column 1 and 2 in Fig. 4 show examples of being of ‘Good’, and the column 3 and 4 show examples of being of ‘intermediate’. The proposed skull stripping algorithm was also compared to two BET based methods developed in [22] using fractional intensity (FI) thresholds of 0.01 and 0.1, denoted as Method-0.01 and Method-0.1, and an atlas registration based method, the SwissSkullStripper method [4] in the software of Slicer [12]. Note that the Method-0.01 and Method-0.1 methods were implemented by thresholding and smoothing images before running the BET algorithm as suggested by [22].
3.3. Implementation details The proposed skull stripping algorithm was implemented in both CPU and GPU versions in Matlab (Natick, MA). The GPU based algorithm was developed using parallel computing architecture (CUDA, NVIDIA Corp., Santa Clara, CA). A non-optimized CPUimplemented Matlab code was also developed. All the experiments were conducted on a Windows desktop with an Intel Core i7-6700 CPU (2.6 GHz) and a graphics card of NVDIA Geforce 970M. It is worthy of noting that the CPU implementation is publicly available.1 All the parameters used in the algorithm were tuned using additional 3 NCCT and 3 CTA images with manual segmentations, which were excluded and different from 40 images used for evaluations. During the evaluation stage, all these parameters were kept constant. 4. Results 4.1. Quantitative evaluations For NCCT images, the quantitative analysis in Table 2 for the proposed algorithm reported mean Dice and Jaccard values of 0.965 ± 0.014 and 0.935 ± 0.025, respectively. This was followed by Method-0.1, Method-0.01, and lastly 3D Slicer that had the lowest Dice and Jaccard values. Fig. 4 shows an example of brain extraction of an NCCT image using these four different methods. The proposed algorithm produced the fastest computational time for brain extraction of NCCT images in terms of CPU computation. In particular, the GPU segmentation time was the fastest of all automated programs, with a mean time of 3.273 ± 0.611 s. This was followed by the CPU computation time of our method, and then Method-0.1, Method-0.01, and 3D Slicer being the slowest by a substantial amount. For CTA images, the quantitative analysis in Table 3 shows that all four methods generated reasonable Dice and Jaccard values. The proposed algorithm reported mean Dice and Jaccard values of 0.965 ± 0.018 and 0.932 ± 0, 065 respectively. The Method-0.1 closely followed this, followed by Method-0.01 and 3D Slicer that had the lowest Dice and Jaccard values. The obtained Jaccard and Dice values for the CTA images were much more similar than those obtained for the NCCT images. Method-0.1 produced the fastest time for automated brain extraction of CTA images (19.603 s). This was followed closely by Method-0.01 and our GPU-implemented version (22.924 s), then the 3D Slicer. Our CPU-implemented version and 3D slicer were substantially slower than the rest. 4.2. Qualitative evaluations Visual check performed on 1736 NCCT and CTA images in Table 4 show that the proposed brain extraction method is able to generate good segmentation results for ≈ 98% of NCCT and CTA images. The failure rate is only 0.53% for 1331 NCCT images, 1
https://github.com/WuChanada/StripSkullCT.git.
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
5
Fig. 4. Brain extraction results of an NCCT image, using our proposed method, Method-0.1, Method-0.01, and Slicer with DSCs in parentheses. Row 1–4: axial, sagittal, coronal, and volume rendering views. Red: exact overlap between manual and algorithm segmentations; blue: over-estimation by algorithm; green: under-estimation by algorithm. Table 2 Quantitative evaluation results using 20 NCCT images acquired from 20 AIS patients, in terms of Dice, Jaccard, and computational time (s) when compared to the ground truth manual segmentation. All values are represented as mean ± standard deviation. In particular, the proposed method has a CPU time and a GPU time (the bottom value in bracket) recorded.
Dice Jaccard Time (s) (CPU and GPU)
Our method
Method-0.1 [22]
Method-0.01 [22]
3D Slicer
0.965 ± 0.014 0.935 ± 0.025 7.807 ± 1.847 (3.273 ± 0.611)
0.899 ± 0.034 0.818 ± 0.06 7.913 ± 0.581
0.877 ± 0.038 0.784 ± 0.06 8.158 ± 0.976
0.859 ± 0.037 0.761 ± 0.057 39.24 ± 4.653
Table 3 Quantitative evaluation results using 20 CTA images acquired from 20 AIS patients, in terms of Dice, Jaccard, and computational time (s) when compared to the ground truth manual segmentation. All values are represented as mean ± standard deviation. In particular, the proposed method has a CPU time and a GPU time (the bottom value in bracket) recorded.
Dice Jaccard Time (s) (CPU and GPU)
Our method
Method-0.1 [22]
Method-0.01 [22]
3D Slicer
0.965 ± 0.018 0.932 ± 0.032 48.795 ± 27.921 (22.924 ± 12.749)
0.947 ± 0.024 0.913 ± 0.043 19.60 ± 6.657
0.916 ± 0.05 0.83 ± 0.08 20.51 ± 6.971
0.861 ± 0.036 0.755 ± 0.056 45.095 ± 6.082
Table 4 Qualitative evaluation results of using 1331 NCCT images and 405 CTA images.
Good quality Intermediate quality Failure
NCCT (%)
CTA (%)
In all (%)
97.74(1301/1331) 1.73(23/1331) 0.53(7/1331)
98.52(399/405) 1.48(6/405) 0
97.93(1700/1736) 1.67(29/1736) 0.40(7/1736)
and it worked successfully for all 405 CTA images. For the 250 images randomly selected out of 1736 images, the review results performed by the post doc stroke fellow are consistent to the results done by the image scientist. As multiple scanners were used in this qualitative evaluation, we also analyzed the rate of failure and intermediate scans for each vendor in Table 5. There
6
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8 Table 5 Rates of failure and intermediate quality scans for each vendor of brain extraction of the 1736 scans. Vendor
Rate of failure and intermediate scans (%)
GE Siemens Philips Toshiba Canon
0.020 (14/694) 0.014 (6/417) 0.025 (7/278) 0.022 (5/232) 0.035 (4/115)
was no significant difference observed ( p = 0.89) across the five different vendors by a post-hoc pairwise Fisher’s exact test. In addition, in order to increase the reliability of qualitative evaluation, another 40 images randomly selected from these 1736 patient images were manually segmented as ground truth to quantitatively evaluate the performance of the proposed method regarding the metric of Dice further. This additional experimental results show that the proposed method is capable of generating a Dice of 0.957 ± 0.024 for these 40 images, consistent to the quantitative results in Tables 2 and 3. 4.3. Comparisons with other methods Although brain extraction has been well established for MRI, leading to some successful tools, such as BET, there is no software that is extensively validated on head CT and CTA images. The BET method was adapted on extracting brains from NCCT images, and validated on 1095 scans from 126 patients in [22]. The results in [22] showed that the adapted BET method was able to generate great accuracy for their datasets with certain FI values. However, two BET adapted methods: Method-0.1 and Method-0.01 developed in the paper worked inconsistently for NCCT and CTA images in our experiments, only generated Dice values of < 0.90 for the 20 NCCT images in our experiments. When the Method-0.1 and Method-0.01 methods were applied on CTA images with thinner slice thickness, the generated Dice of > 0.9 was decent, but lower than that generated by the proposed method. Another brain extraction method, the SwissSkullStripper method [4] in the Slicer software, was also evaluated using the same data. This method yielded consistent results for both NCCT and CTA images with Dice values of ≈ 0.86. For the compared three methods, most segmentation errors exist at the base region of the brain, as shown in the first row in Fig. 4. The big difference of slice thickness between NCCT and CTA images, ranging from 0.5 mm to 5 mm in Table 1, might be the reason resulting in segmentation accuracy discrepancy for the Method-0.1 and Method-0.01 methods. In contrast to these two methods, the proposed approach is insensitive to slice thickness in the range of [0.5, 5] mm, performing consistently well for both modality images. We additionally performed another sensitivity analysis by changing the slice thickness of 5 CTA images randomly selected from 20 AIS patients used for quantitative evaluation, through resampling images in z axis with thickness of 5, 10, 15, 20, 25 mm. It was found that the average DSC was consistently >= 90% when slice thickness was <= 15 mm. However, the DSC decreased to <= 75% when slice thickness was >= 20 mm. 5. Discussion A brain extraction method for NCCT and CTA images is proposed in this study, which is based on a novel contour evolution technique with a shape propagation strategy. The proposed approach was quantitatively evaluated using 40 NCCT and CTA images of 20 stroke patients (one NCCT and one CTA scan for each patient), which were acquired from four different vendors (Toshiba, GE, Philips, and Siemens) with five NCCT and five CTA
Fig. 5. Four challenging cases successfully segmented by the proposed method.
scans from each vendor. Compared to manual segmentations as ground truth, the proposed brain extraction method was able to generate high Dice and Jaccard values for both NCCT and CTA images, outperforming the benchmark methods, two BET adapted methods: Method - 0.1 and Method - 0.01, and 3D Slicer. Although the scans used for the quantitative validation were obtained from different patients, different centers, and different scanners, the sample size was relatively small, which does not allow for exposure to many different image types, such as brain scans with, such as hemorrhage stroke, craniotomy, surgery, as well as scans for neonates and children. Further evaluation on 1736 NCCT and CTA images acquired from 1331 stroke patients was performed. These images were acquired from three multi-national and multi-centric clinical trials representing a diverse range of clinically relevant images with different levels of noise, artifacts, motion, and scanning parameters. The high degree of segmentation accuracy for both NCCT and CTA images in qualitative analysis indicates that the proposed algorithm can maintain accurate and specific results in a larger dataset. Fig. 5 shows a few cases, which were successfully processed by the proposed method, but not by other three methods. In particular, cut-off scans and the scans with irregular brain shape are quite challenging for the BET and Slicer methods. Additionally, after the qualitative evaluations, the processed 1331 NCCT images excluding 7 failure cases have been subsequently used to develop and validate an automated Alberta Stroke Program Early CT Score (ASPECTS) system [17], and the processed 405 CTA images have been used to develop a computer aided diagnosis system for characterizing intracranial thrombus predicting early recanalization with intravenous alteplase [27]. We do believe that through both quantitative and qualitative evaluations, the proposed brain extraction method is shown to be accurate and reliable for both NCCT and CTA scans. Another advantage of the proposed method is that it is a stand-alone algorithm without any dependency on any other third-party library or atlas, such as the BET toolbox. Thus, it can serve as a fundamental image processing step, therefore being integrated into any other algorithms seamlessly, such as spatial normalization [16]. More importantly, the CPU implementation of the developed algorithm is publicly accessible to benefit more for the community. The developed algorithm was the fastest automated method of brain extraction for NCCT images even with the CPU imple-
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
7
bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patentlicensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Acknowledgements
Fig. 6. Two examples of failure cases.
mentation. Our GPU implementation sped up calculation for NCCT images by at least two times, which is much faster than the other three methods. However, the computational time of the proposed method increased substantially for CTA images. Our CPU implementation for CTA images is the slowest compared to the other three methods. However, the GPU version is comparable to the BET method and faster than the Slicer. The Slicer software provided the slowest solution for both modality images as it relies heavily on image registration. Since our proposed method extracts brains in a slice by slice fashion, the greater slice number in CTA images dramatically slows down the calculation. However, the computation speed can be potentially accelerated. For instance, the computation of the data cost functions Ds,t (x ) in (8) and (9) and edge indicator g (x ) in (7) were calcuated using serial Matlab code and run on the CPU. This could be sped up by parallelizing Matlab code or converting it to C language. In addition, based on the segmentation results on NCCT images with greater slice thickness, it is reasonable to believe that the computation time for CTA images could be reduced greatly by enlarging the propagation interval between slices, while maintaining similar segmentation accuracy. However, finding the optimal interval and interpolating segmentations between slices are required to be further explored. We acknowledge that the proposed brain extraction method has some limitations. For example, poor image quality caused by excessive movement and incomplete acquisition of brain, as shown in Fig. 6 affects the proposed method.In addition, the images used in this study were without craniotomy and surgery, which would potentially increase the failure rate of the proposed method. Another limitation of the proposed method is that it cannot deal with images with large gantry tilt directly. For these cases, linear registration onto an atlas is required before processing to spatially recover brain orientation, which increases computational time accordingly.
Ethics approval was obtained from the local institutional review boards. Dr. Menon holds the Heart and Stroke/University of Calgary Professorship in Stroke Imaging. Dr. Goyal has a patent pending on systems of stroke diagnosis using multiphase computed tomographic angiography and a licensing agreement with GE healthcare for the same. Dr. Goyal also provided consulting services to Covidien for design and conduct of SWIFT PRIME trial and for teaching engagements. Dr. Hill and Dr. Goyal have a research grant with Medtronic (Covidien) to conduct clinical trials. The other authors report no conflicts. Appendix A. Pseudo code for the developed algorithm
Algorithm 1 Brain extraction. Require: Input image: I, slice number: S Ith = threshold (I ) for k = 1, A(k)=Cal cul ateArea(Ith (k )), while k < S [, Nmid ] = max(A(: )) UNmid = Local izeInitial Sl ice(Ith (Nmid )) while N = 0 do if N == Nmid then U prev = UNmid else U prev = UN+1 end if UN = Ext ract Brain2D(I (N ), U prev ) N ←N−1 end while while N < S do if N == Nmid then U prev = UNmid else U prev = UN−1 end if UN = Ext ract Brain2D(I (N ), U prev ) N ←N+1 end while U = Hol eF il l 3D(U )
6. Conclusion To summarize, an accurate, robust, and efficient method of brain extraction from NCCT and CTA images is proposed in this study. This method has demonstrated a high degree of accuracy when compared to manual segmentation. Qualitative evaluations of the method show that the specificity of the proposed algorithm can be maintained with a larger dataset. The proposed method has the potential of being used in different clinical settings, thus serving to improve current CT and CTA based image analysis in the field of neuroimaging.
Appendix B. Continuous Max-flow model and algorithm A flow maximization configuration is set up in a spatially continuous setting as shown in [25,26]. Let ps (x) and pt (x) be source and sink flows to and from pixel x to the source and sink terminals. The spatial flow q(x) exists around the neighborhood of pixel x. Then, the continuous max-flow model that maximizes the flow can be formulated as follows:
max P ( ps , pt , q ) =
ps ,pt ,q
ps dx
(10)
Conflicts of interest
subject to the following flow capacity and conservation constraints:
All authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’
| q ( x ) | ≤ g( x ) , p s ( x ) ≤ D s ( x ) , (divq − ps + pt )(x ) = 0.
pt (x ) ≤ Dt (x ),
8
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
The continuous max-flow model (10) can be solved using the augmented Lagrangian method [31]:
ps dx + u, div q − ps + pt −
c div q − ps + pt 2 2
(11)
leading to an algorithm as following iterations:
q
k+1
2
c uk := arg max − div q + pkt − pks − , 2 c q∞ ≤α
pkt +1
pks +1 := arg max ps
u
k+1
2
c uk k+1 k , := arg max − p + div q − p − t s 2 c pt (x )≤ρ (l,x )
ps dx −
2
k
c ps − pkt +1 + div qk+1 + u , 2 c
= u − c div qk+1 − pks +1 + pkt +1 . k
Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cmpb.2019.04.030. References [1] Z. Akkus, P.M. Kostandy, K.A. Philbrick, B.J. Erickson, Extraction of brain tissue from CT head images using fully convolutional neural networks, in: Medical Imaging 2018: Image Processing, 10574, International Society for Optics and Photonics, 2018, p. 1057420. [2] F. Al Ajlan, A. Al Sultan, P. Minhas, Z. Assis, M. de Miquel, M. Millán, L. San Román, A. Tomassello, A. Demchuk, T. Jovin, et al., Posttreatment infarct volumes when compared with 24-hour and 90-day clinical outcomes: insights from the revascat randomized controlled trial, Am. J. Neuroradiol. (2018). [3] F.S. Al Ajlan, M. Goyal, A.M. Demchuk, P. Minhas, F. Sabiq, Z. Assis, R. Willinsky, W.J. Montanera, J.L. Rempel, A. Shuaib, et al., Intra-arterial therapy and post-treatment infarct volumes, Stroke 47 (3) (2016) 777–781. [4] S. Bauer, T. Fejes, M. Reyes, A skull-stripping filter for ITK, Insight J. 2012 (2013). [5] A. Boers, I. Zijlstra, C. Gathier, R. van den Berg, C.H. Slump, H. Marquering, C. Majoie, Automatic quantification of subarachnoid hemorrhage on noncontrast CT, Am. J. Neuroradiol. 35 (12) (2014) 2279–2286. [6] Y. Boykov, V. Kolmogorov, An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision, IEEE Trans. Pattern Anal. Mach. Intell. 26 (2001) 359–374. [7] Y. Boykov, V. Kolmogorov, D. Cremers, A. Delong, An integral solution to surface evolution PDES via geo-cuts, in: In ECCV, 2006, pp. 409–422. [8] Y. Boykov, O. Veksler, R. Zabih, Fast approximate energy minimization via graph cuts, IEEE Trans. Pattern Anal. Mach. Intell. 23 (11) (2001) 1222–1239. ¯ , M. Nikolova, Algorithms for finding global minimizers of [9] T.F. Chan, S. Esedoglu image segmentation and denoising models, SIAM J. Appl. Math. 66 (5) (2006) 1632–1648(electronic). [10] T.F. Chan, L.A. Vese, An efficient variational multiphase motion for the Mumford–Shah segmentation model, in: Proc. Conf Signals, Systems and Computers Record of the Thirty-Fourth Asilomar Conf, 1, 20 0 0, pp. 490–494. [11] A. Criminisi, T. Sharp, A. Blake, Geos: Geodesic Image Segmentation, in: Computer Vision–ECCV 2008, Springer, 2008, pp. 99–112. [12] A. Fedorov, R. Beichel, J. Kalpathy-Cramer, J. Finet, J.-C. Fillion-Robin, S. Pujol, C. Bauer, D. Jennings, F. Fennessy, M. Sonka, et al., 3d slicer as an image computing platform for the quantitative imaging network, Magn. Reson. Imaging 30 (9) (2012) 1323–1341. [13] A.M. Frölich, S.L. Wolff, M.N. Psychogios, E. Klotz, R. Schramm, K. Wasser, M. Knauth, P. Schramm, Time-resolved assessment of collateral flow using 4d CT angiography in large-vessel occlusion stroke, Eur. Radiol. 24 (2) (2014) 390–396. [14] M. Goyal, A.M. Demchuk, B.K. Menon, M. Eesa, J.L. Rempel, J. Thornton, D. Roy, T.G. Jovin, R.A. Willinsky, B.L. Sapkota, et al., Randomized assessment of rapid endovascular treatment of ischemic stroke, N. Engl. J. Med. 372 (11) (2015) 1019–1030. [15] P. Jaccard, The distribution of the flora in the alpine zone., New Phytol. 11 (2) (1912) 37–50. [16] J.S. Kim, H. Cho, J.Y. Choi, S.H. Lee, Y.H. Ryu, C.H. Lyoo, M.S. Lee, Feasibility of computed tomography-guided methods for spatial normalization of dopamine transporter positron emission tomography image, PLoS One 10 (7) (2015) e0132585. [17] H. Kuang, M. Najm, D. Chakraborty, N. Maraj, S. Sohn, M. Goyal, M. Hill, A. Demchuk, B. Menon, W. Qiu, Automated aspects on noncontrast ct scans in patients with acute ischemic stroke using machine learning, Am. J. Neuroradiol. 40 (1) (2019) 33–38.
[18] Z. Ma, J.M.R. Tavares, R.N. Jorge, T. Mascarenhas, A review of algorithms for medical image segmentation and their applications to the female pelvic cavity, Comput. Methods Biomech. Biomed. Eng. 13 (2) (2010) 235–246. [19] J.G. Mandell, J.W. Langelaan, A.G. Webb, S.J. Schiff, Volumetric brain analysis in neurosurgery: part 1. particle filter segmentation of brain and cerebrospinal fluid growth dynamics from mri and CT images, J. Neurosurg. 15 (2) (2015) 113–124. [20] B.K. Menon, C.D. d Esterre, E.M. Qazi, M. Almekhlafi, L. Hahn, A.M. Demchuk, M. Goyal, Multiphase ct angiography: a new tool for the imaging triage of patients with acute ischemic stroke, Radiology 275 (2) (2015) 510–520. [21] B.K. Menon, M. Najm, F. Al-Ajlan, J.P. Alcantara, D. Dowlatshahi, A. Calleja, S.-I. Sohn, S.H. Ahn, A. Poppe, R. Mikulik, et al., IV TPA recanalization rates by site of occlusion and time after TPA bolus-main results of the interrsect multinational multicenter prospective cohort study, 2017. [22] J. Muschelli, N.L. Ullman, W.A. Mould, P. Vespa, D.F. Hanley, C.M. Crainiceanu, Validated automatic brain extraction of head CT images, Neuroimage 114 (2015) 379–385. [23] D.L. Pham, C. Xu, J.L. Prince, Current methods in medical image segmentation, Annu. Rev. Biomed. Eng. 2 (1) (20 0 0) 315–337. [24] T. Pock, A. Chambolle, H. Bischof, D. Cremers, A convex relaxation approach for computing minimal partitions, in: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2009. Miami, Florida [25] W. Qiu, Y. Chen, J. Kishimoto, S. de Ribaupierre, B. Chiu, A. Fenster, B. Menon, J. Yuan, Longitudinal analysis of pre-term neonatal cerebral ventricles from 3D ultrasound images using spatial-temporal deformable registration, IEEE Trans. Med. Imaging (2017). [26] W. Qiu, Y. Chen, J. Kishimoto, S.d. Ribaupierre, B. Chiu, A. Fenster, J. Yuan, Automatic segmentation approach to extracting neonatal cerebral ventricles from 3D ultrasound images, Med. Image Anal. 35 (2017) 181–191. [27] W. Qiu, H. Kuang, J. Nair, Z. Assis, M. Najm, C. McDougall, B. McDougall, K. Chung, A. Wilson, M. Goyal, et al., Radiomics-based intracranial thrombus features on CTand CTA predict recanalization with intravenous alteplase in patients with acute ischemic stroke, Am. J. Neuroradiol. 40 (2019) 39–44. [28] W. Qiu, J. Yuan, M. Rajchl, J. Kishimoto, Y. Chen, S. de Ribaupierre, B. Chiu, A. Fenster, 3D MR Ventricle segmentation in pre-term infants with post-hemorrhagic ventricle dilatation (phvd) using multi-phase geodesic level-sets, Neuroimage 118 (2015) 13–25. [29] W. Qiu, J. Yuan, E. Ukwatta, A. Fenster, Rotationally resliced 3D prostate trus segmentation using convex optimization with shape priors, Med. Phys. 42 (2) (2015) 877–891. [30] W. Qiu, J. Yuan, E. Ukwatta, Y. Sun, M. Rajchl, A. Fenster, Prostate segmentation: an efficient convex optimization approach with axial symmetry using 3D TRUS and MR images, IEEE Trans. Med. Imaging 33 (4) (2014) 1–14. [31] R.T. Rockafellar, Augmented Lagrangians and applications of the proximal point algorithm in convex programming, Math.Oper.Res. 1 (2) (1976) 97–116. [32] C. Rorden, L. Bonilha, J. Fridriksson, B. Bender, H.-O. Karnath, Age-specific CT and MRI templates for spatial normalization, Neuroimage 61 (4) (2012) 957–965. [33] C. Rorden, M. Brett, Stereotaxic display of brain lesions, Behav. Neurol. 12 (4) (20 0 0) 191–20 0. [34] D. Rueckert, L.I. Sonoda, C. Hayes, D.L. Hill, M.O. Leach, D.J. Hawkes, Nonrigid registration using free-form deformations: application to breast MR images, Med. Imaging, IEEE Trans. 18 (8) (1999) 712–721. [35] P.W. Schaefer, L. Souza, S. Kamalian, J.A. Hirsch, A.J. Yoo, S. Kamalian, R.G. Gonzalez, M.H. Lev, Limited reliability of computed tomographic perfusion acute infarct volume measurements compared with diffusion-weighted imaging in anterior circulation stroke, Stroke 46 (2) (2015) 419–424. [36] D. Shen, G. Wu, H.-I. Suk, Deep learning in medical image analysis, Annu. Rev. Biomed. Eng. 19 (2017) 221–248. [37] S.M. Smith, Fast robust automated brain extraction, Hum. Brain Mapp. 17 (3) (2002) 143–155. [38] J. Solomon, V. Raymont, A. Braun, J.A. Butman, J. Grafman, User-friendly software for the analysis of brain lesions (able), Comput. Methods Programs Biomed. 86 (3) (2007) 245–254. [39] P.J. Toivanen, New geodesic distance transforms for gray-scale images, Pattern Recognit. Lett. 17 (5) (1996) 437–450. [40] Z. Wang, K.K. Bhatia, B. Glocker, A. Marvao, T. Dawes, K. Misawa, K. Mori, D. Rueckert, Geodesic Patch-based Segmentation, in: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2014, Springer, 2014, pp. 666–673. [41] J. Yuan, E. Bae, X. Tai, A study on continuous max-flow and min-cut approaches, CVPR, 2010. [42] J. Yuan, E. Bae, X.-C. Tai, Y. Boykov, A continuous max-flow approach to Potts model, ECCV, 2010. [43] P.A. Yushkevich, J. Piven, H. Cody Hazlett, R. Gimpel Smith, S. Ho, J.C. Gee, G. Gerig, User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability, Neuroimage 31 (3) (2006) 1116–1128. [44] Y. Zhuge, J.K. Udupa, Intensity standardization simplifies brain MR image segmentation, Comput. Vision Image Underst. 113 (10) (2009) 1095–1103.
Contents lists available at ScienceDirect
Computer Methods and Programs in Biomedicine journal homepage: www.elsevier.com/locate/cmpb
Automated brain extraction from head CT and CTA images using convex optimization with shape propagation Mohamed Najm, Hulin Kuang, Alyssa Federico, Uzair Jogiat, Mayank Goyal, Michael D. Hill, Andrew Demchuk, Bijoy K. Menon, Wu Qiu∗ Department of Clinical Neuroscience, University of Calgary, Calgary, Alberta, Canada
a r t i c l e
i n f o
Article history: Received 26 March 2019 Revised 20 April 2019 Accepted 28 April 2019
Keywords: Brain extraction Convex optimization Geodesic level-set Non-contrast computer tomograph Computer tomograph angiogram
a b s t r a c t Background and objective: Non-Contrast Computer Tomography (NCCT) and CT angiography (CTA) are the most used and widely acceptable imaging modalities in clinical practice for the diagnosis and treatment of acute ischemic stroke (AIS) patients. Brain extraction of CT/CTA images plays an essential role in stroke imaging research. There is no robust automated brain extraction method in the literature that is well established for both NCCT and CTA images. Thus, a validated and automated brain extraction tool for CT imaging would be of great value for both research and clinical practice. Methods: The proposed brain extraction method is based on the contour evolution technique, which extracts brain tissues from acquired NCCT and CTA images in a slice-by-slice fashion. Specifically, the proposed approach makes use of a novel propagation framework, which is initialized by a localized slice with the largest brain section in axial views, followed by a geodesic level-set evolution for automatically extracting the brain section in each slice. In particular, the segmented contour propagated from the previous slice is reused to penalize the defined object function for contour evolution to enforce the shape continuity between any two adjacent contours. We show that the defined contour evolution function can be solved iteratively by globally optimal convex optimization. Results: The proposed brain extraction approach is quantitatively evaluated using 40 NCCT and CTA images acquired from 20 AIS patients and drawn from 4 different vendors, compared to manual segmentations using Dice and Jaccard coefficient metrics. The quantitative results show that the proposed segmentation algorithm is consistently accurate for both NCCT and CTA images using Dice metric. The proposed method is further validated on 1736 NCCT and CTA images of 1331 AIS patients acquired from three multi-national multi-centric clinical trials. A visual check performed on these data demonstrates a low failure rate of 0.4% for 1331 NCCT images and a zero-failure rate for 405 CTA images. Conclusions: Both quantitative and qualitative evaluation suggest that the proposed brain extraction approach for NCCT and CTA images can be used for different clinical imaging settings, thus serving to improve current image analysis in the field of neuroimaging. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Automated skull stripping or brain extraction is the process of isolating brain tissue on radiologic images. This technique has become a valuable tool for image processing in neuroimaging research [37]. The literature on brain extraction tools for magnetic resonance imaging (MRI) is well established, and the brain extraction tool (BET) from the FMRIB Software Library (FSL) is a robust and validated tool for brain extraction [37]. However, there is a gap in research pertaining to computer tomography (CT) brain extraction [22]. This is especially poignant given that CT is ∗
Corresponding author. E-mail address: [email protected] (W. Qiu).
https://doi.org/10.1016/j.cmpb.2019.04.030 0169-2607/© 2019 Elsevier B.V. All rights reserved.
more widely used, compared to MR imaging, in clinical practice for the diagnosis and treatment of stroke patients. Post processing of CT images is a necessary step in ischemic stroke lesion and haemorrhage segmentation on non contrast CT (NCCT) [2,3,5,35], collateral grading on CT angiography (CTA) [13], automated Alberta Stroke Program Early CT Score (ASPECTS) scoring [17] etc. Thus, a validated and automated brain extraction tool for CT imaging would be of great value for both research and clinical practice. There have been adaptations to the BET in order to utilize CT images, however these methods have not been validated against the gold standard of manual segmentation [32,38]. Mandell et al. [19] developed a brain extraction method for CT images that was validated; however, this method was not completely automated. Muschelli et al. [22] developed an automated and validated
2
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
method of brain extraction for CT images based on the BET with successful results. However, the algorithm created for this method was found to be parameter dependent (fraction intensity (FI)), and was only validated on NCCT, not on CTA images. Recently, deep learning technique based on convolutional neural network (CNN) has been applied on this application, demonstrating promising results [1]. This method however, was only validated on 22 NCCT images, and no results were reported on CTA images. In addition, like any other deep learning based segmentation methods [36], it requires a large amount of training images labeled appropriately, which is time and labour intensive. In this study, we propose a contour evolution based segmentation approach [18,23] for extracting brain tissues from acquired NCCT and CTA images in a slice by slice fashion. The proposed approach is based on a novel propagation framework, which is initialized by a localized slice with the largest brain section in axial views, followed by a variational time-implicit geodesic levelset evolution for automatically extracting the brain in each slice. The proposed contour evolution approach incorporates intensity probability density functions (PDFs) of brain tissue as well as a shape constraint propagated from the slice previously segmented. Specifically, the segmented contour propagated from the previous slices, reused as a shape constraint, is used to penalize the convex optimization method to assist segmenting the other slices, making use of the shape continuity between any two adjacent slices. We show that the defined contour evolution function can be iteratively solved by global optimal convex optimization. Moreover, the algorithm of convex optimization can be parallelized to obtain a high computational performance. The proposed skull stripping approach fills a gap in the established literature with extensive validation using clinically used NCCT and CTA images. The purpose of this paper is to provide an algorithm that is automated and validated against manual segmentation, and has high accuracy and robustness in brain extraction for both CTA and NCCT images.
Fig. 1. Approach overview.
2. Method 2.1. Algorithm overview Fig. 1 shows the pipeline of the proposed brain extraction method, which consists of the following steps: (1) The input 3D head CT or CTA axial images are first thresholded by a pre-defined thresholding values, followed by a component analysis in each slice along the z-axis. The generated component profile could help to localize an initial 2D slice with approximately the largest brain cross section. (2) 2D morphological operations and component analysis are used to segment a brain section from the localized initial slice. (3) The obtained segmentation result in the initial slice is then propagated to the head apex and base, respectively, for segmenting its two neighbouring slices, in which the segmented result in the initial slice is reused as both initial contour and the propagation shape constraint for the proposed convex optimization segmentation algorithm. This step is repeated until all 2D slices in a 3D image are segmented. (4) Post processing, including hole filling and small island removal, is finally used to refine the segmented brain.
Fig. 2. Component analysis along the z-axis.
The cases with a large oblique angle compared to a NCCT atlas in MNI space [32] were reoriented by registering them onto the atlas using affine transformation in the toolbox of NiftyReg [34]. Note that the affine transformation might change intracranial volume, so the developed technique might not be suitable for brain volume measurement [22]. Each NCCT/CTA image was thresholded using a thresholding Hounsfield unit (HU) value of 100 for NCCT images and 400 for CTA images. The thresholded images were then smoothed using Gaussian filter with a kernel size of σ = 1 mm3 . Intensity normalization typically performed for preprocessing MR images was not used in this application. [44] 2.3. Initial slice localization and segmentation
2.2. Preprocessing Although orientation correction is not necessary for most of cases, it is required for the cases with large gantry tilt. The cases with a large gantry tilt were recovered using dcm2nii tool [33].
Component analysis was performed on each slice after thresholding, therefore leading to the largest component in each slice reserved, which was assumed as the approximated brain cross sections. Fig. 2 shows an example of component analysis, in
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
which, the area of the largest component in each slice as well as the area difference of the largest component between any two adjacent slices are plotted. The slice with the largest area of the brain section was localized as the initial slice for the propagation strategy. After a 2D erosion with a disk-structuring element in the initial section, 2D connected component analysis was applied on the initial slice to extract the brain section, therefore generating two intensity probability density functions (PDFs) for foreground (brain) and background (skin, skull, muscle, etc.). These two PDFs serve as intensity priors to segment the remained 2D slices, which will be used in the contour evolution algorithm described below. 2.4. Propagation based segmentation using geodesic contour evolution with shape constraint 2.4.1. Time-implicit level-set evolution Recent studies [7,28] demonstrated a fully time-implicit levelset scheme using global optimization, which shows advantages over the classical level-set approach [10] in terms of both implementation and computation. Given an outer force f and meancurvature k, the level-set problem driven by mean-curvature [7]:
∂t C = −k + f
min C
C C t
1 d(x, ∂ C t ) dx + ds − f dx , h ∂C C
(2)
for each discrete time frame from t to t + h, where h > 0 is constant, and d(. ) indicates the signed distance of x to ∂ Ct , and C C t is the symmetric difference between the contour C and C t . The function (1) or (2) can be globally minimized by the continuous min-cut algorithm [9,24,42] or graph-cuts [6,8]. This implies that the level-set function can be evolved to a globally optimal position during each discrete time-frame. For image segmentation problem, the evolution of level-set is typically driven by the distance function as well as image feature. Generally, the cost functions c− (x ) and c+ (x ) associated with contour region changes are defined as the combination of image feature cost and distance function, such that
c (x ) = c (x ) = −
+
ω1 d(x, ∂ C )/h + ω2 f (x ) t
min C
c (x ) dx + −
C−
C+
c (x ) dx +
+
∂C
g(s ) ds
(4)
2.4.2. Propagation shape constraint for contour evolution In this study, in order to enforce the consistency between the current evolved brain contour and the propagated brain contour from previous slice, the difference between these two contours is also taken into account during the contour evolution. Specifically, an additional energy term is proposed to penalize the difference between the evolved contour and the propagated contour to account for the consistency of brain contours along the z-axis of 3D NCCT/CTA image. In particular, the distance between the voxel on the current contour ∂ C t and the contour C∗ propagated from the previous slice is minimized, assisting the evolution of the current contour, for which, the optimized energy function can be formulated as follows:
1
min C
C C t
x to the contour ∂ C t . Geodesic distance function measures the shortest path along the image intensities while considering gradient information, therefore giving more adaptive spatial contexts compared to Euclidean distance [11,28,30,39,40]. The cost function (3) can be extended as:
c− (x ) = c+ (x ) = − ω1 log f (x ) +
1 gd(x, ∂ C t ) + gd(x, ∂ C ∗ ) dx + ds h1 h2 ∂C
(5)
Geodesic distance function gd(x, C ) is used in this study, instead of Euclidean distance function d(. ), measuring the distance from
ω2
+ ω3
1 gd(x, ∂ C t ) h1 1 gd(x, ∂ C ∗ ) h2
(6)
∀x ∈ C t
in which the parameters ω1 , ω2 , ω3 > 0 (ω1 + ω2 + ω3 = 1) weight the contributions from three aspects: image intensity, geodesic distance, and shape constraint. f (x ) is denoted as intensity PDFs derived from the localized initial slice in Section 2.3. If L(x ) ∈ {0, 1} is the indicator function of the region C, the minimization of (5) can be equivalently formulated as a spatially continuous min-cut problem:
min
L(x )∈{0,1}
1 − L, Ds + L, Dt +
g(x )|∇ L| dx ,
(7)
where Ds and Dt are cost functions for the object and background pixels, respectively, which are associated with the current contour evolution [25,26,28–30], and defined as:
Ds :=
(3)
where f (x ) denotes the specified image feature, such as intensity distributions. Thus, the corresponding optimization function can be formulated as:
Fig. 3. Contour evolution with propagation constraint. (a) Propagation direction; (b) Contour evolution with shape constraint in single slice.
(1)
can be iteratively solved by minimizing the following energy function:
3
Dt :=
c− ( x ) , 0,
where x ∈ C t otherwise
(8)
c+ ( x ) , 0,
where x ∈ / Ct . otherwise
(9)
An efficient continuous max-flow algorithm is introduced to minimize the nonlinear energy function (7) [28,30]. If the binary constraint Li (x ) ∈ {0, 1} in (7) is relaxed by its convex relaxation Li (x ) ∈ [0, 1], optimizing the energy function (7) is converted into a continuous min-cut problem, which can be globally minimized [41]. The continuous min-cut problem indicates that the given contour can be evolved to its globally optimal position at each time step. In addition, the evolution of the contour at each step is time-implicitly calculated, which is beneficial for algorithm parallelization. We refer to Yuan et al. [41] for the configuration of the continuous max-flow used in this study. (Fig. 3) 3. Experiments 3.1. Image acquisition Forty NCCT and CTA images of 20 AIS patients at admission (median age: 73 years; Interquartile Range (IQR): 67–81 years; 60% male), one NCCT and one CTA image for each patient, were acquired for quantitative evaluations. The National Institutes of Health Stroke Scale (NIHSS) score at admission for the 20 patients is a median value of 17 (IQR, 12–19). Manual segmentations of the
4
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8 Table 1 Image information of 40 NCCT and CTA images acquired from 20 AIS patients used for quantitative evaluation. NCCT
CTA
Vendor
Image size
Image spacing (mm3 )
Image size
Image spacing (mm3 )
Siemens GE Philips Toshiba
512 × 512 × (18–33) 512 × 512 × (28–32) 512 × 512 × (30–49) 512 × 512 × (29–35)
0.41 × 0.41 × 3.1 − 0.47 × 0.47 × 5 0.45 × 0.45 × 5 0.37 × 0.37 × 3.1 − 0.49 × 0.49 × 5 0.43 × 0.43 × 5
512 × 512 × (37–109) 512 × 512 × (177–293) 512 × 512 × (268–378) 512 × 512 × (143–160)
0.36 × 0.36 × 1.5 − 0.40 × 0.40 × 4 0.45 × 0.45 × 0.63 − 0.49 × 0.49 × 0.63 0.41 × 0.41 × 0.5 − 0.49 × 0.49 × 5 0.43 × 0.43 × 1
40 images served as ground truth. The detailed image information is listed in Table 1. The small sample size used for quantitative comparison of the manual and automatic segmentations did not allow for exposure to many different image types. Thus, we opted to test the proposed algorithm on a larger dataset with a greater variety within the images. We then analyzed 1736 CT/CTA head (without neck) images, comprising of 1331 NCCT and 405 CTA baseline images. These images were acquired from 1331 stroke patients with different levels of stroke severity at admission. The patients were drawn from three multi-national multi-center hospital-based cohort studies that focus on the diagnosis and treatment of acute ischemic stroke patients. Those trials were Prove-IT (Precise and Rapid Assessment of Collaterals Using Multi-phase CTA in the Triage of Patients with Acute Ischemic Stroke for IA Therapy) [20], INTERRSECT (Identifying Nonocclusive and occlusive Thrombus characteristics associated with recanalization and ideal parenchymal and arterial imaging findings predictive of Early neurological Recovery after Recanalization using Serial CT angiography) [21], and ESCAPE (Endovascular treatment for Small Core and Anterior circulation Proximal occlusion with Emphasis on minimizing CT to recanalization times) [14]. These acquired images were scanned across more than 5 vendors in 22 centers in 5 countries, representing most of image types currently clinically used. 3.2. Quantitative evaluation metrics Our data sets were manually segmented by consensus readings of two experts. The brain were manually contoured through the axial planes on parallel slices using ITK-SNAP [43]. 2R ∩R We used region-based metrics, Dice = R M+R A [28,30] and JacR ∩R
M
A
card = RM ∪RA [15] coefficients, to evaluate the performance of the M A proposed segmentation approach, where RM and RA define the regions enclosed by the manual surface M and algorithm-segmented surface A. In addition, a visual check done by an imaging scientist with more than 10 years experience was performed for these 1736 images to qualitatively determine if the results generated by the proposed automatic segmentation was a ‘failure’, of ‘intermediate’, or of ‘good’ quality. ‘Good’ quality was defined as > 90% accuracy, ‘intermediate’ was defined as 70 − 90% accuracy, and ‘failure’ was defined as < 70% accuracy. In order to increase the reliability of visual check, a post doc stroke fellow was asked to review brain extraction results of 250 cases randomly selected from 1736 images. The column 1 and 2 in Fig. 4 show examples of being of ‘Good’, and the column 3 and 4 show examples of being of ‘intermediate’. The proposed skull stripping algorithm was also compared to two BET based methods developed in [22] using fractional intensity (FI) thresholds of 0.01 and 0.1, denoted as Method-0.01 and Method-0.1, and an atlas registration based method, the SwissSkullStripper method [4] in the software of Slicer [12]. Note that the Method-0.01 and Method-0.1 methods were implemented by thresholding and smoothing images before running the BET algorithm as suggested by [22].
3.3. Implementation details The proposed skull stripping algorithm was implemented in both CPU and GPU versions in Matlab (Natick, MA). The GPU based algorithm was developed using parallel computing architecture (CUDA, NVIDIA Corp., Santa Clara, CA). A non-optimized CPUimplemented Matlab code was also developed. All the experiments were conducted on a Windows desktop with an Intel Core i7-6700 CPU (2.6 GHz) and a graphics card of NVDIA Geforce 970M. It is worthy of noting that the CPU implementation is publicly available.1 All the parameters used in the algorithm were tuned using additional 3 NCCT and 3 CTA images with manual segmentations, which were excluded and different from 40 images used for evaluations. During the evaluation stage, all these parameters were kept constant. 4. Results 4.1. Quantitative evaluations For NCCT images, the quantitative analysis in Table 2 for the proposed algorithm reported mean Dice and Jaccard values of 0.965 ± 0.014 and 0.935 ± 0.025, respectively. This was followed by Method-0.1, Method-0.01, and lastly 3D Slicer that had the lowest Dice and Jaccard values. Fig. 4 shows an example of brain extraction of an NCCT image using these four different methods. The proposed algorithm produced the fastest computational time for brain extraction of NCCT images in terms of CPU computation. In particular, the GPU segmentation time was the fastest of all automated programs, with a mean time of 3.273 ± 0.611 s. This was followed by the CPU computation time of our method, and then Method-0.1, Method-0.01, and 3D Slicer being the slowest by a substantial amount. For CTA images, the quantitative analysis in Table 3 shows that all four methods generated reasonable Dice and Jaccard values. The proposed algorithm reported mean Dice and Jaccard values of 0.965 ± 0.018 and 0.932 ± 0, 065 respectively. The Method-0.1 closely followed this, followed by Method-0.01 and 3D Slicer that had the lowest Dice and Jaccard values. The obtained Jaccard and Dice values for the CTA images were much more similar than those obtained for the NCCT images. Method-0.1 produced the fastest time for automated brain extraction of CTA images (19.603 s). This was followed closely by Method-0.01 and our GPU-implemented version (22.924 s), then the 3D Slicer. Our CPU-implemented version and 3D slicer were substantially slower than the rest. 4.2. Qualitative evaluations Visual check performed on 1736 NCCT and CTA images in Table 4 show that the proposed brain extraction method is able to generate good segmentation results for ≈ 98% of NCCT and CTA images. The failure rate is only 0.53% for 1331 NCCT images, 1
https://github.com/WuChanada/StripSkullCT.git.
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
5
Fig. 4. Brain extraction results of an NCCT image, using our proposed method, Method-0.1, Method-0.01, and Slicer with DSCs in parentheses. Row 1–4: axial, sagittal, coronal, and volume rendering views. Red: exact overlap between manual and algorithm segmentations; blue: over-estimation by algorithm; green: under-estimation by algorithm. Table 2 Quantitative evaluation results using 20 NCCT images acquired from 20 AIS patients, in terms of Dice, Jaccard, and computational time (s) when compared to the ground truth manual segmentation. All values are represented as mean ± standard deviation. In particular, the proposed method has a CPU time and a GPU time (the bottom value in bracket) recorded.
Dice Jaccard Time (s) (CPU and GPU)
Our method
Method-0.1 [22]
Method-0.01 [22]
3D Slicer
0.965 ± 0.014 0.935 ± 0.025 7.807 ± 1.847 (3.273 ± 0.611)
0.899 ± 0.034 0.818 ± 0.06 7.913 ± 0.581
0.877 ± 0.038 0.784 ± 0.06 8.158 ± 0.976
0.859 ± 0.037 0.761 ± 0.057 39.24 ± 4.653
Table 3 Quantitative evaluation results using 20 CTA images acquired from 20 AIS patients, in terms of Dice, Jaccard, and computational time (s) when compared to the ground truth manual segmentation. All values are represented as mean ± standard deviation. In particular, the proposed method has a CPU time and a GPU time (the bottom value in bracket) recorded.
Dice Jaccard Time (s) (CPU and GPU)
Our method
Method-0.1 [22]
Method-0.01 [22]
3D Slicer
0.965 ± 0.018 0.932 ± 0.032 48.795 ± 27.921 (22.924 ± 12.749)
0.947 ± 0.024 0.913 ± 0.043 19.60 ± 6.657
0.916 ± 0.05 0.83 ± 0.08 20.51 ± 6.971
0.861 ± 0.036 0.755 ± 0.056 45.095 ± 6.082
Table 4 Qualitative evaluation results of using 1331 NCCT images and 405 CTA images.
Good quality Intermediate quality Failure
NCCT (%)
CTA (%)
In all (%)
97.74(1301/1331) 1.73(23/1331) 0.53(7/1331)
98.52(399/405) 1.48(6/405) 0
97.93(1700/1736) 1.67(29/1736) 0.40(7/1736)
and it worked successfully for all 405 CTA images. For the 250 images randomly selected out of 1736 images, the review results performed by the post doc stroke fellow are consistent to the results done by the image scientist. As multiple scanners were used in this qualitative evaluation, we also analyzed the rate of failure and intermediate scans for each vendor in Table 5. There
6
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8 Table 5 Rates of failure and intermediate quality scans for each vendor of brain extraction of the 1736 scans. Vendor
Rate of failure and intermediate scans (%)
GE Siemens Philips Toshiba Canon
0.020 (14/694) 0.014 (6/417) 0.025 (7/278) 0.022 (5/232) 0.035 (4/115)
was no significant difference observed ( p = 0.89) across the five different vendors by a post-hoc pairwise Fisher’s exact test. In addition, in order to increase the reliability of qualitative evaluation, another 40 images randomly selected from these 1736 patient images were manually segmented as ground truth to quantitatively evaluate the performance of the proposed method regarding the metric of Dice further. This additional experimental results show that the proposed method is capable of generating a Dice of 0.957 ± 0.024 for these 40 images, consistent to the quantitative results in Tables 2 and 3. 4.3. Comparisons with other methods Although brain extraction has been well established for MRI, leading to some successful tools, such as BET, there is no software that is extensively validated on head CT and CTA images. The BET method was adapted on extracting brains from NCCT images, and validated on 1095 scans from 126 patients in [22]. The results in [22] showed that the adapted BET method was able to generate great accuracy for their datasets with certain FI values. However, two BET adapted methods: Method-0.1 and Method-0.01 developed in the paper worked inconsistently for NCCT and CTA images in our experiments, only generated Dice values of < 0.90 for the 20 NCCT images in our experiments. When the Method-0.1 and Method-0.01 methods were applied on CTA images with thinner slice thickness, the generated Dice of > 0.9 was decent, but lower than that generated by the proposed method. Another brain extraction method, the SwissSkullStripper method [4] in the Slicer software, was also evaluated using the same data. This method yielded consistent results for both NCCT and CTA images with Dice values of ≈ 0.86. For the compared three methods, most segmentation errors exist at the base region of the brain, as shown in the first row in Fig. 4. The big difference of slice thickness between NCCT and CTA images, ranging from 0.5 mm to 5 mm in Table 1, might be the reason resulting in segmentation accuracy discrepancy for the Method-0.1 and Method-0.01 methods. In contrast to these two methods, the proposed approach is insensitive to slice thickness in the range of [0.5, 5] mm, performing consistently well for both modality images. We additionally performed another sensitivity analysis by changing the slice thickness of 5 CTA images randomly selected from 20 AIS patients used for quantitative evaluation, through resampling images in z axis with thickness of 5, 10, 15, 20, 25 mm. It was found that the average DSC was consistently >= 90% when slice thickness was <= 15 mm. However, the DSC decreased to <= 75% when slice thickness was >= 20 mm. 5. Discussion A brain extraction method for NCCT and CTA images is proposed in this study, which is based on a novel contour evolution technique with a shape propagation strategy. The proposed approach was quantitatively evaluated using 40 NCCT and CTA images of 20 stroke patients (one NCCT and one CTA scan for each patient), which were acquired from four different vendors (Toshiba, GE, Philips, and Siemens) with five NCCT and five CTA
Fig. 5. Four challenging cases successfully segmented by the proposed method.
scans from each vendor. Compared to manual segmentations as ground truth, the proposed brain extraction method was able to generate high Dice and Jaccard values for both NCCT and CTA images, outperforming the benchmark methods, two BET adapted methods: Method - 0.1 and Method - 0.01, and 3D Slicer. Although the scans used for the quantitative validation were obtained from different patients, different centers, and different scanners, the sample size was relatively small, which does not allow for exposure to many different image types, such as brain scans with, such as hemorrhage stroke, craniotomy, surgery, as well as scans for neonates and children. Further evaluation on 1736 NCCT and CTA images acquired from 1331 stroke patients was performed. These images were acquired from three multi-national and multi-centric clinical trials representing a diverse range of clinically relevant images with different levels of noise, artifacts, motion, and scanning parameters. The high degree of segmentation accuracy for both NCCT and CTA images in qualitative analysis indicates that the proposed algorithm can maintain accurate and specific results in a larger dataset. Fig. 5 shows a few cases, which were successfully processed by the proposed method, but not by other three methods. In particular, cut-off scans and the scans with irregular brain shape are quite challenging for the BET and Slicer methods. Additionally, after the qualitative evaluations, the processed 1331 NCCT images excluding 7 failure cases have been subsequently used to develop and validate an automated Alberta Stroke Program Early CT Score (ASPECTS) system [17], and the processed 405 CTA images have been used to develop a computer aided diagnosis system for characterizing intracranial thrombus predicting early recanalization with intravenous alteplase [27]. We do believe that through both quantitative and qualitative evaluations, the proposed brain extraction method is shown to be accurate and reliable for both NCCT and CTA scans. Another advantage of the proposed method is that it is a stand-alone algorithm without any dependency on any other third-party library or atlas, such as the BET toolbox. Thus, it can serve as a fundamental image processing step, therefore being integrated into any other algorithms seamlessly, such as spatial normalization [16]. More importantly, the CPU implementation of the developed algorithm is publicly accessible to benefit more for the community. The developed algorithm was the fastest automated method of brain extraction for NCCT images even with the CPU imple-
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
7
bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patentlicensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Acknowledgements
Fig. 6. Two examples of failure cases.
mentation. Our GPU implementation sped up calculation for NCCT images by at least two times, which is much faster than the other three methods. However, the computational time of the proposed method increased substantially for CTA images. Our CPU implementation for CTA images is the slowest compared to the other three methods. However, the GPU version is comparable to the BET method and faster than the Slicer. The Slicer software provided the slowest solution for both modality images as it relies heavily on image registration. Since our proposed method extracts brains in a slice by slice fashion, the greater slice number in CTA images dramatically slows down the calculation. However, the computation speed can be potentially accelerated. For instance, the computation of the data cost functions Ds,t (x ) in (8) and (9) and edge indicator g (x ) in (7) were calcuated using serial Matlab code and run on the CPU. This could be sped up by parallelizing Matlab code or converting it to C language. In addition, based on the segmentation results on NCCT images with greater slice thickness, it is reasonable to believe that the computation time for CTA images could be reduced greatly by enlarging the propagation interval between slices, while maintaining similar segmentation accuracy. However, finding the optimal interval and interpolating segmentations between slices are required to be further explored. We acknowledge that the proposed brain extraction method has some limitations. For example, poor image quality caused by excessive movement and incomplete acquisition of brain, as shown in Fig. 6 affects the proposed method.In addition, the images used in this study were without craniotomy and surgery, which would potentially increase the failure rate of the proposed method. Another limitation of the proposed method is that it cannot deal with images with large gantry tilt directly. For these cases, linear registration onto an atlas is required before processing to spatially recover brain orientation, which increases computational time accordingly.
Ethics approval was obtained from the local institutional review boards. Dr. Menon holds the Heart and Stroke/University of Calgary Professorship in Stroke Imaging. Dr. Goyal has a patent pending on systems of stroke diagnosis using multiphase computed tomographic angiography and a licensing agreement with GE healthcare for the same. Dr. Goyal also provided consulting services to Covidien for design and conduct of SWIFT PRIME trial and for teaching engagements. Dr. Hill and Dr. Goyal have a research grant with Medtronic (Covidien) to conduct clinical trials. The other authors report no conflicts. Appendix A. Pseudo code for the developed algorithm
Algorithm 1 Brain extraction. Require: Input image: I, slice number: S Ith = threshold (I ) for k = 1, A(k)=Cal cul ateArea(Ith (k )), while k < S [, Nmid ] = max(A(: )) UNmid = Local izeInitial Sl ice(Ith (Nmid )) while N = 0 do if N == Nmid then U prev = UNmid else U prev = UN+1 end if UN = Ext ract Brain2D(I (N ), U prev ) N ←N−1 end while while N < S do if N == Nmid then U prev = UNmid else U prev = UN−1 end if UN = Ext ract Brain2D(I (N ), U prev ) N ←N+1 end while U = Hol eF il l 3D(U )
6. Conclusion To summarize, an accurate, robust, and efficient method of brain extraction from NCCT and CTA images is proposed in this study. This method has demonstrated a high degree of accuracy when compared to manual segmentation. Qualitative evaluations of the method show that the specificity of the proposed algorithm can be maintained with a larger dataset. The proposed method has the potential of being used in different clinical settings, thus serving to improve current CT and CTA based image analysis in the field of neuroimaging.
Appendix B. Continuous Max-flow model and algorithm A flow maximization configuration is set up in a spatially continuous setting as shown in [25,26]. Let ps (x) and pt (x) be source and sink flows to and from pixel x to the source and sink terminals. The spatial flow q(x) exists around the neighborhood of pixel x. Then, the continuous max-flow model that maximizes the flow can be formulated as follows:
max P ( ps , pt , q ) =
ps ,pt ,q
ps dx
(10)
Conflicts of interest
subject to the following flow capacity and conservation constraints:
All authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’
| q ( x ) | ≤ g( x ) , p s ( x ) ≤ D s ( x ) , (divq − ps + pt )(x ) = 0.
pt (x ) ≤ Dt (x ),
8
M. Najm, H. Kuang and A. Federico et al. / Computer Methods and Programs in Biomedicine 176 (2019) 1–8
The continuous max-flow model (10) can be solved using the augmented Lagrangian method [31]:
ps dx + u, div q − ps + pt −
c div q − ps + pt 2 2
(11)
leading to an algorithm as following iterations:
q
k+1
2
c uk := arg max − div q + pkt − pks − , 2 c q∞ ≤α
pkt +1
pks +1 := arg max ps
u
k+1
2
c uk k+1 k , := arg max − p + div q − p − t s 2 c pt (x )≤ρ (l,x )
ps dx −
2
k
c ps − pkt +1 + div qk+1 + u , 2 c
= u − c div qk+1 − pks +1 + pkt +1 . k
Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cmpb.2019.04.030. References [1] Z. Akkus, P.M. Kostandy, K.A. Philbrick, B.J. Erickson, Extraction of brain tissue from CT head images using fully convolutional neural networks, in: Medical Imaging 2018: Image Processing, 10574, International Society for Optics and Photonics, 2018, p. 1057420. [2] F. Al Ajlan, A. Al Sultan, P. Minhas, Z. Assis, M. de Miquel, M. Millán, L. San Román, A. Tomassello, A. Demchuk, T. Jovin, et al., Posttreatment infarct volumes when compared with 24-hour and 90-day clinical outcomes: insights from the revascat randomized controlled trial, Am. J. Neuroradiol. (2018). [3] F.S. Al Ajlan, M. Goyal, A.M. Demchuk, P. Minhas, F. Sabiq, Z. Assis, R. Willinsky, W.J. Montanera, J.L. Rempel, A. Shuaib, et al., Intra-arterial therapy and post-treatment infarct volumes, Stroke 47 (3) (2016) 777–781. [4] S. Bauer, T. Fejes, M. Reyes, A skull-stripping filter for ITK, Insight J. 2012 (2013). [5] A. Boers, I. Zijlstra, C. Gathier, R. van den Berg, C.H. Slump, H. Marquering, C. Majoie, Automatic quantification of subarachnoid hemorrhage on noncontrast CT, Am. J. Neuroradiol. 35 (12) (2014) 2279–2286. [6] Y. Boykov, V. Kolmogorov, An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision, IEEE Trans. Pattern Anal. Mach. Intell. 26 (2001) 359–374. [7] Y. Boykov, V. Kolmogorov, D. Cremers, A. Delong, An integral solution to surface evolution PDES via geo-cuts, in: In ECCV, 2006, pp. 409–422. [8] Y. Boykov, O. Veksler, R. Zabih, Fast approximate energy minimization via graph cuts, IEEE Trans. Pattern Anal. Mach. Intell. 23 (11) (2001) 1222–1239. ¯ , M. Nikolova, Algorithms for finding global minimizers of [9] T.F. Chan, S. Esedoglu image segmentation and denoising models, SIAM J. Appl. Math. 66 (5) (2006) 1632–1648(electronic). [10] T.F. Chan, L.A. Vese, An efficient variational multiphase motion for the Mumford–Shah segmentation model, in: Proc. Conf Signals, Systems and Computers Record of the Thirty-Fourth Asilomar Conf, 1, 20 0 0, pp. 490–494. [11] A. Criminisi, T. Sharp, A. Blake, Geos: Geodesic Image Segmentation, in: Computer Vision–ECCV 2008, Springer, 2008, pp. 99–112. [12] A. Fedorov, R. Beichel, J. Kalpathy-Cramer, J. Finet, J.-C. Fillion-Robin, S. Pujol, C. Bauer, D. Jennings, F. Fennessy, M. Sonka, et al., 3d slicer as an image computing platform for the quantitative imaging network, Magn. Reson. Imaging 30 (9) (2012) 1323–1341. [13] A.M. Frölich, S.L. Wolff, M.N. Psychogios, E. Klotz, R. Schramm, K. Wasser, M. Knauth, P. Schramm, Time-resolved assessment of collateral flow using 4d CT angiography in large-vessel occlusion stroke, Eur. Radiol. 24 (2) (2014) 390–396. [14] M. Goyal, A.M. Demchuk, B.K. Menon, M. Eesa, J.L. Rempel, J. Thornton, D. Roy, T.G. Jovin, R.A. Willinsky, B.L. Sapkota, et al., Randomized assessment of rapid endovascular treatment of ischemic stroke, N. Engl. J. Med. 372 (11) (2015) 1019–1030. [15] P. Jaccard, The distribution of the flora in the alpine zone., New Phytol. 11 (2) (1912) 37–50. [16] J.S. Kim, H. Cho, J.Y. Choi, S.H. Lee, Y.H. Ryu, C.H. Lyoo, M.S. Lee, Feasibility of computed tomography-guided methods for spatial normalization of dopamine transporter positron emission tomography image, PLoS One 10 (7) (2015) e0132585. [17] H. Kuang, M. Najm, D. Chakraborty, N. Maraj, S. Sohn, M. Goyal, M. Hill, A. Demchuk, B. Menon, W. Qiu, Automated aspects on noncontrast ct scans in patients with acute ischemic stroke using machine learning, Am. J. Neuroradiol. 40 (1) (2019) 33–38.
[18] Z. Ma, J.M.R. Tavares, R.N. Jorge, T. Mascarenhas, A review of algorithms for medical image segmentation and their applications to the female pelvic cavity, Comput. Methods Biomech. Biomed. Eng. 13 (2) (2010) 235–246. [19] J.G. Mandell, J.W. Langelaan, A.G. Webb, S.J. Schiff, Volumetric brain analysis in neurosurgery: part 1. particle filter segmentation of brain and cerebrospinal fluid growth dynamics from mri and CT images, J. Neurosurg. 15 (2) (2015) 113–124. [20] B.K. Menon, C.D. d Esterre, E.M. Qazi, M. Almekhlafi, L. Hahn, A.M. Demchuk, M. Goyal, Multiphase ct angiography: a new tool for the imaging triage of patients with acute ischemic stroke, Radiology 275 (2) (2015) 510–520. [21] B.K. Menon, M. Najm, F. Al-Ajlan, J.P. Alcantara, D. Dowlatshahi, A. Calleja, S.-I. Sohn, S.H. Ahn, A. Poppe, R. Mikulik, et al., IV TPA recanalization rates by site of occlusion and time after TPA bolus-main results of the interrsect multinational multicenter prospective cohort study, 2017. [22] J. Muschelli, N.L. Ullman, W.A. Mould, P. Vespa, D.F. Hanley, C.M. Crainiceanu, Validated automatic brain extraction of head CT images, Neuroimage 114 (2015) 379–385. [23] D.L. Pham, C. Xu, J.L. Prince, Current methods in medical image segmentation, Annu. Rev. Biomed. Eng. 2 (1) (20 0 0) 315–337. [24] T. Pock, A. Chambolle, H. Bischof, D. Cremers, A convex relaxation approach for computing minimal partitions, in: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2009. Miami, Florida [25] W. Qiu, Y. Chen, J. Kishimoto, S. de Ribaupierre, B. Chiu, A. Fenster, B. Menon, J. Yuan, Longitudinal analysis of pre-term neonatal cerebral ventricles from 3D ultrasound images using spatial-temporal deformable registration, IEEE Trans. Med. Imaging (2017). [26] W. Qiu, Y. Chen, J. Kishimoto, S.d. Ribaupierre, B. Chiu, A. Fenster, J. Yuan, Automatic segmentation approach to extracting neonatal cerebral ventricles from 3D ultrasound images, Med. Image Anal. 35 (2017) 181–191. [27] W. Qiu, H. Kuang, J. Nair, Z. Assis, M. Najm, C. McDougall, B. McDougall, K. Chung, A. Wilson, M. Goyal, et al., Radiomics-based intracranial thrombus features on CTand CTA predict recanalization with intravenous alteplase in patients with acute ischemic stroke, Am. J. Neuroradiol. 40 (2019) 39–44. [28] W. Qiu, J. Yuan, M. Rajchl, J. Kishimoto, Y. Chen, S. de Ribaupierre, B. Chiu, A. Fenster, 3D MR Ventricle segmentation in pre-term infants with post-hemorrhagic ventricle dilatation (phvd) using multi-phase geodesic level-sets, Neuroimage 118 (2015) 13–25. [29] W. Qiu, J. Yuan, E. Ukwatta, A. Fenster, Rotationally resliced 3D prostate trus segmentation using convex optimization with shape priors, Med. Phys. 42 (2) (2015) 877–891. [30] W. Qiu, J. Yuan, E. Ukwatta, Y. Sun, M. Rajchl, A. Fenster, Prostate segmentation: an efficient convex optimization approach with axial symmetry using 3D TRUS and MR images, IEEE Trans. Med. Imaging 33 (4) (2014) 1–14. [31] R.T. Rockafellar, Augmented Lagrangians and applications of the proximal point algorithm in convex programming, Math.Oper.Res. 1 (2) (1976) 97–116. [32] C. Rorden, L. Bonilha, J. Fridriksson, B. Bender, H.-O. Karnath, Age-specific CT and MRI templates for spatial normalization, Neuroimage 61 (4) (2012) 957–965. [33] C. Rorden, M. Brett, Stereotaxic display of brain lesions, Behav. Neurol. 12 (4) (20 0 0) 191–20 0. [34] D. Rueckert, L.I. Sonoda, C. Hayes, D.L. Hill, M.O. Leach, D.J. Hawkes, Nonrigid registration using free-form deformations: application to breast MR images, Med. Imaging, IEEE Trans. 18 (8) (1999) 712–721. [35] P.W. Schaefer, L. Souza, S. Kamalian, J.A. Hirsch, A.J. Yoo, S. Kamalian, R.G. Gonzalez, M.H. Lev, Limited reliability of computed tomographic perfusion acute infarct volume measurements compared with diffusion-weighted imaging in anterior circulation stroke, Stroke 46 (2) (2015) 419–424. [36] D. Shen, G. Wu, H.-I. Suk, Deep learning in medical image analysis, Annu. Rev. Biomed. Eng. 19 (2017) 221–248. [37] S.M. Smith, Fast robust automated brain extraction, Hum. Brain Mapp. 17 (3) (2002) 143–155. [38] J. Solomon, V. Raymont, A. Braun, J.A. Butman, J. Grafman, User-friendly software for the analysis of brain lesions (able), Comput. Methods Programs Biomed. 86 (3) (2007) 245–254. [39] P.J. Toivanen, New geodesic distance transforms for gray-scale images, Pattern Recognit. Lett. 17 (5) (1996) 437–450. [40] Z. Wang, K.K. Bhatia, B. Glocker, A. Marvao, T. Dawes, K. Misawa, K. Mori, D. Rueckert, Geodesic Patch-based Segmentation, in: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2014, Springer, 2014, pp. 666–673. [41] J. Yuan, E. Bae, X. Tai, A study on continuous max-flow and min-cut approaches, CVPR, 2010. [42] J. Yuan, E. Bae, X.-C. Tai, Y. Boykov, A continuous max-flow approach to Potts model, ECCV, 2010. [43] P.A. Yushkevich, J. Piven, H. Cody Hazlett, R. Gimpel Smith, S. Ho, J.C. Gee, G. Gerig, User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability, Neuroimage 31 (3) (2006) 1116–1128. [44] Y. Zhuge, J.K. Udupa, Intensity standardization simplifies brain MR image segmentation, Comput. Vision Image Underst. 113 (10) (2009) 1095–1103.