Segmentation and Automatic Identification of Vasculature in Coronary Angiograms.

Coronary angiography is the "gold standard" for the diagnosis of coronary heart disease, of which vessel segmentation and identification technologies are paid much attention to. However, because of the characteristics of coronary angiograms, such as the complex and variable morphology of coronary artery structure and the noise caused by various factors, there are many difficulties in these studies. To conquer these problems, we design a preprocessing scheme including block-matching and 3D filtering, unsharp masking, contrast-limited adaptive histogram equalization, and multiscale image enhancement to improve the quality of the image and enhance the vascular structure. To achieve vessel segmentation, we use the C-V model to extract the vascular contour. Finally, we propose an improved adaptive tracking algorithm to realize automatic identification of the vascular skeleton. According to our experiments, the vascular structures can be successfully highlighted and the background is restrained by the preprocessing scheme, the continuous contour of the vessel is extracted accurately by the C-V model, and it is verified that the proposed tracking method has higher accuracy and stronger robustness compared with the existing adaptive tracking method.

1. Introduction

Cardiovascular disease is currently recognized as one of the most important chronic diseases leading to human death in the world. In recent years, morbidity and mortality from cardiovascular diseases continue to increase, ranking first among various diseases. Coronary angiography (CA) is a common and effective method for diagnosing coronary heart disease. It is regarded as the “gold standard” for the diagnosis of coronary heart disease and is widely used in clinical diagnosis [1].

Normally, human arteries and vessels are invisible under X-rays. However, by injecting X-ray impervious substances into the coronary arteries and then irradiating the coronary artery area with X-rays, the arteries and vessels can be visualized. To decide the treatment plan, doctors need to find the location and degree of coronary artery stenosis based on the image by themselves. Nevertheless, in this way, a large amount of repetitive work and subjective errors are inevitable. Thus, it is of great benefit to invent technologies to segment and identify vessels in angiograms. For this reason, many scholars have proposed various methods.

For many years, image segmentation is one of the focuses of image processing. Up to now, many segmentation technologies for vessels have been proposed. Based on the two characteristics of discontinuity between regions and similarity within regions, we can divide vessel segmentation technologies into three categories: boundary-based segmentation technologies [2-8], region-based segmentation technologies [9-11], and technologies combined with specific theories and tool segmentation [12-15]. Sahoo et al. [16] adopted the maximum entropy method and the gray threshold that maximizes entropy corresponded to the optimal segmentation threshold. Sato et al. [17] constructed a multiparameter similarity function for enhancing vessels by analyzing the properties of the eigenvalues of the Hessian matrix of spherical, tubular, and sheet-like structures at a certain scale. Based on the simplified Mumford-Shah model and the level set idea, Chan and Vese [18] proposed a new active contour C-V model to evolve the curve through the minimization of the energy function. Most recently, deep learning methods have also been widely used in the field of vessel segmentation. For example, Chen et al. [19] trained the 3D U-Net to perform three-dimensional vessel segmentation and achieved high segmentation accuracy.

Moreover, people have studied a variety of methods for vascular identification, such as multiscale-based methods [20-24] and tracking-based methods [25-29]. Among these methods, the tracking-based method has been proved to be very effective. It can detect coronary information based on the local response of angiogram without scanning the entire image. In the process of coronary artery extraction, the extraction result is unstable due to the manual setting of seed points. Aiming at this problem, Xiao et al. [30] proposed an automatic seed point acquisition method based on ridge point detection. These ridge points serve as seed points for adaptive tracking of the centerline of the coronary artery. Aylward and Bullitt [26] proposed a multiscale spatial centerline tracking algorithm based on ridge detection, which uses the eigenvalue decomposition of the Hessian matrix to extract the ridge. However, due to limitations in algorithm design and the effects of low image quality, noise, etc., the accuracy and robustness of these methods still have room for improvement.

Our main work and contributions are as follows: first, we designed a preprocessing scheme to increase the quality of the image and enhance the vascular structure. Then, we used the C-V model to achieve vessel segmentation. Finally, we proposed an improved adaptive tracking algorithm to realize automatic identification of the vascular skeleton, which achieved better effects than former methods according to our experiments.

This paper is organized as follows. In Section 2, we introduce our scheme of image preprocessing. Section 3 describes the active contour model to extract the vascular contour. Section 4 describes the details of our proposed improved adaptive tracking method. Section 5 presents the analysis and experimental results of testing the robustness and accuracy of our methods. Finally, conclusions are drawn in Section 6.

2. Image Preprocessing

The complex and varied configuration of the coronary artery structure, noise caused by various factors, artifact caused by the beating of the heart, and low contrast of terminal vessels make precise segmentation very challenging. Therefore, before the extraction of coronary artery structure, coronary angiograms should be preprocessed to enhance the vascular structure and suppress the background noise. In this paper, block-matching and 3D filtering (BM3D) [31] is used to effectively filter out noise. Unsharp masking (UM) [32], contrast-limited adaptive histogram equalization (CLAHE), [33] and multiscale image enhancement [34] are used to improve image contrast and highlight the vascular structure.

2.1. Block-Matching and 3D Filtering

BM3D is a 3D block-matching algorithm used primarily for noise reduction in images. Firstly, by the grouping technique of block-matching, image fragments are grouped based on similarity and are integrated into a three-dimensional matrix. Then, filtering is done on every fragment group. At last, the image is transformed back into its two-dimensional form and all overlapping image fragments are weight-averaged to ensure that they are filtered for noise yet retain their distinct signal. This algorithm can effectively remove image noise.

2.2. Unsharp Masking

The main procedures of UM algorithm are as follows: first, a passivated fuzzy image is generated after low-pass filtering of the original image. Then, the image with high-frequency components is obtained by subtraction of the original image and the fuzzy image. Finally, the high-frequency image is enlarged with a parameter and superimposed with the original image; that is, an image with enhanced edges is generated. The specific algorithm steps are as follows: where I(x, y) represents the gray of the pixel (x, y), represents the gray of the pixel (x, y) after low-pass filtering, and gmask(∙) generates the high-frequency component of the image where k is the enlarge coefficient and g(∙) generates the image with enhanced edges

Generate the smoothing result:

(2) Add the passivation template to the original image with a certain proportion:

2.3. Contrast-Limited Adaptive Histogram Equalization

As a variant of adaptive histogram equalization, the CLAHE method limits the contrast amplification to reduce excessive amplification of noise. In CLAHE, the contrast amplification in the vicinity of a given pixel is given by the slope of the transformation function, which is proportional to the slope of the neighborhood cumulative distribution function (CDF) and therefore to the value of the histogram. CLAHE limits the amplification by clipping the histogram at a predefined value before computing the CDF. This limits the slope of the CDF and therefore of the transformation function. It is advantageous not to discard the part of the histogram that exceeds the clip limit but to redistribute it equally among all histogram bins. The process of clipping the histogram is shown in Figure 1.

Figure 1

The clipping process of CLAHE.

2.4. Multiscale Image Enhancement

Frangi et al. [34] proposed the multiscale enhancement method based on the Hessian matrix of the image. In this method, the relationship among the eigenvalues, eigenvectors of the Hessian matrix, and the orientation of vascular structure are utilized, combined with the multiscale theory. Then, the vascular structure in the coronary angiogram is detected by constructing an appropriate vascular similarity function. At present, the method has become one of the most commonly used methods of multiscale enhancement.

The vascular similarity function is established as follows: where λ1 and λ2 are two eigenvalues of the Hessian matrix, ∣λ1 | ≤|λ2|, P is an arbitrary point in the image, σ is the scale parameter, R = |λ1|/|λ2|,, and R is the vascular structure enhancement factor, which is used to distinguish the globular structures from the tubular structures; S is the norm of the Hessian matrix; and β, c, and m control the overall smoothness of linear objects.

When the scale factor σ is consistent with the width of the tubular structure, the filtering result V(P; σ) gets the maximum value. By iterating the scale parameters σ, the values V(P; σ) under different scales are obtained, and the maximum value is taken as the actual output of the point P: where σmin and σmax are the minimum and maximum sizes of the vascular structure, respectively.

3. Vessel Segmentation

In this section, we will introduce the active contour model to extract the vascular contour of coronary angiograms.

Kass et al. [35] proposed the active contour model (ACM). This method converts the image segmentation problem into solving an energy minimization problem. The contour curve is the edge of the blood vessel when the energy function reaches the minimum. The active contour model is mainly divided into edge-based and region-based according to the different construction methods of the energy function. The most prominent advantage of the ACM is its resistance against strong noise.

The C-V model [18, 36, 37] is a representative region-based active contour model. The specific algorithm steps are as follows: where c0, c represent the average gray levels of the outside and inside areas of the curve C, respectively; L(C) represents the length of the closed curve C; S(C) represents the area of the inner area of C; and u, v, λ0, and λ represent the weights of items in the energy function.

Put forward the energy function:

(2) Introduce the level set method, set w(x, y) as a sign distance function with positive, negative, zero representing inside, outside, and right on the curve C, respectively:

Introduce the following H and δ functions:

Rewrite the energy function as a level set equation: where λ and λ0 are the iterative parameters in the C-V model, and their values affect the evolution rate of the curve C. When the curve C contains the segmentation target, the internal homogeneity of the curve C is low; thus, it is necessary to enlarge λ0 to accelerate the evolution of the curve C to the target, and vice versa.

(3) The energy minimization problem can be solved by minimizing the level set equation iteratively

The C-V model minimizes the energy function to obtain the evolution curve that approaches the edge of the blood vessels and finally segments the target. Compared with other methods, it has better effects on the continuous gradient.

4. Improved Adaptive Tracking

In this section, we will propose an improved adaptive tracking method, which is more robust and has fewer misjudgments in the tracking process, to automatically extract the skeleton of the coronary blood vessels.

4.1. Ridge Point Detection

Ridge point detection is important for seed point selection, blood vessel tracking, and bifurcation point detection. Ridge point is the local gray maximum point of the two-dimensional image. After multiscale enhancement, a ridge point of the blood vessel is usually located at the maximum point perpendicular to the direction of the blood vessel. The gradient of the local maximum point in the image is zero, and its Hessian matrix is negative [38]. Since the coordinates of image pixels are all integers, according to the principle of linear interpolation, if the point (ε, η) (x < ε < x + 1, y < η < y + 1) satisfies the following conditions: where ∇(x, y) is the gray gradient of point (x, y) and λ(x, y) are the eigenvalues of the Hessian matrix of point (x, y); then, (ε, η) can be considered as a local maximum point, and the pixel (x, y), as its approximate solution, is defined as a ridge point.

The ridge points may be misjudged due to the image noise caused by the uneven distribution of contrast agents and other factors. Thus, a gray threshold is used to screen out those misjudged ridge points. This method can effectively remove most of the ridge points outside the blood vessel.

4.2. Tracking Process

The tracking algorithm starts from a seed point and gradually tracks to the end of the vessel, extracting the blood vessel skeleton. We randomly select the seed point from the detected ridge points.

The initial tracking direction can be calculated from the gray information around the seed point. According to [38], take the seed point as the center and search for the gray maximum point P+ on the circle with radius d. P+ is the first point of forward tracking, the forward initial tracking direction u+ and angle θ+ can be expressed as

After obtaining the forward tracking direction, we search for the local maximum point P− on arc l(2π − θ+ − Δθ, 2π − θ+ + Δθ) centered in the opposite direction (2π − θ+) of the forward tracking angle θ+. The search area is shown in Figure 2.

Figure 2

Initial direction detection.

The backward direction of the initial trace u− can be calculated as

Tracking from the current point forward to the next point is the main step of this algorithm. The current tracking direction is determined by the direction from the previous point P to the current point P:

After determining the tracking direction, we search for the local maximum point P on arc l(θ − Δθ, θ + Δθ) and the following conditions should be met: where I(P) is the gray of P, N(P) is the number of tracking points around P, and I0 and τ are two thresholds. The first condition is to prevent the overtracking beyond the vessel area, while the second condition can avoid repeatedly tracking the vessel and being trapped in a local endless loop. If both conditions are satisfied, we continue to track from P. Otherwise, P is the endpoint of the vessel. We illustrate the tracking process in Figure 3.

Figure 3

Forward tracking.

Due to the noise and other issues mentioned above, a few tracking points may deviate from the center of the vessel. The tracking point can be adjusted to the center by center adjustment, which combines the blood vessel contour and tracking direction information. The specific steps are as follows: get the normal line of the vessel through the vertical direction of the current tracking direction and find the intersection points G1, G2 of the normal line and the blood vessel contour;, then tracking point P can be adjusted to meanwhile, change the tracking direction u to

The adjustment process is illustrated in Figure 4.

Figure 4

Centerline adjustment.

Bifurcation detection is another important process of the tracking algorithm. Ideally, we only need to distinguish two different vessel branches at the vessel bifurcation. However, in the actual tracking process, the accurate positions of vessel bifurcations are usually unknown. Thus, bifurcation detection is required at each point of the tracking process. We propose a robust bifurcation detection method. It includes two main steps: first, obtain one branch point (tracking point) P by the tracking method, and second, search in the fan ring area between angle (θ − Δθ′,θ + Δθ′) and radius (r1, r2) to find a ridge point that satisfies the following conditions: where P and u(P) represent the detected ridge point of the new branch (the branch point) and its direction, respectively; N(P) is the number of bifurcations around P; and τ is a threshold. The first three conditions mean that when u(P) significantly differs from u(P) and u(P) and the distance between P and P is large enough, the new branch has a large gap with the former branch. The last condition indicates that N(P) should be smaller than τ to avoid duplication with existing tracking. If all the conditions are satisfied, P is detected as a bifurcation point and we keep tracking the branch vessels. The schematic diagram of bifurcation detection is shown in Figure 5.

Figure 5

Vascular branch detection.

In addition, before tracking, the ridge image can be preprocessed to remove the scattered and distributed ridge points, it can also reduce the misjudgments of bifurcations. The specific steps areas follows: set a threshold τ, count the number of surrounding ridge points for each ridge point N(P), and then remove this ridge point if N(P) < τ and keep it otherwise.

5. Results and Analysis

In this section, we will conduct several experiments to justify the effectiveness of the proposed method. All the images are captured from the video data of coronary angiograms provided by Qilu Hospital (Qingdao). The experiments are implemented on an Intel Core i5-8300H and 8 GB of RAM processor using MATLAB software of version 2019b.

We carefully selected the parameters used. In multiscale enhancement, we set β, c, and σ to 0.5, 20, and [1 : 10]. In the proposed tracking method, we set the radius to 5 pixels and ∆θ to 45° in forward tracking. For bifurcation detection, it needs a larger area for searching; thus, we set the radius (r1, r2) to (7, 12) pixels and Δθ′ to 135° which can avoid backward tracking. Note that we set other thresholds I0, τ1, τ2, d, τ, τ, and τ to 10, 45°, 30°, 5, 4, 2, and 3.

Three images with the different vascular structures were selected for independent experiments, which are shown in Figure 6. We applied our methods for these images, and the results are shown in Figure 7.

Figure 6

Three selected original images.

Figure 7

Experimental results of three original images obtained by applying the proposed method. (a) Images preprocessed. (b) Vascular contour segmentation. (c) Improved adaptive tracking (red dots are bifurcation points, green dots are normal tracking points).

As can be seen from Figure 7, even though the vascular structures in the images are very different, the proposed method still has a nice experimental effect. From the images preprocessed (a), we can find that after the image preprocessing, the vascular structures were successfully highlighted and the background was restrained. The extraction of vessel contours (b) obtained vascular contours accurately and completely. The improved adaptive tracking method (c) is a core part of our work: compared with the original adaptive tracking method of [38], one of the main improvements in this approach is the bifurcation point detection part. We changed the originally fixed search radius to a proper search scope, which enhanced the capacity of the retrieval of bifurcation, and we used four conditions in Equation (17) to judge bifurcation point instead of only using the first condition, which greatly improved the detection accuracy and reduced the misjudgment. The results can be seen in Figure 8.

Figure 8

Comparison of the tracking effect between our proposed method and the method of [38]. (a) Results of [38]. (b) Results of the proposed method.

To achieve the completely automatic identification of vessels, we need to test the robustness of the proposed tracking method for randomly selected seed points from the detected ridge points. Taking Figure 6(c) as an example, three seed points were selected from different positions. The experimental results are shown in Figure 9. It can be seen that the results of the proposed method have strong robustness; that is, our method is generally applicable for automatically selected seed points. Meanwhile, our method is more accurate than the method of [38], which is clear in Figure 9 that the points of different types we identified are more approaching to the real vessel.

Figure 9

Experimental results of different seed points (blue). (a) The method of [38]. (b) The proposed method.

Even though our method has an improvement in accuracy and robustness compared to the former one, it still has some shortcomings. For example, the image preprocessing method is not effective enough for some images with complex vascular structures. Although the detection of bifurcation points has been improved compared with the method in [38], there are still a few misjudgments. This phenomenon can be seen in Figure 8(b). In the case of a more complex vascular structure, the tracking effect varies with the selection of seed points, and some vessel segments may be lost, as shown in Figure 7(c2).

6. Conclusion

In this paper, we designed a scheme of image preprocessing, used the C-V model, and proposed an improved adaptive tracking method, with which we can realize segmentation and automatic identification of vessels in coronary angiograms. Among these methods, the improved adaptive tracking method contains our major innovations that can enhance the capability of identifying vessels. Besides, we did many experiments to test our proposed method and the results turned out that our method is more robust and accurate than the former method.

Due to the complexity of coronary angiograms described above, traditional image processing methods are not effective enough. Hence, in the following work, we will continue to optimize the tracking algorithm and carry out image process research on deep learning to achieve a better effect.