1. Introduction

Automatic segmentation of coronary arteries in X-ray angiography images can help cardiologists for diagnosing and treating vessel abnormalities. Hence, the development of efficient methods for automatic vessel segmentation has become essential for computer-aided diagnosis (CAD) systems. The main disadvantages of X-ray coronary angiograms are the nonuniform illumination and the low contrast between blood vessels and image background which are illustrated in Figure 1. Since the presence of several shades along blood vessels generates multimodal histograms, the vessel segmentation problem has been commonly addressed in two different steps. The first step is vessel enhancement, where the methods try to remove noise from the image and to enhance the vessel-like structures. After vessel enhancement, the second step consists in the classification of vessel and nonvessel pixels by using different strategies such as the selection of an optimal threshold value.

Figure 1 First row: three coronary angiograms. Second row: hand-labeled ground-truth of angiograms of the first row, as drawn by an expert. Third row: histograms of the angiograms. 

In recent years, several techniques have been reported for this purpose in different types of clinical studies. Some of the proposed methods are based on mathematical morphology. The method of (^{Eiho and Qian, 1997}) represents one of the most commonly used interactive vessel segmentation methods, because its efficiency and ease of implementation. The method uses the single-scale top-hat operator (^{Serra, 1982}) to enhance vessel-like structures followed by different morphological operators such as erosion, thinning, and watershed transformation. This method was improved in (^{Qian et al., 1998}), by introducing the multiscale top-hat operator using the linear combination of fixed scales of the top-hat operator in order to capture different diameters of blood vessels. Another morphological-based method was proposed in (^{Maglaveras et al., 2001}), which it uses the skeleton, thresholding and connected-component operator for automatic segmentation of coronary arteries in angiograms. (^{Bouraoui et al., 2008}) introduced an automatic segmentation method for coronary angiograms consisting in the application of the hit-or-miss transform and the region-growing strategy. (^{Lara et al., 2009}) proposed an interactive segmentation method based on the region-growing strategy and a differential-geometry approach for coronary arteries in angiograms. In general, these morphological-based methods provide satisfactory results when vessels have high contrast; however, they do not perform well when vessels have different diameters and low contrast.

Due to the shortcomings of the morphological-based methods for automatic segmentation of coronary angiograms, some strategies using spatial filtering have been introduced. The most commonly vessel enhancement technique is the Gaussian matched filters (GMF) method proposed by (^{Chaudhuri et al., 1989}) and improved by (^{Al-Rawi et al., 2007}). The fundamental idea of this method is to approximate the shape of blood vessels through a Gaussian curve as the matching template in the spatial image domain. The Gaussian template is rotated at different orientations and convolved with the original angiogram. From the filter bank of oriented responses, the maximum response at each pixel is taken to obtain the final enhanced image. The GMF method is used as a vessel enhancement procedure in different automatic segmentation methods. (^{Chanwimaluang and Fan, 2003}) used the GMF method to enhance retinal blood vessels, and then applied the entropy thresholding method proposed by (^{Pal and Pal, 1989}) to obtain the final vessel segmentation result. This automatic segmentation method was used successfully in (^{Chanwimaluang et al., 2006}) for a retinal vessel registration system. (^{Kang et al., 2009}) proposed a fusion enhancement strategy for coronary arteries where the top-hat operator and the GMF are applied separately on the original angiogram. Both resulting enhanced images are thresholded by the maximizing entropy thresholding method and the segmentation result is acquired by pixel-wise multiplication of the two binary images. This method was also used in (^{Kang et al., 2013}) just replacing the thresholding method by the degree segmentation method. These template matching methods present low efficiency when images contain occlusions or vessels with different diameters to be detected.

More recently, different vessel enhancement methods based on the properties of the Hessian matrix have been proposed (^{Lorenz et al., 1997}; ^{Frangi et al., 1998}; ^{Jin et al., 2013}; ^{Xiao et al., 2013}). The Hessian matrix is computed by convolving the second-order derivatives of a Gaussian kernel with the original angiogram. The method of (^{Frangi et al., 1998}) uses the properties of the eigenvalues of the Hessian matrix to compute a vesselness measure. This measure is used to classify vessel and nonvessel pixels from the resulting image of dominant eigenvalues at different vessel width scales. Taking into account the properties of the Hessian matrix, different vesselness measures have been introduced (^{Salem and Nandi, 2008}; ^{Tsai et al., 2013}; ^{Mhiri et al., 2013}). (^{Li et al., 2012}) and (^{Wang et al., 2012}) proposed different vesselness measures, which are segmented by using similar procedures of the region-growing operator. These methods have the ability of detecting vessels of different caliber although they are highly sensitive to noise because of the second-order derivative.

Moreover, from the revised morphological-based methods, the top-hat operator presents low efficiency since it is not able to capture vessels of different caliber, which is an advantage of the methods based on the properties of the Hessian matrix. In the present paper, a novel method based on two stages for automatic enhancement and segmentation of coronary arteries in angiograms is proposed. In the enhancement stage, a new multiscale top-hat operator based on the eigenvalues of the Hessian matrix is introduced in order to capture most of the vessel diameters in the angiograms. The method is compared with the multiscale top-hat operator proposed by (^{Qian et al., 1998}), Gaussian matched filters (^{Al-Rawi et al., 2007}), and the Hessian-based method proposed by (^{Wang et al., 2012}) using the area A_z under the ROC curve. In the segmentation stage, the multiscale top-hat response is thresholded by a new global thresholding technique based on multiobjective optimization using the weighted sum method (^{Marler and Arora, 2010}). In this stage, seven state-of-the-art thresholding methods are compared with our multiobjective method using the measures of sensitivity, specificity, and accuracy. In addition, the segmentation results of the proposed method are compared with those obtained with five state-of-the-art vessel segmentation methods.

The remainder of this paper is organized as follows. In Section 2, the database of coronary angiograms and the fundamentals of the proposed vessel segmentation method along with a set of evaluation metrics are introduced. The experimental results are discussed in Section 3, and conclusions are given in Section 4.

2. Materials and Methods

2.2 Proposed Vessel Segmentation Method

Since coronary arteries present lower reflectance compared with the background angiogram, the proposed coronary artery segmentation method is composed of two steps. In this section, the fundamentals of the proposed multiscale top-hat operator for vessel enhancement and the thresholding strategy based on multiobjective optimization approach are described in detail.

2.2.1 Vessel Enhancement

The coronary arteries in angiograms can be described as tubular structures with a nonuniform illumination, different diameters and diverse orientations. Due to these factors, a vessel enhancement procedure is commonly used as a preprocessing step before applying a specific segmentation method.

The morphological top-hat operator has been widely used to suppress local noise and to enhance vessel-like structures in angiograms (^{Eiho and Qian, 1997}; ^{Qian et al., 1998}; ^{Kang et al., 2009}; ^{Kang et al., 2013}). The top-hat operator can be defined as follows:

(1)

where Io is the original angiogram, Ie is the resulting enhanced angiogram, and γ(Io, A) represents the grayscale morphological opening operator applied on the original angiogram with a predefined structuring element (SE). In order to obtain the best performance by using the top-hat operator, the parameters of size and shape of the SE have to be tuned. In general, the size parameter is chosen slightly larger than the maximum caliber of the vessel and a disk-shaped SE is commonly used. The main disadvantage of some previously mentioned vessel enhancement methods on applying the conventional top-hat operator, is the fact that by using only one structuring element with a fixed size, cannot be properly enhanced all of the vessel diameters in the angiogram (^{Qian et al., 1998}). In our work, we propose the use of the eigenvalues computed from the Hessian matrix, in order to select the best scales for the multiscale top-hat operator.

The Hessian matrix is computed from the second-order derivatives of a Gaussian kernel at each pixel of the intensity angiogram L(x, y) as follows:

(2)

where Lxx, Lxy, Lyx, and Lyy denote the second-order derivatives of the intensity image as follows:

(3)

(4)

(5)

where the symbol “ * ” represents a convolution operator, σ is the the length scale factor, and Gxx, Gyy, and Gxy are the second-order derivatives of the Gaussian kernel G at different scales of σ, which is defined as follows:

(6)

Due to the Hessian matrix is symmetric; the large eigenvalue λ₁ can be computed as follows:

(7)

(8)

The largest eigenvalue λ_max over all the scales of the σ value at each image pixel is calculated as following

(9)

The final step of the enhancement stage consists in computing the quartiles method of descriptive statistics from the λ_max image of vessel width scales. Figure 2 illustrates the distribution of the vessel scales of the 40 training angiograms, which is acquired by the number of pixels for each scale in the corresponding angiogram. In our experiments and based on the ROC curve analysis using the training set of angiograms, the second and third quartiles contain the optimal scales for the image to be processed. These vessel scales are used as the size parameter for the multiscale top-hat operator. For each vessel scale, a disk-shaped structuring element of fixed size (scale) is applied on the angiogram.

Figure 2 Distribution of the vessel width scales according to the number of pixels using the training set of 40 angiograms. 

Moreover, to detect most of the vessels at different diameters, the final enhanced image is obtained by preserving the maximum response for all the enhanced pixels at different scales. Since the number of optimal scales for each angiogram is not constant, the proposed enhancement method can automatically detect the appropriate vessel scales for a particular angiogram, which is one of the main advantages regarding morphological-based methods.

2.2.2 Thresholding based on Multiobjective Optimization Approach

Multiobjective optimization (MOO) is the process of optimizing systematically and simultaneously a number of different objective functions (^{Marler and Arora, 2010}). One of the most widely used strategies for solving MOO problems is the weighted sum method (WSM). The WSM consists in the arrangement of different objective functions and the selection of a particular set of weights to form a single objective function as follows:

(10)

where n is the number of objective functions of the problem to be optimized, w_i is the weighting factor for each function, and F_i is the particular objective function. Generally, the weight factors are positives, established by the preference of the user, and their cumulative sum has to be equal to 1.

Although a multiobjective optimization problem has multiple solutions, the weighted sum method provides a single optimal solution using a set of fixed weights, which represents an ideal case to classify vessel and nonvessel pixels for our coronary vessel segmentation application. The proposed thresholding method uses three objective functions, as it is shown in the following composite objective function:

(11)

where w₁, w₂, and w₃ are the assigned weights. In our work, these weights were statistically determined as w₁ = 0.3, w₂ = 0.4, and w₃ = 0.3, respectively, using the average of the sensitivity, specificity and accuracy measures on the training set of 40 angiograms. F₁ is the between-class variance criterion based on the method of (^{Otsu, 1979}). F₂ and F₃ represent entropy criteria based on the methods of (^{Kapur et al., 1985}) and (^{Pal and Pal, 1989}), respectively, which are described below.

The between-class variance criterion is applied on the histogram of intensity levels L of an image. This criterion uses the mean intensity of the given image defined as μ, the probability distribution W_j and the mean intensity for each class of pixels μ_j, where n = 2, since only two classes of pixels are required (vessel and nonvessel pixels). The intensity level with the maximum between-class variance gives the optimal threshold value wich is acquired by maximizing the following:

(12)

The entropy criterion of (^{Kapur et al., 1985}) is given by the entropies of the corresponding classes A and B of pixels in which the image will be segmented. The entropy of class A and class B is computed using the accumulated probabilities P_A and P_B as follows:

(13)

(14)

The optimal threshold value of the intensity levels L is acquired by maximizing the total entropy as follows:

(15)

The two-dimensional entropy criterion of (^{Pal and Pal, 1989}), is based on the co-occurrence matrix of the image and the second-order entropy. The second-order entropy corresponding to the foreground and background image can be defined as follows:

(16)

(17)

where s represents a threshold value, L is the number of intensity levels of the input image, and P^A and P^B represent the probability distributions of the two classes of pixels through the co-occurrence matrix. The total second-order entropy of the image has to be computed for all the intensity levels, and the optimal threshold is the value that maximizes the entropy by using the following:

(18)

In the proposed thresholding method, each criterion is normalized to the range 0,1, and the optimization process is carried out over all the intensity levels in the multiscale top-hat response. Since a multiobjective optimization problem has multiple solutions (segmentation results), in Figure 3, the Pareto front and the corresponding vectors of weights acquired from the top-hat response of the training set of angiograms are introduced. The advantage of multiobjective optimization following the weighted sum method resides in the fact that the weights can be determined such that the highest trade-off among the three objective functions is achieved.

Figure 3 From left to right, vectors of weights and Pareto front found by the proposed thresholding method through the multiscale top-hat response of the training set of angiograms. 

2.2.3 Evaluation metrics

In order to assess the performance of the proposed vessel segmentation method, the area under the receiver operating characteristic (ROC) curve (A_z) and measures of sensitivity, specificity, and accuracy, have been adopted which are described below.

The ROC curve is a plot of the true-positive rate (TPR) against false-positive rate (FPR) of a classification system. The TPR represents the fraction of the vessel pixels outlined by the specialist that are correctly detected by the method. The TPR is also known as the sensitivity measure which is computed using Eq.(19). FPR represents the fraction of the nonvessel pixels that are incorrectly detected as vessel pixels by the method.

(19)

Specificity measures the true-negative rate (TNR), which represents the fraction of nonvessel pixels that are correctly detected as such by the method (background pixels) as follows:

(20)

Accuracy represents the fraction of vessel and nonvessel pixels correctly detected by the method divided by the total number of pixels in the angiogram as follows:

(21)

In this work, the area under the ROC curve is used to determine the most suitable parameters of the vessel enhancement methods and also to evaluate their performance using the training set. The measures of sensitivity, specificity, and accuracy were used to assess the segmentation results between the regions obtained by computational methods and the regions outlined by the specialist. In these measures when the regions are completely superimposed the obtained result is 1, and 0, when these regions are completely different.

In Section 3, the segmentation results obtained from the proposed method on X-ray angiographic images are presented and analyzed by the evaluation metrics.

3. Results and Discussion

The computational experiments were tested on a computer with an Intel Core i3, 2.13 GHZ procesor, and 4GB of RAM using the Matlab software.

In the vessel detection stage, the proposed multiscale top-hat operator, the multiscale top-hat method (^{Qian et al., 1998}), the Hessian-based method (^{Wang et al., 2012}), and the Gaussian matched filters (^{Al-rawi et al., 2007}) were compared against each other using the training set of angiograms via ROC curve analysis.

The parameters of the multiscale top-hat method (^{Qian et al., 1998}) was obtained by varying the size of the structuring element and taking into account three different scales. The best set of scales was determined to be S₁ = 19, S₂ = 5, and S₃ = 3 pixels. The multiscale Hessian-based method (^{Wang et al., 2012}), was tested with different values of σ in the range [1, 20] pixels, and with step sizes (s) in the range [0.5, 4]. The most suitable parameters for vessel width σ, and for the step size were determined to be σ = [1, 10], and s = 0.5. Since the Gaussian matched filters represent the enhancement step of different vessel segmentation methods, the GMF was tuned with the same set of parameters for the following methods (^{Chanwimaluang and Fan, 2003}; ^{Kang et al., 2009}; ^{Kang et al., 2013}). The best Gaussian matched filters parameters for length of the vessel segment, width of the vessel segment, and number of oriented filters, were statistically determined to be L = 11 , σ = 1.9, T = 8, and κ = 12, respectively.

From the best set of parameters for each enhancement method, Figure 4 illustrates the area (A_z) under the ROC curves obtained using the training set of angiograms. The results of vessel enhancement show that the proposed multiscale top-hat operator provides a higher level of classification between vessel and nonvessel pixels than the comparative methods. The most suitable range of vessel width scale for the Gaussian kernel of the proposed enhancement method was set as σ = 1, 25, with a step size of s = 3. The proposed multiscale top-hat operator with the best set of parameters provided A_z = 0.965 with the test set of angiograms. In comparison, using the best set of parameters of the method of (^{Qian et al., 1998}), (^{Wang et al., 2012}), and Gaussian matched filters, they provided A_z = 0.938, A_z = 0.941, and A_z = 0.923, respectively, using the test set of angiograms. Additionally, Figure 5 shows a subset of angiograms with the corresponding ground-truth, and the vessel enhancement results of the proposed multiscale top-hat operator.

Figure 4 Comparative analysis of ROC curves for vessel enhancement with the training set of angiograms, using the proposed multiscale top-hat operator and the methods of (^{Qian et al., 1998}), (^{Wang et al., 2012}), and the Gaussian matched filters (^{Al-rawi et al., 2007}). 

Figure 5 First row: subset of original angiograms. Second row: ground-truth of images of the first row. Last row: results of applying the proposed multiscale top-hat operator on the angiograms of the first row. 

On the other hand, the proposed thresholding strategy based on the weighted sum method was compared with seven state-of-the-art thresholding methods using the multiscale top-hat response of the training set of angiograms. In Table 1, the segmentation performance is presented. In this experiment, the measures of sensitivity, specificity and accuracy were used for evaluation and comparison. These three measures have been widely used in vessel segmentation as the main measures for evaluation. The values of these measures were obtained by concatenating the forty segmented angiograms of the training set in order to form just one large binary image. However, since vessel pixels often occupy less than 15% of an angiogram, the measures of specificity and accuracy are always high. Consequently, the average of these three measures has been considered as the main performance metric. The obtained results show that the proposed multiobjective thresholding method provides the best segmentation performance using the training set of angiograms and hence it was selected for further analysis. Due to the thresholding based on multiobjective optimization uses the intensity values in the image histogram, a low computational time is required to perform the optimization task. The proposed thresholding strategy obtained an average execution time of 0.329 seconds for each angiographic image.

Table 1 Comparative analysis of the proposed multiobjective strategy with seven thresholding methods (^{Pal and Pal, 1989}; ^{Kapur et al., 1985}; ^{Kittler et al., 1985}; ^{Ridler and Calvard, 1978}; ^{Otsu, 1979}; ^{Rosenfeld and De La Torre, 1983}; ^{Tsai, 1985}) using the multiscale top-hat response of the training set of angiograms. 

Moreover, the proposed vessel segmentation method consisting of the steps of multiscale top-hat operator for vessel enhancement, and multiobjective thresholding for segmentation, was compared with five state-of-the-art vessel segmentation methods. The methods of (^{Qian et al., 1998}) and (^{Wang et al., 2012}) were tuned taking into account the aforementioned parameters. The methods of (^{Chanwimaluang and Fan, 2003}), and the fusion methods of (^{Kang et al., 2009}; ^{Kang et al., 2013}) were tuned according to the best parameters of the Gaussian matched filters. Since the two fusion methods of Kang et al., also use the classical top-hat operator as enhancement method, the size of the structuring element for the top-hat operator was set as S = 19, and it was obtained by using the area under the ROC curve with the training set of angiograms. In Table 2, the comparative analysis of the proposed method with the five vessel segmentation methods using the test set of angiograms is presented. In this experiment, the average of the sensitivity, specificity and accuracy measures is presented. The values of these measures were obtained by concatenating the forty segmented angiograms of the test set to form just one large binary image. The obtained results show that the proposed method provides the highest segmentation performance using the test set of angiograms.

Table 2 Comparative analysis of the proposed method with five state-of-the-art vessel segmentation methods using the test set of angiograms. 

Figure 6 introduces a subset of X-ray coronary angiographic images with the corresponding ground-truth for each image. The segmentation results of the method of (^{Qian et al., 1998}) shows difficulties to detect different diameters of vessels obtaining a high rate of false-positive pixels in most of the low contrast angiograms. The fusion method of (^{Kang et al., 2009}) based on the entropy of the enhanced images presents a low rate of false-positive pixels illustrating broken vessels, which reduces the rate of accuracy segmentation. The method of (^{Chanwimaluang and Fan, 2003}) can detect a high rate of true-positive and false-positive pixels, which increase the sensitivity rate and decrease the specificity and accuracy values. The results of the method of (^{Wang et al., 2012}) shows a uniform segmentation over the main tubular structures at different scales while it presents a low rate of true-positive pixels and broken vessels. The method of (^{Kang et al., 2013}) based on the degree segmentation method shows low performance in angiograms with an uneven illumination and high rate of false-positive pixels in most of the angiograms. The proposed method presents the highest average performance in the test set of angiograms. The method presents an appropriate rate of true-positive pixels, detecting vessels with different diameters and obtaining a low rate of broken vessels and false-positive pixels.

Figure 6 First row: subset of original angiographic images from the test set. Second row: ground-truth of the images of the first row. The remaining six rows present the segmentation results of the methods of (^{Qian et al., 1998}), (^{Kang et al., 2009}), (^{Chanwimaluang and Fan, 2003}), (^{Wang et al., 2012}), (^{Kang et al., 2013}), and the proposed method, respectively, applied on the angiograms of the first row. 

Although the state-of-the-art vessel segmentation methods provide an appropriate performance based on the quantitative and qualitative analysis discussed above, the comparative analysis reveals that the proposed method works well in nonuniform illuminated X-ray angiograms providing the highest average rate of vessel enhancement and segmentation. The segmentation results also show that the proposed method is more accurate and robust than the comparative methods considering the vessels outlined by the specialist, which is suitable for computer-aided diagnosis systems.

Conclusions

In this paper, we introduce a novel coronary artery enhancement and segmentation method in X-ray angiographic images. In the first stage, a new multiscale top-hat operator based on the eigenvalues of the Hessian matrix has demonstrated to be more efficient than three enhancement methods, achieving A_z = 0.942 with a training set of 40 angiograms and A_z = 0.965 with a test set of 40 angiograms. In the second stage, a new thresholding strategy based on a multiobjective optimization approach has obtained the highest average rate performance (0.911) regarding seven thresholding methods. Experimental results have proven that the proposed method consisting of multiscale top-hat operator and a thresholding strategy produces the highest performance rate (0.923) using the test set of angiograms compared with five state-of-the-art methods for automatic vessel segmentation. In addition, the experimental results have also shown that, based on the angiograms hand-labeled by a specialist, the performance in vessel pixel detection obtained from the proposed method is highly suitable for computer-aided diagnosis systems of cardiac abnormalities.