Automatic Vasculature Identification in Coronary Angiograms by Adaptive Geometrical Tracking

As the uneven distribution of contrast agents and the perspective projection principle of X-ray, the vasculatures in angiographic image are with low contrast and are generally superposed with other organic tissues; therefore, it is very difficult to identify the vasculature and quantitatively estimate the blood flow directly from angiographic images. In this paper, we propose a fully automatic algorithm named adaptive geometrical vessel tracking (AGVT) for coronary artery identification in X-ray angiograms. Initially, the ridge enhancement (RE) image is obtained utilizing multiscale Hessian information. Then, automatic initialization procedures including seed points detection, and initial directions determination are performed on the RE image. The extracted ridge points can be adjusted to the geometrical centerline points adaptively through diameter estimation. Bifurcations are identified by discriminating connecting relationship of the tracked ridge points. Finally, all the tracked centerlines are merged and smoothed by classifying the connecting components on the vascular structures. Synthetic angiographic images and clinical angiograms are used to evaluate the performance of the proposed algorithm. The proposed algorithm is compared with other two vascular tracking techniques in terms of the efficiency and accuracy, which demonstrate successful applications of the proposed segmentation and extraction scheme in vasculature identification.


Introduction
World Health Organization's survey of "The top ten causes of death" acknowledged the fact that coronary artery diseases (CADs) are the leading cause of human deaths worldwide. CADs were responsible for 7.25 million deaths in 2008, which accounted for 12.8% of the total deaths worldwide, and this number of deaths has been increasing ever since [1]. The Xray angiography is an effective technique for imaging of the coronary artery and is considered as the "golden standard" for clinical observation of coronary anatomy and identification of vascular stenosis [2]. Therefore, it is widely used in clinical diagnosis and monitoring of disease. The coronary artery obtained from the X-ray angiograms can provide useful parameters for quantitative assessment and diagnosis of cardiovascular disease. Furthermore, the extraction of coronary arteries from the sequence of angiographic images is an important basis for heart motion analysis [3][4][5] and 3D vascular reconstruction [6][7][8][9]. However, fully automatic, robust, and accurate extraction of coronary artery from angiograms is still a challenging task so far. The main difficulties for the accurate vascular structure extraction or identification in angiograms are as follows: (1) the irregular gray level distribution of blood vessel due to the uneven perfusion of the contrast agent; (2) the nonvascular structures background interference including bones, catheters, and soft tissues; (3) the diversity of direction and width of the vessels; and (4) the commonly existing motion artifacts due to the heart motion and the presence of some pathological lesions.
In the past two decades, a few methods have been studied for the extraction of the blood vessels in angiographic images, such as morphology based methods [10], tracking based methods [11][12][13][14][15][16], multiscale based, methods [17,18], edge detection methods, and registration based methods [19,20]. Among them, the tracking-based methods proved to be very effective, which can detect coronary artery according to the local response of the angiograms and do not need to scan the whole image. Also, the methods extract the parameters including centerlines, diameter, or bifurcations using an adjustable pattern element to fit incremental section tracking procedures.
For conventional tracking-based algorithm, adapted diameter measurement [11] often suffers from artifacts due to nonuniform contrast distribution of the contrast materials. Diameter measurement approaches have been applied to the single segment and the full vascular network and can produce acceptable quality in the coronary artery extraction [21]. And also, certain new template-based techniques have been introduced [14][15][16] to fit the coronary artery structure in the tracking procedure, such as rectangular or circular templates [15], Gabor filters [14], and Gaussian kernels [16]. The templatebased model construction is known to be complicated and time consuming. Most of them need to create massive templates according to the variation of diameter and direction of coronary artery.
To reduce the calculation complexity, an alternative solution of intensity ridge detection is utilized to approximate the medial axes of the coronary artery. Aylward and Bullitt [12] proposed a centerline tracking algorithm based on ridge detection in multi-scale space, which is constructed by extracting the ridge using eigen decomposition of the Hessian matrix. Zamani Boroujeni et al. [13] employed a ridge scanning scheme for reliable identification of the vessel points and calculating the magnitude by adaptive look-ahead distance method. However, due to the effects of low image quality and noise, ridgeline cannot be exactly located at the geometrical center of coronary artery.
In this paper, an adaptive geometrical vessel tracking (AGVT) algorithm is proposed for the automatic extraction of coronary artery from angiograms. There are two main contributions of this study. First, an automatic initialization algorithm including location and direction of the seed point is proposed. Second, ridge points are recursively detected in consecutive scanlines, which are then adjusted to the geometrical centers according to the estimated diameters. The proposed method does not require any human intervention beforehand and can simultaneously estimate the parameters for quantitative analysis of vasculature including centerlines, diameter, and bifurcations.

Methods
The proposed AGVT algorithm is composed of three main steps: ridge enhancement, seeds determination, and adaptive tracking. The brief description of the calculation procedures is as follows.
(1) Ridge enhancement: the angiogram is first convoluted with a Gaussian kernel with different standard deviation. For each scale of convolution, the blood vessel is enhanced by combination of eigenvalue response of the Hessian matrix, and the maximum responses for each level of scale space are extracted as the RE image. The enhancement procedure guarantees that the background is suppressed and the vascular structures are highlighted in the RE image.
(2) Seeds determination: the seed points are detected by scanning of local maximum, for which the intensity is brighter than other points in the local area. The forward and backward tracking directions of the seeds are designed to detect the neighboring ridge points.
(3) Adaptive tracking: rough centerlines are extracted through sequentially detecting ridge point in RE image, which are then further refined adaptively to the geometric center according to the distance between the vascular boundaries and current calculating point. If the bifurcation is identified by discriminating connecting relationship of the tracked ridge points, the tracking will be divided into two paths for each branch. Termination criteria are assigned to ensure that the tracking is within the region of coronary artery without repetition. Finally, all the detected centerlines are merged and smoothed, and also false tracking results are removed.

Ridge Enhancement.
In this paper, we use a multiscale vascular enhancement filter (MVEF) [22] to obtain the RE image. It enhances all the tubular targets along the centerline and fades out the background. Suppose the ( ) is the Hessian matrix of a point on image ; then, the Hessian matrix of the each pixel can be calculated as follows: where subscripts and denote the second order derivate along or direction. The eigenvalues of the matrix are denoted by 1 and 2 ( 1 ≤ 2 ). According to the analysis of [22], different eigenvalues of the Hessian matrix are corresponding to different types of structure, such as plate-like, blob-like, tubular, and noises. In order to detect different size of vessels, Gaussian kernel with different standard deviation is usually convolved with the angiogram [23,24]. For a certain scale , the intensity differential of the point is expressed as follows: where is the normalization coefficient as defined in [24] and (0, ) denotes a Gaussian function with the mean of 0 and standard deviation of . The of the Gaussian kernel is designed as a variant value to enhance different scale of vessels. In implementation, is usually designed as a value between the maximum and minimum size of the vascular diameter to be enhanced. The enhancement response of each pixel in scale space is computed, and the maximum response is then utilized as the final enhancement result. The enhancement response ( ) of pixel can be calculated as follows: where and are control parameters, while [ min , max ] is the size of the scale space.

Seed Points.
After the ridge enhancement, the local maximum points in the gray level space of the RE image are located on the ridgelines, and they are detected as the initial positions for the tracking. According to [25], if a point ( , ) is a local maximum, then its gradient is equal to zero and its Hessian matrix is negative. However, the points satisfied with these two conditions are usually with float coordinates. Therefore, the seed points need to be interpolated according to their neighboring coordinates. For any pixel ( , ) and its neighbor pixels ( + 1, ), ( , + 1), ( + 1, + 1), if there is a point ( , ) ( < < +1, < < +1) that meets the conditions that ∇( , ) = 0 and 1 ( , ) < 0, 2 ( , ) < 0, then ( , ) is a local maximum on the image. According to the bilinear interpolation equation we have The solutions of (5) are implicit; hence, the approximate solutions can be solved by calculating the continuity of the image and its differential information. If a point satisfied the following equations: there will be a point ( , ) ( < < + 1, < < + 1), which satisfies the conditions ∇( , ) = 0 and 1 ( , ) < 0, 2 ( , ) < 0. Due to the continuity of image, ( , ) is an approximate solution of ( , ).
Since there are still numerous noises with low gray value in RE image, the extracted seed points are refined by an intensity threshold . If the intensity value of a candidate seed point is below the predefined threshold, it will be discarded.

Tracking Directions.
In this section, the initial tracking directions of the detected seed points are obtained in RE image. Suppose that is a detected seed point, while + and − are its neighboring forward (with greater abscissa value) and backward (with smaller abscissa value) points. Then As shown in Figure 1, for a seed point with coordinate ( , ), we define a circle with radius centered at , and the points ( ) ( ( ) = ( , )) on the circle can be expressed as parametric equations as follows: As there are two intersections between the circle and the centerline, the intersection points ( ) and ( ) hence can be detected by calculating gray level difference along the path of the circle. Here, the point ( ) can be detected by finding the maximum on the circle as follows: Obviously, ( ) is one of the points + or − , while ( ) can be detected by finding the local maximum on the local arc defined by : where and represent the forward and backward tracking angle of pixel . [2 − − Δ , 2 − − Δ ] represents the searching interval. The condition of cos( ) ≥ 0 means that ( ) is located at the right side of . Since noise is widely scattered over RE image, it is necessary to refine the seed points from the noise. Here, a rectangle with size of × is constructed on the center of each seed point. Then, if the average intensity within the defined rectangle is below a predefined threshold , this seed point is removed.  Figure 2 illustrates the process of ridge tracking. The location and tracking angle of current centerline point are denoted by and ; then, the relationship between the tracking direction u and tracking angle is defined as follows: To detect the ridge point̃+ 1 , the intensities of the pixels located on arc ( , , − Δ , + Δ ) are enumerated and compared. And the local maximum with intensity equal to the ridge point̃+ 1 can be found as follows: where 2Δ is the search scope and is the calculating step size. Then, the initial tracking directionũ +1 can be defined by and̃+ 1 as follows: Let̃+ 1 and̂+ 1 represent the detected ridge point and its corresponding geometrical centerline point. As shown in Figure 3, the intensity profilẽ+ 1 ( ) is defined by the scanline of̃+ 1 , which is perpendicular toũ +1 . And is denoted by the curve parameter of̃+ 1 ( ).
Hence, two estimated edge points̃+ 1 ( + ) and̃+ 1 ( − ) can be detected according to the gradient information of +1 ( ): Computational and Mathematical Methods in Medicine where the parameter is the defined as the search size, while + and − denote the upward and downward direction of̃+ 1 , respectively.
Then, we count the average gray value of the vessel points and background to obtain the two corresponding edge points +1 (̃+) and̃+ 1 (̃−). The average gray valuẽ+ 1 of the vessel points in the profilẽ+ 1 ( ) is calculated as follows: And the average gray value of background̃+ 1 iñ+ 1 ( ) is calculated as follows: The two edge points̃+ 1 (̂+) and̃+ 1 (̂−) with gray value less than (1/2)(̃+ 1 +̃+ 1 ) on either side of̃+ 1 are identified. Then, the geometrical centerline point̂+ 1 corresponding tõ+ 1 can be calculated as follows: According to (13), the adjusted tracking directionû +1 of +1 can be then updated aŝ Once the positions of the two edge points̃+ 1 (̂+) and +1 (̂−) are determined, the diameter̂+ 1 of the centerline point̂+ 1 can be measured aŝ (b) Bifurcations Identification. As shown in Figure 4, −1 and are two detected centerline points in the previous tracking process. If is a bifurcation, then two ridge points 1 +1 and̃2 +1 can be found on RE image at each branch of the vessel due to the process of MVEF. According to the adaptive tracking described in the previous section, we are able to get two geometrical centerline points 1 +1 and 2 +1 . Thereafter, the tracking process is proceeded in the two directions of (c) Termination Criteria. The tracking procedure starts from a randomly selected seed point and then iteratively extends to the other seed points. Three termination rules are defined to prevent the tracking path from going out of the image range, or being repeatedly calculated.
(1) If the detected centerline point goes beyond the scope of the image.
(2) If the detected centerline point goes beyond the scope of the vessel.
(3) If the detected centerline point intersects an extracted centerline.
Condition (1) is designed to guarantee that both abscissa value and ordinate value of current calculating point are within the border of the image while condition (2) is designed to ensure that the gray value of current calculating pixel is within the range of the vascular boundaries. Condition (3) is devised to determine if current calculating pixel is being detected or not.
A seed point tracking process is stopped immediately if any of the above three conditions occurs, and a new seed point tracking procedure will be launched. By setting these three conditions, it can be guaranteed that the tracking procedure will be carried out in the region of the blood vessels without repetition.
(d) Vasculature Refinement. After all discrete incremental sections are obtained from the angiographic image, vasculature refinement needs to be done to remove the false tracking noise, such as hair-like noise and discrete small edges. All the discrete centerlines are merged into several connected centerline sets, and the short false tracking results are removed by a predefined length threshold. Then, the vessel structure with the largest connected components is reserved as the tracking result. To avoid "jagged" tracking phenomenon, the cardinal spline interpolation [26] is used to preserve the smoothness of the tracked paths.

Results and Discussion
To validate the performance of the proposed vascular extraction method, series of experiments are designed and tested on both synthetic data and real clinical angiograms. And the comparative evaluation of the results demonstrates the efficiency of the proposed method over the several existing methods. Our tracking algorithm is implemented in Microsoft Visual Studio 2010 on an Intel Core i7 PC (with CPU 3.5 G and 16 G memory), and all the simulated angiographic images are with the resolution of 512 × 512.

Synthetic Data.
In order to quantitatively evaluate the performance of the proposed method, a series of angiograms with defined vascular structures are designed and simulated. The angiograms are simulated by projecting 3D synthetic cylindrical segments onto the image plane according to the perspective projection model [27] of the angiography system as developed in [13,28,29]. In order to simplify the simulation procedure, the background image is acquired prior to the injection of contrasting substance, so there is no visible vessel or catheter in the obtained image.
The projection intensity of vessel segment model adopted for the simulation comes from [27], which can be defined as follows: where is the distance between current point and tentative centerline point, is the radius of vessel, is the linear attenuation coefficient, and ( ) is the intensity of the projected point.
To test the proposed algorithm over noise interference, 25 angiograms containing a single vessel segment with different additive Gaussian noise are generated. The proposed algorithm then proceeds on the simulated image, and the extracted centerline and diameter are compared with the predefined vessel, and the error distributions over different noise scale for the centerline and diameter estimation are shown in Figures 5(a) and 5(b), respectively. From Figure 5 it can be seen that with the adding of noise power, the mean error of the centerline and diameter estimations increase gradually. Specifically, the mean error for centerline estimation ranges from 0.13 pixels to 0.33 pixels for noise power between 1 and 20; while for the diameter, the estimated errors range from 0.39 pixels to 0.62 pixels for noise power between 1 and 20. For noise power of 25, the image is considerably corrupted; however, the tracking centerline and diameter estimation errors are less than 0.4 and 0.9 pixels. It is obvious that the proposed method is every effective and robust even under terrible image qualities.
In order to simulate the vascular stenosis, the diameter is adjusted by a designed exponentially decaying function as formulated as follows [29]: where 0 is an initial diameter and is taper coefficient. and are Gaussian function parameters. For this calculation, at the stenosis part ( = ), 100% stenosis can be simulated by setting = 1. Stenosis rate of vessel can computed as follows: According to (21) and (22), nine types of arterial segments with stenosis are generated, as shown in Figure 6, and Computational and Mathematical Methods in Medicine the performance of the proposed algorithm is compared with the other two existing centerline tracking algorithms. To simplify the following description, nine types of arterial segments named V1 to V9 are designed as follows.
The tracking results of the three algorithms are listed in Table 1, which contains the maximum, minimum, mean, and RMS (root mean square) of the estimated centerline error. The mean errors are 0.34, 1.27, and 0.48 pixels for the proposed algorithm, Aylward, and Boroujeni algorithms, respectively. The proposed method represents 73.22% and 29.17% improvement over Aylward and Boroujeni algorithms. For the most complex structure V9, the mean tracking error for the proposed algorithm is 0.41 pixels, while it is 1.61 pixels and 0.87 pixels for Aylward and Boroujeni algorithms. The proposed method represents 74.53% and 39.87% improvement over Aylward algorithm and Boroujeni algorithm. Obviously, the proposed algorithm provides a more accurate scheme for coronary artery extraction and identification in X-ray angiograms.
The above experiments are repeatedly tested on 45 groups of synthetic images, for which each type of the vessel structure is randomly generated for five times. For these experiments, the average computational times are 3.66 s, 3.82 s, and 3.85 s for Aylward, Boroujeni, and the proposed algorithms, respectively. It shows that all these three algorithms take roughly the same computational time for the vessel extraction in angiograms. The computation efficiency can be further improved by integrating GPU based implementations in the future.

Clinical Angiograms.
In order to investigate the applicability and efficiency of the proposed algorithm, the proposed method is also applied to 20 coronary angiograms which were randomly selected from a database routinely acquired from the Philips Digital Imaging System at the Beijing Chaoyang Red-Cross Hospital during invasive catheterization procedures. Here, the experimental result of a random angiogram (Figure 7(a)) is selected to illustrate each step of the calculation. Figure 7 Figure 7(f) shows the boundaries of the extracted vessels. It can be seen from the images that through the process of MVEF, vessel is highlighted and the background is suppressed in RE image, and it serves as useful auxiliary information for the followed tracking steps. Through the seed determination of the proposed algorithm, it provides a set of seed points on the vascular centerline, and their directions are roughly parallel to the path of the vessels. From the figures, it can be seen that the proposed AGVT algorithm can extract main structure of coronary artery. Figure 8 shows four representative tracking results from these 20 angiograms, and we can see that the main structure of the coronary arteries can also be extracted by the proposed algorithm. Figure 9(a) shows a section of curved vessel for which the gray distribution is irregular. Figures 9(b), 9(c), and 9(d) show the tracking results of Aylward, Boroujeni, and the proposed algorithms, respectively. It can be seen that the tracking result of Aylward algorithm is more likely to change with the variation of gray level, which hence is more dependent on the gray level distribution of the angiogram. For Boroujeni algorithm, the detected centerline jitters significantly because it tracks the ridgeline rather than the geometrical center of the vascular structures. The proposed algorithm works better than the other algorithms as it adaptively extracts the geometrical centerline of the coronary artery and is more robust in processing irregular vascular segments.
Despite the advances of the proposed algorithm demonstrated in this paper, Figure 10 shows two challenging  angiograms from the set of 20 angiograms on which the proposed algorithm is only partially successful. It can be seen from Figure 10 that the proposed algorithm is failed for areas with overflowing (a) or incomplete filling (b) of the contrast agents in the angiograms. Up to know, segmentation of the vessels in such an angiogram is still a very challenge task.

Conclusions
The major contribution of this proposed algorithm is providing an automatic AGVT framework for coronary artery extraction by putting forward three major key steps named ridge enhancement, seed determination, and adaptive tracking. The proposed algorithm makes three key advances. First, considerable validated seed points and their tracking directions are automatically detected from angiograms, and it ensures the tracking affected without the interference of human factors. Second, it tracks coronary artery from angiograms without constructing templates, and, hence, the calculation complexity is greatly reduced. Third, the extracted ridge points are adaptively adjusted to the corresponding geometrical center, and it has been shown to work well with angiograms of low image quality. From the experiments on synthetic data and real angiograms, it can be stated that the proposed AGVT algorithm is more accurate and stable than the current widely used algorithms.