Abstract
Of interest is the accurate and robust delineation of vessel center-lines for complete arterial tree structure in coronary angiograms which is an imperative step towards 3D reconstruction of coronary tree and feature-based registration of multiple view angiograms. Most existing center-line tracking methods encounter limitations in coping with abrupt variations in local artery direction and sudden changes of lumen diameter that occur in the vicinity of arterial lesions. This paper presents an improved center-line tracing algorithm for automatic extraction of coronary arterial tree based on robust local features. The algorithm employs an improved scanning schema based on eigenvalues of Hessian matrix for reliable identification of true vessel points as well as an adaptive look-ahead distance schema for calculating the magnitude of scanning profile. In addition to a huge variety of clinical examples, a well-established vessel simulation tool was used to create several synthetic angiograms for objective comparison and performance evaluation. The experimental results on the accuracy and robustness of the proposed algorithm and its counterparts under difficult situations such as poor image quality and complicated vessel geometry are presented.
1. Introduction
In recent years, the utilization of computerized technologies in cardiovascular examinations has introduced a great deal of precision and speed to the diagnosis of coronary artery disease (CAD). In most of the vascular analysis applications, fast and accurate delineation of the arterial center-lines is a major prerequisite which provides a basis for subsequent image analysis steps. The aim of this study is to develop an efficient algorithm for producing a skeleton representation of whole coronary arterial tree. This is typically performed by either a pixel-based segmentation method followed by a skeletonization of the segmented image or direct exploratory center-line extraction in which coronary arterial segments are extracted through a recursive artery tracking algorithm. The objective of the first approach is to produce a separable representation of the foreground and background that entails a broad range of vessel enhancement or feature extraction methods such as matched/nonlinear filtering [1, 2], morphological filtering [3], eigenvalues of Hessian matrix [4], hysteresis thresholding [5], pixel classification methods [6], and many others. Unfortunately, most of these methods produce a large number of unconnected clusters of pixels instead of a single connected arterial tree, especially when the images contain nonuniform illumination.
As opposed to the tedious and error-prone pixel based segmentation approach, exploratory tracing methods directly extract the features of interest, circumventing low level processing of every pixel in the image. These algorithms are based on sequential searches through examining a small number of pixels that are close to the vasculature which results in efficient extraction of pixels located on the medial axis of the arterial segments. Several properties such as producing useful partial results, upon the occurrence of a deadline, and their computational efficiency make them attractive for real-time, live, and high resolution processing. Figure 1 depicts a general schema of the artery tracing approach followed by most existing methods in the literature. This schema is an iterative process to delineate a sequence of center-line points for each vessel segment. The steps of the process can be described as follows.(i)Initialization. The initialization step provides preliminary information for starting the tracing procedure. The information typically includes the location of the starting point, the initial direction of the vessel, and distances between the starting point and the left and right boundaries. In manual or semiautomatic algorithms, the location of the seed point and the initial direction are manually defined by the user. In case of fully automatic center-line extraction, however, a valid seed point is selected from a set of validated seed points provided by an automatic seed point detection algorithm. (ii)Estimation. At the initial point with directions , the location of the next center-line point is estimated simply by linear extrapolation in equation of the form: where parameter is a step size and determines the distance between the current point and the next point. Since an exploratory search produces a sequence of connected pixels, this distance is filled by a set of intermediary points to connect the two end points using a straight line drawing algorithm [7]. In some algorithms the step size is determined in a self-adaptive manner which yields more accurate results in comparison with the case that the tracing algorithm uses a fixed value for the step size, for example, [8].(iii)Updating the Tracing Direction. The extrapolation equation uses the vessel direction calculated only at the current center-line point and does not take into account the vessel geometry at the estimated next point. Hence, the estimated tracing direction is not always accurate and needs to be updated according to the local geometric and intensity-based information at the site of estimated next point. The idea of updating the preliminary estimate for the vessel direction has been suggested by many known tracing methods [8–10]. The idea is to refine the first estimate of the vessel direction by adjusting the location of the estimated next point based on its distances from the nearest boundary points.(iv)Tracking. In this step the tracing algorithm verifies the stopping conditions and repetitious traces and collects information about vessel intersections and branching or crossover points. If the estimated point is still located inside the vessel and the current vessel segment has not already been traced, it returns to the “Estimation” step to collect the subsequent center-line points; otherwise the tracking is terminated and a new trace is started from the next validated seed point.(v)Filtering the Small or Disconnected Segments. According to the application and required output, the tracing sequences which are smaller than a certain threshold and/or do not have any intersections with other segments are considered as false traces and are removed from the tracing result. In some applications, the operator can be given the ability to complete the tracing results by adding new seed points or deleting false traces.(vi)Curve Smoothing. Many well-known tracing algorithms produce jaggy or indented center-line due to their coarse angular quantization [9, 11–13]. In the sequel, they require a curve smoothing step to enhance the accuracy of the tracing results at the cost of more computation.
Different tracing algorithms vary primarily in their mechanism to estimate the next center-line point and local orientation of the arterial segment. Bolson et al. [14] proposed a method based on geometric properties of the vessel structures in the image. By manually defining the starting point and an initial direction, the algorithm estimates a new center-line point position and orientation by using a T-shaped structure which is rotated to find the best location for the next center-line point. Another representative algorithm for geometric-based vessel tracing approaches is proposed by Sun [10] that is the basis of many successful QCA systems. Instead of using directional filters or T-shape structure, their method was based on recursive sequential tracking of the vessel’s center-line with the assumption of geometric and densitometric continuity of the arteries in each incremental section. Their algorithm employs an iterative process with two extrapolation-update stages. In the first stage, an initial guess for the location of the next center-line point is made, assuming rectangular pattern for intensity profile defined perpendicular to the initial vessel direction. Then, the previous estimation is corrected by defining another profile at the new center-line point to find the point between two detected edges. Despite its accuracy and robustness, there are some major drawbacks with the Sun’s algorithm. Specifically, the algorithm uses matched filtering mechanism for identification of new center-line points in vessel cross-sectional density profiles which performs inefficiently when arterial segments with nonuniform intensity distribution are encountered. Moreover, their edge detection method is based on identification of roll-off points using pure intensity values; this can cause their method to have difficulties in coping with situations such as sudden changes in path-line orientation and vessel diameters [15].
In another study, Haris et al. [16] proposed a recursive tracking algorithm based on circular template analysis and appropriate model of the vessels. Although they showed that their method is very robust and outperforms its predecessors in terms of handling bifurcations and vessel crossings, it heavily depends on the vessel center-line and contour points detected at the artery tree approximation stage which may fail to detect all of the arterial segments in poor quality angiogram images. These drawbacks are stemmed from calculating the vessel direction merely based on the local geometric features extracted from the arteries and neglecting inherent intensity and contrast variations between the two corresponding edge points at the same coronary segment. Xu et al. [8] proposed a method to improve the Sun’s algorithm by combining it with a ridge based method proposed by Aylward and Bullitt [17]. In their method, the vessel direction is calculated based on a weighted combination of geometrical topology information obtained from Sun’s algorithm and intensity distribution information obtained from Hessian matrix calculation in Aylward’s method. They also proposed a self-adaptive look-ahead distance schema to increase the accuracy of the algorithm for extracting highly curved segments, and a dynamic size search window to cope with situations where two arteries are overlapped. Yet, some of the problems originated from Sun’s algorithm are still remained unsolved, causing Xu’s algorithm to deal with deviations at the site of severe stenoses. Also, there are many other recently published artery tracing algorithms in the literature. However, most of them have been developed for different applications or image modalities such as ophthalmic artery images [18] and CT angiography [19, 20].
The above mentioned limitations have motivated us to propose an improved algorithm which incorporates a semicircular vesselness profile for robust identification of next center-line point in the sequential tracking process. It uses reliable features to discriminate between the true vessel points and the points that do not coincide with arterial segments in the angiogram. In fact, instead of using pure intensity values to identify the true vessel points, it takes advantage of a feature image based on the eigenvalues of the Hessian matrix. Each pixel in the feature image represents a vesselness measure which associates the likelihood of being a vessel point to the corresponding pixel in the original image, allowing the tracing algorithm to robustly identify the center-line path along the arteries. In addition, an adaptive schema for magnitude of the search profile is incorporated to avoid divergence and premature termination of the tracing process.
The remainder of this paper is organized as follows: a complete explanation of our proposed method is presented in Section 2. In Section 3, the experiments conducted for parameter tuning are described and the results are presented. It also includes the experimental results obtained from comparative performance evaluation. Finally, Section 4 contains the conclusion and our decision for the future work.
2. The Proposed Algorithm
The focus of this study is to propose a robust and accurate algorithm for automatic extraction of complete coronary arterial tree from angiogram images. Toward this aim, a two step solution for fully automatic vessel center-line tracing algorithm is employed which is comprised of two main steps: automatic seed point detection and center-line extraction. In this study, we present those aspects of the tracing approach necessary to extract the center-line of the arterial segments; we do not discuss the details of seed point detection and validation steps as they are previously described in the literature [9, 13, 21].
2.1. Seed Point Detection
In the fully automatic tracing algorithms, the final output of the algorithm highly depends on the initial points that are provided for the tracing algorithm to start its process. In this work, we used our previously published method for fully automatic seed point detection [21] due to its capability to provide optimal balance between the accuracy of the validation procedure and the computational efficiency. To avoid pixel by pixel processing, the seed point detection algorithm samples the image by defining a sparse grid over the image and searching for the edge pixels along the horizontal and vertical lines. The number of grid lines determines the number of edge pixels where the grid lines cut across the vessel segments. The searching process involves identification of candidate boundary points by convolving the profile produced by each grid line with the first derivative of 1D Gaussian low-pass filter. Due to computational efficiency considerations, a 1D kernel of the form is considered as a discrete approximation of the continuous filter. By convolving this kernel with the intensity values along each grid line, local peak values of the filter response are identified and collected within a small neighborhood distance. Figure 2 illustrates the result of boundary point collection process in an example angiogram.
(a)
(b)
In the validation procedure, the collected boundary points are verified against a set of rules which are defined to discriminate between the actual seed points and misdetections. The algorithm incorporates the symmetry of geometric properties of the gradient vectors calculated at each candidate boundary point and its neighbors. After the false candidates are removed, the locations of center-line seed points are estimated by calculating the mid-point between the validated seed points and their corresponding point on the opposite edge. Also, the initial direction of the vessel can be represented by a unit vector perpendicular to the gradient vector which is calculated at the validated boundary point. It should be noted that the space between the grid lines is chosen to be much smaller than the longitude of the smallest arterial segments. Hence, the seed point detection algorithm may produce more than one seed point for each single segment which results in repetitious traces. To avoid this problem, an efficient mechanism is applied before validating each candidate seed point to prevent the assignment of a seed point to a previously traced segment; see Section 2.5.
2.2. Estimation
Once a reliable seed point is selected, the next step (so-called the estimation step) is to determine the exact location of the next center-line point. This necessitates the identification of true vessel points in a certain neighborhood of the selected seed point. However, X-ray images suffer from a high level of noise, nonuniform background, and existence of many image artifacts which make the task more difficult to perform. Among many vessel enhancement (and extraction) approaches proposed to overcome these difficulties, the outputs of vessel enhancement methods that are based on eigenvalues of Hessian matrix exhibit more attractive properties for our case. The first reason is that, instead of producing a logical value for each pixel (to show whether or not it is located on a vessel), these methods assign a continuous vesselness value to each pixel allowing the algorithm to identify true vessel point by finding the maximum vesselness value from a set of candidate next points in a small area surrounding the pixel. The second reason is that the innate computational characteristics of the eigenvalues of Hessian matrix allows for calculating the vesselness value based on local intensity information at each individual pixel in the image. This eliminates the need for low-level pixel processing for the whole image. The Hessian matrix at a given point is represented by: where denotes the second-order spatial derivative of the image at point , calculated by convolving the input image with the second-order derivative of the 2-D Gaussian function at a certain scale . According to [4], the eigenvalues and eigenvectors of Hessian matrix can be used to extract the principal directions of the local second-order variations at the vessel points. The two-dimensional Hessian matrix has two eigenvalues and and their corresponding eigenvectors and . The eigenvalues are assumed to be ordered such that: In coronary angiograms, vessels appear darker than the background. Thus, if we consider the input image as a 3D curvature surface, the vessel center-lines are represented by intensity valleys. Hence, for a given center-line point the eigenvector which corresponds to the stronger eigenvalue, that is, , reflects the direction of the stronger curvature within the small neighborhood around the center-line which is perpendicular to the vessel’s long axis. Since the eigenvectors are orthogonal, the second eigenvector, that is, , is parallel with the direction of the vessel. Based on the above considerations, the vessel points can be identified by examining the eigenvalues and as follows: However, this condition may also be met for some non-vessel points due to noise or line-like background structures. To obtain more deterministic criteria, Frangi, et al. [4] developed a multi-scale vesselness function which provides a value between zero and one for each point at a certain scale. We are interested in the 2-D version of their function which combines the measures of the curvature strength and the ratio of and to a single value measure as follows: where is the second-order structureness which accounts for the strength of the local contrast and is the blobness measure which differentiates between the tube-like and blob-like structures. The parameters and determine the sensitivity of the filter to the measures and , respectively. In this study, the optimal values of and are obtained based on the image characteristics.
2.2.1. Estimating the Vessel Direction
Since there is a wide range of blood vessel diameters in each angiogram, it is required to calculate the vessel resemblance values at various scales and combine them to obtain a single-valued metric. The combination process simply selects the scale which yields the maximum value of function at a given point : The selected scale is then used for selecting the resemblance value and the best estimate of the vessel direction as follows: Due to its computational complexity, this multi-scale calculation is considered as a major drawback of this vesselness function, since the core function only calculates the measure of vessel resemblance at a single scale. In our case, however, the vessel diameter, estimated for the current center-line point, can be used to calculate an appropriate range of discrete scales for the vesselness function, obviating the need for time-consuming calculations for various scales [22].
2.2.2. Semicircular Scanning Profile
In order to find a reliable next point, the value of the vesselness function is calculated for the pixels in a small neighborhood around the current known point where superscript indicates the th iteration of the algorithm. For this purpose, a semicircular scanning profile is defined which samples the vesselness measure for neighboring pixels around point within a certain radius . The semicircular scanning profile () is mathematically described as: where is the angle between the current tracing direction and the -axis and is a unit vector with direction . The radius , that is, the look-ahead distance, is adapted to the current vessel’s half width (the adaptation schema is described later in this section). Specifically, the value of the vesselness function (as well as the vessel direction ) is calculated at each point on the semicircular search area.
This scanning profile schema has been adopted by many other methods including square scanning profiles [23–25] and complete circumferential profile functions [16, 26] which proved to be useful in providing a uniform look-ahead distance in all directions. Nevertheless, the proposed method employs a semicircular scanning profile instead of a complete circle, employed by the previous methods, to avoid unwanted backward tracing and to maintain the computational efficiency. The proposed scanning profile schema is similar to the method proposed by Schrijver [26]. The main difference between the two methods is that the Schrijver’s method suggests the use of a single seed point as a starting point for recursive artery tracing algorithm. He described several conditions in which his seed point detection technique fails to provide a reliable seed point and the user intervention is required for manual selection of the initial seed point. As a consequence of using a single seed point, the scanning process requires relying on confidence scores and many threshold values to choose the potential tracing directions from several candidates for tracing the whole arterial tree. These difficulties make their algorithm unsuitable for real-time tracing applications.
2.2.3. Selecting the Correct Vessel Point
Figure 3 shows the scanning profile drawn at current point to find the next point at distance from point . In the process of finding the next point , the following situations can be distinguished.
Correct Vessel Point
If a given point is on the vessel, it can be recognized by a local maximum of the scanning profile. However, since some local maximum points may not coincide with the arteries in the image, another criterion should also be verified. Specifically, if the point is on the vessel, the direction of the vessel segment between points and is parallel to the direction field estimated by the second eigenvector of the Hessian matrix calculated at point, that is, in (8). Furthermore, the direction field calculated at point involves vectors with similar directions.
Nonvessel Point
A large part of the scanning profile is occupied by non-vessel points, for example, point , at which the value of the vesselness function is small and their neighboring points constitute a nonuniform direction field.
Points That Belong to a Vessel Branch
The values of scanning profile also attain a local maximum at point . In addition, the direction of the vessel segment between points and is also parallel to the vessel direction . Therefore, the algorithm should select one point between and . The objective is to select the main vessel segment based on the minimum difference between the angle of current tracing direction and the directions calculated by connecting point to each of the two candidate points and . In this situation, the algorithm marks the current point as a bifurcation point and adds point and its initial direction to the list of validated seed points to start a new trace in the next iteration.
Points That Belong to a Neighboring Vessel
Similar to the points that belong to a bifurcation segment, the points that coincide with a neighboring vessel constitute a local maximum in the scanning profile. However, as can be seen in Figure 3, the difference between the vessel direction at point and the direction of the vessel segment between points and indicate that it is very unlikely that the points and lie on the same vessel. In addition, if small values are chosen for radius , the problem of jumping between the vessels can be greatly avoided.
Figures 4(b) and 4(c) show the vesselness values and vessel directions of the pixels on the profile drawn at a bifurcation point in the example angiogram in Figure 4(a). It can be seen that the vesselness graph has two major local maximum for the vessel points located on the main segment and the branching artery. Further, Figure 4(c) illustrates two sets of uniform directions for points 10–26 and 41–50 which correspond to the points of the lobes in Figure 4(b).
(a)
(b)
(c)
It should be noted that the response of the vesselness filter is decreased at the site of branching points. However, since the global shape of the vesselness filter is of main concern and not its exact values, this does not affect the performance of the algorithm. Furthermore, if the current centerline point is on the branching point, the range of scales used to calculate the vesselness values for the branching segments is sufficient to identify the local maximum points on the branches.
2.3. Updating the Tracing Direction
Starting from the current center-line point and its initial values of direction and radius , a semicircular scanning profile is established to find the first estimate of the next vessel point denoted by as described in the previous section. Given the next point , the first estimate of vessel direction is calculated based on the geometric direction of the vector that connects the point to the point as follows: where denotes the magnitude of a vector. In most cases, this direction provides an accurate estimate of the vessel direction. Nevertheless, a new estimate is made by adjusting the location of such that it is located in the middle of local edges. As shown in Figure 5, to find the middle point, two linear density profiles and are drawn at point perpendicular to the direction . Then, two edge points and are detected by any edge detection algorithm such as finding the roll-off point based on signal and background levels of intensity values [10], directional low-pass filters [27], weighted sum of first and second derivatives of gray values [23], and many others.
However, our interest is to find the edges based on contribution of more than one pixel to detect the vessel borders in the original image. Therefore, the edges are identified by finding the maximum value of the local gradient magnitude (contrast) calculated for each point on the profiles and as follows: where is an estimate of gradient magnitude at the th pixel location on the scan profile and is the length of the search profiles which is adapted to the current vessel radius . Since the radius of the semicircular scan profile satisfies our need for defining a search window that sufficiently spans the vessel width, the value of parameter is chosen to be equal to radius . The calculation of right edge point is the same as for , thus its equations are omitted for the sake of brevity. After calculating the location of local edge points and , the location of next center-line point can be updated as follows:
Once the next point is identified, the current vessel radius is updated. Then, the next step is to update the vessel direction vector according to the updated next point as follows:
2.4. The Schema for Adaptation of Step Size
The final position of the next center-line point is determined based on the position of the current point and the value of the step size. An important challenge is to select an appropriate value for. Since thin and small vessels are naturally more flexible and tortuous than the large ones, the tracing algorithm should take smaller steps to describe them with larger number of points. One solution is to take the radius of the scanning profile as the step size to control the distance between the current center-line and estimated next point. As depicted in Figure 6, radius should be greater than the current vessel’s half width because the semicircular profile should cut across the vessel borders at distance from the current center-line point. Therefore, the radius of the semicircular scan profile is calculated adaptively based on the size of vessel half width at the current center-line point: where parameter is a constant factor and and are the vessel’s half width calculated for the current and next center-line points with superscript denoting the iteration number. The term accounts for controlling the magnitude of the step size when a branching point or sudden change in the vessel’s radius, that is, high-grade stenosis, is encountered. As shown in Figure 6, once the tracing algorithm reaches to a severe stenosis, the vessel’s half width calculated at the next point is much smaller than that of the current center-line point. Thus, if this sudden change is not taken into account, the size of subsequent scan profiles would not be large enough to surround the vessel boundaries. In this situation, it is impossible to detect the two edge points and , resulting in premature algorithm termination or divergence.
A large difference between the current and previous estimation of vessel’s radius yields a large scanning profile, allowing the algorithm to overcome the problem of tracing the high-grade stenosis and handling the branch and crossover points. On one hand, should be kept relatively small for very large vessels in order to avoid coincidence with neighboring vessels. On the other hand, large values of should be selected for small vessel widths, for example, less than 5 pixels, because significant rounding errors in calculating result in obtaining values that are equal to estimated vessel half width. The experiment conducted to find an optimal value for parameter will be explained later in this paper.
However, as opposed to look-ahead distance , the values of step size can be less than current vessel radius because it should account for variations of the vessel direction. An alternative solution suggests using the difference between the current and previous vessel directions for adapting the step size . Let be the angular difference between the vessel directions and calculated for previous and current center-line points and , respectively, that is, . The values of represent the change of local curvature along the current vessel segment such that when the tortuousness is small, the value of is small, and vice versa. Therefore, the value of can be used to adopt the step size based on current estimation of the vessel’s half width [8]: where denotes the step size calculated adaptively for th center-line point. As the value of angular difference is between 0 and , the magnitude of the term never exceeds the vessel’s half width. The above equation yields a self-adaptive step size such that the tracing algorithm takes smaller steps over the highly curved arteries. Consequently, the proposed schema provides more accurate tracing results by improving the ability of the tracing algorithm to keep up with the abrupt direction changes and coping with complex vessel geometries. The empirical study for setting the optimal value for will be covered later in this paper.
2.5. Preventing Repetitious Traces
Starting from each validated seed point, the tracing process generates a sequence of vessel center-line points called “center-line segment” in the form of -triplets each of which includes the position of center-line, direction, and local vessel radius as follows:
As mentioned before, validated seed points can be located on any point along the vessel. Hence, starting from each validated seed point, a new center-line segment is created by tracing the vessel once in direction and once along . The results of several traces are stored in a two-dimensional array of integer values that has the same dimensions as the original image to maintain the center-line segments in a single array called “center-line map”.
Initially, the values of all elements in the center-line map are set to zero. When a new segment is traced, a variable called “segment number” is incremented by 1. In order to store a center-line segment, the corresponding pixels in the center-line image are set to the non-zero value of the current segment number. This technique allows the tracing algorithm to prevent repetitious traces. Basically, there are two situations that should be checked to see if the current vessel has already been traced.(i)Before Validating a Seed Point. Before applying the validation rules to a given seed point, the seed point detection algorithm should check the center-line map for the existence of a previously traced segment in a small neighborhood of the candidate seed point (e.g., 5 × 5 neighborhood points). The seed point that is found to be located in the neighborhood of an already traced segment is ignored and a new seed point is selected from the collection of the candidate points.(ii)During the Tracing Process. The center-line image is also used for the detection of intersecting vessel segments. Before calculating the vesselness values for the pixels on the semicircular scanning profile, their corresponding pixels in the center-line map are checked for a non-zero value. Specifically, the pixels of the scanning profile that encounter a non-zero value are collected in a candidate list and the nearest pixel to the current center-line point is considered as the intersection point. In this situation, all the points that are located on the straight line which connects the current center-line point to the intersection point are added to the current center-line segment and the current trace is terminated accordingly.
There are two points that should be noted here: (1) if the length of a traced segment is shorter than a certain threshold (e.g., less than 20 pixels), the segment is discarded and its pixels should not be added to the center-line map; (2) the segments that intersect with themselves are considered as false traces and should be rejected.
2.6. The Stopping Criteria
The traced segments should be limited to the points that belong to the arteries in the image. Accordingly, the tracing algorithm repeats tracing until one of the following criteria is met.(i)One or more of the pixels on the scanning profile does not coincide with the image field.(ii)No valid vessel point is found in the scanning profile .(iii)The current center-line segment intersects a previously traced segment. This condition is checked for all the pixels on the straight line that connects the point to .(iv)The percent dynamic range () of the vesselness values of the points in cross-sectional profiles and is below a certain threshold.
We assume that the vessel segments have continuous densitometric features and the percent dynamic range does not vary significantly along the arteries. Based on this assumption, the percent dynamic range is computed based on the vesselness values of the points of cross-sectional profiles and , and is constantly monitored in the tracing process to check if the stopping condition is met [10]. To calculate the dynamic range of the vesselness measure, the signal level is determined by the average vesselness values between the two edge points and in the cross-sectional profiles and drawn at point : where and denote the offsets of the edge points on the cross sectional profiles and , respectively. Also, the background level is defined as: Based on the above definition, the percent dynamic range of the vesselness measure can be determined by: The parameter is used to detect the situations where the estimated point is located on the background. In normal situations, the blood vessels have greater vesselness values than the background. In case of background tracing, however, the value of the signal level would be very close to the background level, resulting in a significant reduction in the value of . Therefore, the fourth stopping criterion is defined as: where is a threshold value for percent dynamic which is set empirically such that the optimum values for performance measures consistency and discrepancy are achieved.
3. Results and Discussion
The experiments aim at finding optimal settings for parameters used in the proposed algorithm, validating the functionality of the proposed algorithm, and demonstrating the efficiency of the proposed algorithm compared to the conventional methods by conducting comparative performance evaluation. They comprise of two different types of evaluation studies.
3.1. Simulation Study
In this experiment, the synthetic images with known center-line positions and tracing directions are processed by the proposed algorithm. The estimated results are then compared with the optimal results that are generated based on a priori data used in the creation of the synthetic images. This comparison is made to evaluate the ability of the center-line extraction algorithms to keep up with producing satisfactory traces in difficult conditions such as complex vessel geometry and low signal-to-noise ratio. The purpose of the simulation study is to analyze the performance of the proposed center-line extraction algorithm under various geometries of the vessel segment, different vessel contrast, and different values of signal-to-noise ratio. For this purpose, a method for generating a synthetic vessel dataset proposed in Greenspan et al. [15] is adopted. As shown by the authors, the method is able to provide an objective way for comparing different center-line extraction algorithms. Nevertheless, the original method is modified to generate a wider range of geometric features such as symmetric and asymmetric lesions, radial dilation of the vessels, and multiple lesions in a single segment. Figure 7 illustrates sample vessel images in the synthetic dataset. The dataset is composed of the following image groups.(i)19 vessels with zero curvature; zero taper or medium taper value with stenosis.(ii)23 knee-type vessels; no stenosis. (iii)9 vessels with curvature; with stenosis. (iv)9 multiple segment vessels; with stenosis. (v)13 multiple segment vessels with multiple stenosis; zero taper. (vi)13 multiple segment vessels with multiple stenosis; medium taper.
(a)
(b)
(c)
(d)
(e)
(f)
For each image group, four subgroups are generated by adding white Gaussian noise with different variance values to each original image, resulting in 344 synthetic images. The method used to generate synthesized vessels models the coronary angiogram image based on the 2-D geometrical representation of the vessel’s projection using four parameters: vessel taper, percent stenosis, and center-line curvature and curve length.
3.2. Clinical Examples
The performance of the proposed algorithm should also be evaluated by comparing the accuracy of the proposed method with existing methods when applied to real-world images. It is worth noting that this experiment requires executing seed point detection algorithm before the automatic center-line extraction. Since our interest is to compare the performance of individual center-line extraction algorithms and not the combination of seed point detection and center-line extraction algorithms, the starting points are provided by the same seed point detection algorithm for all the experiments performed on the clinical dataset. In this study, the final center-line images are achieved by executing the seed point detection algorithm proposed in [21] with the same parameters settings described in the paper, followed by any one of the opponent center-line extraction algorithms.
To obtain a set of reliable center-line images as ground truth data, a set of 315 angiograms were processed by a modified version of the ground truth estimation method proposed by Al-Kofahi et al. [11]. The images were randomly selected from a database of routinely acquired coronary angiograms with anonymous patient information at UKM Medical Center. It consists of a wide variety of vessels, with different types of coronary lesions (types A, B, and C in AHA classification) and different geometries of vessel segments. The selected images have spatial resolution 512 × 512 and 8-bit quantization acquired by a “GE-Innova 2100-IQ” C-arm system.
This dataset is preprocessed and the vessel center-lines were manually annotated to obtain reference standard center-line images. In the first step, the boundaries of arterial tree in each angiogram are manually traced 5 times by the same person at different times, ignoring small arteries with less than 3 pixels wide. This results in five corresponding edge images for each image in the dataset. Then, the images in a correspondence set were superimposed on each other such that the pixel value in the resulting image is a function of the number of overlapping pixels. This yields average images with unavoidable discontinuities. To remove the holes and discontinuities, a morphological closing operator with a 3 × 3 square identity matrix is used as the structuring element. In the next step, the boundary images are filled (with white color) to obtain binary images which illustrate silhouette of the coronary arterial tree. Finally, the skeletonization algorithm developed by Zhang and Suen [28] was used to estimate the location of the true center-line points. Nevertheless, the resulting center-line images were manually modified when needed. A set of angiograms and their corresponding ground truth images are shown in Figure 8.
(a)
(b)
(c)
(d)
3.3. Performance Measures
Algorithmic evaluation of center-line extraction techniques requires defining a set of performance measures. In this study, the focus is on improving the robustness of automatic centerline feature extraction while maintaining an accuracy level similar to the existing methods. Accordingly, two different types of performance measures are employed: (1) error estimation measures which provide quantitative metrics to evaluate the robustness of the proposed algorithm against difficult morphologies of the arteries, complex lesions and image degradation which are characterized by synthetic images; (2) accuracy measures which are used to assess the ability of the proposed algorithm in generating accurate tracing results in terms of consistency of the results with the ground truth skeleton images in the clinical dataset.
3.3.1. Error Estimation Measures
As mentioned earlier, in order to create a reliable and accurate synthetic dataset, the authors developed a vessel generating tool based on the method described by Greenspan et al. [15]. In their method, the reference center-line is generated by concatenation of semicircle curves with different lengths and constant curvatures. They defined a parametric equation for the semicircle curve as a function of curve arc length as follows:where is a position vector, is the curve arc length variable, and is the total length of the semicircle curve. Further, parameter is the initial position of the semicircle and denotes the tangent at . The next semicircle curve can be attached to the current curve by using the values of and as its initial definition. Based on the above formulation, two error estimates are defined as follows.
Normalized Global Distance Error
This error measure reflects the average radial distance between the points on the reference center-line and their corresponding points on the estimated center-line. In order to assure that the two corresponding points lie on the same curvature radius, the closest point technique is not used to find the correspondence between points on the two center-lines. Instead, for a given point on the reference center-line, a corresponding point is identified on the estimated center-line that lies on the profile which is drawn at point perpendicular to the reference center-line, that is, local tangent vector. Therefore, the normalized global distance error is defined as [15]:
where is the radial distance in pixels between the two corresponding point, is the vessel’s radius at point and is the number of points contributed in calculations. This formula exhibits more emphasis on the distances where the vessel’s diameter is small rather than the distances calculated at vessel areas with large diameter. It is used to evaluate the algorithms against our main objective which is the accurate extraction of vessel center-lines including coronary arterial lesions.
Global Orientation Performance Measure
In addition to the distance error, orientation distance between the corresponding points is also computed as the difference between the corresponding tangents at the selected points. This performance measure is obtained by calculating the mean square orientation error as follows [15]:
where is the orientation error in degrees which reflects the difference between the direction of tangent vectors measured in degrees at the th point of the reference center-line. The traces in which the algorithm fails to cover more than 60% of the ground truth center-line are considered as divergence.
3.3.2. Accuracy Measures
In order to evaluate the accuracy of the proposed algorithm, a validation study is conducted on the accuracy of the tracing results when applied to real world images. In this study, the accuracy is measured with respect to ground truth images obtained from the clinical dataset and is defined based on “discrepancy” and “consistency” measures described in [11]. Discrepancy measures the quality of estimating the true location of the center-line points. It is calculated by computing the average Euclidean distance between the points of the center-line map produced by the algorithm and their corresponding points in the ground truth image.
Let denote a set of centreline points generated by the proposed tracing algorithm and be the set of ground truth points. Let two subsets and be the points of sets and that have a correspondence in another image. The correspondence indicates that for each point in subset , there is a corresponding point such that the Euclidean distance between the points is less than a particular number of pixels . The correspondence can be described by: where notation denotes the Euclidean distance. Similarly, for each there is a corresponding point whose Euclidean distance with is less than and is denoted by . It should be noted that due to the curved nature of the traces, there is no guarantee to find one-to-one correspondence between the points in and .
The spatial discrepancy between the center-line map produced by the proposed algorithm and the center-line image in the ground truth is defined as follows [11]: Consistency measures the ability of the algorithm in detection of true segments characterized by the ground truth and avoiding false traces. The consistency between two trace sets is calculated by finding the percentage of points in one set which has a corresponding point in another set that is within disk radius . Observe that the consistency is a mutual measure with similar definitions as follows:The first definition refers to the ability of the tracing algorithm in preventing false traces, while the second definition indicates the completeness of the tracing output. These measures are equally important to assess the accuracy of the proposed method, thus to compare the performance of different algorithms we calculate a single balancing measure called measure as follows: The values of are calculated for each image in the clinical dataset as a function of disk radius . The average values over all clinical images are used as the basis of our comparisons.
3.4. Parameter Tuning
Before performing experiments for performance evaluation, the optimal values of the algorithm’s parameters should be found. Primarily, two parameters and in (5) are tuned by examining their different value combinations on the performance of the vessel resemblance function. Referring to (5), it is expected that in most cases, the value of is close to 1. This is due to the fact that, on average, the values of and are similar. Therefore, in order to obtain more discrimination between the line-like and blob-like structures, the value of should be in the order of 1. On the other hand, determines the influence of contrast strength in vessel enhancement. By selecting large values for this parameter (e.g., in the order of 10), low-contrast objects are ignored and only vessels with significant contrast are enhanced.
The above conclusions are supported by Figure 9 which shows the effect of selecting different value pairs for and within a particular range of scales, that is, . It can be observed that, high values of increase the response of the vesselness function for vessel structures than the small values at the expense of enhancing more background structures. The small values of the second parameter, for example, , incorporate more noise and background structures in the outcome of the enhancement than the large values. However, the large values, for example, , result in significant reduction of the filter response even for the high contrast vessel areas.
The best values of and are selected based on the experiment conducted to compare the outcome of applying the vesselness algorithm on the images of the dataset with their corresponding ground truth silhouette images. The comparison procedure involves the following steps.(1)Calculating the number of false positives by counting the number of pixels in the vesselness image at which the vesselness value is greater than zero, but their corresponding point in the ground truth image is black.(2)The number of false negatives is also calculated as the number of white pixels in ground truth silhouette image which correspond to a zero vesselness value in the vesselness image.(3)Calculating the normalized sum of false detections which is used as an objective discrepancy measure that quantifies the deviation of vesselness images, obtained by applying different values of parameters and , from the ground truth silhouette images [29]: where is the number of all pixels in the image.
In this study, a set of 12 vesselness images corresponding to three values for parameter , that is, and four values for parameter , that is, were created. The appropriate values of and were chosen from the vesselness image with minimum discrepancy measure . The result of this experiment showed that the values of and are appropriate for the clinical dataset.
The remaining parameters are constant factor in (14) and in (20). To obtain the optimal value for , the proposed algorithm was used to extract the center-lines of clinical images. As mentioned in Section 2.4, a large value of parameter increases the size of the semicircular profile and, accordingly, raises the probability of jumping the tracing point from the current vessel segment to another vessel. Therefore, the the value of this parameter should be determined by measuring the consistency and the discrepancy of the algorithm’s output with the ground truth center-line on the real-world clinical images in which the problem of jumping between the vessels is probable. We did not use synthetic images for this purpose since each synthetic image contains a single arterial segment with no side branches.
For each image, a total of 11 center-line segment were obtained, corresponding to different values of ranging from 1 to 2 in a step of 0.1. Then, the average values are calculated over all images in the clinical dataset within disk radius of 5 pixels. Also, deviations from the ground truth center-line images are measured by calculating average discrepancy as shown in Figure 10.
(a)
(b)
It can be observed that the value of discrepancy rises as larger values are used for . This is due to the fact that the vessel curvature and vessel diameter varies more slowly within a small vicinity than a large distance along the vessel segment. Therefore, the large radius of scan profile has a negative impact on the accuracy of the estimated vessel direction. By comparing different values of measure, the value of is selected as the optimal value.
However, it should be noted that since is an integer value (number of pixels in semicircular scan profile), a fixed value for results in significant rounding errors for the vessels with small diameter (less than 5 pixels). Figure 11 shows the ratio of rounding error relative to the vessel radius for different values of vessel radius where is set to 1.3. Therefore, a constant value for (e.g., ) can be used when the vessel half width is less than 5 pixels.
Another parameter is the threshold for the percent dynamic range of the vesselness measure which should be set based upon the image characteristics. This parameter is tuned by observing the effect of its different values on the performance measures and by applying the tracing algorithm to the images of the clinical dataset. According to Figure 12, the optimal value for is obtained where the performance measure reaches to its peak at . It can be seen that the values of performance measures and almost remain steady when . After this point, the performance measure starts to grow steeply while the performance measure begins to fall more rapidly. By considering the definition of performance measures in (26a) and (26b), the above observation can be related to the reduction in the number of center-line pixels generated by the proposed algorithm. This reduction affects both performance measures and because setting larger values for parameter leads to obtain less false traces as well as less correctly traced segments. However, as the value of threshold increases, more correct traces are lost than the false ones. This implies that the proposed algorithm tends to produce less false traces and more true positives regardless of the value selected for parameter .
3.5. Experimental Results for Algorithm Validation
In this section, the efficiency of the proposed algorithm is assessed to see if it is able to produce satisfactory results. It should be noted that all heuristic schemas which are proposed for step size and look-ahead distance adaptation are employed in all experiments conducted for algorithm validation except otherwise stated. The synthetic images are divided into 7 groups of arterial segments, with different geometries and varying percentage of stenosis, for comprehensive validation of the proposed algorithm. The image groups are listed in Table 1. The first group contains the vessels with zero curvature, zero taper value with different percentage of stenosis. This group is used to evaluate the ability of the algorithm in addressing the problem of algorithm’s divergence at the site of high-grade stenosis. In this experiment, our focus is on cases in which large deviations occur at the site of stenosis on the straight vessels, that is, zero curvature with constant taper segments. As shown in Figure 13, the effect of percent stenosis on the accuracy of the proposed algorithm is markedly low such that for stenoses between 60–90%, the values obtained for both error measures remain in a reasonably low level and almost equal to each other. Furthermore, no divergence was observed for vessels with 95% stenosis.
(a)
(b)
The second group of images differs from the first group by changing the taper value from 0 to 0.00145, while almost the same values of percent stenosis are employed. Figures 14(a) and 14(b) present 7 error samples for the second group corresponding to vessels #13 to #19. By comparing the results in Figures 13(a) and 14(a), it can be seen that the distance errors plotted for different values of percent stenosis in group 2 are slightly higher than that of group 1; while approximately similar range of results (less than 1 degree) are obtained for the orientation error measure in Figures 13(b) and 14(b).
(a)
(b)
The results indicate that the vessel tapering has a trivial effect on the accuracy of the proposed algorithm in terms of estimating the curvature (tangent values) at the center-line points. As expected, for the vessel segments with larger taper values, a slight increase in the distance error measure is observed for all values of percent stenosis. This is due to the sudden change of vessel’s half width before and after the stenosed region which causes jittery behaviour in the estimated center-line.
The third image group is used to assay the behaviour of the proposed algorithm when applied to the vessel segments with different curvatures. To attain this goal, the other contributing factors were removed from the third group, that is, no stenosis and constant tapering value. Two experiments were performed on 24 images (corresponding to vessels #20–#43) with different values of curvature ranging from 0.003 to 0.015 radians. In the first experimental run, no adaptation is used for calculating the step size; while in the second run the step size is adapted based on angular difference between the current and previous estimates of the vessel direction. Figure 15(a) shows the performance results for the first experimental run by applying the proposed algorithm to the knee-type vessels in group 3. The graphs illustrate the distance and the orientation error performance measures as a function of curvature, that is, vessel segments with varying arc-length (from 50 to 175 pixels) are grouped based on their curvature values. The graph of distance error demonstrates an exponential relationship between the performance and the curvature of the vessel segments. In contrast, the values of orientation error in Figure 15(b), that are obtained from the second experimental run, exhibit more gradual increase in orientation error as a function of curvature value. This indicates that without step size adaptation, the proposed algorithm is highly sensitive to the vessel’s curvature in terms of estimating the distance rather than the orientation. In this condition, the algorithm is more accurate in estimating the direction of the vessel segment than estimating the correct position of the center-line points when highly curved segments are encountered. This is due to large errors in the approximation of the local vessel directions by using large values for the step size at the site of curved segments.
(a)
(b)
As explained before, the second adaptation schema is based on choosing values smaller than current vessel radius so as to describe the tortuous vessel center-lines with larger number of pixels. The graphs of error measures in Figure 15(b) show the effectiveness of the proposed step size adaptation schema. In the second experiment, the value of parameter is kept constant and the proposed algorithm utilizes the step size adaptation schema. In contrast to the graph of distance error in Figure 15(a), the values of distance error in Figure 15(b) increase more gradually as the curvature increases from 0.006 to 0.015. This indicates a considerable reduction in the distance error when highly curved vessels are encountered. Although no significant reduction is achieved for values of the orientation error in the second graph, a comparison between the slopes of the corresponding graphs indicates that the step size adaptation schema reduces the sensitivity of the tracing algorithm to the vessel curvature and improves its robustness to geometrically complex structures.
Validating the Step Size Adaptation Schema
However, this improvement is obtained at the cost of more divergence at the site of stenoses that causes rapid changes in the vessel’s diameter. The divergence is defined as a condition where the center-line produced by the proposed algorithm covers less than 60% percent of the ground truth center-line. To evaluate the ability of the tracing algorithms in extracting the center-line at the site of severe stenoses, we measure the algorithm’s success rate as the average number of center-line points in the estimated center-line that have a corresponding point in the ground truth center-line for different values of percent stenosis in vessel groups 4 and 5.
In Figure 16, the success rate of the proposed algorithm is plotted against percent stenosis for two adaptation schemas. In the first schema, only the radius is calculated adaptively and no adaptation is used for the step size , while in the second schema both radius and step size are calculated adaptive to the vessel’s half width. It can be concluded that the drawback of using adaptive schema for calculating the step size emerges in difficulties in coping with abrupt changes of the vessel diameter near to the severe stenoses, that is, vessels with percent stenosis above 90%.
3.6. Experimental Results for Performance Evaluation
To evaluate the accuracy and robustness of the center-line extraction algorithm, two experiments were conducted on the proposed algorithm and its three well-established counterparts: Sun algorithm [10], Aylward algorithm [17], and the algorithm proposed by Xu et al. [8].
In the method proposed by Xu et al., the vessel direction is calculated based on a weighted combination of geometrical topology information obtained from Sun’s algorithm and intensity distribution information obtained from Hessian matrix calculation in Aylward’s method. This combination is achieved by adjusting the weighting factor whose range is. Hence, the tracking direction is determined solely by the geometric direction when and the tracking algorithm becomes very similar to the Sun algorithm; while the tracking direction is determined solely by the intensity direction, that is, the Hessian eigenvector, when and the tracking algorithm becomes somewhat similar to the Aylward algorithm. Xu et al. suggest that more accurate results can be achieved by changing the weighting factor to 0.5. Accordingly, we implemented the algorithm developed by Xu et al. and compared our proposed method with the other three methods by setting parameter to 0, 0.5, and 1. Other parameters are listed in Table 2.
The first experiment aims to evaluate the accuracy of the algorithms in extracting the center-lines of the coronary vessels in clinical images. In the first step, the proposed seed point detection algorithm was used to provide an equal set of validated seed points for all the tracing algorithms. The optimal parameter values are used to setup the seed point detection algorithm. Given the validated seed points, the proposed algorithm and the three existing algorithms were employed to trace the artery center-lines in the images of the clinical dataset. The original image and the tracing outputs using different algorithms in a small area of an example angiogram are shown in Figure 17.
(a)
(b)
(c)
(d)
(e)
It can be clearly seen that the proposed algorithm outperforms the earlier solutions in terms of the accuracy of the tracing output. According to the outputs, the Sun algorithm is significantly distracted by the peripheral image artifacts while the Aylward algorithm fails to trace the highly curved segments. Also, the output obtained from the Xu’s algorithm is considerably accurate. This can be ascribed to taking advantage of the strength of both geometrical and intensity-based approaches for estimating the vessel direction. Yet, it contains small deviations where the minute branches are encountered.
Figure 18 illustrates the figures of average discrepancy measure over all clinical images at different values of disk radius . Observe that the proposed algorithm is superior to its opponents for most values of . For example, the average discrepancy calculated for the proposed algorithm is about 0.88 pixel at , which is approximately 0.1 pixels less than that of Xu’s algorithm at the same disk radius. This can be attributed to the fact that, in contrast to the existing methods, the proposed tracing algorithm uses robust features to identify the location of the next point, and thus it is less sensitive to intensity variations, illumination changes, and image artifacts which increase the discrepancy between the estimated and ground truth center-lines.
In Figure 19, the average of consistency measure over all clinical images is plotted for the proposed algorithm and its opponents. According to the graphs, Aylward has the lowest consistency figures for all values of disk radius. This is due to ignoring the geometric features of the arteries in estimating the location of the next point which results in high estimation error. The plots correspond to the other algorithms run almost closely to each other. However, a more precise observation revealed that, as it was expected, the Sun’s algorithm is outperformed by the algorithm proposed by Xu et al. for all values of disk radius . This superiority comes from utilizing the advantages of estimating the vessel direction based on eigenvalues and eigenvectors of Hessian matrix. Although Xu’s algorithm dominates the competing algorithms in terms of highest values of consistency measure for the first two values of disk radius (26% and 60% resp.), the proposed algorithm exhibits an average of 5% increase in consistency measure where . For instance, an average -measure of 82.1% was obtained for the Xu’s algorithm compared with 88.8% recorded for the proposed algorithm at . By considering the fact that the consistency measure has two contributing factors, namely precision and completeness of the tracing output, this improvement can be related to more complete tracing results or higher recall values obtained for the proposed algorithm at higher values of disk radius.
The second category comprises the experiments to assess the robustness of the algorithm to different amounts of impulse noise. The robustness of the algorithms is measured by calculating the algorithm’s success rate when an increasing amount of noise is added to the synthetic images. The Poisson noise is simulated by White Gaussian noise with a known range of variance. In this experiment, the noise was gradually added to the synthetic images such that the resulting signal-to-noise ratio declined from 20 dB to 10 dB.
As mentioned earlier in this chapter, the algorithm’s success rate refers to the proportion of the true vessel center-line that can be traced by the algorithm without premature termination. Figure 20 shows how different algorithms behave in response to increasing amount of noise which causes the tracing algorithms to diverge and terminate pre-maturely before the vessel segment is completely traced.
As demonstrated in Figure 20, the Aylward algorithm has the lowest success rates for high signal to noise images, that is, SNR ≥ 20 dB. In contrast, the other rivals have almost perfect success rates (at least 94%) within the same range of SNR. By increasing the amount of noise, the success rate of the Sun’s algorithm falls significantly, from about 93% at SNR = 20 dB to approximately 69% at SNR = 16 dB and reaches to the bottom of 9% for the images with lowest signal-to-noise ratio (10 dB). Interestingly, the success rate of the Xu’s algorithm follows the same pattern as the Sun’s algorithm does. However, on average, it exhibits almost 10% improvement in success rate for low quality images. It can be clearly seen that, in images with higher values of additive noise, Aylward algorithm yields more success rates than both methods of Sun and Xu. Also, it can be observed that all existing algorithms diverge if the SNR falls below 16 dB. The figures obtained for the proposed algorithm show the superiority of the proposed method to the existing algorithms in terms of robustness to the inherent noise of angiogram images. The results of this experiment showed that the proposed tracing algorithm obtained about 33% improvement upon the existing methods in terms of algorithm’s success rate in processing low quality images.
4. Conclusion and Future Works
All these observations led us to conclude that estimating the vessel direction based on eigenvalues and eigenvectors of Hessian matrix results in improvement in the robustness of the tracing algorithms. On the other hand, utilizing the geometric features of the arterial segments in estimating the location of the center-line points leads to obtaining more accurate results. It seems possible that the promising results obtained for the proposed algorithm are due to combining the advantages of the above mentioned approaches and avoiding the limitations associated with existing methods in handling highly curved segments and sudden changes of vessel diameter at the site of arterial lesions. The results of comparative performance evaluation showed that, according to expectations, the proposed method achieved a remarkable improvement in the accuracy of the tracing algorithm. Surprisingly, the proposed algorithm was found to be extremely more robust to image noise than existing well-known methods. This makes the proposed algorithm more suitable for feature extraction and quantitative coronary analysis from inherently noisy data in real-world applications. In the future, we plan to conduct a comprehensive study on the effect of utilizing different seed point detection algorithms on the performance of the whole center-line extraction method.
Acknowledgment
The authors wish to thank the medical staff of Cardiology Department of Hospital HUKM Malaysia for their support in the data collection effort.