Open Access

Automatic segmentation of coronary angiograms based on fuzzy inferring and probabilistic tracking

BioMedical Engineering OnLine20109:40

https://doi.org/10.1186/1475-925X-9-40

Received: 16 February 2010

Accepted: 20 August 2010

Published: 20 August 2010

Abstract

Background

Segmentation of the coronary angiogram is important in computer-assisted artery motion analysis or reconstruction of 3D vascular structures from a single-plan or biplane angiographic system. Developing fully automated and accurate vessel segmentation algorithms is highly challenging, especially when extracting vascular structures with large variations in image intensities and noise, as well as with variable cross-sections or vascular lesions.

Methods

This paper presents a novel tracking method for automatic segmentation of the coronary artery tree in X-ray angiographic images, based on probabilistic vessel tracking and fuzzy structure pattern inferring. The method is composed of two main steps: preprocessing and tracking. In preprocessing, multiscale Gabor filtering and Hessian matrix analysis were used to enhance and extract vessel features from the original angiographic image, leading to a vessel feature map as well as a vessel direction map. In tracking, a seed point was first automatically detected by analyzing the vessel feature map. Subsequently, two operators [e.g., a probabilistic tracking operator (PTO) and a vessel structure pattern detector (SPD)] worked together based on the detected seed point to extract vessel segments or branches one at a time. The local structure pattern was inferred by a multi-feature based fuzzy inferring function employed in the SPD. The identified structure pattern, such as crossing or bifurcation, was used to control the tracking process, for example, to keep tracking the current segment or start tracking a new one, depending on the detected pattern.

Results

By appropriate integration of these advanced preprocessing and tracking steps, our tracking algorithm is able to extract both vessel axis lines and edge points, as well as measure the arterial diameters in various complicated cases. For example, it can walk across gaps along the longitudinal vessel direction, manage varying vessel curvatures, and adapt to varying vessel widths in situations with arterial stenoses and aneurysms.

Conclusions

Our algorithm performs well in terms of robustness, automation, adaptability, and applicability. In particular, the successful development of two novel operators, namely, PTO and SPD, ensures the performance of our algorithm in vessel tracking.

Background

Accurate extraction of the coronary artery tree from coronary angiograms is important for the diagnoses, treatment, and clinical study of various coronary artery diseases. In particular, computer-assisted analysis can improve the performance of quantitative evaluation. It can reduce the inter- and intra-observer variations in determining the severity of coronary artery stenosis [1]. An efficient vessel extraction algorithm also enables the detection of coronary artery motion, as well as the reconstruction of 3D vascular structures from a single-plan or biplane angiographic system [2].

Anatomical structures, such as blood vessels, nerves, and bronchi, present themselves as line-liked structures in 2D images or as tubular structures in 3D images. Over the past 10 years, a variety of approaches have been developed for vessel segmentation from 2D and 3D medical images [3], which typically differ in terms of basic strategies or imaging modalities. Common segmentation approaches employ multiscale filters, morphological segmentations, deformable models, front propagation methods, tracking-based methods, and ridge traversal-based methods. Although many of these approaches are promising for vessel segmentation, developing fully automated, faster, robust, and accurate vessel segmentation algorithms remains highly challenging because of the complexity of vascular structures as well as large variations in image intensities and noise.

Tracking-based methods [49] exhibit a natural advantage in extracting arterial axis lines in the angiographic images because of their relative simplicity, as well as their adaptability to variations in vessel diameters. In addition, their ability to capture detailed quantitative descriptions of vessel axes, diameters, and boundaries, leads to high levels of accuracy. Traditional tracking-based methods, however, are guided by simple local features that limit their utility because of difficulties in efficiently integrating advanced features into canonical expressions for guiding the tracking. Therefore, the use of minimum cost functions or multi-scale filters [1013] to construct the tracking operator poses several advantages in the segmentation of coronary angiograms; this approach, however, may also be problematic when meeting bifurcations or vessel crossings because it can follow only the path with the strongest response.

In this work, we propose a fully automatic tracking-based method that can adapt to varying vessel curvatures and diameters resulting from arterial stenoses or aneurysms. It can also walk across intensity gaps along the longitudinal vessel direction, and manage vessel bifurcations and crossings.

Methods

An integrated framework was designed to solve the challenging problem of coronary angiogram segmentation. The main components in this framework are summarized in Figure 1. The automatically detected seed point is generally located near the root of the arterial tree or on the main trunk. The probabilistic tracking operator (PTO) was initialized with the initial seed point to search for a candidate artery element with the most similar vessel geometries (including vessel feature, direction, and diameter) in the neighborhood. The vessel structure pattern detector (SPD) was utilized in parallel to classify each detected artery element as bifurcation or crossing. The artery element classified as bifurcation is placed on the list of source data to enable tracking new vessel segments beginning from this artery element in subsequent tracking. The tracking algorithm stops when all points on the list of source data are tracked.
Figure 1

Overview of the proposed tracking framework for the segmentation of X-ray angiograms with (a) Preprocessing step, (b) Tracking step, and (c) Results.

During the process mentioned above, local features and measures, such as the Gabor response, vessel path, skeleton-line, and directional consistency, were also analyzed to integrate the relational vessel attributes for probabilistic vessel tracking and structure pattern identification. The PTO is designed to incorporate local features and measures to determine candidate artery elements along the same vessel segment. These candidates were selected particularly from the overlapping region between the binarized vessel regions and a sampling disk centered at the current artery element under consideration. The size of the sampling disk corresponds to the current vessel diameter, thus rendering the PTO adaptive to varying vessel diameters in case of arterial stenoses or aneurysms. Simultaneously, the SPD works to infer the vessel structure pattern, such as the distal end, segment, bifurcation, and crossing. To complete structure pattern identification, the vessel feature profile scanned from a circle around the current artery element was analyzed and a fuzzy inferring function was applied. In particular, multiple features such as vessel feature, direction continuity, and Gabor responses at the peaks of the vessel feature profile were integrated into the fuzzy inferring function for pattern identification.

Vessel feature and direction maps

The intensity profile across the longitudinal direction of the vessel may highlight an elevation or even a slight intensity dip at the center, depending on the contrast between vessels and background. On the other hand, the profile along the longitudinal direction of the vessel can be smooth [Figure 2(b)]. Accordingly, the multiscale detection or enhancement of vessel structures can be completed by convolving angiographic image I(x) [where x = (x, y)] with normalized second-order Gaussian derivatives G (x; σ f ) at different scales σ f . The use of a Gaussian filter is important in reducing the influence of noise, properly shaping the line profile and ensuring a large second derivative across the vessel.
Figure 2

Show the angiographic images, feature maps and the cross-sectional profiles of coronary artery. (a1) and (a2) depict two original angiographic images; (a3) and (a4) are their vessel feature maps, respectively. Subfigure (b) depicts the cross-sectional profiles of the vessel in the original image and in the corresponding feature map.

In 2D cases, let 2 I (x, σ f ) be a Hessian matrix combined with Gaussian convolution for angiographic image I(x). It can be written as
2 I ( x , σ f ) = σ f 2 { ( 2 x x 2 x y 2 y x 2 y y ) G ( x , σ f ) } I ( x )
(1)

where G (x; σ f ) denotes a Gaussian filter with standard deviation σ f .

The bright line structure on a dark background can be reflected by a Hessian matrix with a large negative eigenvalue λ1 and a small eigenvalue λ2 of positive or negative sign; that is λ1 < 0 and |λ2| 1|. Various response functions have been proposed in previous studies [1113]. We provide a similar response function to extract the vessel structures by thresholding the eigenvalue map using a low threshold |r λ |, expressed in
Z ( x , σ f ) = { ln ( | λ 1 ( x , σ f ) | ) + c λ 1 < r λ 0 else
(2)
where c is a constant. Note that λ1 is negative in our case, with r λ = 2 π σ f . This response equation provides a vessel discriminant function at scale σ f . To integrate multiscale responses corresponding to different scales at each location x, we select the maximum response across multiple scales, i.e.
Z ( x ) = max ( Z ( x ; σ f ) ) , σ f = σ 1 , ... , σ n
(3)

where σ1 and σ n are the minimum and maximum scales, determined according to the range of vessel diameters in the angiographic images. In this paper, 14 scales were used. For convenience, the resulting vessel features are normalized to 0[1] in this paper. Vessel feature map Z ( x ) (Figure 2 (a3-a4)) may also be binarized for extracting the vessels, and may further be thinned to obtain the skeletal lines representing the vessels.

The vessel direction map is estimated using the best detected scale. First, the maximum of absolute eigenvalue |λ1| is obtained at each location x by | λ o p t ( x ) | = max σ f | λ 1 ( x ; σ f ) | . Assuming that the eigenvector corresponds to λ opt (x) is v opt (x), the direction orthogonal to v opt (x) is then selected as the vessel direction at x. For convenience, we again use v 2 to represent the estimated vessel direction map later in this paper.

Therefore, using this step of multiscale analysis, we obtain three separate maps from the original angiographic images, namely, a vessel feature map, a vessel direction map, and a map of skeleton lines. Next, an initial seed point x 0 is automatically detected from a point within the thinned vessel skeleton-lines, as well as with maximum vessel feature intensity Z ( x 0 ) . The determined initial seed point is generally located near the root of the arterial tree or the main trunk of the arterial tree. The data above are used for the identification of consequential vessel pattern.

Probabilistic tracking of coronary artery

A dynamic mechanism was designed for our tracking framework. First, a list for source points (called source-list) was maintained to keep the initial seed point, as well as the vessel bifurcations that were detected during the vessel tracking procedure. In addition, a second list of crossings (called crossing-list) was designed to keep the possible crossings of two vessels. Then, a sampling disk with about half the size of the vessel diameter was used in the PTO to automatically track the artery segments (or branches) one at a time, beginning from the seed point and stopping at the vessel termination. During the tracking, the proposed SPD functioned simultaneously to report the local vessel pattern as bifurcation, crossing, segment, or termination. All detected vessel bifurcations were stored one at a time at the top of the source-list and were later used as subsequent starting points for tracking their respective artery segments (or branches). Finally, the vessel PTO and the SPD stop tracking and identifying vessels once all vessels have been tracked and registered.

Tracking criterion

The tracking mechanism was designed based on the criterion of continuity properties of luminance, position, and diameter of the general coronary artery in the angiographic images. The appropriate incorporation of these properties into the PTO enables successful management of the cases of vessel gap (caused by nonuniformity of the contrast medium or inappropriate thresholding in the preprocessing step) and arterial lesion (such as stenosis or aneurysm).

Description of PTO

Our PTO performs automatically along the coarse vessel path extracted from the vessel feature map. It begins or stops tracking an artery segment (or branch), depending on the structure pattern detected at current artery element e t . The tracking problem at simple artery positions, with no bifurcations, crossings, or terminations, is to determine the next artery element e t+1in the vessel path. Each artery element e t has four associated attributes, namely position x t , vessel feature Z ( x t ) , direction vector ϕ t = v 2(x t ), and vessel diameter d t . These attributes guide the vessel tracking algorithm. For efficiency, a sampling disk centered at the current artery element e t was used to search for a set of candidate elements { ( x t + 1 i , Z ( x t + 1 i ) , ϕ t + 1 i , d t + 1 i ) | i [ 1 , ... , N s ] } from the overlapping region between the sampling disk and the binarized vessel feature map. From artery element t to t+1, the probabilistic tracking can be formulated by a maximum posterior probability (MAP), expressed as
x t + 1 = arg max i [ 1 , ... , N t ] P ( x t + 1 i | x t )
(4)
The sequent modeling was based on the aforementioned continuity properties of the artery vessel. During tracking vessel length [Figure 3(a)], three local terms were used to constrain the candidates based on vessel feature continuity, longitudinal direction continuity, and diameter continuity. The first term, Z ( x t + 1 i ) , was acquired from the vessel feature map at candidate position x t + 1 i . The second term to some extent, kept the vessel tracking from being interrupted by noise, pseudo contours, and crossing, although more assurances should come from the structure pattern detection result, as described in aforementioned sections. Given vector V t + 1 i pointing from current position x t to candidate point x t + 1 i , a deflection of vessel direction θ t + 1 i = V t + 1 i ϕ t [Figure 3(a)] yields a means of the longitudinal continuity. In the case of vessel crossing or paralleling with each other (especially for thin vessels), the third term of diameter continuity enables the tracking of the current vessel segment/branch to proceed without being affected by neighboring vessels. Given d t + 1 i and d t as the vessel diameters, respectively at candidate point x t + 1 i and current artery element x t , increment d t + 1 i = | d t + 1 i d t | yields the measure of diameter continuity. Altogether, the estimation of the optimal candidate gives emphasis to the continuity properties, which is also the basic idea for probabilistic tracking. In the following paragraph, we provide further discussion to mathematically model.
Figure 3

Show the probabilistic tracking and diameter measurement along vessel path. (a) is the description of the probabilistic tracking. (b) depicts the estimation of vessel diameter using the vessel feature profile of transversal line C t + 1 i .

First, the relationship between the sampling instance and the posterior probability can be described as
( Z ( x t + 1 i ) ; θ t + 1 i d t + 1 i ) | x t P ( x t + 1 i | x t )
(5)
All feature measures in the left of the abovementioned expression are normalized before participation in the estimation of probability. In particular, maximal vessel feature Z max and maximal diameter d max are determined from all N s candidate points in the sampling disk, and are used to normalize the respective measures. Note that in practice, estimated vessel direction ϕ t can be reversed; thus, it should be considered in the normalization of θ t + 1 i . Given normalized measures Q Z ( x t + 1 i ) = | 1 Z ( x t + 1 i ) / Z max | , Q θ ( v t + 1 i ) = 1 π min { | arg θ t + 1 i π | , | arg θ t + 1 i | } , and Q d ( x t + 1 i ) = d t + 1 i / d max , their probabilities can be respectively written as:
q Z ( x t + 1 i ) = 1 2 π σ Z , t exp ( Q Z 2 ( x t + 1 i ) 2 σ Z , t 2 )
(6.1)
q θ ( x t + 1 i ) = 1 2 π σ θ , t exp ( Q θ 2 ( x t + 1 i ) 2 σ θ , t 2 )
(6.2)
q d ( x t + 1 i ) = 1 2 π σ d , t exp ( Q d 2 ( x t + 1 i ) 2 σ d , t 2 )
(6.3)
where σ Z,t , σ θ,t , and σ d,t are the standard deviations, which can be estimated by calculating the variances of the corresponding measures in all N s candidates. For example, σ Z,t is estimated as
σ Z , t = 1 N s i = 1 N S [ Q Z ( x t + 1 i ) Q ¯ Z ] 2
(7)
where Q ¯ Z is the average measure obtained by Q ¯ Z = 1 / N s i = 1 N s Q Z ( x t + 1 i ) . By integrating all distribution functions in Eq. (6), the optimal estimation of artery element e t+1can be obtained by
e t + 1 = x t + 1 = arg max x t + 1 i , i = 1 , ... , N s ( q z ( x t + 1 i ) q θ ( x t + 1 i ) q d ( x t + 1 i ) )
(8)

PTO uses adaptive parameters adjusted according to Eq. (7) based on all candidates in the sampling disk. Thus, it can successfully estimate the artery elements along the vessel path even in cases of gaps and arterial lesions. The gap is generally regarded as the region of a longitudinal narrow crack with a width less than the vessel radius, whereas the arterial lesion presents as a short or long component of the longitudinal region with large variations in vessel diameter and, at times, weak luminance. During the vessel tracking, the candidates running into the gaps can be crossed out by the continuity terms, thus enabling PTO to continue working. On the other hand, in case candidates run into the regions of the artery lesion, the PTO still carries on, even with large variations in vessel diameter.

Artery diameter estimation

Considering that the blood vessels have measurable abnormalities in diameter, intensity, and tortuosity, the artery diameter of candidate x t + 1 i should be estimated adaptively by analyzing its 1-D vessel feature profile orthogonal to the vessel direction. All candidates sitting along the sample profile have the same vessel diameter estimated (as described next). For example, in Figure 3(b), another candidate x t + 1 j sitting at the peak of the profile has the same vessel diameter as x t + 1 i . The profile can be represented by transversal line C t + 1 i with a single pixel width, as shown in Figure 3(b). Candidates x t + 1 i and x t + 1 j have the same profile, and thus the same diameter.

The feature values in the particular profile can be obtained from vessel feature map Z . The vessel diameter on this profile can be measured through the convolution of the normalized second derivative of the 1-D Gaussian function with the profile. Standard deviation σ in the 1-D Gaussian function should be adjusted during the vessel diameter estimation. Standard deviation σ that results in the largest response among all convolution results with the use of different standard deviations is finally selected as the estimated diameter for all candidates sitting on the same profile. Using this method, we can estimate vessel diameter d t + 1 i for candidate x t + 1 i .

This vessel diameter estimation process proceeds automatically without user interaction. The relevant diameter estimation results are demonstrated in the experimental section.

Inferring structure pattern of artery vessel

When the PTO works on current artery element x t , the SPD should also work in parallel to report the structure pattern to guide the PTO. Our SPD performs in the following manner. First, a circular pattern detector is centered on current artery element x t . It then acquires the circular cross-section profile from vessel feature map Z . Subsequently, it captures local vessel measures, including the number of peaks in profile (M t ), the azimuth vector directing from the detector center to every peak (p m ), and the Gabor response on each peak (γ o.m ). Finally, it returns the structure pattern inferred by a fuzzy inferring operator using the fuzzy membership degrees of multi-feature fuzzy subsets.

Detecting the peaks of the vessel-feature profile

The pattern detector samples of vessel features Z along the circular cross-section profile are centered at current artery element x t with radius r t = α·d t /2. Generally, the range of parameter α in α [1,1.5] is preferred for better measurement because in this case, the radius of circular pattern detector is a little bit larger than half of vessel diameter d t of current artery element x t ; thus, the structure around x t can be adaptively and completely captured. The resulting vessel feature profile is called circular template T (x,r,t), which consists of a fixed number (e.g., 150) of vessel features sampled over the circumference of a circle (x t , r t ). The same number of vessel features Z is sampled to make the SPD adaptive to different sizes of vessel diameters. Therefore, even for small vessels, the same number of points is also sampled. Because vessel feature map Z is defined in discrete domain, the bilinear interpolation is used for sampling features from Z when generating circular template T (x,r,t).

The sampling process is demonstrated in Figure 4, in which the vessel feature profiles are sampled by the circular templates (T (x,r,t)) in 4(a2)-4(c2). The typical profiles are shown in 4(a3) - 4(c3). To detect the peaks in each profile, first-order derivative dT (x,r,t)of T (x,r,t)is computed through the circular convolution of T (x,r,t)with the first-order derivative of the 1-D Gaussian filter. Figures 4(a4) - 4(c4) shows the magnitude profile of the resulting first-order derivative |dT (x,r,t)|. The local minima of the derivative magnitude correspond to the maxima or minima of the vessel feature profile along T (x,r,t). The peaks of the profile corresponding to maxima can be detected according to the derivative signs. For example, a peak of the profile corresponds to the sign transition of the first-order derivative from positive to negative. In this manner, we can detect all peaks along with their locations [e.g., z{p m = (x m , y m ), m = 1,...,M t }], where M t is the total number of peaks.
Figure 4

Using the circular pattern detector to characterize the vessel structures, in the cases with bifurcation, paralleling, and crossing patterns. (a1)-(c1) are original images and the respective regions of interest. (a2)-(c2) are the corresponding direction maps overlaping with circular pattern detectors. (a3)-(c3) are the corresponding vessel-feature profiles on the circular template T (x,r,t). (a4)-(c4) are the corresponding magnitude profiles of the first-order derivative dT (x ,r,t).

The resulting peaks provide the SPD vessel with several important properties, such as their positions, vessel features, vessel directions, as well as their relationships with the center of pattern detector (x t ). To better infer the pattern of the vessel structure, more relevant metric features, such as the Gabor features described in the next section, should be extracted.

The Gabor responses at locations of the peaks

Considering that the peaks of the circular profile provide useful information for structure identification, we first discuss how to infer the case in which two points (apart less than a few pixels) are located at the same vessel axis line. This case may also be called "co-vessel." To detect whether peak point p m = (x m , y m ) and detector center (O d ) are located in the same segment of the vessel, the Gabor filter operator [14, 15] was used for local description. As a function of orientation, frequency, and scale, the operator can characterize the dominant attributes at peak point p m = (x m , y m ). Note that this responds strongly to the isolated vessel strip, but exhibits a weaker or no response to the noise.

Given the relationship between the Gabor filter scale and the frequency with 2σ = 2πf = T, several masks by the Gabor operator can be defined by the real component of the Gabor function with a size of T × T. The masks of different scales are shown in Figure 5(b). Thus, to characterize local attributes, such as direction and "co-vessel" at a peak point, we merely need to compute the responses of several Gabor masks, and retain only the related parameters corresponding to the maximum response. Although none of these parameters are invariant to common image transformation, such as scaling and rotation, only one set of parameters corresponds to maximum response g max at the considered position. Given vector V o,m directing from the detector center to the peak position (x m , y m ), the following response measure γ o,m represents the degree with which peak p m = (x m , y m ) belongs to the same vessel axis line of the center point of the pattern detector
Figure 5

The feature integration, Gabor masks, and structure detection near artery bifurcation. (a) illustrates circular pattern detector and feature integration. (b) is an examples of five Gabor masks, and their parameters, where (b1) T 1 = 5, θ 1 = π/2, f 1 = 0.2; (b2) T 2 = 9, θ 2 = π/2, f 2 = 0.111; (b3) T 3 = 13, θ 3 = π/2, f 3 = 0.077; (b4) T 4 = 9, θ 4 = 3π/4, f 4 = 0.111; (b5) T 5 = 5, θ 5 = 3π/4, f 5 = 0.2. (c) illustrates the problem of the uncertainty of the node positions.

γ o , m = 1 g max arg  max σ | G V o , m , σ / π , σ ( x m , y m ) |
(9)

Where 0 < γ o.m < 1. With an image sized 512 × 512, the size of the mask (σ) is equal to the vessel diameter, which ranges from 3 to 14 pixels. Therefore, the resulting measure γ o,m reflects the relationship between peak p m and the detector center. For example, the resulting measure γ o,m at peak p m is lower if mask orientation V o,m does not accord with the local vessel direction, that is, there is no "co-vessel" existing between the peak point and the detector center.

The multi-feature based fuzzy inferring function (FIF)

Four types of patterns can be inferred, namely, vessel termination (i = 1), artery segment (i = 2), bifurcation (i = 3), and crossing point (i = 4), where variable i denotes the index for the ordinal structure pattern of the artery vessel. With effects from noise and pseudo edge, the number of peaks in each of these four types of patterns might be higher than the actual one. For example, by comparing the circular cross-section profile of different cases (Figure 4), we can observe that more than two peaks might occur in the profile of the artery segment, and more than three or four peaks might occur in the profile of the bifurcations or crossings, respectively. Without loss of generality, the circular pattern detector detects the number of peaks from the ordinal structure pattern according to M t ≥ 1, M t ≥ 2, M t ≥ 3 and M t ≥ 4 [Figure 5(a)]. Thus, the number M t can be regarded as a pseudo pattern of the artery structure. On the other hand, each peak p m has two metric features important to pattern identification: vessel feature Z ( p m ) at peak p m , and response measure γ o,m of Gabor mask function G V o , m , σ / π , σ ( x m , y m ) . The relationship between the pattern detector and the metric features for a typical artery bifurcation is shown in Figure 5(a).

In the following section, the fuzzy subsets are defined for all metric features, as well as for number of peaks M t . They will then be integrated for pattern inferring of point x t .

Definition 1: In the region Ω = { p m , Z m , γ o , m , M t | 1 m M t ; Z m , γ o , m [ 0 , 1 ] } , given that fuzzy subsets F m , Z and F m,γ correspond to the metric features of peak p m , and fuzzy subset F M t , x t i corresponds to the i-th pseudo pattern at current position x t , the membership degrees of these fuzzy subsets can be defined by the standard Gaussian function as
μ m , Z = exp ( ( Z max Z m ) 2 2 σ Z 2 ) ; μ m , γ = exp ( ( γ max γ o , m ) 2 2 σ γ 2 ) ; μ M t , x t i = { exp ( ( M t i ) 2 2 σ M t 2 ) M t i 0 else
(40)

where. Z max = max m Z m ., γ max = max m γ o,m , and m = 1,...,M t . The above expressions, μ m , Z and μ m,γ , represent the respective degrees of membership of two metric features at current position x t , which are used to ensure that each peak p m is located at the center of the vessel. Meanwhile, μ M t , x t i provides the membership degree of the pseudo pattern based on the number of peaks. According to the abovementioned membership degrees, the structure patterns of current position x t can be inferred by the following definition.

Definition 2: In the same region Ω, given J i (x t ) the fuzzy subset of the structure pattern belonging to the i-th pattern at current position x t , the relationship between fuzzy set J i (x t ) and the fuzzy sets of the metric features is given by
J i ( x t ) = { I m = 1 M t F m , Z } I { I m = 1 M t F m , γ } I { F M t , x t i }
(11)
The membership degree of J i (x t ) is given by:
μ J i ( x t ) = ( m = 1 M t μ m , Z ) ( m = 1 M t μ m , γ ) μ M t , x t i
(12)
where it is shown that the direct inferring of the local vessel structure can be hardly implemented using μ J i ( x t ) > T α with threshold T α . Therefore, the optimal pattern inference should be based on the principle of maximum membership degree and can be expressed as
μ J i 0 ( x t ) = i = 1 4 [ μ J i ( x t ) ]
(13)

which infers the artery structure at current position x t belonging to the i 0-th pattern.

When the detector center is near the vessel node (e.g., bifurcation or crossing), the resulting node position can be ambiguous; for example, a few of the redundant positions may possibly be the node positions because of the similarity of their metric features to those of the actual node. To solve this problem, a second local maximum operation is performed to obtain a preferable node position. Given n redundant positions near the bifurcation, with the membership degrees { μ J i 3 ( x t 1 ) , ... ,  and  μ J i 3 ( x t n ) } , the preferable one should be located at
x = arg max x t k , k = 1 , ... , n ( μ J i 3 ( x t k ) )
(14)

In Figure 5(c), the numbers 2 and 3 represent the two structure patterns determined in the corresponding locations by the proposed pattern detector. Pattern detection result is used to guide the vessel tracking process.

Manipulation of vessel tracking process and organiztion of tracking results

Appropriate termination conditions for vessel tracking are important to achieving better tracking performance. Our vessel tracking algorithm stops tracking the current artery segment or branch once one of the following conditions is satisfied: (1) a new vessel bifurcation is detected, which should be dealt with appropriately; (2) vessel feature Z ( x t + 1 ) at position x t+1along the vessel path is below a certain threshold, such as when a value is slightly higher than the background (this threshold is set as 0.05 in this paper); (3) the current path overlaps with a vessel path previously detected, except when a crossing point on the vessels has been detected and recorded in advance; (4) the source-list is empty and no new node is found; and (5) the new point is outside the image field. In all cases except for vessel crossing, the last valid vessel point is marked as a termination point, and the algorithm starts tracking the other vessel with the seed point placed at the top of the source-list.

The artery segment is a vessel length with two end-points at two nodes, or at an origin point and a node, whereas the artery branch is a distal vessel with a distal point. Either the artery segment or branch is composed of artery elements. Considering these artery elements arranged at intervals (or with "look-ahead distance") of less than the radius of the sampling disk, polylines are used to connect the artery elements by linear interpolation. The generated polylines are treated as an approximation of vessel centerlines, and have high probability inside the binarized vessel structures.

Results

The proposed algorithm was tested on both simulated images and real image sequences with stenosis vessel and aneurysm. To obtain images with acceptable quality, the cranial and caudal angles, as well as the X-ray dose, were adjusted so that the X-rays travel a moderate distance. Our vessel tracking algorithm was implemented in MATLAB 7.0 on a Pentium IV PC (with CPU 2.8 G and 512 M memory).

Our algorithm was first tested on various simulated artery vessel images to adjust its parameters. Actual images were used to evaluate the actual performance of the algorithm by comparing the algorithm-segmented structure and vessel length of the entire artery network with those produced by expert cardiologists. The tracking algorithm has several distinguishing characteristics: (a) it is robust to the starting position detected; (b) it can detect and identify a larger area of the artery tree without any manual intervention; and (c) it can accurately distinguish the narrows resulting from the gaps in the vessel feature map or those caused by real vessel stenosis.

In our algorithm, total computational time T consisting of two major components T 1, T 2, and one minor component T 3, were spent on preprocessing (36.7 s), tracking (15.7 s), and connecting artery elements (1.4 s), respectively.

General effects of vessel segmentation are listed in the first column of Table 1, where the results of right coronary artery (RCA) and the left coronary artery (LCA) in the real angiographic images are compared separately. In Table 1, L and L ¯ denote the length of algorithm-segmented vessels (in pixels) and the total vessel length delineated by expert cardiologists, respectively. The percentage of the detected artery segments (DAS %) and the extracted vessel length (EVL %), as well as the identification ratio (IR %) of the arterial structure was computed to evaluate the performance of our algorithm. In particular, the IR % is defined as the percentage of the correct parts of vessel length in the entire length of algorithm-segmented vessels. As openly reported in literature, the DAS % produced by Haris et al. [16] reached 90%, while the EVL % produced by Schrijver et al. [4] was less than 78%. On average, our algorithm can reach 94.7% for DAS % and 88.1% for EVL %.
Table 1

Statistical results on artery tracking and structure identification in real angiographic images.

Images, Num.

EVL(%)

IR(%)

L ¯

L

Average identification rate of bifurcation ( N ¯ b / W ¯ b , i = 3 )

Average identification rates of different patterns (%)

     

b = 1

b = 2

b = 3

b = 4

b = 5

Total

i = 1

i = 2

i = 3

i = 4

SAT, 1

99.8

99.6

950

948

1/1

2/2

2/2

3/4

0/0

8/9

100

100

88.9

75

LCA, 128

87.2

98.3

1868

1629

1.00/1.00

2.00/2.00

7.11/7.22

5.13/6.12

0.32/1.34

15.56/17.68

98.9

98.7

88.0

72.8

RCA, 32

91.7

99.1

947

868

1.00/1.00

3.10/3.15

2.65/2.70

0.30/0.33

0.01/0.02

7.05/7.20

98.3

99.5

97.9

77.5

Mean

88.1

98.5

1677

1477

1.00/1.00

(100%)

2.22/2.23

(99.6%)

6.22/6.32

(98.4%)

4.16/4.96

(83.9%)

0.26/1.08

(24.1%)

13.86/15.58

(89.0%)

98.8

98.9

90.0

73.7

Note: L and L ¯ denote the length of algorithm-segmented vessels and the total vessel length, respectively. N b /W b denotes the rate of bifurcation identification., where N b denotes the number of correctly detected bifurcations for the b-th level vessel branch, while W b denotes the total number of the detected bifurcations.

Other detailed experiments and performances are introduced in the next section.

Arterial vessel simulation and algorithm evaluation

Simulated data were created for the arterial segment [Figure 6(b)] and the simulated arterial tree (SAT) [Figure 6(c)]. In particular, various levels of Gaussian smoothing were applied to simulate the spatial blurring in an X-ray projection. The SAT in Figure 6(c) was generated to simulate artery characteristics, such as curving, overlapping, and bifurcating. For these simulated images, all information related to vessel segments, centerline, and local structure are known beforehand. We can analyze the algorithm characteristics given that various levels of white Gaussian are added to these images.
Figure 6

The simulations of arterial vessels: (a) the local coordinate system of vessel simulation; (b) the simulated arterial segments; (c) the SAT.

The performance of our algorithm can first be evaluated using the simulated images, such as vessel segments A to F. Given an initial seed point at the left end of the vessel segment in A, B, E, and F of Figure 6(b), the PTO can automatically track vessels and measure their diameters. The measured diameters along the vessel are consistent with varying morphology, such as stenoses and aneurysms. The resulting mean errors are listed in Table 2, where L x denotes the length of the simulated arterial segment, e 0 and e d denote the mean error of the resulting vessel centerlines and diameters, respectively. Under zero-level Gaussian noise, the mean errors of estimated centerlines and diameters of stenosis or aneurysm vessels are no more than 0.02 and 0.04, respectively. In case of close proximity (such as C) and local gap (such as D), the tracking results depend on distance d P between two neighboring vessels and gap-width w G , respectively. Suppose that r 1 and r 2, and r D denote the radii of two neighboring vessels and the radius of the detector, respectively. To enable the algorithm to operate through the vessel despite influence from the neighboring vessel and the gaps, the limiting distance and gap-width would be d P ≥ (r 1 + r 2)/2 and w G < r D , respectively. The entire length of the skeleton lines of the SAT is 948 pixels, and is entirely involved in the testing.
Table 2

Parameters and mean errors in testing of the phantoms of arterial segments.

 

A

B

C

D

E

F

S V

1

1

1

1

1

1

S B

0

0

0

0

0

0

σ B

3.2

3.2

3.2

3.2

4.0

3.2

L m

13

13

13, 6

13

9

15

L r

13

5

13, 6

13

18

10

L x

120

120

120, 125

120

120

120

σ w

0

0

0

0

0

0

e 0

0.01

0.01

0.01,0.02

0.01

0.02

0.02

e d

0.02

0.04

0.05,0.05

0.02

0.04

0.03

The results of structure identification on the SAT and real angiograms are provided in the following sections.

Experiments on structure pattern inferring

Two experiments were performed on simulated and real images to test the capability of our algorithm to infer the structure patterns in the artery tree. The performances on identification of all four different patterns, as well as separate performances on the detection of bifurcation are reported.

The first experiment was conducted on the SAT image. Based on results of estimated membership degrees μ J i ( X t ) in Eqs. (13) and (14), the attributes of all interesting points in Figure 6(c) can be well characterized by our algorithm. Overall, the captured membership degrees corresponded well to four different structure patterns, such as the membership degree distribution between N1-N4 in Figure 7. For vessel overlapping cases, such as local superposition between the two vessel segments N2-N4 - N5-D and N1-N4-N5-C, the structure ambiguity (between N4 and N5) might be difficult to resolve using only a single-view projection, although it might be resolved by observing vessel motion or radius information from sequence images as done by human observers while reading cine-angiograms. In short, our multi-feature based fuzzy inferring method is capable of accurately estimating the structure patterns from 2D images.
Figure 7

The fuzzy memberships measured along the vessel axis, starting from point N1 to N4 in the SAT.

The second experiment was performed on actual angiogram images selected from the angiogram sequences. The related results are provided in the second and third columns of Table 1. The lower level of the arterial structure can be inferred with higher IR %, and the recognizable vessels can reach up to the fifth level at best. In practice, the cardiovascular surgeon is often interested in the vessels below the fourth level; thus, the performance of our algorithm is promising. The second column of Table 1 also displays the rate of bifurcation identification (N b /W b ). Here, N b denotes the number of correctly detected bifurcations for the b-th level vessel branch, while W b (≥N b ) denotes the total number of the detected bifurcations, including the incorrectly detected ones.

Experiments on vessel tracking and measurement

We also evaluated the respective contributions of the probability terms of the continuity properties given in Eq. (6) to the vessel tracking. By rewriting Eq. (8) as x t + 1 = f ( q θ , q Z , q d ) , three different types of PTO, namely, x t + 1 = f ( q θ ) , x t + 1 = f ( q θ , q Z ) , and x t + 1 = f ( q θ , q Z , q d ) , can be constructed. Their respective results are compared in Figure 8, where the paths traced in the same sub-image were obtained by the same PTO. Figures 8(b1) and 8(b2) show that the direction constraint term in Eq. (6.2) can effectively compel the sampling disk and the pattern detector to move ahead, but it can often be lost in edges, pseudo contours, neighboring vessels, or even noise. By incorporating the vessel feature term, that is, using x t + 1 = f ( q θ , q Z ) , the vessel tracking result can be improved, as shown in Figures 8(c1) and 8(c2). Using all terms together, that is, using x t + 1 = f ( q θ , q Z , q d ) , the vessel tracking result can be further improved; this is reflected as the smooth and accurate tracking paths detected in Figures 8(d1) and 8(d2). In addition, if only the vessel feature or diameter term is used, the tracking model almost stops near the initial seed position. These visual results are not shown here.
Figure 8

Comparison of the results of vessel-tracking using the PTO based on different continuity terms. (a1) & (a2) two original images. (b1) & (b2) using the PTO simply based on q θ . (c1) & (c2) using the PTO based on q θ and q z . (d1) & (d2), using the PTO based on complete terms. (e1) & (e2), measuring the vessel diameters, respectively along the vessel paths in (d1) and (d2).

Figures 8(c1) and 8(c2) also yield a typical case in which the algorithm can successfully manage the close proximity or crossing of vessels, as well as varying vessel curvatures, which validates a number of attractive characteristics of the PTO. For example, the algorithm can run across the local vessel gap, proceed along the current vessel branch in spite of crossing, and successfully manage arterial stenoses and aneurysms.

The diameter estimation errors at the artery elements are less than 0.23 pixels on average. The results of the axis line estimation along the artery segments in the SAT image indicate that more errors appear near the nodes, with 1.36 pixels on average. These results show that the axis line estimation error at the nodes is controlled within a limited extent. For actual angiograms, we compared the automatically extracted arterial axis lines with the standard arterial centerlines delineated by an expert cardiologist. We also obtained good results, with a mean error of 0.17 pixels and the maximum error close to half of the respective vessel radius. These errors are produced by noise, as well as large curvatures or abnormal morphologies at the locations of vessel nodes.

Vessel stenosis

We tested a simple method to determine the position and severity of coronary artery stenosis in the given images, which can be used as an additional evaluation criterion for our tracking algorithm. Our algorithm was able to detect the arterial stenosis. An example in which the diameter of arterial stenosis is tracked and measured is shown in Figure 9. According to stenosis evaluation using the grading scheme of the visual analysis of the American Heart Association [17], %-LENGTH stenosis is estimated as %-LENGTH stenosis = (1-L r /L m ) × 100%, where L r and L m are the estimates of the minimal diameter of the stenotic part of a vessel and the mean diameter of the nonstenotic part, respectively. During tracking and detection, the vessel length of interest begins or terminates when %-LENGTH is 25% higher or lower than the minimum classified value. The vessel length of interest can be described by the rectangular sash in Figure 9.
Figure 9

Tracking and measuring vessel width using the real angiographic images with arterial stenosis. (a) Original angiogram with arterial stenosises; (b) Location and measurement of the stenosises.

In the testing, 35 angiograms with 40 stenosis positions were used to evaluate the correctness of stenosis identification, in which all the standard vessel stenoses, as well as their severities were provided beforehand by the expert cardiologist. As a result, 33 stenosis objects are correctly detected and the number of false negatives and false positives are seven and one, respectively. The average correctness rate of stenosis severity identification reaches 68.3%, in comparison with the given stenosis severity levels.

Discussion

In this paper, we analyzed image data and vessel structures within X-ray projections. A structure distinction exists between the crossing and bifurcation of vessel segments. Two crossing vessel segments have no spatial structure connection, whereas vessel bifurcations (which in general produce three vessel segments) have a structural connection among them. Vessel bifurcation enables low-level vessel segment separation to two high-level ones at the bifurcation point where these vessel segments have different features in their diameters and direction continuity. However, in single-view X-ray projection, two spatial overlapping vessel segments may appear connected to each other at a crossing-point, where the vessel segments maintain their respective features of diameter and direction continuity near the point. The abovementioned vessel continuity of different cases can address, in part, 2D structure ambiguity.

In our algorithm, three parameters can be adjusted to achieve better segmentation results. These include threshold value r λ in the vessel function Eq. (2), the radius of sampling disk r S , and the radius of pattern detector r D . These parameters can be adjusted according to the overall performance integrated from the three percentages, DAS %, EVL %, and IR %. For instance, parameter r λ may be optimized according to
lim r λ r λ ' Δ D A S ( r λ ) + Δ E V L ( r λ ) + Δ I R ( r λ ) Δ r λ ( D A S ( r λ ) + E V L ( r λ ) + I R ( r λ ) ) 0

The optimization for parameters r S and r D can be conducted similarly. Using the simulated images with different levels of noise, we can determine the optimal values for these three parameters. In the absence of noise, r λ is computed as r λ = 2 π σ f

; parameters r S and r D turn out to be equal, that is, r S = r D , and are both equal to current diameter d t of the vessel under tracking.

In the experiments, the time spent in tracking depends on the length of the detected vessel, look-ahead distance s t , and vessel diameter d t . Look-ahead distance s t is determined by the diameter of the sampling disk, which is designed to vary with the diameter of current vessel d t . The ranges of these parameters are 1 ≤ s t ≤ 7 (pixels) and 3 ≤ d t ≤ 15 (pixels), respectively, which agree with the required relationship 1 ≤ s t d t /2. Similarly, more artery elements can result in more computational time.

The rate of bifurcation identification (N b /W b ) in the second column of Table 1 is related to EVL % and independent of IR % because the IR% is calculated only by the detected vessel points. In addition, the IR % of bifurcation and crossing are relatively lower than that of other structure patterns, which is mostly caused by the structure ambiguity in the 2D image acquisition. According to the testing on SAT, the error ratio is correlated with the effects of image noise and the diameter of the detector.

Based on the results on 40 stenosis objects, the application for the diagnosis of stenosis severity levels do not appear to be very reliable because of the use of single-view projection images. The accurate quantification and analysis of stenotic lesions is actually very complicated and has been elaborately discussed by Sato et al. in [18]. The importance of designing an acquisition system for acquiring "good" images, in which stenotic lesions can be displayed clearly, rather than simply analyzing the given images by common clinic systems, is therefore important.

Conclusions

We developed a vessel tracking framework for the segmentation and measurement of the arterial tree in angiographic images, and applied it to actual clinical data. Based on the various evaluation results using both simulated and actual data, our algorithm demonstrates a very impressive performance in tracking and measuring artery trees. The particular advantages include efficient handling of vessel nodes, such as bifurcation or crossing; distinguishing between bifurcations and crossings, and automatically locating their positions; and automatic adaptation to varying vessel diameters in the coronary artery. These advantages were made clear with the introduction of our tracking strategies.

Our algorithms are able to locate stenosis or aneurysm positions in a number of given images. Naturally, the actual 3D morphology cannot be simply inferred through the 2D structure patterns. The use of 2D angiograms can pose ambiguity in determining vessel connectivity, particularly in cases with vessel overlapping. This indicates the importance of extending our algorithm to 3D applications to achieve certain diagnostic functions. Moreover, although our fuzzy inferring function is found relatively effective in integrating spatial characteristics of centerline, orientation, diameter, and density along the coronary blood vessels, the performance of structure pattern identification might be further improved, if we can effectively employ advanced learning-based methods, such as support vector machine, to identify the vessel structure patterns. These ideas will be explored further in our future studies.

Declarations

Acknowledgements

This work was partly supported by the National Basic Research Program of China (Grant No. 2010CB732505) and the National Science Foundation of China (Project No. 60772120, and Project No. 60902103). We would like to thank the Cardiac Center of Beijing Chaoyang Hospital and the Radiology Department of Tianjin General Hospital for clinical assistance and providing coronary angiogram data.

Authors’ Affiliations

(1)
School of Biomedical Engineering, Southern Medical University
(2)
Department of Optic electronics engineering, Beijing institute of technology

References

  1. Kawarada O, Yokoi Y, Morioka N, Nakata S, Higashiue S, Mori T, Iwahashi M, Hatada A: Carotid stenosis and peripheral artery disease in Japanese patients with coronary artery disease undergoing coronary artery bypass grafting. Circ J 2003, 67: 1003–1006. 10.1253/circj.67.1003View ArticleGoogle Scholar
  2. Jian Y, Yongtian W, Yue L, Songyuan T, Wufan C: Novel Approach for 3-D Reconstruction of Coronary Arteries From Two Uncalibrated Angiographic Images. IEEE Trans Image Process 2009, 18(7):1563–1572. 10.1109/TIP.2009.2017363MathSciNetView ArticleGoogle Scholar
  3. Kirbas C, Quek FKH: A review of vessel extraction techniques and algorithms. ACM Comput Surv 2004, 36(2):81–121. 10.1145/1031120.1031121View ArticleGoogle Scholar
  4. Schrijver M, Slump CH: Automatic segmentation of the coronary artery tree in angiographic projections. Proceedings of ProRISC 2002, 449–464.Google Scholar
  5. Behrens T, Rohr K, Stiehl HS: Robust segmentation of tubular structures in 3-D medical images by parametric object detection and tracking. IEEE Transactions on Systems, Man, and Cybernetics - Part B, Cybernetics 2003, 33(4):554–561. 10.1109/TSMCB.2003.814305View ArticleGoogle Scholar
  6. Olabarriaga SD, Breeuwer M, Niessen WJ: Minimum cost path algorithm for coronary artery central axis tracking in CT images. In proceedings of MICCAI 2003, 2879: 687–694. LNCSGoogle Scholar
  7. Al-Kofahi KA, Can A, Lasek S, Szarowski DH, Dowell-Mesfin N, Shain W, Turner JN, Roysam B: Median-based robust algorithms for tracing neurons from noisy confocal microscope images. IEEE Trans Inf Technol Biomed 2003, 7(4):302–17. 10.1109/TITB.2003.816564View ArticleGoogle Scholar
  8. Wink O, Niessen WJ, Viergever MA: Multiscale vessel tracking. IEEE Trans Med Imaging 2004, 23(1):130–133. 10.1109/TMI.2003.819920View ArticleGoogle Scholar
  9. Tolias YA, Panas SM: A fuzzy vessel tracing algorithm for retinal images based on fuzzy clustering. IEEE Trans Med Imaging 1998, 17(2):263–273. 10.1109/42.700738View ArticleGoogle Scholar
  10. Lorenz C, Carlsen I-C, Buzug TM, Fassnacht C, Weese J: Multi-scale line segmentation with automatic estimation of width, contrast and tangential direction in 2-D and 3-D medical images. CVRMED-MRCAS'97 1997, 1205: 233–242. full_textView ArticleGoogle Scholar
  11. Sato Y, Nakajima S, Shiraga N, Atsumi H, Yoshida S, Koller T, Gerig G, Kikinis R: 3D multi-scale line filter for segmentation and visualization of curvilinear structures in medical images. Medical Image Analysis 1998, 2(2):143–168. 10.1016/S1361-8415(98)80009-1View ArticleGoogle Scholar
  12. Frangi AF, Niessen WJ, Vincken KL, Viergever MA: Multiscale vessel enhancement filtering. Medical Image Computing and Computer-Assisted Intervention - MICCAI'98 1998, 1496: 130–137.View ArticleGoogle Scholar
  13. Krissian K, Malandain G, Ayache N, Vaillant R, Trousset Y: Model based detection of tubular structures in 3-D images. Computer Vision and Image Understanding archive 2000, 80: 130–171. 10.1006/cviu.2000.0866View ArticleGoogle Scholar
  14. Moreno P, Bernardino A, Santos-Victor J: Model based selection and classification of local features for recognition using Gabor filters. Lecture Notes in Computer Science, Image Analysis and Recognition 2006, 4142: 181–192. full_textView ArticleGoogle Scholar
  15. Moreno P, Bernardino A, Santos-Victor J: Gabor parameter selection for local feature detection. Lecture Notes in Computer Science, Pattern Recognition and Image Analysis 2005, 3522: 11–19.Google Scholar
  16. Haris K, Efstratiadis SN, Maglaveras N, Pappas C, Gourassas J, Louridas G: Model-based morphological segmentation and labeling of coronary angiograms. IEEE Trans Med Imaging 1999, 18(10):1003–1015. 10.1109/42.811312View ArticleGoogle Scholar
  17. Austen WG, Edwards JE, Frye RL, Gensini GG, Gott VL, Griffith LS, McGoon DC, Murphy ML, Roe BB: A reporting system on patients evaluated for coronary artery disease. Circulation 1975, 51(4):5–40.View ArticleGoogle Scholar
  18. Sato Y, Araki T, Hanayama M, Naito H, Tamura S: A viewpoint determination system for stenosis diagnosis and quantification in coronary angiographic image acquisition. IEEE Trans Med Imaging 1998, 17(1):121–137. 10.1109/42.668703View ArticleGoogle Scholar

Copyright

© Shoujun et al; licensee BioMed Central Ltd. 2010

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Advertisement