A FAST AND ACCURATE METHOD FOR AUTOMATIC CORONARY ARTERIAL TREE EXTRACTION IN ANGIOGRAMS

Coronary arterial tree extraction in angiograms is an essential component of each cardiac image proces sing system. Once physicians decide to check up coronary arteries from x-ray angiograms, extraction must be done precisely, fast, automatically and including w hole arterial tree to help diagnosis or treatment d uring the cardiac surgical operation. This application is ver y helpful for the surgeon on deciding the target ve ssels prior to coronary artery bypass graft surgery. Some techniques and algorithms are proposed for extract ing coronary arteries in angiograms. However, most of t hem suffer from some disadvantages such as time complexity, low accuracy, extracting only parts of main arteries instead of the full coronary arterial tree, need manual segmentation, appearance of artifacts a nd so forth. This study presents a new method for extracting whole coronary arterial tree in angiogra phy images using Starlet wavelet transform. To this end, firstly we remove noise from raw angiograms and the n sharpen the coronary arteries. Then coronary arte ri l tree is extracted by applying a modified Starlet wa velet transform and afterwards the residual noises and artifacts are cleaned. For evaluation, we measure p roposed method performance on our created data set from 4932 Left Coronary Artery (LCA) and Right Coro nary Artery (RCA) angiograms and compared with some state-of-the-art approaches. The proposed meth od s ows much higher accuracy 96% for LCA and 97% for RCA, higher sensitivity 86% for LCA and 89% for RCA, higher specificity 98% for LCA and 99% for RCA and also higher precision 87% for LCA and 9 3% for RCA angiograms.


INTRODUCTION
Coronary arteries are one of the vital parts of heart which mostly rely on the heart surface and carry the blood to the heart muscle. Coronary arteries are categorized into two main types: Right Coronary Artery (RCA) and Left Coronary Artery (LCA). RCA stems from the right aortic sinus and LCA stems from the left aortic sinus. One of the heart disorders is Coronary Artery Disease (CAD) which cause of death nowadays in the most of the countries and it is related to blockage and narrowing the left or right coronary arteries. Therefore, it is important to choose the precise method to visualize the coronary arteries. There are several medical imaging techniques which can be used for this purpose; X-Ray Angiography, Cardiac CT Angiography (CTA), MR Angiography (MRA), Cardiac PET and SPECT. Among these techniques, X-ray Angiography, x-ray arteriography or briefly angiography, is the best one. Angiography can be well used to visualize the coronary Science Publications JCS arteries and diagnose the blockage and stenosis in realtime and most of the physicians prefer to use this modality instead of others to diagnose and treat the cardiac coronary artery diseases. Coronary artery angiography is done by injecting the radio opaque contrast agent into the coronary arteries and imaging using X-ray based techniques such as fluoroscopy. Series of radiographs of blood vessels is called an angiograph or angiogram.
Precise computer-assisted coronary artery analysis in angiograms is very consequential in later diagnosis and treatment with the cardiologists and cardiac surgeons. Coronary artery extraction is one of the crucial steps which should be done precisely. Because other processes, such as the stenosis measurement, 3D reconstruction of coronary arterial tree and blood flow analysis will be done based on the extraction's result. There are some difficulties in coronary artery extraction from angiography images. Firstly, angiograms have poor signal to noise ratio and the second one would be artifacts by existing other organs in angiography images, such as backbone, ribs or other cardiac parts. As presented by Tayebi et al. (2013), a number of methods have proposed to extract coronary arteries in angiograms. One of the defect of them would be they extract only small parts of the coronary arteries, instead of the whole coronary arterial tree. Also, some of them only extract centerline or counter of arteries and other ones extract main arteries and refuse to recognize sub arteries. Xu et al. (2007) proposed an algorithm for coronary artery centerline tracking in angiograms by using matched filter on the eigenvalues. In another study, Hernández-Vela et al. (2012) introduced an accurate coronary centerline extraction in angiograms. In these methods whole coronary arterial tree is not extracted and physicians cannot visualize and follow the sub-arteries, especially for diagnostic procedure in sub arteries.
Several semiautomatic algorithms have been proposed regarding to the coronary artery extraction. One of the disadvantages of proposed algorithms is involving user in defining seed points to locate the coronary arteries. Wang et al. (2005) also proposed a method for coronary artery extraction in angiograms. This method requires seed points to start extraction procedure and only extracts centerline and counter of main vessels instead of the whole coronary arteries. Furthermore it works only on Digital Subtracted Angiography (DSA) images while our proposed approach works on the plain angiography images. It is worth noting that extraction from plain angiograms is much more difficult than DSA. Considering the fact which the grayscale values between the coronary arteries and background in plain angiograms are mostly similar and also there are a lot of noise and reflection of x-ray in this type of image.
Another drawback would be, the proposed algorithms cannot extract coronary arteries precisely. Sun et al. (2005) extracted coronary arteries with fuzzy morphological method. However, their proposed algorithm have some drawbacks; some parts of the coronary arteries, especially in the thinner vessels do not have the same size as the original ones and the vessels are appeared thicker as well.
Other algorithms suffer from high computational complexity. Zhou et al. (2008) proposed an automatic approach for segmenting coronary artery in angiography images. One of the disadvantages of this study is the computational time. For Left Coronary Artery (LCA) extraction, the average computational time is about 32s. Another one is related to non-precise extracting arteries. In fact, when 3D structures are superposing or crossing together, especially with their similar diameters, this algorithm cannot correctly determine the corresponding 2D projection structures.
Some other algorithms have proposed segmentation for all types of the vessels. Using them on the coronary artery in angiograms will appear the artifacts in results. For example, Li et al. (2008) proposed a region based active contour model for vessel segmentation. After running their algorithm on angiograms, disorders are appeared; some parts of background are extracted as a coronary artery and some of the arteries are missed.
Mentioned limitations motivated us to propose a fast, accurate and automatic method for coronary arterial tree extraction which visualizes all parts of the coronary arteries precisely. The proposed method consists of three major phases: Preprocessing, processing and postprocessing. In the first one, after removing the noise from the raw angiograms with discrete wavelet transform, we sharpen the coronary arteries with Starlet wavelet transform. In the processing phase, coronary arterial tree is extracted by applying the modified Starlet wavelet transform and in the last phase, post-processing phase, the residual noises and artifacts will be cleaned.

THE PROPOSED METHOD
This research concentrates to propose a fast and accurate method for coronary artery extraction in angiograms automatically. To this end, three phases are defined: Preprocessing, processing and post-processing.

JCS
The details of our proposed method are described in the following sections.

Preprocessing Phase
Each type of x-ray angiography has a moving DICOM format result. Dealing with this type of imaging is too difficult. At first, the angiograms should be converted to 2D bitmap frames {f 1 , f 2 ,….f n } on each DICOM movie angle, such that varies from 40 to 70 frames. One optimum frame (f optimum ) must be selected from these frames and it is defined as a well x-ray injected frame that the whole of coronary arterial tree is contrasted with that. After choosing an appropriate frame, it should be converted from RGB format to grayscale; we call it as the original input image f for next steps. We define two main stages here for preprocessing. At first, we remove the noise from original input image by using Discrete Wavelet Transform (DWT) and then sharpen the arteries by applying Starlet Wavelet Transform (SWT).

Noise Removal with Discrete Wavelet Transform
Angiograms still suffer from a large amount of noise. Aim of noise removal process is to eliminate all the noises while keeping the quality of images. The quality hints to retain the arteries and particularly subarteries in angiograms here. When we use the traditional algorithms for removing noise in angiogram, such as smoothing methods included mean or median filters, we see some arteries disappear, especially thinner ones and sub-arteries. Another technique for removing noise is defined by converting images into the transform domain such as wavelet and then comparing the transform coefficients with a threshold. In this way, we hope the arteries' structure will be kept. So, we prefer to use Discrete Wavelet Transform (DWT) for removing noise from angiograms in this research.
The wavelet family comes from an original function which is called a mother wavelet. Let x is a real variable. We call function ψ(x) as a mother wavelet on the condition that it oscillates and averaging to zero This function is well localized; it means that quickly declines to 0 when x tends to ±∞. In addition, this function is normalized Equation 1: The Continuous Wavelet Transform (CWT) of a signal f(x) at the time and scale is defined as Equation 2: where, 〈f, ψ u,s 〉 denotes the scalar product of f and ψ u,s and also ψ*(x) represents the complex conjugate of ψ(x). W f (u, s) describes the fluctuations of the signal f(x) in position u and scale s. In fact, W f (u,s) is two dimensional representation ((u,s)-plane) plane) of a signal f(x) which is in 1D and so it is extremely redundant. To solve this problem, the Discrete Wavelet Transforms (DWT) is defined to reduce this redundancy by utilizing one of the countable infinite wavelets {ψ j,k }, where s = 2 j and u = To reduce redundancy, the family {ψ j,k } must establish an orthonormal basis of 2 2 ( )( ( ) L L ℝ ℝ is finite energy function). These conditions are related to multiresolution analysis MRA (or multiscale approximation) design method that proposed by Mallat (1989). The MRA of a signal is the representation of a signal in the finite sequential approximations. In this representation, each approximation is a smoothed type of the prior one and with the different resolutions. This concept means the signal f(x) will be divided into set spaces V j , in order to as j increases the detail information will be lost. So, in each resolution, the signal f(x) belong to space V j is divided into the subspace V j+1 orthogonal to the subspace W j+1 , so that the detailed information is retained in W j+1 (Ouahabi, 2013).
The MRA technique can be applied for some important applications, such as de-noising, compression, smoothing. The wavelet coefficients retain the difference between each two sequential resolutions and information in each two subspaces (Candes and Donoho, 2000) (Mallat, 1989). The DWT use the pair orthogonal high-pass ( ) g and low-pass ( ) h filters to decompose the signal f(x) into w j

JCS
(wavelet coefficient) and c j (scale coefficient) components in each resolution, respectively. In Fig. 1, implementation of MRA of signal f(x) in 1D case is illustrated. We can consider input signal f(x) is decomposed iteratively by high-pass and low-pass filters and then down-sampled by factor 2, ↓2. The down-sampled (or decimated), means to retain every other sample, instead of every one. The signal f(x) can be reconstructed by inverse Discrete Wavelet Transform (IDWT), as illustrated in Fig. 2.
Reconstruction includes up-sampled ↑2 by inserting zero between wavelet coefficients and scale coefficients to double their length and then followed by h and g filters, dual to h and g . Because for an image (2D), DWT is applied in each coordinates, the 2D signal (image) will be divided into four regions; (1) LL: Achieved by using low-pass filter on both coordinates, (2) LH and (3) HL: Achieved by using low-pass filter on one coordinate and high-pass filter on another one and (4) HH: Achieved by using high-pass filter on both coordinates. Similar in 1D signal, LL can be decomposed recursively to four regions in the next resolution (level) by using high-pass and low-pass filters. This scenario is illustrated in Fig. 3.
We can consider MRA for de-noising angiograms so that noisy signal is retained in HL, LH and HH regions and smoothed image is retained in the LL region, in every level. In this research, we used MRA with three levels. So, we can consider here region LL3 as an image version which probably does not contain any noise. Since the information in other three regions in each level will have useful information, such as arteries edge and boundary, we should keep them in reconstructing phase. For this reason, we define a threshold for these three regions to remove noise and reconstruct de-noise image version. There are two types of thresholding, soft and hard. We use soft thresholding here, which is defined by Equation 4: Which j denotes the levels (resolution) and k denotes the regions; value 1 for HL, value 2 for LH and value 3 for HH region in each level. Thresholds ( ) k j τ must be considered carefully for retaining useful information and processing noise removal. We calculated lots of thresholds and found for the type of angiograms which using Haar mother wavelet transform, the best Peak Signal to Noise Ratio (PSNR) value is obtained by using these thresholds: 4.3 for wavelet coefficients in level 3 and 4.7 for wavelet coefficients in level 2 and 4.9 for wavelet coefficients in level 1. These values are chosen based on this fact that if we consider higher values, the thinner arteries are disappeared and also PSNR is decreased; and if we consider lower values, noise still remains and PSNR decreased as well.

Sharpening Arteries with Starlet Wavelet Transform
In this stage, we intend to sharpen the edges of arteries and also cleaning the surrounding of arteries as background. The traditional methods for sharpening are not applicable here; because they sharpen not only the arteries, but also background of the angiogram and it ends to increase the noise in whole image and also may some artifacts appear in the segmentation phase.
Here, we proposed a new method for sharpening angiograms by focusing only on large angiogram's features; in this place the coronary arteries. We employed Starlet Wavelet Transform (SWT) on angiograms which causes to show the arteries and their borders more clear; especially when two separate arteries are closed together. SWT or Isotropic Undecimated Wavelet Transform (IUWT) is well known in the field of biology (Genovesio et al., 2006), astronomy (Starck and Murtagh, 2006) and nowadays in medical (Bankhead et al., 2012) applications. SWT decomposes an image l into w j , as a wavelet coefficient and c j as a scale coefficient in each iteration j (Starck et al., 2007). The sketch of SWT algorithm that proposed by Starck et al. (2010) is presented in Fig. 4.
As it is shown in this algorithm, for each scale in SWT, we have only one wavelet coefficient and not three as other wavelet types which introduced earlier.
It worth noting that, the value of J defines the maximum level of using the SWT, which possibly is different in various applications.
In this section, we modified the SWT algorithm to obtain sharped angiograms and in the next step we will aim to modify this algorithm again and use as a segmentation tool for coronary artery extraction. For the artery structure in angiograms, we have applied different filters with various levels and finally found out using two different filters which are We have chosen several levels because boundary of arteries will be shown by using smaller wavelet levels and surrounding background of the arteries will be cleaned by using the larger ones; and also we have chosen two filters because they have better result in angiograms in compare with other filters, such as [1,4,6,4,1]/16 which is mentioned in the original SWT algorithm. We will discuss more about choosing filters in result and discussion part, section 3. The proposed algorithm for sharpening coronary arteries in angiogram is illustrated in Fig. 5.
After running this algorithm, the sharpened image is defined by summation W and I. In Fig. 6, the output of this algorithm is shown. As it is demonstrated in this figure, coronary arteries are shown sharpened and also their surroundings are cleaned as the background. This process will be useful, when we want to extract coronary arteries. We will use this output as input for the next section, which will be segmentation part.

Processing Phase
SWT is a type of wavelet transform which is fast in computation and can be used for objects' segmentation. In this phase, we propose a new method for coronary artery extraction in angiograms based on applying SWT. Since in segmentation with SWT, intensity values of objects should be higher than background, at first we invert angiogram images. Then we use modified SWT for extracting coronary arteries. As mentioned before, using different filters and also levels in SWT, will give us different results. In this case, thinner arteries are shown in the smaller levels and thicker ones in the larger levels. In addition using different filters, return back different qualities of coronary arteries.
An example of using filter [1,3,3,1] with some different levels on the inverted version of angiogram are shown in Fig. 7 respectively on the inverted version of preprocessed angiogram, will be the best way of extracting the whole coronary arteries. Our proposed algorithm for processing phase is shown in Fig. 8. It worth nothing that above compound levels and filters should be used consecutively.
As it is shown in proposed algorithm, at first we run SWT by filter 0 1 [1,2,1] / 4 h = and levels {2,3} on the input image. Then, we add this output to input image to show arteries clearly; note that adding two images means to add intensity values of each two correspond pixels. Afterwards we use the result for next step with the same filter but different levels {2,3,4,5} and then add to previous output. In the second iteration of main outer loop, the same procedure run, but by a new filter 0 1 [1,3,3,1] / 8 h = . Finally, we use output w 2 second inner loop, as output of our algorithm in this phase. The obtained results from each step are shown in Fig. 9.
As it is shown in this Fig. 9, coronary arteries will be extracted after using the proposed algorithm. In next section, we are going to do post-processing on the result for thresholding and removing some residual noises, such as backbone, ribs and small segmented objects which do not belong to coronary arteries and also filling holes in the arteries.

Post-Processing Phase
In post-processing phase, final coronary arterial tree is extracted by thresholding followed by the length refinement and then filling holes methods. For thresholding, we examined some threshold values (T E ) and found out the best value for coronary artery extraction in angiogram which is 28. Result of using this threshold is shown in Fig. 10.
As it is shown in Fig. 10, there are some residual noises, such as ribs and small segmented objects and also some holes in along the arteries. The traditional method for removing these types of noises is morphological approach which uses dilation and erosion operators. But by using this method, width and diameter of the coronary arteries will be changed after using these operators and subsequently accuracy will be decreased. Therefore, we defined another method, as follow. At first, we applied length refinement procedure. Because coronary arteries have a connected vascular structure, we supposed objects which have connected pixels more than a specific threshold (T L ) are belong to the coronary arterial tree and other objects are assumed as ribs and small segmented objects which must be removed. In the length refinement procedure, we experimentally examined different threshold values for size of coronary arterial tree parts and finally we found the best one which is 1500 pixels with image size 512×512. For angiograms with larger size, it must be changed to bigger values. We also used another procedure for filling the holes. As it is shown in Fig. 10, there are some holes which appear after thresholding in along the arteries. Addressing this problem, we defined a threshold for minimum size of holes. It means, smaller holes will be filled after filling procedure. We experimentally examined some threshold values (T F ) and found out the best one is 55 pixels. In result and discussion, we will discuss more about length refinement's threshold and minimum holes' threshold. The final result of coronary artery extraction after applying length refinement and filling holes procedures is shown in Fig. 11.

EXPERIMENTAL RESULTS
In this section, we intend to find optimum values for parameters used in proposed algorithm and afterwards evaluate the efficiency of the algorithm in compare with other methods and finally compare running time of the proposed algorithm with the latest methods. We obtained the database of angiograms from the UiTM Medical Center. The database includes 4932 angiograms from 9 different patients and for every image from 8 angles. The format of images is 24-bit greyscale BMP in 512×512 pixels size acquired by PHILIPS angiography system. We randomly selected 20 angiograms (10 LCA and 10 RCA) from this database as our data set in this study and manually created ground truth by labeling the coronary artery pixels in each of them. It is worth noting that the ground truth was validated accurately by the cardiologist and cardiac surgeon.
We applied the quantitative methods to determine the performance, such as accuracy, precision, sensitivity and specificity values (Han et al., 2006) In this research, the TP, TN, FP and FN metrics are defined as follow: TP: The number of coronary artery pixels detected correctly TN: The number of non-coronary artery pixels detected correctly FP: The number of non-coronary artery pixels detected as coronary artery FN: The number of coronary artery pixels not detected

Parameters Optimization
In this section, we evaluated the effect of different values for parameters, such as filters 0 1 (h ) , wavelet levels (Levels), Threshold for Extraction (T E ), Threshold for Length refinement (T L ) and Threshold for Filling the holes (T F ), on accuracy of our proposed coronary artery extraction method. We examined these parameters on 20 different angiograms from data set and finally calculated the average of each of them for doing optimization.
To obtain the best result in coronary artery extraction in this research, we applied filters [1,2,1]/4, [1,3,3,1]/8 and [1,4,6,4,1] on data set; firstly one by one and afterwards all together and then calculated accuracy for each of them. Results are shown in Fig. 12.
As mentioned before, thinner arteries are detected by smaller levels and thicker ones by bigger levels. We applied the range of wavelet levels 1,2,3,4,5 for this research and examined the effect of them one by one and all together. The obtained results are shown in Fig. 13.
As it is demonstrated in Fig. 13, the best result is obtained by using the combination of wavelet levels {2,3} and {2,3,4,5}. In fact, the thinner arteries are presented by levels {2,3} and thicker ones presented by levels {2,3,4,5}. In post-processing phase, we introduced some parameters for the coronary artery extraction, length refinement and filling the holes. Here, we discuss more about these parameters. The first parameter which should considered is threshold value for coronary artery extraction T E . We empirically examined range of 20 to 30. The obtained result is shown in Fig. 14. As it is shown in this figure, the best result for coronary artery extraction is related to threshold value 28.
Another parameter is length refinement threshold T L . For this parameter, we examined values between 1000 and 2000. The reason for choosing this range is the minimum value for cleaning the segmented objects should be 1000 and the maximum value for avoiding cleaning the parts of arteries should be less than 2000. Therefore, we examined coronary artery extraction in that range and we found out the best value for this threshold is 1500. The obtained result is shown in Fig. 15.
The last parameter is threshold for filling the holes T F . The result of applying different values for T F is shown in Fig. 16. We examined the values between 30 and 70. The reason of choosing this range is the minimum value for filling the holes which appears in coronary artery extraction thresholding is 30 and the maximum value for avoiding filling the holes between two close crossed arteries is 70. As it is shown in Fig.  16, we analyzed this range on the angiograms in data set and reached to this conclusion that the best value is 55.

Coronary Artery Extraction Evaluation
In this section, we evaluate capability of coronary artery extraction in proposed method on data set. At first, we run proposed algorithm on 20 angiograms on data set; includes 10 LCA and 10 RCA. Then, we calculated averages of accuracy, precision, sensitivity and specificity values for each of them.
As it is shown in Table 1, the average accuracy, precision, sensitivity and specificity values of the proposed method in LCA angiograms are 0.96934, 0.87797, 0.86014 and 0.98439, respectively. And as it is shown in Table 2, the average accuracy, precision sensitivity and specificity values of the proposed method in RCA angiograms are 0.97425, 0.93021, 0.89962 and 0.99587, respectively.
To emphasis on capability of proposed method, we compare our algorithm with some state-of-the-art coronary artery extraction methods on the same data set. To this end, proposed methods by Khaleel et al. (2012), Li et al. (2008), Frangi et al. (1998) and Bankhead et al. (2012) were used for comparison. The result is summed up in Table 1 and 2. As it is shown in these tables, all performance metrics such as, accuracy, precision, sensitivity and specificity of our method are much higher than other ones for both LCA and RCA angiograms.
Since Bankhead et al. (2012) applied Starlet wavelet transform method for vessel extraction in retina images, we discuss about their method here also. They used original Starlet wavelet transform algorithm by using filter [1,4,6,4,1]/16, but as mentioned before we have modified this algorithm, especially in filters and wavelet levels for extracting coronary artery in angiograms.
For visual observation and comparison, we showed in Fig. 17 and 18 obtained results from our proposed method and above methods on LCA and RCA angiograms, respectively. As it is illustrated in these figures, there are so many non-coronary artery pixels which detected as coronary artery in Khaleel et al. (2012), Li et al. (2008) and Frangi et al. (1998) methods. This issue increased the FP rate in these methods and consequently decreased performance as mentioned in Tables 1 and 2.
Also, using the original Starlet wavelet transform algorithm, as used in Bankhead et al. (2012), demonstrated arteries thicker than actual size and also increased FN rate; because it could not detect coronary artery pixels in some parts and as a result it decreased performance especially in sensitivity value as shown in Tables 1 and 2.

Running Time Evaluation
We have implemented our proposed method in MATLAB R2013a. In the implementation, we only used original functions in Image Processing Toolbox and because we aimed to decrease the running time, the additional compiled MEX code is not used here. In Table 3, we compared proposed method running time with some state-of-the-art methods. In this table we calculated running time of Khaleel et al. (2012), Frangi et al. (1998), Li et al. (2008) and also our proposed method on same PC and conditions, where the running time of the Zhou et al. (2008) method was obtained from their paper. As it is shown in Table 3, proposed method needs low computation time, less than 1 second and it is a positive point of using this method in real-time systems. Therefore, we can consider it as a fast method for coronary artery extraction in angiograms.

CONCLUSION
In this research, we proposed a fast and accurate approach for automatically coronary artery extraction in angiograms. The proposed method is based on Starlet wavelet transform and discrete wavelet transform. To this end, we applied discrete wavelet transform for de-noising angiograms and then improved original Starlet wavelet transform algorithm for applying on coronary artery extraction aim. Finally, we proposed length refinement and filling holes techniques in post-processing step to eliminate the residual noises. The average accuracy, sensitivity, specificity and precision values of proposed method in LCA angiograms from data set are 0.96934, 0.86014, 0.98439 and 0.87797 and in RCA angiograms are 0.97425, 0.89962, 0.99587 and 0.93021, respectively. Obviously, the proposed method is robust in all of performance metrics. Also comparing results shows that proposed method has minimum artifacts and it also extracts whole parts of the coronary arteries. Furthermore, proposed method running time is much better than other methods; the whole process is done in less than 1 second.