Robust Segmentation of the Full Cerebral Vasculature in 4D CT of Suspected Stroke Patients

A robust method is presented for the segmentation of the full cerebral vasculature in 4-dimensional (4D) computed tomography (CT). The method consists of candidate vessel selection, feature extraction, random forest classification and postprocessing. Image features include among others the weighted temporal variance image and parameters, including entropy, of an intensity histogram in a local region at different scales. These histogram parameters revealed to be a strong feature in the detection of vessels regardless of shape and size. The method was trained and tested on a large database of 264 patients with suspicion of acute ischemia who underwent 4D CT in our hospital in the period January 2014 to December 2015. Five subvolumes representing different regions of the cerebral vasculature were annotated in each image in the training set by medical assistants. The evaluation was done on 242 patients. A total of 16 (<8%) patients showed severe under or over segmentation and were reported as failures. One out of five subvolumes was randomly annotated in 159 patients and was used for quantitative evaluation. Quantitative evaluation showed a Dice coefficient of 0.91 ± 0.07 and a modified Hausdorff distance of 0.23 ± 0.22 mm. Therefore, robust vessel segmentation in 4D CT is feasible with good accuracy.

Another problem encountered in segmentation of the cerebral vasculature is the area of the skull base. In this region the ICA, the jugular bulb and soft bone tissue are close to each other with similar intensities between these structures. Vessel pathologies (e.g. calcifications and aneurysms) often appear in the skull region compared to other parts in the brain 5 . Region based approaches that find the vessel boundaries directly without first centerline tracking 6 , may have difficulties in this area as well.
CT is the primary imaging modality in an emergency setting because of its short acquisition time and because of the fewer contraindications compared to magnetic resonance (MR). Furthermore, CT has a high diagnostic value in the acute phase of stroke, due to its high sensitivity to hemorrhages and vascular pathology [7][8][9] . Modern CT scanners are capable of acquiring contrast-enhanced 4D CT images with high spatial and temporal resolution at large coverage, enabling visualization of the contrast dynamics at submillimeter resolution of the full brain. A 4D CT acquisition is less critical concerning timing of the peak arterial enhancement compared to conventional 3D CTA.
In this work we consider vessel segmentation in 4D CT. Only one related vessel segmentation method was found in 4D CT 10 , which was based on fixed thresholding of a global intensity histogram. Another approach towards vessel segmentation is a method based on pattern recognition. This way local differences of vessel shape and intensity can be learned, which makes it less prone to discontinuities of the vessels e.g. due to the presence of pathologies or foreign objects. In addition, a pattern recognition method requires no explicit seed points for initialization, as opposed to tracking based methods.
Pattern recognition based approaches have hardly been used for vessel segmentation. One reason may be the difficulty of establishing a reliable reference standard. A precise manual delineation of the tortuous small vessels is difficult and time-consuming. Two related works have been found based on pattern recognition: on vessel segmentation in MR angiography (MRA) 11 and on vessel segmentation in lung CT 12 . The latter achieved a top rank in the VESSEL12 challenge 13 demonstrating the potential of such an approach.
In this work, we present a pattern recognition approach for segmentation of the complete cerebral vasculature in 4D CT. The proposed method emphasizes robustness and should be capable of handling patient data typically seen in everyday clinical practice. The method also emphasizes completeness over accuracy in defining the vessel boundaries by including the segmentation of the smaller vessels more distally in the vascular tree. Obtaining accurate boundary segmentations is easier starting from the complete cerebral vasculature than the other way around. Furthermore, completeness of the segmentation is more important, because the desired accuracy of the lumen boundary is dependent on the subsequent application.

Results
Observer Variability. Table 2 shows the quantitative evaluation of the inter-observer variability, intra-observer variability and the performance of the method in comparison to each observer for a subset of five 4D CT images. DSC is Dice similarity coefficient 14 , HD is Hausdorff distance 15 , MHD is modified HD 16 , 95% HD is 95 th percentile HD, CMD is contour mean distance 17 and AVD is absolute volume difference. The DSC for the inter-observer variability was 0.75 ± 0.06. The method scored higher, although this difference is only significant in observer 1. A boxplot graphical overview is shown in Fig. 3 where the center line indicates the median value, Qualitative Evaluation. Visual inspection of all segmentations (n = 242) revealed a total of 16 failures (<8%). In 8 patients this failure was the result of severe motion artifacts in multiple volume acquisitions unaccounted for by the registration of the 4D CT or due to beam hardening artifacts during the acquisition. This resulted in a high intensity variation over time leading to erroneous segmentation. An example of this is shown in Fig. 4. These 16 failures were not considered for quantitative evaluation. Smaller movement artifacts outside the cranial cavity (e.g. eye or jaw movement) posed no problem for the method. The vast majority of the vessel segmentations were evaluated as excellent (Fig. 5). Presence of foreign objects or severe trauma proved no difficulty for the method nor resulted in over segmentation. Figure 6 shows a patient with foreign objects implanted as a result of surgical intervention. The skull base region is known to be the most difficult region for vessel segmentation. This region holds the ICA and a network of small veins intertwined surrounding the ICA, making distinguishment difficult. This region of interest in cerebral vessel segmentation was segmented in detail (Fig. 7).
Quantitative Evaluation. The segmentations were compared with the reference standard provided by the manual annotations for quantitative evaluation. The same metrics as described in Section Observer Variability have been used. The results of the quantitative evaluation for 143 4D CT images are summarized in Table 3. The overall results and results per subvolume are presented. Subvolume 1 included the ICA, the CoW and the middle  cerebral artery (MCA), subvolume 2 was defined in the frontal area including the anterior cerebral artery (ACA), subvolume 3 and 4 were defined left and right in the brain including the M2 and M3 vessel segments, and finally subvolume 5 at the parieto-occipital included the transverse sinus and the superior sagittal sinus.
A boxplot graphical summary of the DSC score and the MHD is shown in Figs 8 and 9 respectively. Overall, the sensitivity, specificity, and accuracy of cerebral vasculature segmentation reported for all 143 patients were 0.938, 0.997 and 0.995 respectively.

Discussion
In this work we have presented a method for full cerebral vessel segmentation in 4D CT using a pattern recognition framework with minimal pre-and postprocessing. The results show that the method correctly handles foreign objects without losing segmentation accuracy. Furthermore, our method is able to segment vessels of different sizes at different locations.
The presented method requires no tracking or active contours as these may be unable to handle the discontinuities in the vessels due to steno-occlusive disease, vessel wall pathology, calcification, or artifacts resulting from foreign objects. In addition, natural intensity differences throughout the vascular system in the brain may pose problems for these approaches. Intensity variations limit the application of global parameters and may require extra user input if the method depends on a manual chosen starting point.
The novelty of our method lies in two distinct features. First, the WTV image, which represents the intensity variance in temporal direction. This feature reduces the dimensionality of the image data and is highly sensitive to contrast variation over time. It is therefore very suitable for the visualization and detection of vessels in 4D CT. In addition, in contrast to approximated measures as the area under the curve, the WTV is parameter free and independent of the temporal distance between the volumetric acquisitions of the 4D CT protocol. The second feature is the use of the parameters (mean, standard deviation, mode and entropy) of a local intensity histogram to train the classifier. By analyzing the contrast variation in a histogram in a subvolume, each vessel, regardless of size, can be detected by the profile represented in the parameters of the histogram. The local histogram features are preferable to global thresholds as smaller vessel segments may be under segmented or missed as they fall below the global threshold dominated by the large vessels. By comparing vessel intensity to its immediate Inter-observer (n = 5) Intra-observer (n = 5)  Table 2. Quantitative evaluation of inter-observer variability, intra-observer variability and method versus each observer independently (n = 5), reported as mean ± standard deviation. p-values were computed with a paired sample t-test between the inter-observer variability and proposed method for each evaluation measure. † Indicates a significant difference (p-value < 0.05, i.e. less than 5% chance that the samples have identical averages). surrounding this can be avoided, as even the smallest vessels containing contrast agent have a higher intensity than the environment. The method used pre-and postprocessing steps. The candidate vessel selection preprocessing step prior to the classification removes voxels which do not show any contrast variation over time. Removing these voxels reduces the computation time without losing possible vessel voxels.
One postprocessing step removed small objects from the segmentation result. The component size constraint of 25 voxels translates to 2.3 μL. These small objects mainly occurred outside the cranial cavity around the facial area of the patient, due to noise, or eye or jaw movement during the acquisition.
Another postprocessing step was applied to compensate for under segmentations in large vessels. Under segmentation can occur in the main supplying arteries and main draining veins (e.g. the ICA and the transversus sinus). A reason for this under segmentation is a high blood flow causing high contrast variations over time. These high flows may not be captured by the 19 patients in the training set, resulting in holes in the segmentation. Although the classifier was trained with only 19 partially annotated patients, the method already performed well with this limited representation of the full population. Extra training sets can be added from patients with various arterial and venous flow speeds to capture a wider range in contrast variation to further increase the robustness of the method.
The inter-observer variability demonstrates the difficulty in the definition of the vessel border. The low inter-observer variability DSC of 0.75 ± 0.06 can be explained by the small size of the distal intracranial vessels. A vessel in this region can amount to merely two or three voxels in diameter. With a low component size the DSC is then easily affected. A 95% HD of 7.57 ± 3.07 mm is reported, which translates to a few voxels on both sides of the vessel, indicating that there is discrepancy on the definition of the exact border of the vessel. An exact definition of the vessel boundary is difficult because of the partial volume voxels. Also, what appears as vessel boundaries  in the image is influenced by other factors including the injection rate of the contrast agent, the CT acquisition protocol, cardiovascular condition of the patient, and window width and level settings.
The proposed method has proven to be robust by evaluating on a large database comprised of 242 patients with foreign objects and image artifacts found in everyday clinical practice. The failures found were a result of insufficient intra-patient registration or as a result of acquisition artifacts resulting in loss of diagnostic value. Image artifacts caused by surgical intervention (e.g. clip, coil or drain) posed no problem for the method and did not result in under or over segmentation. Only severe movement between temporal acquisitions that has not been mitigated by the intra-patient registration showed to be a problem for the method. Smaller movement outside the cranial cavity only proved to be a problem in a small set of patients. The intra-patient registration registers bone tissue in each time point to the first time point. Soft tissues as the eyes and nose are not registered and can cause motion artifacts after the registration. Although our method handles minor soft tissue movements without any problems, severe movement of the jaw and nose may result in over segmentation.  Accurate and complete vessel segmentation in different brain regions is shown by the results in Table 3. The different subvolumes contained vessels of varying sizes, both arterial and venous. The DSC score showed similar results regardless of different subvolume evaluation region, which further supports the claim of independence towards vessel size or location. The MHD showed a distance similar to or smaller than the voxel spacing of 0.43 mm for all subvolumes combined as well for individual subvolumes, indicating a segmentation accurate to or below a single voxel. The method was not explicitly evaluated for topological consistencies, instead, we used similar evaluation measures as were used in some of the larger vessel segmentation challenges [18][19][20] .   Table 3. Quantitative evaluation of the segmentation results overall (n = 143) and per subvolume, reported as mean ± standard deviation. DSC is Dice similarity coefficient, HD is Hausdorff distance, MHD is modified HD, 95% HD is 95 th percentile HD, CMD is contour mean distance and AVD is absolute volume difference.  Temporal information, in combination with contrast agent, is important for vessel segmentation as is reflected by the WTV feature. The added value of 4D CT with improved evaluation of intracranial hemodynamics comes at a cost, as a 4D CT protocol is associated with a higher radiation dose. Although 4D CT imaging is not common practice, applications of 4D CT are expanding 21 . We expect 4D CT to become a single acquisition for stroke workup as it contains both noncontrast CT and CTA information [22][23][24] . These modalities might be reconstructed from a 4D CT acquisition, resulting in a reduction of acquisitions and radiation dose. In addition, studies suggest that 4D CT can be acquired at half the dose of standard clinical protocol 25 , further reducing the radiation dose for the patient.
Although the images for this study were acquired with a high-end 320-row detector scanner, the voxel based classification method is not constrained to a specific image size nor images from scanners with a different number of detectors.
To the best of our knowledge this is the largest quantitative and qualitative evaluation reported for vessel segmentation in general. Robust cerebral vessel segmentation in 4D CT has not been performed on a similar scale.
In conclusion, we presented a method using candidate selection combined with a RF classifier trained with localized histogram features to segment a full cerebral vasculature. The method provides accurate and robust vessel segmentation on large amounts of everyday clinical data and forms the cornerstone of many subsequent applications.

Methods
We propose a three step approach: first, coarse candidate vessel detection, then application of a random forest (RF) classifier 26 using local and global image features and finally postprocessing to ensure full and smooth connectivity. The RF classifier is a supervised learning method that enables the learning of the background and the properties of vessels with different sizes and intensities. In addition, RF enables the learning of foreign objects or image artifacts which may be present in the data through a number of decision trees. The method was developed using the Insight Segmentation and Registration (ITK) 17 library.
Candidate Selection. The candidate vessel voxels were selected from the temporal variance image weighted to the exposures of the individual time points to maximize the signal to noise ratio. This WTV image is sensitive  Weighted Temporal Variance. Similar as to the WTA, the temporal variance can be calculated with a weighted quadratic difference to the WTA. The WTV holds the largest contrast between tissue containing contrast agent and background tissue (Fig. 10) and is given by: Local Histogram Features. To capture local image properties such as varying vessel intensity values, an intensity histogram is computed within a 5 × 5 × 5 and 9 × 9 × 9 neighborhood in the WTV image. The following histogram parameters were included: the mean, standard deviation, mode and the entropy given by , with p i the probability of intensity i, and ε = 0.0001 to prevent zero-logarithm calculation.
Distance to Intracranial Cavity. An intracranial cavity mask is calculated and used as distance feature. The cranial cavity contains all soft tissues and cerebrospinal fluid, including the meninges, cerebrum and ventricles, cerebellum and brain stem. The cavity mask was created by a multi-atlas registration with label fusion followed by a geodesic active contour levelset refinement of the segmentation 28 . The internal carotid arteries and jugular veins are located outside of the intracranial cavity mask. These will be missed if the cavity mask is added as a binary feature, therefore the Euclidean distance to the border of the intracranial cavity mask was computed. Distances from inside the mask to the border were denoted as negative distances, distances from outside the mask to the border were denoted as positive distances. Hessian. The eigenvalues of the Hessian system are calculated on the WTV at four different scales on an evenly distributed range from 0.5 to 2.0 mm. T 0 . The 4D CT acquisition contains 19 volumetric acquisitions. The first time point volume was selected to represent the tissue with minimal contrast. This volume was used in the feature set as T 0 .
Classifier. An RF classifier is applied to the candidate vessel voxels. The RF classifier was trained with 100 trees and a maximum tree depth of 30. Postprocessing. After the RF classification two postprocessing steps were applied. First, a connected components (26-neighborhood) analysis was carried out, discarding components smaller than 25 voxels to remove background noise. Second, morphological hole filling with a 3 × 3 × 3 kernel size was carried out for 10 iterations to fill the center of some larger vessels such as the ICA in the event of under segmentation by the classifier.

Materials
Patient Data. This retrospective study was approved by our local institutional review board and the requirement for informed consent was waived. In total, 264 consecutive patients (age: 68y ± 15.0, male: 149 (55.6%)) with suspicion on acute ischemic stroke who were admitted to the emergency department of our academic hospital, between January 1st 2014 and December 31st 2015, were retrospectively collected for this study. Exclusion criteria were severe motion artifacts resulting in loss of diagnostic value (n = 3). This was determined from the radiology report provided by the radiologist on call at time of admission. Image acquisition was done on a 320row Toshiba Aquilion ONE CT scanner (Toshiba Medical Systems Corporation, Japan). The protocol consisted of 19 whole brain volume acquisitions starting with a single high dose acquisition at 200 mAs 5 seconds after contrast injection, followed by 13 scans every 2 seconds at 100 mAs, followed by 5 scans every 5 seconds at 75 mAs. Each volumetric acquisition was made at 80 kV at 0.5 seconds rotation time with 16 cm z-coverage. A dose of 80 ml Iomeron is injected in the antecubital vein at the start of the first acquisition. Image reconstruction was done with a smooth reconstruction kernel (FC41), resulting in image sizes of 512 × 512 × 320 voxels with voxel sizes of 0.43 × 0.43 × 0.5 mm. All temporal volumes were rigidly registered to the first time point to resolve any patient movement during the acquisition.

Volumetric Cluster Annotation and Segmentation Tool (VCAST).
Annotating vessels in 3D is a laborious and difficult task because of their small sizes and the presence of partial volume voxels, and because vessels meander through the different orthogonal planes. To facilitate the annotation process, a dedicated 3D annotation tool called volumetric cluster annotation tool was developed 29 . VCAST supports interactive region growing and super voxel grids overlaid on the input image to let the user assign or reject vessel labels to clusters of voxels at multiple scales. Additional functionalities include the ability to define subregions confining the region growing algorithm, a 3D rendering of the topological information of the annotations and a colour map defining connected regions with the same colour. Annotations were done on the WTA because of its maximal signal to noise ratio and on the WTV because of its high sensitivity to contrast voxels. Observers were able to switch between these reconstructions for optimal interaction. Reference Standard. Annotating the full cerebral vasculature is too labour intensive and not required, because many vessel segments have similar appearances and intensities. Therefore, annotations were carried out in five subvolumes representing different parts of the vasculature, as shown in Fig. 11.
Observers were asked to annotate high confidence vessel voxels in these subvolumes, thereby focusing more on the completeness of the vasculature and less on the accuracy of the vessel boundaries at a single voxel level. In the training dataset (n = 19) all five subvolumes were annotated in each patient. In the testing dataset (n = 159) one of the five subvolumes was randomly selected per patient for annotation. These annotations were carried out by two medical assistants, supervised by an experienced neuroradiologist (FJAM).

Evaluation Measures.
The segmentations were compared with the reference standard provided by the annotations using the following measures: sensitivity, specificity, accuracy, DSC, HD, MHD, 95% HD, CMD, and AVD. The MHD is defined as the maximum value of the mean distance calculated between two objects. The CMD is defined as the mean distance between the segmentations along their boundary voxels. The sensitivity, specificity, and accuracy were calculated within the annotated subvolume on a voxel basis. The mean and standard deviation was taken over all patients for the other measurements. These evaluation measures are similar to the measures used in some vessel segmentation challenges [18][19][20] .
The goal of our method is to robustly segment a complete cerebral vasculature. To that end, the focus of the evaluation will not be around the border of vessel lumen. Therefore, the lumen boundaries are excluded in the final evaluation by dilating the manual annotations by 1 voxel (3 × 3 × 3 kernel). Applying this boundary mask will only boost the specificity of the algorithm in regions where exact definitions of the vessels border is difficult, but it will not affect the sensitivity.

Experiments
A schematic overview of the acquired imaging data and the experiments described below is shown in Fig. 12. Training and Validation. The annotation of the training set provided positive samples from 5 subvolumes within 19 patients. A random subset of 10% of the annotated vessel voxels was used as positive samples to train the classifier. A dataset of the same size as the positive data set was randomly sampled from all nonannotated voxels within the subvolumes. These voxels served as background. The RF classifier was optimized using leave-one-out crossvalidation within the training set. The accuracy of the receiver operator characteristics (ROC) analysis and visual inspection of the output was used to determine the optimal operating point. Observer Variability. A subset of 5 patients randomly selected from the training set were annotated by both trained observers. The same subset was annotated by one observer once again after a waiting period of two weeks. All evaluation measures described in Section Evaluation Measures were reported for the inter-and intra-observer variability. Quantitative evaluation of observer variability was performed on the full annotation without excluding boundary voxels. Because the annotation sets consisted of five nonconnected subvolumes the distance evaluation measures were first calculated per subvolume, then combined between subvolumes according to their measurement. That is, the mean for CMD and the maximum for HD, MHD, and for 95% HD. Finally, the mean and standard deviation per evaluation measure were taken over all patients and reported.
Qualitative Evaluation. The robustness was verified by imposing minimal exclusion criteria on the test data set, thus resembling clinical practice. All vessel segmentations were initially inspected to assess failures. Failure was reported in cases of severe under or over segmentation. In addition, the entire branch network of the cerebral structure was investigated for completeness.
Quantitative Evaluation. The completeness of vasculature was assessed by quantitatively evaluating different parts of the cerebral structure in a large number of patients. The quantitative evaluation was performed on one of five randomly annotated subvolumes. These subvolumes were defined in a similar region as with the training data, described in Section Results under Quantitative Evaluation.