Vasculature analysis based on CT images is indispensable in liver surgical planning and diagnosis of liver disease. We propose a novel framework and an efficient workflow with minimal user interaction to analyze liver vasculature in multi-phase CT images. To ensure segmentation quality and efficiency, a set of semi-automatic algorithms are applied to initially segment different vascular structures at different phases. A fully automatic vessel separation procedure runs parallel to separately connected hepatic and portal veins. In addition, an interactive editing method is integrated into the framework to refine the segmentation of each individual structure. Quantitative evaluations of segmented vessels conducted for 60 test data sets are demonstrated.
The processing pipelines for different selected phases also vary. Figure 1 illustrates the processing pipelines for HV/PV, HA, and BD phases.
HV/PV segmentation: Adequate delineation of the HV and PV is important for liver surgical planning. During acquisition of multiphase CT images, both of these venous vessels are enhanced and acquired, either in a common phase or in two individual phases. The presence of large lesions and inadequate enhancement of venous vessels, which often occur during imperfect image acquisition in oncological patients, hinder the segmentation task. To overcome these difficulties, we employ a semi-automatic segmentation framework requiring minimal user interactions for quick and robust segmentation of hepatic vein and portal vein presented in common or individual phases. First, in the preprocessing stage, we segment the liver on the selected venous phase and analyze the histogram of intensities to preclude regions with extremely low intensities.
Hessian-based filters have been widely employed to enhance tube-like structures in 3D images. Various vesselness filters published in previous works [4] [5]. Because the contrast level of venous vessels in different input images are quite heterogeneous, we choose the multi-scale, Hessian-based vesselness filter introduced by Frangi et al. to enhance vessels in CT volumes within the liver mask [4]. Major benefits of Frangi’s vesselness filter include incorporation of information from all three Eigen values and independence from their absolute values, which are associated with absolute contrast level of vessels; its output is scaled in range of 0 and 1. Assuming the eigenvalues of Hessian matrix are sorted in order: |λ1| ≤ |λ2| ≤ |λ3|, the Frangi’s vesselness filter is defined as: where | |, | | and . The ratio is designed to differentiate vessels from sheet-li|ke| structure|s, w|hereas is used to distinguish vessels from blob-like structures. The term aims to suppress noise structures. The scale parameter σ indicates the size of Gaussian kernel used for calculating the Hessian matrix. The parameters used in this filter, are set to 0.5, 0.5, and 10. Considering the radius range of HV and PV, we choose three optimized scales based on our experimental tests for σ: 1.5 mm, 2.25 mm, and 3 mm. These can capture vessels with thin, medium, and thick radii (see Fig. 2(a)). The ultimate vesselness response is obtained by extracting the maximum across all scales. Normally, calculating the Hessian matrix is quite expensive when increasing the Gaussian kernel size. To accelerate the vesselness filter, we adopt an inverse multi-scale strategy which keeps the kernel size constant at 1.5mm but downsamples input volumes to larger voxel spacing.
Based on vesselness response, we initially segment the HV and PV using an automatic region-growing algorithm. With histogram analysis, all voxels with a vesselness value between 90th and 99th percentiles are taken as seed points. The lower and upper thresholds of the region growing are chosen as the 75th percentile and maximal response, respectively. If the HV and PV possess a high contrast compared to their surroundings, initial segmentation normally yields satisfactory results (see Fig. 2(b)). Nevertheless, in cases for which HV and PV exhibit extremely low contrast or where large hypodense lesion areas with internal hyperdense structures are present, many irrelevant structures will be enhanced by the vesselness filter and captured in the region-growing step. Therefore, we first transform the segmented vessel tree structure into a graph representation [6] and validate the graph tree using two criteria: 1) The volume of each independent graph tree should exceed a lower limit Tmin.
2) The edge length of branches should be longer than a minimum threshold Emin .
These two graph filters can rule out most irrelevant tree structures and prune disturbing branches. Other issues lie in the areas of inferior vena cava and entrance of PV, which are normally not enclosed by the segmented liver mask. Moreover, the inferior vena cava can be hypodense or hyperdense depending on the contrast agent density. We propose a robust interactive solution to recover the segmentation in both regions. Two markers are placed in the regions where main branches of HV and PV join (see Fig. 2(c)). Then, a fast marching algorithm is applied in the predefined neighborhoods of these two markers [7]. The efficient computation of the algorithm permits instant display of segmented regions when marker positions are adjusted until optimal positions are reached.
HV/PV separation: The HV and PV will be simultaneously segmented when they present in a common phase. For liver surgical planning, it is demanding to separate and analyze them individually. The separation process is triggered when placing HV and PV markers are placed. To meet computation time constraints, we employ an interactive watershed transform algorithm that takes the placed HV and PV markers as seed points of two different classes [8]. The cost images of watershed transforms are original intensity images smoothed with a diffusion filter. The entire volume is classified into two classes associated with HV and PV markers, respectively. Consequently, the segmented HV and PV can be separated into two classes when we overlay the class labels over the segmentation result (see Fig. 2(d)).
HA and BD segmentation: Compared to HV and PV, segmentations of HA or BD are relatively less challenging because HA and BD normally show very high contrast against other structures in HA and BD phases. However, one issue that we observed in some of HA phases is that other vascular structures, such as HV and PV, are slightly enhanced and that contrast enhancement of HA is heterogeneous. Hence, we use a single-scale vesselness filter again to enhance HA in the preprocessing step. The scale parameter of the kernel size is set to 1.0 mm. Afterwards, a semi-automatic region-growing scheme is conducted which requests a single marker placed in the HA or BD volume. The cost images used in region growing are vesselness response image and intensity image for HA and BD, respectively.
To allow modification of the initial segmentation and separation results by the user, we provide a set of interactive editing methods. The methods are based on a graph structure (consisting of a set of edges and nodes) calculated from the existing segmentation mask. As an input for each editing step, only a single click by the user is required. The following correction modes are available:
Deleting branches
Relabeling branches (required for incorrect vessel separation) Adding branches is based primarily on Dijkstra's algorithm, which uses the underlying CT image as a cost image. Given a selected click point, the algorithm calculates the cheapest path from this point to the segmentation mask. This path builds the centerline for the new branch. Subsequently, vessel radii along the path are approximated using a ray-casting algorithm that evaluates grey values around the centerline. Finally, the new branch is added to the existing graph structure.
Whereas adding branches is only applicable to 2D slices, the other two editing modes can also be applied by interacting with the 3D vessel models in the 3D viewer. Analogous to the adding mode, only a single click point by the user is required. For deleting and re-labeling of branches, the user controls the influence of the correction by selecting either a local or a global influence. For the local influence mode, only the edge selected by the user is modified. For the global influence mode, all parent edges are modified, i.e., removed or relabeled. These three correction modes provide the user with the necessary tools to refine the results of vascular segmentation and separation steps describe in section 2.1. 3
To evaluate the performance of the presented segmentation techniques, a test data set enclosing 60 multiphase hepatic CT scans acquired from 30 liver donors and 30 oncological patients was collected. HV/PV phases were taken for all scans. Depending on diagnostic circumstances, HA and BD phases were not always present in the original scans. Segmenting HA and BD is somewhat less challenging compared to HV/PV segmentation, and the results were normally satisfactory after a few manual correction steps. Therefore, in this work, we quantitatively evaluated the quality of HV/PV segmentation. The vasculatures of HV/PV were manually segmented by experienced radiologists and were used as references for the quantitative evaluation. In addition, the markers used to refine the segmentation in the inferior vena cava and entrance of the PV were placed manually and saved prior to automatic testing. The image resolution of HV/PV phases in test images ranged from 512×512×189 to 512×512×310. The calculation of multi-scale vesselness filter required 15 to 40 seconds, depending on the number of vasculatures in the scans. The initial segmentation and subsequent refinement of HV/PV lasted a maximum of 3 seconds using a 3.07 GHz Intel CPU.
The Dice coefficient (DC) is computed to measure the overlap between segmented HV/PV and the ground truth. It reflects both misalignment of branches and errors of the radii. We introduce an additional metric to further review the alignment between segmented and reference vessel trees. The percentage of the skeleton of segmented vessels covered by the reference vessel mask is calculated, and the same percentage is computed inversely for the skeleton of reference vessels. The average of two percentages depicts the branching consistency of both vessel trees. Table 1 lists the statistical results of DC and branching alignment (BA) for all test cases. The branching alignment BA1 is larger than BA2 on average because the reference vessel tree tends to delineate the complete vessel structures, including very thin vasculatures, whose radii are less than 1mm and are of less interest for liver surgery planning.
The multi-scale vesselness filter is sensitive in detecting vascular structures in different contrast levels. Increasing the number of scales will deliver better enhancement response for the vessels with varying radii. Considering the balance between performance and computational expense, we choose three scales in this work, which requires acceptable computation time. Further experiments are required to investigate how the performance relates to the scale numbers. In practice, the disturbing structures resembling vasculatures locally will be enhanced as well, especially in patients with oncological lesions. The graph filters and manual editing tools in the post-processing step are capable of ruling out these false positives. The optimal parameter settings for the graph filters used in post-processing were obtained through a brute-force test iterating through all possible combinations. The success of watershed transform algorithm used to separate HV and PV depends on the marker positions. We search for the local maxima near two placed seed markers and extend the number of seeds by adding neighboring points. The correction of seed markers helps improve the robustness and reliability of the separation method. Generally, segmenting HA and BD using region-growing methods yields reasonably good results due to high contrast. Unsatisfactory segmentation can be efficiently corrected with manual editing. Due to time constraints, we did not adopt advanced segmentation techniques such as level-set or graph-based methods, which could be part of future work. 5 [1] [2] [3] [4] [5]