Unsupervised Three-Dimensional Tubular Structure Segmentation via Filter Combination

Tubular structure enhancement plays an utmost role in medical image segmentation as a pre-processing technique. In this work, an unsupervised 3D tubular structure segmentation technique is developed, which is mainly inspired by the idea of filter combination. Three well-known vessel filters, Frangi’s filter, the modified Frangi’s filter and the Multiscale Fractional Anisotropic Tensor (MFAT) filter, separately enhance the original images. Next, the enhanced images obtained using three different filters are combined. Different categories of vessel filters have the ability of complementarity, which is the main motivation of combining these three advanced filters. The combination of them ensures a high diversity of the enhancing results. Weighted mean and median ranking methods are used to conduct the operation of filter combination. Based on the optimized weights for all the three individual filters, fuzzy C-means method is then applied to segment the tubular structures. The proposed technique is tested on the public DRIVE and STARE datasets, the public synthetic vascular models (2011 and 2013 VascuSynth Sample), and real-patient Coronary Computed Tomography Angiography (CCTA) datasets. Experimental results demonstrate that the proposed technique outperforms the state-of-the-art filter combination-based segmentation methods. Moreover, our proposed method is able to yield better tubular structure segmentation results than that of each individual filter, which exhibits the superiority of the proposed method. In conclusion, the proposed method can be further used to facilitate vessel segmentation in medical practice.


Introduction
In medical image processing and segmentation, tubular structure enhancement plays an utmost role as a medical image pre-processing technique. For example, efficient tubular structure enhancement in clinical practice results in fast and accurate coronary arteries' segmentation, which significantly contributes to the cardiovascular disease screening [1]. Traditional medical image enhancement is usually actualized by improving data acquisition techniques, and performing model-based algorithms [2]. A wide range of tubular structure segmentation approaches are established in the literature [3]. Generally, these approaches are divided into two main types: supervised methods and unsupervised methods [4].
Supervised methods [5][6][7][8] are usually easy to be implemented, however, they are dependent on numerous ground truth images with labels. As described by Hoover et al. [9] and Fraz et al. [10], different expert observers may have significant disagreement in identifying the vessel structure regions. Second, the intensity difference between the vessel structures and the background near the boundary of a vessel is so low. This kind of uncertainty leads to that some problems cannot be clearly labeled. Third, vessel structures near the pathology, e.g., coronary artery stenosis lesion, are usually difficult to be labeled [9].
On the other hand, unsupervised techniques [9,[11][12][13] are normally capable of revealing patterns, correlations and clusters, which can be performed without hand-labeled ground truth images. Many tubular structure enhancement techniques have been developed based on the Hessian matrix [14][15][16][17][18]. Frangi's filter [14] is the pioneer work to capture vessel structures by exploiting the eigensystem of the Hessian matrix, which calculates the second-order Gaussian derivatives for multiple values of standard deviation. This kind of multiscale operation helps to detect vessels of different sizes. A novel enhancement filter is proposed in [15] , which utilizes the ratio of the Hessian eigenvalues in a multiscale manner. Moreover, a modified Frangi's filter (mFrangi) is developed in Cui et al. [17], which multiplies the original Frangi's filter by the novel exponential term, for the purpose of tackling the problem of non-smoothness at the origin. Nevertheless, only local geometric features are explored by the Hessian matrix. Other kinds of vascular shapes may also be enhanced during the filtering process. Therefore, it increases the difficulty in identifying the target tubular structures from the surrounding tissues.
Another category of unsupervised tubular structure segmentation approaches rely on the use of diffusion tensors [2]. For instance, the Fractional Anisotropic Tensor (FAT) [19] captures different shapes of tubular structures by measuring the variance of anisotropy. Due to its strong ability of diffusion tensor regularization, FAT holds great potential in tubular structure enhancement. To better return uniform vessel responses and detect junctions, Alhasson et al. [2] developed a Multiscale Fractional Anisotropic Tensor (MFAT). However, the introduction of too many threshold values and parameters makes the MFAT method difficult to be duplicated for other medical datasets. Additionally, the principle limitation of diffusion tensors is that they are not able to simulate other complex image structures.
Oliveira et al. [4] proposed an unsupervised approach, which combined matched filter (MF) [20], Frangi's filter (FR) and Gabor Wavelet (GW) filter [7] to automatically segment vessel structures from the retinal images. These three filters are based on different concepts, and each individual of them produces relatively good result. The combination of them ensures a high diversity of the enhancing results. The segmentation results enhanced by combining several individual filters are better than that of each single filter, which is the necessary condition for filter combination. However, the work by Oliveira et al. [4] is only applicable to 2D retinal images.
In this work, an unsupervised 3D tubular structure segmentation technique is developed, which is mainly inspired by the idea of filter combination. First, the Hessian matrix of the input images is computed. The second step, i.e., image enhancement, is the major novelty of this work. Given the obtained eigenvalues, Frangi's filter, the modified Frangi's filter and the MFAT filter separately enhance the original images. Next, the enhanced images obtained using three different filters are combined. Third, weighted mean and median ranking method are used to conduct the operation of filter combination. Based on the optimized weights for all the three individual filters, fuzzy C-means method is then applied to segment the tubular structures. The proposed technique is completely unsupervised, and yields better tubular structure segmentation results than that of each individual filter. Figure 1 shows the flowchart of the proposed vascular structure segmentation technique. Our main contributions are as follows: Fig. 1 Flowchart of the developed vascular structure segmentation technique. First, the Hessian matrix of the input images are computed. The second step, i.e., image enhancement, is the major novelty of this work. Given the obtained eigenvalues, Frangi's filter, the modified Frangi's filter and the MFAT filter separately enhance the original images. Next, the enhanced images obtained using three different filters are combined. Third, weighted mean and median ranking method are used to conduct the operation of filter combination. Based on the optimized weights for all the three individual filters, fuzzy C-means method is then applied to segment the tubular structures • A framework to segment tubular structures in medical images using an unsupervised method is presented. Three well-known vessel filters, Frangi's filter, the modified Frangi's filter and the Multiscale Fractional Anisotropic Tensor (MFAT) method, are combined. • Different categories of vessel filters have the ability of complementarity, and the combination of them is shown to ensure a high diversity of the enhancing results. • Weighted mean and median ranking method are used to conduct the operation of filter combination. The superiority of the combination strategy is exhibited on a wide range of validation datasets.
The remaining part of this work is organized as follows. The developed technique for vascular structure segmentation based on filter combination is described in Sect. 2, which can be divided into three parts: (1) filter combination; (2) weight optimization; (3) tubular structure segmentation. Section 3 gives the experimental results on the public DRIVE and STARE datasets, the public synthetic vascular models (2011 and 2013 VascuSynth Sample), and real-patient CCTA datasets. Moreover, the filter combination-based segmentation technique for tubular structure segmentation is discussed in Sect. 4. Finally, conclusions are given in Sect. 5.

Materials and Methods
The developed technique for vascular structure segmentation based on filter combination can be divided into three parts: (1) filter combination; (2) weight optimization; (3) tubular structure segmentation. First, the Hessian matrix of the input images are computed. The second step, i.e., image enhancement, is the major novelty of this work. Given the obtained eigenvalues, Frangi's filter, the modified Frangi's filter and the MFAT filter separately enhance the original images. Next, the enhanced images obtained using three different filters are combined. Third, weighted mean and median ranking methods are used to conduct the operation of filter combination. Based on the optimized weights for all the three individual filters, fuzzy C-means method is then applied to segment the tubular structures.

Frangi's Filter
The eigensystem of the Hessian matrix, which represents the second-order image information (curvature), is widely used for geometrical interpretation of tubular structures. The curvature extreme values are stored in the eigenvalues. Original Frangi's vesselness filter is defined as follows [14]: where | 1 | ≤ | 2 | ≤ | 3 | are the eigenvalues, , and are parameters. Moreover, A, B and S are used to capture different vascular structures:

Modified Frangi's Filter
Unfortunately, Frangi's filter is not smooth at the origin, which needs, therefore, to be modified. Consider the secondorder partial derivative of 1 [21]: Both V F and V (n) F have the following terms: where M, N represent polynomial functions in , and U, V are polynomial functions in without any constant terms. U, V will remain the same and M, N vary when taking any Therefore, T ′ is well defined, which approaches 0 when → 0 . For all the exponential terms of V F , it can be shown that | 2 | 2 3 is the smallest common denominator. Therefore, modified Frangi's filter is multiplied by the term e − 2c 2 3 . The small value of parameter c ensures that it has influence only when → 0 . The new vesselness function V m can be proposed [22]: The multiscale approach is performed by computing the Hessian H by taking the 2nd Gaussian derivatives and normalizing by 2 at multiple scales. Then, the maximum vesselness response is

MFAT
The Fractional Anisotropic Tensor (FAT) [19] captures different shapes of tubular structures by measuring the variance of anisotropy. The FAT value is usually computed as The value of FAT ∈ [0, 1] . D indicates the mean diffusivity and is given by D = T r ∕3 . T r is the trace of diffusion tensor, which is defined as Figure 2 depicts different diffusion ellipsoids spanned by the eigenvectors. However, the main limitation of eigenvaluebased measures is that they can be described as only partial representations of tensor information. Recently, Pardos et al. [23] use anisotropic diffusion tensors in probabilistic form, which is called the Probabilistic Fractional Anisotropic Tensor (PFAT). PFAT is proved to have better performance of capturing tubular structures, which is given by where p i represents relative importance and is defined as and D p is set to be 1 3 . In this work, we combine both forms of fractional anisotropic tensor and make several modifications: (i) absolute values are added to account for differently signed eigenvalues, which results in a more uniform response. (ii) 1 is eliminated to get more normalized results. (iii) 3 is regularized at each scale based on cutoff threshold value . The new filter function is defined as where threshold p corresponds to p , another cutoff threshold value v is used to regulate 3 at any point in each scale to fulfill the following condition: Here, p and v are obtained from: Therefore, by introducing the above eigenvalues regularization, Equation (13) can be rewritten as follows: where the coefficients of relative importance are defined as The inverted responses of Equation (18) at vessel positions give positive values. Therefore, the vessel response R ,p can be defined as The junction positions are also reconstructed by magnitude regularization, and a maximized co-addition of response at each scale is computed. Therefore, the final filter function MFAT is defined as where indicates the current scale, − 1 is the previous scale, and represents the step size.

Filter Combination
Different categories of vessel filters have the ability of complementarity, which is the main motivation of filter combination. The combination of three filters, Frangi's filter, the modified Frangi's filter and the MFAT filter, leads to the improvements in the overall segmentation results, compared with each single filter. In this section, weighted mean method and median ranking method are described in detail for filter combination.
Let G denote the enhanced images and W denote the set of non negative weights. The weighted mean F is normalized between 0 and 1, which is given by with g ∈ G and w ∈ W , and the number of filtered images N. The weights should be adjusted for the purpose of enhancing vessel structures of small sizes and suppressing background noises. The weights are constrained to: Genetic algorithm (GA) [24] is widely used as an optimization technology. The basic idea is that it evolves through selection and reproduction from a population, which comprises the set of individuals that represent possible solutions. The final set of parameters minimizing the object function represents the desired solution. For each element in the set of all possible filter combinations , the above algorithm is executed. Assume the number of training images is denoted as M. For every combination ∈ , there are a total of N filters. Then, the weights W ∈ W are chosen by Genetic Algorithm at each iteration. Let I i denote a training set, where i = 1, 2, … ., M . The filtered images subset G i ∈ G i are computed by the proposed algorithm. Next the combined image F i is obtained based on Eq. (23). Fuzzy C-means method is then used to segment the combined image F i . The weights W are optimized by the algorithm according to the fitness function (Eq. 25). With the number of iterations increases, the value of the fitness function f declines. When the minimum vale of f is reached, the maximum number of iterations is chosen and defined as the stopping criterion. In this work, the stopping criterion is experimentally selected as 50 iterations. Once the stopping criterion is reached, the whole algorithm terminates. Finally, the proposed technique returns the optimized weights W ′ that minimizes f: where A i is the accuracy of the image segmentation i defined by with TP for True Positive, TN for True Negative, FP for False Positive and FN for False Negative. As described in [7,9], the accuracy is calculated based on the hand-labeled ground truth images.
It is described by Belkin et al. [25] that when multiple systems have incompatible scores, ranked outputs-based methods are more suitable for the combination. This characteristic is also applicable for the enhanced images using different filter methods. The enhanced images are first normalized. The median ranking is then applied to combine the normalized images. For each filtered image, all the voxels are ordered from high to low based on their intensities. The Pearson's correlation coefficients [26] for the filtered images obtained by three different filters are given in Table 1.
For the pixel p in the enhanced image obtained using filter k, k = 1, 2, … , N , where N denotes the total number of filters, the position in the rank is denoted as r K p . The pixel with higher intensity of image enhanced by filter k is assigned rank by r K p = top , the runner-up, r K p = top − 1 , and so on, where top is the total number of pixels. The median ranking of p is thus computed by

Fuzzy C-Means Segmentation
In the literature about segmentation and pattern recognition, the operation of clustering usually assembles a collection of objects into groups. The objects that belong to the same group are more similar to each other than the objects in other groups. In this part, the collection of enhanced images obtained by the filter combination are segmented by applying the fuzzy C-means algorithm. As a widely used technique for clustering problems, the fuzzy C-means is improved by Bezdek [27]. One of the main advantages of the fuzzy C-means algorithm is that it enables that one object can be assigned to one or more groups to a certain degree. In this work, it is observed that great variations in grayscale values exist for the non-vascular structures. Therefore, c groups are selected, in which the vessel voxels are grouped to one cluster and the voxels of non-vessel structures are grouped to c − 1 clusters.

Results
Three different types of datasets, i.e., retinal images (DRIVE [28] and STARE [9]), 3D synthetic tubular structures and CCTA datasets from the hospital, are used to demonstrate the accuracy of the proposed segmentation approach. For the modified Frangi's filter, the parameter c is determined experimentally as 10 −5 . For the DRIVE dataset, 20 training images are provided. The Genetic Algorithm is performed over the training set. On the other hand, only 20 testing images are provided for the STARE dataset. In this work, only 6 retinal images are selected as training images, the remaining images are treated as testing images. The number of 6 images represents 30% of the total dataset, which is a compromise between a good representation of the training set and a sufficient number of testing images.
The parameters are selected experimentally by configuring the algorithm. For the DRIVE dataset, the Genetic Algorithm is able to select the best weights of 0.1046, 0.3837 and 0.5117 for the Frangi's filter, the modified Frangi's filter and the MFAT respectively. While for the STARE dataset, the best weights are found to be 0.0945, 0.3782, and 0.5273 for the Frangi's, the modified Frangi's and the MFAT, respectively. The fuzzy C-means segmentation method requires the number of classes C to be defined. In this work, two classes are used for both DRIVE and STARE dataset. For the median ranking method, a threshold value of 0.8611 is used to segment the retinal vessels for the DRIVE dataset, while a threshold value of 0.8807 is used for the STARE dataset. The two threshold values are optimized by configuring the Genetic Algorithm on the 20 training images of the DRIVE dataset and the 6 training images of the STARE dataset.

Retinal Images
In the literature, two public datasets DRIVE [28] and STARE [9], are widely used to validate the accuracy of the vessel structure segmentation approaches. The Area Under the Curve (AUC) value is used in this part as the evaluation metric. The AUC value normally varies between 0 and 1. When the AUC value is higher, the segmentation result is closer to the ground truth. The value of 1 means that the segmentation result is identical to the ground truth.
The comparison of the AUC values for different kinds of segmentation approaches over the DRIVE and STARE datasets is summarized in Table 2. Frangi's filter can only achieve the AUC values of 0.8880 and 0.8980 for the DRIVE and STARE dataset. The filtering results are poor, which may be explained by the fact that it dose not perform well for the vessels of small sizes. The GW filter achieves poor segmentation result on the STARE dataset. The MF filter can achieve relatively good results for both datasets. Oliveira's method   [14] 0.8880 0.8980 GW [7] 0.9306 0.8117 MF [20] 0.9309 0.9422 Oliveira's [4] 0.9379 0.9392 mFrangi [17] 0.9381 0.9472 MFAT [2] 0.9400 0.9500 Proposed 0.9634 0.9657 is able to improve the AUC values up to 0.9379 and 0.9392, which is found to outperform each single filter. Unfortunately, it can be seen that the improvement in the segmentation accuracy is marginal. Moreover, the MFAT method and the modified Frangi's filter both slightly outperform Oliveira's method. Finally, the proposed method, which combines Frangi's filter, modified Frangi's filter and MFAT, is shown to achieve the AUC value of 0.9634 and 0.9657 for the public DRIVE and STARE dataset respectively. Regarding retinal vessel structure segmentation, the proposed method outperforms Oliveira's method, which is the state-of-the-art filter combination-based segmentation methods. Furthermore, the proposed filter combination method outperforms the modified Frangi's method and the MFAT method, which reveals that the combination strategy in this work can significantly improve the segmentation results. Figure 3 presents two examples of segmentation results on DRIVE dataset using the modified Frangi's filter, the MFAT method and the proposed method based on filter combination. The first column (Fig. 3a, e) gives the original retinal images in the DRIVE dataset. The second column (Fig. 3b, f) gives the segmentation results using the modified Frangi's filter. Third, the segmented retinal vessels obtained by the MFAT method are presented in Fig. 3c, g. Finally, the segmentation results by the proposed method based on filter combination are shown in Fig. 3d, h. Our proposed method is able to detect more vessel structures of small sizes, which outperforms the modified Frangi's filter and the MFAT method.

Synthetic Vascular Models
The proposed filter approach was evaluated using the public March 2013 VascuSynth Sample (10 datasets) [29,30]. Each dataset represents a 3D synthetic vascular tree with a volume of 100 × 100 × 100 voxels [29]. Three synthetic tubular structure models in Group 4 (data3, data4 and data5) depicted in Fig. 4, are used as the ground truth to demonstrate the accuracy of the proposed filter method. Four types of evaluation metrics, i.e., true positives (TP), false negatives (FN), false positives (FP) and overlap measure (OM) [30] between the ground truth tubular models and the computed vessel structures, are used in this part for accuracy validation. The average TP, FN, FP, and OM measures of the proposed segmentation technique, for all the ten datasets of 2013 VascuSynth Sample, are summarized in Table 3. The average values with standard deviations of TP, FN, FP, and OM measures turn out to be 97.61 ± 0.56, 2.39 ± 0.56, 5.03 ± 1.18 and 97.70 ± 0.57, respectively.
Furthermore, the 2011 VascuSynth Sample Data [29,30] was used to demonstrate the proposed method's ability of detecting complete vascular structures from noisy images. Low level of Gaussian noise ( 2 = 20 ) was added to the 2011 VascuSynth Sample Data. The proposed method was tested on all the three datasets and compared with the MFAT method. Table 4 gives the comparison results. The proposed method achieves higher values for TP and OM, and lower

CCTA Datasets
The proposed tubular structure segmentation technique was further tested on three real-patient CCTA datasets, that is ccta 1, ccta 2 and ccta 3. The volume sizes of these three datasets were 320 × 320 × 460 , 410 × 410 × 380 , and 400 × 400 × 470 , respectively. Complete coronary artery regions, i.e., left anterior descending artery (LAD), left circumflex artery (LCX) and right coronary artery (RCA), were manually labeled by a professional cardiologist as the ground truth. Table 5 shows the comparison of the evaluation metrics (Dice, Precision and Sensitivity) obtained by different kinds of methods for the real-patient CCTA datasets. The proposed segmentation approach based on filter combination is capable of achieving the highest measure of Dice, Precision and Sensitivity. Among all the CCTA datasets, the proposed method is able to gain the average value of 0.906, 0.927 and 0.917 for Dice, Precision and Sensitivity, respectively.
Moreover, visual inspection is also conducted to demonstrate the effectiveness and accuracy of the proposed segmentation technique. Figure 5 depicts a comparison of  coronary artery lumen region segmentation results obtained using the MFAT method (yellow) and the proposed method (blue) with the ground truth regions (red). The lumen shapes produced by our proposed technique are observed to be closer to the ground truth regions. Furthermore, Fig. 6 presents a detailed comparison of the segmented coronary artery regions by the proposed method with the ground truth regions is presented in three views, for ccta 1, ccta 2 and ccta 3. The segmented coronary artery regions of our method are marked in blue, while the hand-labeled ground truth regions are marked in red. Figure 6a-f is for ccta 1, Fig. 6g-l is for ccta 2 and Fig. 6m-r is for ccta 3. The first two columns give the axial view, the second two columns give the sagittal view and the third two columns show the coronary view. The segmentation results produced by our proposed technique are highly close to the hand-labeled ground truth regions, which exhibits the superiority of our proposed method. Figure 7 shows the rendering of the segmented coronary arteries.

Discussion
As described in [4] On the other hand, the sample size of hand-labeled ground truth images is relatively small, especially the CCTA datasets. The main reasons can be explained as follows. First, different professional cardiologists have significant disagreement in identifying the coronary artery vessels. Second, the intensity difference between the vessel structures and the background near the boundary of a vessel is so low. This kind of uncertainty makes it difficult to correctly label the vessel structures. Third, vessel structures near the pathology, e.g., coronary artery stenosis lesion, cannot be clearly labeled. In our future work, the sample size can be increased by performing the operation of data augmentation. In addition, multiscale-and multi-view-based techniques can be used. Finally, the idea of transfer learning will also be considered.

Conclusion
This work presents a framework to segment tubular structures in medical images using an unsupervised method. Three basic vessel filters, Frangi's filter, the modified Frangi's filter and the MFAT method, are combined. These three filters are based on different concepts, and each individual of them produces relatively good result. Different categories of vessel filters have the ability of complementarity, which is the main motivation of combining these three advanced filters. The combination of them ensures a high diversity of the enhancing results. Weighted mean method and median ranking method are performed to combine the basic filters. Compared with each single filter, the final combined filter benefits from multiple intensity values to the same voxel, which yields better responses. The combination of three filters leads to the improvements in the overall performance. The proposed method is shown to achieve the AUC value of 0.9634 and 0.9657 for the public DRIVE and STARE dataset, respectively. Regarding retinal vessel structure segmentation, the proposed method outperforms Oliveira's method, which is the state-of-the-art filter combination-based segmentation methods. Furthermore, the proposed filter combination method outperforms the modified Frangi's method and the MFAT method, which reveals that the combination strategy in this work can significantly improve the segmentation results. Moreover, on the 2013 VascuSynth Sample, the average values with standard deviations of TP, FN, FP, and OM measures turn out to be 97.61 ± 0.56, 2.39 ± 0.56, 5.03 ± 1.18 and 97.70 ± 0.57, respectively. The proposed method achieves higher values for TP and OM, and lower values for FN and FP. With the presence of low level noise, the proposed segmentation method is able to outperform the MFAT method. Among all the three CCTA datasets, the proposed method gains the average value of 0.906, 0.927 and 0.917 for Dice, Precision and Sensitivity, respectively. The segmentation results produced by our proposed technique are highly close to the handlabeled ground truth regions, which exhibits the superiority of our proposed method. The proposed method is observed to outperform the state-of-the-art filter combination-based segmentation methods. Furthermore, the proposed filter combination method outperforms each individual filter, which reveals that the combination strategy in this work can significantly improve the segmentation results. In conclusion, the proposed method can be used as a medical images pre-processing tool and facilitate vessel segmentation.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.