Automatic Segmentation of Coronary Arteries in X-ray Angiograms using Multiscale Analysis and Artificial Neural Networks

This paper presents a novel method for the automatic segmentation of coronary arteries in X-ray angiograms, based on multiscale analysis and neural networks. The multiscale analysis is performed by using Gaussian filters in the spatial domain and Gabor filters in the frequency domain, which are used as inputs by a multilayer perceptron (MLP) for the enhancement of vessel-like structures. The optimal design of the MLP is selected following a statistical comparative analysis, using a training set of 100 angiograms, and the area under the ROC curve (Az) for assessment of the detection performance. The detection results of the proposed method are compared with eleven state-of-the-art blood vessel enhancement methods, obtaining the highest performance of Az = 0.9775, with a test set of 30 angiograms. The database of 130 X-ray coronary angiograms has been outlined by a specialist and approved by a medical ethics committee. On the other hand, the vessel extraction technique was selected from fourteen binary classification algorithms applied to the multiscale filter response. Finally, the proposed segmentation method is compared with twelve state-of-the-art vessel segmentation methods in terms of six binary evaluation metrics, where the proposed method provided the most accurate coronary arteries segmentation with a classification rate of 0.9698 and Dice coefficient of 0.6857, using the test set of angiograms. In addition to the experimental results, the performance in the detection and segmentation steps of the proposed method have also shown that it can be highly suitable for systems that perform computer-aided diagnosis in X-ray imaging.


Introduction
The task of automatic segmentation of coronary arteries in X-ray angiograms can be integrated in a computer-aided diagnosis system (CAD) to support the analysis of coronary heart disease. In cardiology, the X-ray coronary angiogram is the gold standard employed to evaluate and monitor abnormalities of the coronary arteries. However, the detection of blood vessels in coronary angiograms is a challenging task for the current CAD systems. The most common problems that are present in coronary angiograms involve the non-uniform illumination, weak contrast between arterial features and image background, and occlusion caused by different organs or medical devices in the acquired images.
In the literature, the challenge of automatic detection of blood vessels has been widely addressed by different methods. A common approach is the use of operators and transformations from the mathematical morphology framework. The most used strategy to enhance vessel-like structures is the Top-hat operator, which was applied by Eiho and Qian [1] in coronary angiograms. The Top-hat response is combined with an interactive procedure to segment the main artery using another morphological operator. Later, a multiscale approach composed for the Top-hat filter responses of a number of structuring elements at different size was proposed by Qian et al. [2]. This method allows the automatic detection of vessels with different diameters in a single medical image.
Other common methodology is based on the information of the Hessian matrix of the medical image. In this approach, the gradient of an image is defined as the convolution of itself and the first order derivative of a Gaussian curve. Consequently, the Hessian matrix of the image is the second order derivative of a Gaussian curve convolved with the image. Frangi et al. [3] proposed the enhancement of tubular structures, such is the case of blood vessels, by using a vesselness measure defined from the eigenvalues of the Hessian matrix. For the particular problem of detection of coronary arteries in X-ray angiograms, multiple vesselness measures have been recently proposed. Wang et al. [4] used the eigenvalues of the Hessian matrix to construct a vessel probability map; in a similar manner, Li et al. [5] proposed the generation of a vessel feature map also from the Hessian matrix using a different vesselness measure. These two works use the vessel map as starting point for a region growing method to segment the coronary artery tree.
On the other hand, methods based on matching templates have considered the shape of a blood vessel in the spatial domain to propose a suitable model. Chaudhuri et al. [6] introduced a two dimensional Gaussian Matched Filter (GMF) for the enhancement of blood vessels in retinal fundus images. GMF is founded under the fact that in the spatial domain the cross-section of a blood vessel is Gaussian shaped. In addition, GMF assumes that vessels can be modeled in segments that share the same orientation and diameter. To detect vessel segments at different orientations, a directional filter bank is built by rotating the GMF template into a number of evenly spaced directions in the range [−π/2, π/2] to be later convolved with the input image. Methods based on the GMF have modified the parameters that govern the shape of the matching template and the number of orientations. Most of them have selected the parameters empirically, such is the case of the method proposed by Cinsdikici and Aydın [7], where the number of directional templates is increased in order to improve the detection performance. In the method proposed by Kang et al. [8], the parameters were empirically modified to allow the detection of arteries in coronary angiograms. Moreover, some methods have addressed the selection of the parameters by a training step prior to its application. A first approach is the proposed by Al-Rawi et al. [9], where a full-search of the parameters was employed. Later, Al-Rawi and Karajeh [10] proposed a method from the Evolutionary Computation algorithms (EC) to reduce the computational time consumed by the exhaustive search. More recently, Cruz et al. [11] performed a study to determine the most suitable EC method for the automatic selection of parameters for the GMF applied to the coronary artery segmentation problem.
An alternative to the analysis of the image in the spatial domain is the use of Gabor filters in the frequency domain. Gabor filters for blood vessel enhancement in medical images was proposed by Soares et al. [12] based on its capability of detecting oriented features. Similarly to GMF, Gabor filters can be rotated into different directions in order to detect vessel-like structures in multiple orientations. An important modification to the Gabor filter was proposed by Rangayyan et al. [13], where a simplified definition of the filter is presented, in addition to the use of only the real part of the filter, discarding the imaginary part which is commonly used to detect edge features. For Gabor filters, recent methods have been proposed to select the most appropriate parameters for the detection of blood vessels in image data. The method proposed by Rangayyan et al. [13,14] performs a full-search for the parameters of elongation, average width, and number of oriented filters before selecting the most suitable values for the detection of blood vessels in retinal fundus images. Later, Cervantes et al. [15] proposed the use of a Boltzmann based EC algorithm for the automatic selection of the Gabor filter parameters for the coronary artery segmentation problem.
Considering that the filter responses obtained by the blood vessel detection methods can be treated as multilevel classifiers, the Receiver Operating Characteristic curve analysis (ROC) has been frequently used to evaluate and compare their vessel detection efficiency.
Gabor filters have shown the best detection performance in coronary arteries compared to the state-of-the-art approaches [15,16], followed in efficiency by the Gaussian matched filters [11]. Because these two methods have obtained superior performance to carry out the automatic detection of blood vessels in medical images, in this work they are thoroughly studied. However, it is important to point out that coronary arteries present a natural variant thickness along their span; consequently, a multiscale approach is convenient in order to enhance vessel-like structures properly.
In this paper a novel method that uses the multiscale filter responses of the GMF and Gabor filters to improve the detection of coronary arteries in X-ray angiograms is presented. The proposed method addresses the limitation that single-scale versions of the GMF and Gabor filters present when enhancing vessel-like structures of a unique thickness at a time. In the method, a multilayer perceptron network (MLP) is fed with the multiscale filter responses of the GMF and Gabor filters methods in order to obtain an improved detection response using information from the spatial and frequency domains of the X-ray image.
In addition, to select the optimal MLP configuration for detecting coronary arteries, a hypothesis test for stochastic domination between multiple possible architectures has been performed. To evaluate the detection step of the proposed method, eleven methods for blood vessel enhancement are compared with the proposed method in terms of the area under the ROC curve. Finally, in the segmentation step of the proposed method, twelve specialized blood vessel segmentation methods are used for comparative analysis in terms of six binary classification metrics.
The organization of this paper is the following. In Section 2, the basis of the Gaussian matched filters, Gabor filters and Artificial Neural Networks, along with the proposed method including multiscale analysis, are described in detail. Computational experiments including numerical and qualitative results are discussed in Section 3, and finally, the concluding remarks can be found in Section 4.

Methods
In this section, the single-scale and multiscale versions of the Gaussian matched filters and Gabor filters are described first in Sections 2.1 and 2.2. Later, an introduction to Artificial Neural Networks that are used by the proposed method to perform the multiscale analysis of the coronary angiograms is given in Section 2.3. Additionally, the proposed method and the metrics used to evaluate its detection and segmentation performance can be found in Sections 2.4 and 2.5, respectively.

Multiscale Gaussian Matched Filters
Gaussian matched filters were proposed for vessel enhancement in medical images of the retina [6]. The original GMF method follows the assumption that a vessel cross-section can be approximated by a single-scaled Gaussian curve. Moreover, this method considers that blood vessels can be studied as individual segments that share a similar orientation and diameter, and consequently can be detected with matching templates. In the detection process of the GMF, the vessel structures are enhanced by the convolution of the original image with a gray-scale template defined under a Gaussian distribution. After filtering, the vessel detection response shows an increased contrast between vessel-like structures and the image background. The main idea behind the GMF is supported by the single-scale Gaussian template, which is defined as follows: where x ∈ [−T/2, T/2]. The parameter T delimits the Gaussian curve trails and produce a filter of finite width. The value of L corresponds to the length of a vessel segment which shares the same cross-section profile and orientation. The average vessel diameter to be detected is specified by σ as the standard deviation of the Gaussian distribution. In addition, because vessel structures can be found in multiple orientations, the filter template is transformed in order to detect segments at different directions with the rotation matrix as follows: Commonly, K directions are considered by defining θ to be evenly spaced in the range [−π/2, π/2].
Observing the similarity to blood vessels in the retina, the GMF method has been applied to the detection of coronary arteries in X-ray angiograms [11]. However, due to the original GMF works with a single scale of σ for the whole image, only vessels with the same diameter are correctly enhanced. This fact represents a disadvantage because coronary arteries show a natural thickness variation over their span, and ramifications of vessels with minor diameters are commonly present along the coronary artery tree. To overcome this drawback, a multiscale approach was proposed by Cruz et al. [17] to increase the vessel detection capability of the GMF method.
In order to illustrate the detection process carried out by the GMF, in Figure 1, Gaussian templates with different parameter values are presented. The variation between single-scale and multiscale filtering is the number of templates defined at different scales of σ that produce the directional filter bank, which is shown in Figure 1d,g. Each template is intended to enhance vessels of different diameters, providing a multiscale response for the same angiogram, which can be appreciated in Figure 1f,i. The multiscale analysis of the GMF presents a trade-off between the gain in vessel detection efficiency and the increment of computational cost because the high number of convolution operations.

Multiscale Gabor Filters
The Gabor filter is a Gaussian curve modulated by a sinusoid function which presents an anisotropic property capable of detecting directional features [18]. Soares et al. [12] proposed the use of Gabor filters for the detection of blood vessels of the retina. The Gabor wavelet is employed to detect the edges of blood vessels in the space domain by identifying singularities in the frequency domain. This method was originally proposed to work in a single scale, which increase the contrast of features that share a similar diameter and orientation. The real part of the single-scale Gabor filter can be defined as follows: where σ x and σ y control the standard deviation of a two dimensional Gaussian curve in the x and y axis, respectively. The value of f 0 determines the frequency at which the Gaussian curve is modulated. The filter design was simplified by Rangayyan et al. [13,14] in the sense of using only two parameters to define the Gabor template. The simplification states that the spreads σ x and σ y can be defined as σ x = τ/2 √ 2 ln 2 and σ y = lσ y . In addition, the frequency f 0 is set to modulate the curve under a period of τ. The reformulated Gabor filter depends only on two parameters: τ and l. The diameter of the features to be detected by the Gabor filter is controlled by τ, while l governs the elongation of a segment that shares a similar orientation. For the detection of vessel-like structures in multiple directions, the Gabor filter is rotated into κ evenly spaced angles in the range of θ ∈ [−π/2, π/2]. Under the same principle that the single-scale Gabor filter was proposed to detect blood vessels in the retina, this method has been successfully applied to the detection of coronary arteries in X-ray angiograms [16]. Although single-scale Gabor filters are capable of detecting the main vessels in coronary angiograms, the response obtained tends to ignore vessels of diameters lower than τ. To avoid this problem, a multiscale approach was proposed by Rangayyan et al. [14] for the detection of blood vessels in the retina, and later Cruz et al. [16] applied this method to segment coronary arteries.
In Figure 2, two templates defined with different width of the Gabor filter are presented in order to illustrate the detection process carried out by this method. The main difference between single-scale and multiscale Gabor filters is the number of templates at different width employed for the directional filter bank construction, which is illustrated in Figure 2d,g. The variation of the filter width allows to enhance blood vessels of different diameters in the coronary angiogram, as shown in Figure 2f,i. The multiscale analysis with Gabor filters provides an improved vessel detection quality at expenses of increasing the computational time, which is mostly because of the number of image transformations between the spatial and frequency domains.

Artificial Neural Networks
An Artificial Neural Network (ANN) is a machine learning technique consisting on the assembly of multiple processing units called neurons. ANN are capable of replicate the behavior of unknown functions as universal function approximators [19]. A common function to replicate by an ANN is the judgment of an expert during image labeling, which could be considered as a classification problem.
In classification problems, the prediction process of an artificial neural network could be interpreted as a sequence of projections into new feature spaces, where the linear separation of patterns represents a simpler task. In the projection sequence, each neuron performs a weighted sum of its corresponding inputs; in addition, a bias term is frequently employed to center the neuron response around a decision boundary.
In order to simulate the biological activation of the neuron, a function is applied to the resulting neuron response. The asymptotic behavior of the activation of a neuron is commonly approximated by a nonlinear step function, such as hyperbolic tangent and sigmoid functions. Moreover, differentiable activation functions are selected when gradient-based methods are used for the optimization of the artificial neural network parameters.
The operation executed by a single neuron can be expressed as where x is an input pattern, w i is the associated weight corresponding to the i-th variable of the pattern, b is the bias corresponding to the neuron, n is the number of input features, and g is a nonlinear activation function.
The classification performance of an ANN can be increased by adding more neurons to the network; however, the selection of the proper configuration represents a computational challenge given the large possible neuron arrangements. Layered configurations are a particular arrangement for the ANN neurons, where the computing units apply the same transfer function and share the input data. Moreover, layers of neurons can be connected to other layers, being the Multilayer Perceptron (MLP) illustrated in Figure 3, a frequently used model. The design of a MLP network to carry out a specific task relies on the selection of the appropriate number of hidden layers, and the number of neurons that compose them. Additionally, a training stage is essential in order to find the most appropriate weight values for each neuron of the network to obtain high approximation performance.

Proposed Method
To enhance and segment coronary arteries in X-ray angiograms, the proposed method consists of four different steps, which are illustrated in Figure 4.

Gabor filters
Vessel segmentation

Detection response Original angiogram
Ground-truth Binary classification Segmentation response Step 2 Step 4 The first step enhance vessel-like structures in the original angiogram by using multiscale Gaussian matched filters. The parameters of the GMF are directly related to spatial features found in coronary angiograms, such as length, diameter, and orientation of vessel segments. In order to detect any coronary artery efficiently, the parameters of the GMF have to be selected from a search space properly defined. The combination of the parameters σ and T control the coronary arteries diameter that can be detected, for this reason the value of T is selected from the range [8, 9, . . . , 15] pixels, and σ in the range [1.5, 1.6, . . . , 2.5]. The different vessel segments frequently present lengths in the range [8, 9, . . . , 15] pixels, which is governed in the GMF by the parameter L. The orientation of the vessel segments is controlled by the parameter K, which is defined to cover the angular range [−π/2, π/2] and is commonly fixed to sample 12 evenly spaced orientations.
To perform the multiscale analysis in the spatial domain, only σ is varied within its corresponding search space in order to enhance blood vessels of different diameters, while T, L, and K are kept fixed. The multiscale GMF response is conformed by 11 detection response images, corresponding to the different scales of σ, and can be considered as a three-dimensional image of size 300 × 300 × 11 pixels.
The second step consists of the analysis of the coronary angiograms using multiscale Gabor filters. For the simplified parameterization of Gabor filters, the parameter l governs the elongation of the vessel-like structures detector in a search space of [1, 1.1, 1.2, . . . , 20]. The parameter τ controls the width of the detector in a range of [1, 3, . . . , 19]. The value of κ determines the number of directional filter templates, which is commonly fixed to 45 evenly spaced orientations in the range [−π/2, π/2]. The multiscale analysis of the coronary angiograms in the frequency domain is performed by varying the value of τ within its corresponding search space, in order to enhance vessel-like structures of different diameters, while the parameters l and κ are kept unchanged. The resulting multiscale Gabor filter response is conformed by 10 different images, one per each scale of τ, that can be considered as a three-dimensional image of size 300 × 300 × 10 pixels.
In the third step, a four-layered perceptron network provides the coronary artery detection. The multiscale filter responses obtained from the first and second steps are concatenated to form a three-dimensional image of size 300 × 300 × 21. The MLP performs the vessel detection by using each 1 × 1 × 21 pixel-wise inputs from the concatenated multiscale detection response, forming the input layer of the artificial neural network with 21 neurons. The second and third layers of the MLP are hidden layers, where the number of neurons in each is selected by following a statistical hypothesis test for stochastic domination between a set of candidate architectures. The sigmoid function is selected as activation function for both hidden layers to provide a nonlinear transformation of the neuron responses. The set of all possible architectures is large; for that reason, in this work the search space has been limited to i, j = 1, 2, . . . , 10, where i and j are the number of neurons in the first and second hidden layers, respectively. The fourth layer is the output layer, which provides the final coronary artery detection response in an unbounded continuous scale by using a linear activation function. The detailed procedure followed for the selection of the best MLP architecture, along with the selection of the parameters employed by the multiscale Gaussian matched filters and Gabor filters, is presented in Section 3.2.
Finally, the fourth step consists of the binary classification of the pixels in the detection response obtained from the third step. Each pixel is discriminated between two classes: vessel features, and image background. In literature, multiple classification methods exist for the binarization of images that can be employed over the arbitrary values of a filter response image. However, none of them has been able to obtain an efficient use of the high detection performance presented in the previous steps of the proposed method. This reason leads to the definition of a fixed threshold value in the detection response scale that separates the image pixels into the two different classes. The selection process of the binary classification method can be found in Section 3.4.

Performance Evaluation
The area under the Receiver Operating Characteristic (ROC) curve has been widely used to evaluate the capacity that a method has to detect features in medical images. The ROC curve computes the True Positive Fraction (TPF) and False Positive Fraction (FPF) by applying a sliding threshold to a given filter response, which is in general a gray-scale image, and comparing the number of resulting binary images with the corresponding ground-truth angiogram. The area under the ROC curve (A z ) can be computed using the trapezoidal rule, and it is in the range [0, 1]. Being A z = 1 a perfect detection of the features of interest, and A z = 0 the opposite meaning.
To evaluate the performance of a binary segmentation method, the accuracy metric has been the most commonly used in literature. However, because accuracy metric is low sensitive to imbalanced classes, in this work the metrics of mean class accuracy [20], sensitivity, specificity, positive predictive value, and Dice similarity coefficient are also employed to evaluate the segmentation performance of the proposed method. The six metrics are in the range [0, 1]. An accuracy of 1 means a perfect extraction of the features of interest, and 0 otherwise. The mean class accuracy provides a metric that is not biased by the difference of sizes between the classes of blood vessel pixels and image background that are present in a coronary angiogram. The percentage of feature pixels correctly classified, and background pixels correctly classified, are measured by the sensitivity and specificity metrics, respectively. The positive predictive value corresponds to the percentage of actual feature pixels classified by the feature of interest. The Dice coefficient measures the similarity between the features of interest extracted by a segmentation method and the ground-truth, ignoring the background pixels correctly classified.
The computation of the six metrics for binary classification performance require to obtain the total number of feature pixels correctly classified as features (TP), feature pixels incorrectly classified as background (FP), background pixels correctly classified as background (TN), and background pixels incorrectly classified as features (FN). The classification accuracy is then defined as following: Sensitivity and specificity are defined respectively as follows: When only two classes are present in an image, mean class accuracy can be computed from sensitivity and specificity as follows: The positive predictive value is computed as follows: The Dice similarity coefficient is computed as follows: In the context of artificial neural networks, the approximation performance is generally measured using the Mean Squared Error (MSE). In addition, a scheme of cross validation with random sub-sampling is integrated to the training stage in order to prevent the parameters overfitting. The cross validation scheme allows to employ images from the training set repeatedly without invalidating the results obtained during the training stage [21]. This scheme separates the training set into two subsets, one to train the network, and the second one to validate the generalization capability of the neural network. In this work, for each training trial, the validation set is conformed by a random sample of 25% of the angiograms of the training set, while the remaining X-ray images of this set are used to train the artificial neural network.

Experimental Results
Since the proposed method consists of detection and segmentation steps, the results of each step are presented in separated sections. In the detection step, the methodology employed to select the most appropriate configuration for the MLP is presented. In addition, a comparison with eleven vessel enhancement methods is carried out by using the area under the ROC curve (A z ) as comparative metric. On the segmentation step, the selection process to define a specific method that complements the proposed detection step as an integral vessel extraction approach is introduced. Moreover, this section presents the comparative analysis of the proposed coronary arteries segmentation method with respect to twelve state-of-the-art blood vessel extraction methods in terms of six binary classification metrics.
The computational experiments were implemented on a computer with an Intel Core i7, 4 GB of RAM, and 3.50 GHz processor using the version 2016b of Matlab.

Database of Coronary Angiograms
The database used in this work consists of 130 X-ray coronary angiograms, and their corresponding ground-truth image outlined by an expert cardiologist. Each angiogram is a 300 × 300 pixels, gray-scale image in pgm (Portable gray map) format. The Cardiology Department of the Mexican Social Security Institute, UMAE T1-León has provided the whole image database, and the ethics approval for its use in the present research, under reference R − 2019 − 1001 − 078.
To evaluate the detection and segmentation performance of the comparative methods presented in this paper, the database of X-ray images was divided into two subsets: training and testing. The training set consists of 100 images, and the remaining 30 angiograms form the test set.
The complete database of X-ray images, along with the ground-truth images employed in this research will be available via web page (http://personal.cimat.mx:8181/~ivan.cruz/DB_Angiograms. html). This database represents the first set of angiograms and ground-truth images made accessible to the scientific community for research and comparison purposes.

MLP Architecture Selection
In order to determine the optimal Multilayer perceptron architecture, the number of layers plays an essential role. In general, the input and output layers are directly related with the nature of the problem to be solved. The main challenge is given to establish the number of hidden layers along with the number of hidden neurons for each of them. To establish the number of hidden layers in the proposed Multilayer perceptron, the range of hidden layers h = [1, 2, 3, 4] was explored. MLP architectures with more than two hidden layers present low performance in terms of the area under ROC curve and high computational time using the present training set. Consequently, MLP architectures with one and two hidden layers were selected to be analyzed in detail. According to the statistical analysis presented below, the optimal performance of the MLP architecture can be obtained using two hidden layers; therefore, a Multilayer perceptron with a four layers architecture was proposed in the present paper. In addition, the number of hidden neurons for each hidden layer represents one of the main contributions of the present paper, since it has been determined by using a stochastic dominance approach, which is also discussed below.
Since the Gaussian matched filter is governed by 4 parameters (K, L, T, σ), and Gabor filters are governed by 3 parameters (τ, l, κ), these seven parameters have to be determined previously to be used in the artificial neural network. The parameter that govern the directional filters is fixed as K = 12 for the GMF, and κ = 45 for the Gabor filters, because the number of oriented filters is independent from the training set of images, and they are the most common values in literature. Moreover, to select the optimal set of values for the remaining 5 parameters (T, L, and σ for Gaussian matched filter and l, and τ for Gabor filters), the training set of 100 X-ray images, and the area under the ROC curve were used to evaluate the filter responses (gray-scale images) for each parameter combination applying an exhaustive full-search procedure over the predefined search space (See Section 2.4).
In this procedure, each possible combination of the Gaussian matched filter parameters are taken from the search space and used to perform a single-scale filtering for each image of the training set, which provides a detection response with the same dimension as the input image (300 × 300 pixels). The detection responses obtained with each parameter combination are evaluated in terms of the area under the ROC curve with respect to the ground-truth images. The optimal parameter combination is the one that provides the highest A z value from the single-scale detection responses.
The optimal values of L and T for the single-scale Gaussian filtering are used later with each scale of σ when the multiscale analysis is performed in the spatial domain. Moreover, the same procedure was performed in order to determine the optimal l (elongation) parameter for the Gabor filters. In this case, the optimal elongation parameter found for the single-scale filtering is combined with each scale of τ during the multiscale analysis in the frequency domain.
After the exhaustive search, the most appropriate set of values for the parameters of the GMF are defined to be T = 15 and L = 13 pixels, and for the Gabor filters the elongation is defined to be l = 3. The multiscale analysis was performed in the range σ = [1.5, 1.6, . . . , 2.5] for Gaussian matched filters (11 scales) and using the range τ = [1, 3, 5, . . . , 19] for Gabor filters (10 scales).
On the other hand, the selection of the MLP architecture requires a procedure based on a statistical hypothesis test. The main reason is the stochastic nature of the training process of an artificial neural network, which is introduced by the random initialization of the weights for each neuron, and the sampling process used to form the validation and training subsets. For the training stage of each possible architecture of the MLP, the back-propagation with the Levenberg-Marquardt optimization algorithm [22] are employed for the weights fitting, using the training set of 100 angiograms. After the training stage, the generalization performance of the artificial neural network is evaluated using the test set of 30 angiograms, and the area under the ROC curve as metric. The ROC analysis is used since the final detection response of the proposed method is provided in an unbounded continuous scale, which can be considered as a multiple level classifier for the pixels of the coronary angiograms.
To illustrate the stochastic behavior of the generalization performance of the tested MLP architectures, Figure 5 presents the distribution of the performance achieved using 30 trials by each artificial neural network configuration. In addition, the Shapiro-Wilk test is applied to determine whether the behavior of the generalization performance of each tested configuration follows a normal distribution. The results of the Shapiro-Wilk test are presented in Table 1, where the analysis suggests that normality cannot be assumed for all MLP configurations with a significance level of 5%. Consequently, a non-parametric comparison approach is required to select the most appropriate MLP architecture.  In literature, the Wilcoxon-Mann-Whitney test is the most common non-parametric technique used to compare the performance of two methods applied to independent samples. This test is used to determine the most appropriate MLP configuration for coronary arteries detection in X-ray angiograms by pairwise comparison of all candidate architectures using the testing set of images. The hypothesis test determines whether the detection performance between two configurations are similarly distributed, by using a significance level of 5%. On the other hand, when the hypothesis of similarly distributed detection performance is rejected, the sum of ranks is employed to determine whether the left or right hand configuration is stochastically dominant. For the best 10 MLP configurations, Table 2 illustrates the hypothesis testing process by reporting the p-value and the corresponding left and right hand sum of ranks, denoted respectively by − and +. The stochastic dominance comparison is summarized in Figure 6, where a high detection performance can be observed for configurations with more than 4 neurons in each hidden layer. In this figure, the configuration [9,9] is highlighted with a dotted line, corresponding to the architecture that stochastically dominates most of the other configurations of the comparative analysis. Figure 6. Distribution of the number of architectures dominated by each configuration according to the ranking procedure using the testing set of images. The configuration [9,9] that dominates most of the other architectures is highlighted with a dotted line.
The stochastic dominance comparison results are ranked according to the number of other architectures stochastically dominated, in order to determine the most appropriate configuration to be used in the proposed method. Table 3 shows the comparative results of the top 10 architectures from the 100 comparative architectures tested. According to comparative analysis, the configuration that stochastically dominates the most of other architectures is the [9,9] architecture, being superior to other 90 configurations. In addition, the [9,9] architecture obtains a detection performance stochastically similar to only other 9 configurations, and it is not stochastically dominated by any other configuration. Because the MLP [9,9] configuration presents the highest stochastic dominance results, it was selected as the most appropriate architecture to provide the final detection response of the proposed method. Table 2. P-values, left-hand (−), and right-hand (+) sum of ranks from the Wilcoxon-Mann-Whitney test for stochastic difference of the best 10 MLP configurations. The configurations are presented as [i, j], where i is the number of neurons in the first hidden layer, and j the number of neurons in the second hidden layer. [8,10] [9, 5] [9,8] [9,9] [9, 10] [10,6] [10, 8] [10,9] [10, 10] [8,8]

Detection Results
To evaluate the performance of the detection step of the proposed method, the comparative analysis with respect blood vessel enhancement methods from the state of the art is reported in this section. The proposed method is compared with four artificial neural network based methods, and seven state-of-the-art methods for blood vessel enhancement, using the test set of 30 images. The detection performance of each blood vessel enhancement method is evaluated in terms of the area under the ROC curve.
The methods based on neural networks were trained using the training set of 100 images, and were selected from several candidates that were previously tested, and according to their relevant high blood vessel detection performance. The multiscale Gaussian matched filter [17], which employs a multilayer perceptron network with 3 neurons in a first hidden layer, and 8 neurons in a second hidden layer. This method obtains a vessel detection performance of A z = 0.9362, using a database of 100 angiograms, and the set of parameters l = 13, T = 15, K = 12 and σ = [1.5, 1.6, . . . , 2.5] for the multiscale GMF. The multiscale Gabor filter [16], that uses a multilayer perceptron with two hidden layers with 5 neurons each. This coronary artery detection method achieves a performance of A z = 0.9520 using a database of 80 angiograms, whereas the parameters for the multiscale Gabor filters are set as l = 3, τ = [2, 3, . . . , 20] and κ = 45. The following two methods are convolutional neural networks (CNN) that were originally proposed to perform object recognition in images; however, in this work those architectures are adapted to perform the detection of blood vessels in coronary angiograms. The architecture proposed by Lecun et al. [23] which is applied to the multiscale responses of the Gaussian matched filters and Gabor filters, using the same parameters of the proposed method for both multiscale vessel enhancement methods. Finally, the architecture proposed by Simonyan and Zisserman [24] is employed using 32 × 32 pixel sized patches from the original angiogram without additional processing.
The following seven blood vessel enhancement methods are the most referenced in literature. The method of Cruz et al. [11] employs the single-scale GMF, where the filter parameters are selected with an algorithm inspired by natural evolution. Cinsdikici and Aydin [7] proposed to use the single-scale GMF extending the number of orientations to K = 18, keeping the other parameters in the original values proposed by Chaudhuri et al. [6]. Nguyen et al. [25] addressed the detection of blood vessels by using a line detector and varying its scale in a range of 1 to 15 pixels. The method of Al-Rawi and Karajeh [10] is based on the single-scale GMF, where the parameters T and L that govern the filter are selected by a genetic algorithm from the evolutionary computation family, preserving the orientations in K = 12. The mathematical morphology based method of Eiho and Qian [1] applies the Top-hat operator with a disk-shaped structuring element at multiple scales over the original image. The single-scale Gabor filter proposed by Rangayyan et al. [13] defines the parameters with a full-search, resulting on a filter with parameters τ = 8, l = 2.9, and κ = 180 orientations. The method of Frangi et al. [3] is based on a vesselness measure computed from the Hessian matrix of the original angiogram using multiple scales in the range [1, 1.25, 1.5, . . . , 21] in order to detect vessels of different diameters.
The comparative analysis between the proposed method, and the eleven blood vessel enhancement methods described above, is illustrated in Figure 7 using the test set of 30 images of the present database. The comparative analysis reveals that the proposed method presents the highest vessel detection performance with A z = 0.9775, compared with the eleven methods from the state of the art. To visualize the vessel detection results, Figure 8 presents the responses obtained by applying the different methods over a subset of angiograms of the testing set.

Segmentation Results
This section presents the process to determine the most suitable binary classification technique to be used in the fourth step of the proposed method. Later, the evaluation of the segmentation performance of the proposed method is presented in a comparative analysis approach with other blood-vessel segmentation methods from the state of the art. For the evaluation of the blood vessels classification efficiency, the accuracy metric is employed.
To select the most suitable technique for the segmentation step of the proposed method, fourteen binary classification tools are compared using the detection response of the multiscale analysis performed in the previous steps. Thirteen methods of the comparative set are automatic binary classification techniques frequently used in vessel-segmentation literature. The remaining method is a soft computing approach which selects the threshold value from the detection response scale that maximizes the accuracy metric. The soft computing procedure followed to define the most appropriate threshold value is illustrated in Figure 9, where the best value found is 0.87 in the detection response scale using the training set of 100 images and a full search strategy.    In addition, Table 4 presents the comparative analysis of the fourteen automatic binary classification tools in terms of six different binary evaluation metrics, using the test set of 30 angiograms. The analysis is performed by applying each of the compared segmentation tools to the multiscale response of the proposed method, which is illustrated in Figure 10. After visual examination of Figure 10, it can be noticed that the soft computing approach provides the lowest incorrect classification rate of background pixels as vessel features, reflected in its high specificity. On the other hand, most of the other classification tools have presented high sensitivity at the cost of a considerable high number of false positives that results in a low positive predictive value, whereas the fixed threshold performs the correct classification of vessel pixels with a significant rate of positive predictive value. Moreover, the Dice coefficient achieved by the fixed threshold is significantly higher than the obtained by the other binary classification methods, being the soft computing approach the one that provides the most similar vessel segmentation to the ground-truth. Because of its high overall segmentation performance presented in the comparative analysis, the best classification technique to segment the multiscale detection response is the binarization under a threshold value of 0.87 obtained with a soft computing procedure.  The following rows are the segmentation results of (c) the fixed threshold (0.87) method, (d) the method of Otsu [26], (e) the method of Ridler and Calvard [27], (f) the method of Niblack [28], (g) the Moment-preserving [29] method, (h) the RATS [30] method, (i) the Vessel repair [5] method, (j) the method of Sauvola [31], (k) the Entropy maximization [32] method, (l) the Background unification [1] method, (m) the Degree-based [33] method, (n) the method of White [34], (o) the Histogram concavity [35] method, and (p) the Local entropy [36] method.
Moreover, in order to evaluate the results of the proposed method, twelve specialized vessel segmentation methods of the state of the art have been selected to perform a comparative analysis. From the set of comparative methods, the method introduced by Eiho and Qian [1] uses the Top-hat operator for coronary arteries enhancement in angiograms, which combined with a background unification algorithm, extracts vessel features from coronary angiograms. Qian et al. [2] used the background unification algorithm on the multiscale response of the Top-hat operator, applied at three different scales over the original angiogram. Chanwimaluang and Fan [37] addressed the vessel segmentation from the single-scale GMF response. The binary segmentation is obtained by selecting the threshold value that maximizes the local entropy computed from the co-occurrence matrix. Kang et al. [8] applies two different enhancement methods to the original image: the Top-hat operator, and the single-scale GMF. The segmentation result is obtained with a threshold level that maximizes the entropy of the gray-scale intensity of each detection response. The method proposed by Kang et al. [38] introduces a segmentation procedure based on a similarity degree between pixels of the gray-scale detection response. The approach of Nguyen et al. [25] segments vessel pixels from the multiline detector response by using a fixed threshold value of 0.56. The value is selected as the binarization level that maximizes the classification accuracy using a training set of images. The method of Li et al. [5] defines a vessel probability map from a coronary angiogram using a Hessian matrix based method. The vessel segmentation is performed by a region growing algorithm and a subsequent vessel detail repairing procedure. The inter-class variance maximization method of Otsu [26] has been an efficient separation tool employed by the following methods. The method of Cervantes et al. [15] employs the single-scale Gabor filter, which parameters are selected using a nature inspired algorithm for numerical optimization. The U-Net [39] is a convolutional neural network that has been widely used for biomedical image segmentation tasks in recent years. For the comparative analysis, the U-Net [39] was trained using the training set of 100 images, by taking patches of size 32 × 32 pixels from the zero-padded angiograms as inputs, and using the cross entropy function as loss metric. Finally, the method proposed by Cruz et al. [17] introduced the multiscale GMF that improves the detection performance of its single-scale version.
The result of the comparative analysis of the twelve specialized vessel segmentation methods, and the proposed method, is presented in Table 5. The test set of 30 angiograms is used to evaluate the segmentation methods in terms of binary classification accuracy. According to Table 5, the proposed method obtains the highest segmentation performance with respect to the twelve methods from the state of the art by reaching a classification accuracy of 0.9698. The resulting segmentation responses of the comparative segmentation methods are illustrated in Figure 11. By visual inspection of Figure 11, it can be noticed that the method of GMF/Local Entropy [37] presents responses with a high rate of false positive pixels, obtaining the lowest segmentation accuracy of the comparison, this problem is also reflected in the low positive predictive value. The method of GMF/Degree-based [38] can obtain a relative high number of false negative pixels, affecting the segmentation accuracy and sensitivity of the method. The method based on single-scale Gabor filters [15] presents a frequent problem on the resulting image edges reducing the sensitivity of the method. The high false positive rate in those areas is associated to the strong detection intensities obtained by using only the average vessel width as single scale. In the method based on multiscale Gaussian templates [17], the accuracy and positive predictive value are affected by the incorrect classification of the pixels that are closer to the main coronary artery. This problem is because of the strong detection responses obtained by the multiscale Gaussian matched filter for small tubular structures in the background. The methods Hessian Matrix-based/Vessel Repair [5], and Multiline Detector/Threshold (0.56) [25] achieve similar classification accuracy. However, their segmentation performance is affected by the presence of irregular edges, holes, and broken vessels, being the method of Hessian Matrix-based/Vessel Repair [5] the most affected between both, reducing its segmentation sensitivity. The methods Top-hat/Background Unification [1], and Multiscale Top-hat/Background Unification [2] achieve high classification accuracy, only limited by the frequent presence of irregular borders and broken vessels in the segmentation result which reduces their positive predictive value and Dice similarity coefficient with the ground-truth. The U-Net (Ronneberger et al.) [39] method presents the highest specificity and positive predictive value; however, the segmentation result show a high rate of false negative pixels in several images, resulting in broken vessels and a reduction of its overall vessel extraction performance. The method of GMF/Entropy Maximization [8] presents high segmentation results in overall, being affected only by the presence of holes, broken vessels, and false positives in image edges, reducing its positive predictive value. On the other hand, the experimental results obtained by the proposed method achieved the highest segmentation accuracy, specificity, positive predictive value, and Dice similarity coefficient, using the testing set of the present database. Because of the suitable performance of the proposed method, the segmented coronary arteries present a low rate of false positives along with a low presence of holes and isolated vessels.
The comparison of the average execution time of the proposed method, and twelve vessel segmentation methods from the state of the art, is presented in Table 6. The methods of GMF/Entropy Maximization [8] and GMF/Degree-based [38] present the lowest computational times, followed by the methods Top-hat/Background Unification [1] and GMF/Local Entropy [37], which low computational time can be attributed to the single-scale analysis performed in the detection step. The methods of Multiscale Top-hat/Background Unification [2] and Multiline Detector/Threshold (0.56) [25] employ a multiscale analysis during the detection step, which increases the execution time with respect to single-scale approaches. The Single-scale Gabor filters/Otsu [15] involve transformations between the space and frequency domain that affects the overall execution time. The Multiscale GMF/Otsu [17] method performs a multiscale analysis in the spatial domain, which demand a greater execution time than other single-scale GMF methods. The Hessian Matrix-based/Vessel Repair [5] method employs a multiscale approach during the detection step, and a vessel repair process in the segmentation step. This method requires a greater execution time than other multiscale methods and is reflected on its average execution time. The proposed method presents the highest execution time in the comparative as a trade-off for its high segmentation accuracy. The computational time demanded by the proposed method is directly related to the number of multiple scales of the Gaussian filter applied in the spatial domain and the Gabor filters performed in the frequency domain. m n o Figure 11. (a) Ten angiograms from the testing set of 30 images, and (b) their respective ground-truth. The following rows are the segmentation response obtained by (c) the proposed method, (d) the the GMF/Entropy Maximization [8] method, (e) the U-Net (Ronneberger et al.) [39] method, (f) the Multiscale Top-hat/Background Unification [2] method, (g) the Top-hat/Background Unification [1] method, (h) the Multiline Detector/Threshold (0.56) [25] method, (i) the Hessian Matrix-based/Vessel Repair [5] method, (j) the Multiscale GMF/Otsu [17] method, (k) the Single-scale Gabor filters/Otsu [15] method, (l) the CNN (Simonyan and Zisserman) [24] method, (m) the GMF/Degree-based [38] method, (n) the CNN (Lecun et al.) [23], and (o) the GMF/Local Entropy [37] method. Table 6. Average execution time of the proposed method and the twelve vessel segmentation methods from the state of the art presented in Table 5, using the testing set of 30 images.
Finally, it is important to point out that the proposed method for the automatic segmentation of coronary arteries achieved the most accurate results using the present database and different evaluation metrics for each of the vessel detection and segmentation steps. The high segmentation accuracy can be attributed to the multiscale analysis performed in both spatial and frequency domains along with the statistically determined artificial neural network. In addition, given the efficiency obtained by the proposed method and the competitive computational time, it could be integrated into a computer-aided diagnosis system to support cardiologist during the decision-making process.

Conclusions
In this paper, a new method for automatic segmentation of coronary arteries, based on multiscale analysis of X-ray angiograms, has been introduced. The proposed method employs the multiscale Gaussian matched filters, applied in the image domain, and the multiscale Gabor filter, applied in the frequency domain, in order to improve the contrast between vessel features and background, suppress noise in the image, and homogenize the illumination of the angiogram. Both multiscale responses are used as input layer in a four layered perceptron network. A statistical analysis has shown that the Multilayer perceptron with 9 neurons in each of the two hidden layers is the appropriate using a comparative set of configurations. Moreover, the proposed method based on an artificial neural network with multiscale analysis has performed better than eleven vessel detection methods from the state of the art by reaching an A z = 0.9775 when applied to 30 test images. The use of a fixed threshold value of 0.87, over the multiscale response of the proposed method, has proved to be the most appropriate from fourteen segmentation techniques. Finally, according to comparative experiments with twelve state-of-the-art vessel segmentation methods, the proposed method has provided the most accurate vessel segmentation of the test set of 30 angiograms with a classification accuracy of 0.9698, and Dice coefficient of 0.6857. In addition to the experimental results, the high performance and robustness in the detection and segmentation steps have shown that the proposed method could be a potential method to be integrated in a computer decision support system to aid diagnosis in medical practice.

Conflicts of Interest:
The authors declare no conflict of interest.