Keywords

1 Introduction

Cardiovascular diseases are the leading cause of death in the developed countries essentially due to coronary atherosclerosis [1]. The medical imaging modality and the minor procedures most currently used to diagnosis the coronary diseases are the X-ray angiography. As angiographic images are subject to noise and radiological contrast agent that is widely heterogeneous, it is difficult to identify the arteries compared to the background. Segmentation is therefore necessary to extract the coronary arteries by eliminating artifacts contained in the background. Several methods were proposed for the coronary artery segmentation. In fact, Vessel segmentation algorithms are the fundamental component of automated radiological diagnostic systems such as diagnosis of the vessels (e.g. stenosis or malformations) and registration of patient images obtained at different times. Two categories of vessel segmentation algorithms are distinguished. The first one is skeleton based technique which aim first to extract median blood vessel and then connect these centerlines and estimate the vessel width to develop the vessel tree. The main idea of non skeleton methods is to directly extract blood vessels based mostly to the pixels intensity. In this context, several works were proposed to solve the vessel segmentation issue including basic ones such as thresholding and morphological operator. In [2], the authors propose a fuzzy clustering where each data point can belong to more than one cluster and clusters are determined via similarity measures, such as distance, connectivity, or intensity. [3, 4] use thresholding to extract vessels from angiographic images using the intensity information to assign pixels into categories. Mathematical morphology was proposed in [5], the authors provide systematic approach to analyze geometric Characteristic of artery by applying structuring elements :Dilation and erosion (SE) to images. In [6] the authors present a region growing procedure that segment images by incrementally selecting pixels into the artery region based on similarity and spatial proximity criteria. In [7], authors define a new segmentation method using an iteration region growing and the Sato vesselness function by merging both vesselness and direction information for vessel tracking. This method was extended in [8] where authors change sato filter with Frangi filter because of the sensitivity of the Sato filter which is higher than the Frangi filter [9].

More sophisticated techniques were also used to detect vessels such as deformable models which find the vessel contours using parametric curves that deform under the influence of internal and external forces [10].

In [14], the authors propose a method for segmenting coronary arteries based on graph cuts using multi-scale vesselness measure. In [11], the authors adapt the graph cut method for segmentation of thin blood vessels since the original method suffers from the shrinking bias drawback [12] which makes it not suitable for segmentation of thin objects like blood vessels.

In the latter work, the authors offer an energy functional GC which takes into account the probability of tabular structures, local connectivity of the vessel region and measuring edgeness.

In this paper we propose a new multi-scale Graph cut method for coronary artery segmentation in 2D X-ray angiograms by combining Vesselness measure to visualize the local vessel, geodesic paths to emphasize the local connectivity of another vessel region, a multi-scale edgeness map using the adaptive Canny detector and the directional information to guide the segmentation of blood vessels. The remainder of the paper is organized as follows: In Sect. 2, we remind the Graph cut method, Sect. 3 describes our proposed method. Experiments are presented and discussed in Sect. 4. And finally we conclude the paper by giving some details on future work.

2 Graph Cut Method

In this section, we overview graph cut (GC) method presented in [11]. The main idea of this latter paper is to use an energy formulation which takes into account the local vessel appearance using a vesselness measure, the local connectivity to other vessel regions using geodesic paths and a measure of edgeness based on a multi-scale adaptive Canny detector.

Image is seen as non-oriented graph where V is the set of nodes corresponding to the pixels and L is the set of arcs referring to the adjacency relations between pixels with two additional nodes the source s and the sink t. The graph arcs are composed of:

  • s-links: arcs connecting a node i to the source s;

  • t-links: arcs connecting a node i to the sink t;

  • n-links: arcs connecting neighbor pixels i and j.

A cut is a partition \(\left\{ S, T \right\} \) of nodes into two disjoint sets S and T where the source s belongs to S and the sink t belongs to T. The aim of the method is to find the optimal cut C with the lowest cost. The cost of the cut can be defined as:

$$\begin{aligned} C(S,T)= \sum _{i\in S,\ j\in T} C\left( i,j\right) . \end{aligned}$$
(1)

Where \( C\left( i,j\right) \) is the capacity of the arc connecting nodes i and j.

The t-links and s-links between a pixel i and the two nodes S and T terminals are weighted by a term named unary term U related to regional properties of the segmentation. This term can be considered as the probability that the pixel i belongs to the class \(\omega \) (i.e. ). Every pixel is connected with its N-neighborhood (N = 4,8) by a segment called n-link. The n-links are weighted by a regularization term designed to ensure the spatial coherence in a pixel neighborhood named boundary term \( B \left\{ i, j\right\} \).

Any \( B \left\{ i, j\right\} \ge 0 \) can be interpreted as the presence of a continuity between i and j. The total energy of a Cut C in the graph is defined by:

$$\begin{aligned} E(C)= \sum _{i\in V} U_i\left( w\right) + \sum _{i\in N ,j \in N} B \left\{ i, j\right\} \times \delta \left( i,j\right) . \end{aligned}$$
(2)

Where V is the set of vertices (nodes), \( \delta \left( i,j\right) =0 \) if \( \omega _i = \omega _j \;and \; 1,\) otherwise.

In [11] authors define the unary terms by combining vesselness map \(V \left( i \right) \) and the geodesic distance map Di, we note that the measure of vesselness is not totally effective in the detection of tabular structures since there are some vessel regions, which can have low vessel probabilities, so authors defeated this limitation by introducing a vesselness-geodesic measure VG which employ the geodesic paths among vessel seeds and compute geodesic distance among vessel seeds. In fact, the geodesic distance map includes the distance of each pixel to a set of centroid of a k-means clustering on the Cartesian coordinates of the vessels seeds. The map VG for pixel i is defined by the authors as:

$$\begin{aligned} VG_i=max \left( V \left( i \right) \ ,\frac{\left( \frac{1}{\left( D_i+ \mu \left( D \right) \right) }\right) }{\left( max \left( \frac{1}{\left( D+ \mu \left( D \right) \right) }\right) \right) } \right) \end{aligned}$$
(3)

where D correspond to the geodesic distance map and \(\mu \left( D \right) \) is its mean and \(V \left( i \right) \) is the measure of vesselness for the pixel i.

This measure \(V \left( i, \sigma \right) \) with i: indice, \( \sigma \): scale is calculated as follows:

$$\begin{aligned} V(i,\sigma ) = \left\{ \begin{array}{cc} 0 &{} if \; \lambda _2 > 0 \\ \exp \left( {\frac{ - R_\beta ^2 }{2 \times \beta ^2 }} \right) \times \left( 1 - exp\left( {\frac{ - \varDelta ^2 }{2 \times C^2 }}\right) \right) &{} otherwise \end{array} \right. \end{aligned}$$
(4)

where \(R_\beta =\frac{\lambda _1 }{\lambda _2 }\), \( \varDelta =\lambda _1 ^2 +\lambda _2 ^2\), \(\beta \) and c are thresholds which control the filter’s sensitivity to \( R_\beta \) and \( \varDelta \) respectively, \( \lambda _1(i,\sigma ) \; and \; \lambda _2(i,\sigma )\) where \( |\lambda _1(i,\sigma )|<= |\lambda _2(i,\sigma )| \) are the eigenvalues of the Hessian matrix of image I computed at scale \(\sigma \) and location i.

\( R_\beta \) is the blob-like structure measure, this measure is used to distinguish between structures which form lines and flat structures.

\(\varDelta \) is the Frobenius norm of the Hessian matrix, this value is small in the background with individually small eigenvalues by lack of contrast.

This Eq. (4) can be evaluated at different scales Q, \(\sigma = \left\{ \sigma ^{1} ,.., \sigma ^{Q} \right\} \), and the V for each pixel is given by \(V (i) = max_{\sigma \in \left\{ \sigma ^{1} ,.., \sigma ^{Q} \right\} }V(i,\sigma ) \).

The unary terms for each pixel i can therefore be written as:

(5)
(6)

Where is the probability of a pixel belonging to a vessel region, is the opposite probability.

In order to initialize the boundary term, authors use the Canny edge detector algorithm on the image at different threshold levels.

3 The Proposed Method

In this work, we propose a GC new energy functional adapted to the segmentation of vessel problem taking into consideration direction information given by the first Hessian eigenvector in order to extract thin vessels in addition to the vesselness measure, the geodesic path and the multi-scale edgeness measure. The direction and the normal direction of the potential linear structure describe the local shape of the intensity and can be determined respectively by computing the eigenvectors \(e_1\) and \(e_2\) of the hessian matrix. The direction information \(D_{\sigma }\) of a pixel i is given by the following expression:

$$\begin{aligned} D_{\sigma }(i) = \left\{ \begin{array}{l} 0 \;if \;\lambda _2 < 0 \\ e_1 \; otherwise \cdot \end{array} \right. \end{aligned}$$
(7)

Where \(\lambda _2 \) is the eigen value of the Hessian matrix of image I computed at scale \( \sigma \) and location i.

The unary term is computed identically to the described graph cut method by combining the vesselness map and the geodesic map. For boundary terms, the key idea is to compute them in such a way they guarantee the continuity of the segmentation along the artery direction.

The initial boundary potential over the multiscale edgeness map is computed as \( B\left\{ i,j\right\} =J_i^*\).

\(J_i^*\) is computed as follows:

$$\begin{aligned} J_i^*=min_f \frac{1}{n } \sum _{e=1}^{n} J_{i,\theta _e,\alpha ^f } \end{aligned}$$
(8)

where \(\theta _e\) is the threshold, \(\alpha ^f\) is a scale and \( J_{i,\theta _e,\alpha ^f }\) is the binary edge map for pixel i.

Scales should be adjusted between the approximate width of the smallest and largest vessel to be detected. For different scales \(\sigma \), the direction \((D_{\sigma })\) is computed and the measure \(B\left\{ i,j\right\} \) is increased only if the pixels i and j have close directions even if their the edge probability values are lower than a given threshold.

The final boundary potential is computed as:

If \( B\left\{ i,j\right\} >\eta \left( 1- \varPhi _\sigma \left( i,j\right) \right) \) then \( B_{current} \left\{ i,j\right\} =B_{ancient} \left\{ i,j \right\} + J_i^*\).

Where \(\eta \) is a parameter is a basic threshold of the edgenesse value and \(\varPhi _\sigma \left( i,j\right) \) is a correlation index between vessel orientations at pixels i and j defined as:

$$\begin{aligned} \varPhi _\sigma \left( i,j\right) )= \frac{D_\sigma \left( i \right) \times D_\sigma \left( j \right) }{|| D_\sigma \left( i \right) || \times ||D_\sigma \left( j \right) ||} \end{aligned}$$
(9)

\(\varPhi _\sigma \left( i,j\right) \) equal to 1 if the two local direction vectors are parallel, 0 if they are orthogonal, and a value in ]0..1[ otherwise.

Finally, min cut algorithm [13] is applied to find the segmentation with minimum energy and we keep the largest connected component in the final segmentation.

4 Experimental Result

We conduct experiments on the three datasets DS1, DS2 and DS3 defined in [8] with in total 91 images from 60 patients. Each dataset offers a challenging difficulty to evaluate the proposed algorithm. Experiments on Ds1 aim to evaluate the general performance of the algorithm while Ds2 contains simple vessel illumination cases and stenosis degrees. The images in DS3 present more complicated cases such as Severe stenosis, stent(s) and low signal to noise ratio. Expert hand segmentation is provided for all images of the datasets. In all experiments, \( \sigma \in \left\{ width \;of\; smallest\; vessel ,...,width\; of\; largest\; vessel \right\} \), in our dataset, the approximate width of vessels vary from less than 1 pixel to more than 10 pixels, \( \eta =60 \% \). In Figs. 1, 2 and 3, we present original sample images respectively from DS1, DS2 and DS3 and their corresponding segmentation results with the AQCA method and the proposed method.

Fig. 1.
figure 1

Segmentation results: Top to bottom: real angiograms images from DS1; ground truth; AQCA and our method

Full size image
Fig. 2.
figure 2

Segmentation results: Top to bottom: real angiograms images from DS2; ground truth; AQCA and our method

Full size image
Fig. 3.
figure 3

Segmentation results: Top to bottom: real angiograms images from DS3; ground truth; AQCA and our method

Full size image
Fig. 4.
figure 4

Dice similarity measure for: (a): DS1, (b): DS2 and (c): DS3

Full size image

We can visibly notice that the proposed method can detect more vessels and is much closer to the ground truth than the AQCA method for the three datasets. In fact, the results given by the AQCA method describe the main vessels in the image and ignore others relevant ones, while our method successfully detect most of them. In order to quantitatively assess the quality of the proposed algorithm, we compute the false positive rate RFP, the sensibility and Dice similarity coefficient. DSC indicates the spatial cover between the segmented vessel and the ground truth and defined as:

$$\begin{aligned} DSC=\left( 2\times |GT \cap Seg|\right) /(|GT |+|Seg|) \end{aligned}$$
(10)

Where |GT| and |Seg| are the number of pixels classified as vessels respectively in the ground truth image and segmented image. \(|GT \cap Seg|\) gives the number of pixels in the intersected region between these two images.

The sensitivity is the percentage of pixels belonging to vessels which correctly segmented and defined as:

$$\begin{aligned} Sensibility=\left( TP\right) /(TP + FN) \end{aligned}$$
(11)

Where TP is the number of true positive cases and FN is the number of false negative cases.

The false positive rate RFP is the percentage of pixels belonging to vessels which incorrectly segmented as belonging to the background and defined as:

$$\begin{aligned} RFP=\left( FP \right) /(FP + TN) \end{aligned}$$
(12)

where FP is number of false positives and TN is number of true negatives.

The obtained results are presented as boxplots that indicate the first, second and third quartile with whiskers from minimum to maximum. In Fig. 1, Boxplots of DSC measures are drawn for both AQCA and the proposed method and for each of the three datasets.

Figure 4 shows that the DSC measures given by the proposed method exceed 0.7 for Ds1 and Ds2 which is considerate as good index. The measures are less satisfying in the case of Ds3 which could be explained by the quality of the images with a very non homogeneous background. The results of the proposed method outperform the ones given by the AQCA method. This difference is more remarkable for the two first datasets but less significant for the third one (Fig. 5).

Fig. 5.
figure 5

Sensitivity measure for: (a): DS1, (b): DS2 and (c): DS3

Full size image

We can clearly notice that the proposed method is more accurate than the AQCA method in fact for the three datasets, the sensitivity values of our method exceed those of the AQCA method.

Fig. 6.
figure 6

RFP measure for: (a): DS1, (b): DS2 and (c): DS3

Full size image

From Fig. 6, we can see that the RFP values of our algorithm are higher than the RFP values of AQCA due to the appearance of false positives frequently on thin vessels.

5 Conclusion

In this paper, we have presented a new coronary artery segmentation method based on multiscale analysis and the graph cut algorithm by merging the multiscale edgeness with Hessian geometrical features. The main idea is to introduce the multiscale boundary term that ensure the continuity of the segmentation along the artery direction thanks to the direction information and detect the vessel boundaries using the multi-scale edgeness measure. A comparison was made with the AQCA method and showed that the results given by the proposed method are very encouraging allowing to detect more vessels than the AQCA method.

Despite the robustness of the approach, the results obtained by our algorithm show the appearance of false positives frequently on thin vessels. In order to overcome this limitation, we will try, as a future work, to locally adapt the scale to avoid the overemphasis of these vessels by processing them at a scale less or equal than their sizes.