Asma Kerkeni
email: asma.kerkeni@gmail.com

Asma Benabdallah

Antoine Manzanera

Mohamed Hedi

Mohamed Hedi Bedoui

A Coronary Artery Segmentation Method Based On MultiScale Analysis And Region Growing

Keywords: Coronary Artery, Hessian, Segmentation, Region Growing, Multiscale, Vesselness, Direction

A coronary artery segmentation method based on multiscale analysis and region growing

Introduction

X-ray Coronary Angiography (XCA) is the gold standard for the assessment of clinically significant coronary artery diseases (CAD) [START_REF] Li | A robust coronary artery identification and centerline extraction method in angiographies[END_REF]. The angiograms obtained by the XCA enable to reveal the initial CAD symptoms by the morphological features of the coronary arteries such as diameter, length, branching angle, and tortuosity. However, complex vessel structure, image noise, poor contrast and non-uniform illumination make vessel tracking a tedious task.

Accurate coronary vessels detection is a fundamental step in various medical imaging applications such as stenosis detection [START_REF] Sato | A viewpoint determination system for stenosis diagnosis and quantification in coronary angiographic image acquisition[END_REF], 3D reconstruction [START_REF] Liao | 3-d reconstruction of the coronary artery tree from multiple views of a rotational x-ray angiography[END_REF] and cardiac dynamics assessing [START_REF] Zheng | Assessing cardiac dynamics based on x-ray coronary angiograms[END_REF]. Vessel detection is generally related to two important tasks which are vascular features enhancement and blood vessel segmentation.

The enhancement step aims to improve the vessels delineation while reducing background artefacts. So far, a variety of vessel enhancement methods have been proposed [START_REF] Koller | Multiscale detection of curvilinear structures in 2D and 3D image data[END_REF][START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF][START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF][START_REF] Lorenz | Multi-scale line segmentation with automatic estimation of width, contrast and tangential direction in 2d and 3d medical images[END_REF][START_REF] Krissian | Model-based detection of tubular structures in 3d images[END_REF][START_REF] Erdt | Automatic hepatic vessel segmentation using graphics hardware[END_REF][START_REF] Kroon | Coherence filtering to enhance the mandibular canal in cone-beam ct data[END_REF][START_REF] Xiao | Multiscale bi-gaussian filter for adjacent curvilinear structures detection with application to vasculature images[END_REF][START_REF] Dufour | Filtering and segmentation of 3d angiographic data: Advances based on mathematical morphology[END_REF]. Truc et al. [START_REF] Truc | Vessel enhancement filter using directional filter bank[END_REF] distinguish three method classes : linear filtering [START_REF] Zana | Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation, Image Processing[END_REF], non-linear anisotropic filtering [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF][START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF] and Hessian-based multiscale filtering [START_REF] Koller | Multiscale detection of curvilinear structures in 2D and 3D image data[END_REF][START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF][START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF][START_REF] Lorenz | Multi-scale line segmentation with automatic estimation of width, contrast and tangential direction in 2d and 3d medical images[END_REF][START_REF] Krissian | Model-based detection of tubular structures in 3d images[END_REF][START_REF] Erdt | Automatic hepatic vessel segmentation using graphics hardware[END_REF].

Linear enhancement methods generally use Gaussian kernels or Gabor filters in order to denoise images. They are inappropriate to complex vascular structures since they blur vessel borders as well as thin vessels.

Unlike the linear smoothing filters, non-linear anisotropic diffusion filtering adjusts the filter for local variations by acting mainly along the preferred structure direction [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF]. Hence, important features are better preserved during the smoothing process. Those filters are widely applied in many image processing tasks such as fingerprint image filtering [START_REF] Greenberg | Fingerprint image enhancement using filtering techniques[END_REF], Optical Coherence Tomography image denoising [START_REF] Cabrera Fernández | Automated detection of retinal layer structures on optical coherence tomography images[END_REF] and cell membrane enhancement [START_REF] Pop | Image filtering using anisotropic structure tensor for cell membrane enhancement in 3d microscopy[END_REF]. The main drawback of diffusion based methods is that they usually act at a fixed scale and are therefore unable to detect vessels within a wide range of size.

Most of the research on vessel enhancement falls into falls into Hessian-based methods. These methods utilize the second order derivatives matrix (Hessian matrix ) of the image intensity in order to detect tubular-like structures.

For that, the Hessian matrix at each pixel is computed by convolving the initial image with second-order Gaussian derivatives. Then, the Hessian eigenvalues and/or eigenvectors are often used to propose a vesselness measure which can be interpreted as the probability of belonging to a blood vessel of a pixel. The vesselness generates maximum response at the scale which matches the diameter of the vessel to be detected. The main advantage of methods in this category is that can perform in a multiscale fashion, thus detecting different sized objects.

In [START_REF] Kerkeni | Coronary artery multiscale enhancement methods: A comparative study[END_REF], the authors compared four different Hessian-based multiscale filters, namely those proposed by Koller [START_REF] Koller | Multiscale detection of curvilinear structures in 2D and 3D image data[END_REF],

Sato [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF], Frangi [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] and the Vessel Enhancing Diffusion (VED) filter. VED is a diffusion filter built upon the Frangi filter. It can be seen as a generalization of anisotropic diffusion to iteratively smooth the image, while preserving the vessel structure [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF]. This filter was initially proposed by [START_REF] Caero | Vesselness enhancement diffusion[END_REF] and extended by Manniesing et al. [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF] who added a smoothness constraint and an enhanced scheme. Results of this study showed that the two latter outperform the others with a far better background suppression performance and vascular structure enhancement.

The second important step in vascular analysis is vessel segmentation. Existing segmentation methods can be divided into two general categories [START_REF] Suri | A review on mr vascular image processing: skeleton versus nonskeleton approaches: part ii[END_REF]: skeleton based, and non-skeleton based. The principle of skeleton-based techniques is the segmentation of vessels by first detecting the centerlines, then estimating the vessel width. A typical problem with this class is that it can fail to approximate the vessel segments which present stenosis and bifurcations.

Non-skeleton-based segmentation techniques, on the other hand, are those that extract the vessels directly. Methods in this category vary from thresholding [START_REF] Condurache | Vessel segmentation in angiograms using hysteresis thresholding[END_REF][START_REF] Blondel | Reconstruction of coronary arteries from a single rotational x-ray projection sequence[END_REF][START_REF] Andriotis | A new method of three-dimensional coronary artery reconstruction from x-ray angiography: Validation against a virtual phantom and multislice computed tomography[END_REF][START_REF] Liao | 3-d reconstruction of the coronary artery tree from multiple views of a rotational x-ray angiography[END_REF], region growing [START_REF] O'brien | Automated segmentation of coronary vessels in angiographic image sequences utilizing temporal, spatial, and structural constraints[END_REF][START_REF] Eiho | Branch-based region growing method for blood vessel segmentation[END_REF][START_REF] Sekiguchi | Blood vessel segmentation for head mra using branch-based region growing[END_REF][START_REF] Martinez-Perez | Segmentation of blood vessels from red-free and fluorescein retinal images[END_REF][START_REF] Li | A novel method of vessel segmentation for x-ray coronary angiography images[END_REF][START_REF] Orkisz | Segmentation of the pulmonary vascular trees in 3d ct images using variational region-growing[END_REF][START_REF] Zhao | Retinal vessels segmentation based on level set and region growing[END_REF], fuzzy clustering [START_REF] Shoujun | Automatic segmentation of coronary angiograms based on fuzzy inferring and probabilistic tracking[END_REF], mathematical morphology [START_REF] Ko | A fully automated identification of coronary borders from the tree structure of coronary angiograms[END_REF][START_REF] Sun | Extraction of vascular tree on angiogram with fuzzy morphological method[END_REF][START_REF] Dufour | Filtering and segmentation of 3d angiographic data: Advances based on mathematical morphology[END_REF], deformable models such as active contour and level set [START_REF] Brieva | A level set method for vessel segmentation in coronary angiography[END_REF][START_REF] Sum | Vessel extraction under non-uniform illumination: a level set approach[END_REF][START_REF] Tagizaheh | Segmentation of coronary vessels by combining the detection of centerlines and active contour model[END_REF][START_REF] Chen | A semi-automatic method for vascular image segmentation[END_REF] to graph cuts [START_REF] Schaap | Coronary lumen segmentation using graph cuts and robust kernel regression[END_REF][START_REF] Esneault | Liver vessels segmentation using a hybrid geometrical moments/graph cuts method[END_REF][START_REF] Hernandez-Vela | Accurate coronary centerline extraction, caliber estimation, and catheter detection in angiographies[END_REF]. A more general review on vessel segmentation can be found in [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF] and [START_REF] Lesage | A review of 3d vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF]. Hereafter, we especially focus on region-based vessel segmentation which is used in our work. Region-growing has been widely used for image segmentation [START_REF] Zucker | Region growing: Childhood and adolescence[END_REF], and in particular medical image applications such as vessel extraction. The rationale for this popularity is that these methods are based on a connectivity assumption, which is naturally suited to the case of the vascular trees [START_REF] Yim | Gray-scale skeletonization of small vessels in magnetic resonance angiography[END_REF][START_REF] Hu | Segmenting linear parts using layered region growing[END_REF]. Region-growing incrementally segment an object by recruiting neighboring pixel starting from seed points or regions located inside a vessel based on some inclusion rule. Defining a robust region growing rule for X-ray vessel segmentation may encounter several difficulties such as image artefacts, lesions, noise and very low contrast between the vessels and the background particularly in thin vessels location. Therefore, classical region growing methods based on grey level values and/or spatial proximity are often sensitive to noise, and inhomogeneous contrast, which often leads to false negative (holes), false positive (leakage) and unwanted stop of region-growing [START_REF] Lesage | A review of 3d vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF].

In coronary artery angiograms, several factors such as image artefacts, stenosis and noise may introduce discontinuities, hence region growing may result in holes and over-segmentation. To avoid such problems, the growing process was performed in many frames of the same sequence in [START_REF] O'brien | Automated segmentation of coronary vessels in angiographic image sequences utilizing temporal, spatial, and structural constraints[END_REF]. Although discontinuities may be avoided by temporal tracking, other difficulties such as centerline matching and user interaction in each step are necessary in this approach.

Instead of using the intensity feature in the growing process, authors in [START_REF] Lara | A semi-automatic method for segmentation of the coronary artery tree from angiography, Computer Graphics and Image Processing[END_REF] incorporate the Frangi vesselness filter [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] in order to introduce more seed points when the growing process stops.

In [START_REF] Li | A novel method of vessel segmentation for x-ray coronary angiography images[END_REF], the authors proposed a novel vesselness function and performed vessels segmentation in two steps. First, large vessels are extracted from the maximum vesselness response by region growing. However, this final vesselness response given by the maximum response among the scale space may present low values for thin vessels, junctions and stenosis location. Thus discontinuities may occur if the growing process is applied directly on this final vesselness image. That's why a detail repairing process is launched in order to extract thin vessels using a direction information given by the first directional derivatives i.e. image intensity gradient.

Combining the vesselness function with a directional information for vessel tracking is certainly a good idea but the image gradient isn't, in our opinion, the appropriate direction information descriptor. A more precise direction descriptor would be given by the Hessian matrix eigenvectors. This idea was investigated in [START_REF] Nimura | Blood vessel segmentation using line-direction vector based on hessian analysis[END_REF] where the authors combined the Sato vesselness function [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF] with the line direction information given by the first Hessian eigenvector to guide the segmentation of 3D blood vessels. As discussed in [START_REF] Kerkeni | Coronary artery multiscale enhancement methods: A comparative study[END_REF], the Sato filter is more sensitive to noise than the Frangi filter when applied to 2D X-ray angiograms. Fig. 1 shows that background noise is better removed by the Frangi filter than the Sato one.

In this paper, a multiscale vessel segmentation method for X-ray coronary angiography images is presented. This method takes advantages from multiscale Hessian analysis strength and aims to overcome the two major problems encountered in region growing techniques which are the difficulty of detection of poorly enhanced vessel segments and the occurrence of false positives [START_REF] Nimura | Blood vessel segmentation using line-direction vector based on hessian analysis[END_REF]. The major contributions of this paper is threefold: (1)From methodological point of view, we combine both Hessain geometrical features which are eigenvalues and eigenvectors to define a robust new region growing criterion tailored to the coronary artery segmentation problem. This criterion is integrated into a mulitscale region growing algorithm ensuring the detection of different sized vessels. [START_REF] Sato | A viewpoint determination system for stenosis diagnosis and quantification in coronary angiographic image acquisition[END_REF] We propose an evaluation database containing two datasets with different challenging degrees, which allow the quantitative evaluation of the method in terms of region overlap, sensitivity and precision. To the best of our knowledge, there is no standard database available for this task. The only database available in the literature was proposed by [START_REF] Hernandez-Vela | Accurate coronary centerline extraction, caliber estimation, and catheter detection in angiographies[END_REF] and was designed to evaluate centerline extraction methods rather than segmentation evaluation. (3) We evaluate our method from two aspects which are the enhancement filter and the segmentation methods. The remainder of this paper is organised as follows. In section 2, background concepts related to our approach are reviewed. Section 3 illustrates the proposed method in detail. Experiments are described in section 4. Some illustrative results are given in section5. Section 6 presents discussion and concluding remarks.

Background on Multiscale Analysis

The idea of multiscale image analysis is to add a new dimension to the analysis which is the image scale. The image is transformed into a set of blurred images, each representing the original image, but at a different scale [START_REF] Lindeberg | Principles for automatic scale selection[END_REF].

These blurred images are obtained by convolving the initial image I 0 (p) = I 0 (x, y) with a Gaussian kernel to represent the information at a certain scale.

I σ (p) = I σ (x, y) = I 0 (x, y) ⊗ G σ (x, y) (1)
where I σ is an image of the scale space, p = (x, y) is a pixel location, ⊗ represents the convolution operation and G σ (x, y) is the 2D Gaussian kernel with standard deviation σ defined as:

G σ (p) = G σ (x, y) = 1 √ 2πσ 2 exp - (x 2 + y 2 ) 2σ 2 (2)
where σ ∈ Σ = {σ min , .., σ max }, σ min and σ max are set according to the approximate width of the smallest and largest vessel to be detected [START_REF] Lindeberg | Principles for automatic scale selection[END_REF].

In the scale space framework, differentiation is defined as a convolution with derivatives of Gaussian:

( ∂ n1+n2 I ∂ x n1 ∂ y n2 ) σ = I ⊗ ∂ n1+n2 G σ ∂ x n1 ∂ y n2 (3)
Thus, the Hessian matrix which is the second order derivative of an intensity image in the scale space can be obtained at each point by Eq.4:

H σ (p) =   ∂ 2 I σ (x,y) ∂ x 2 ∂ 2 I σ (x,y) ∂ y∂ x ∂ 2 I σ (x,y) ∂ x∂ y ∂ 2 I σ (x,y) ∂ y 2   . (4)
The Hessian matrix is symmetric, therefore has two real eigenvalues λ 1 and λ 2 and two associated eigenvectors e 1 and e 2 . The signs as well as the relative and absolute magnitudes of eigenvalues characterize the local shape of the intensity in the image. If we assume that we seek to characterize black vessels on a white background, which is the case of the angiography images, a pixel belonging to a vessel region will be given by small (ideally 0) λ 1 and a high positive value of λ 2 (with λ 1 < λ 2 ).

Frangi Vesselness Filter

Frangi et al. [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] note that the Hessian eigenvalues ordered by magnitude (i.e |λ 1 | < |λ 2 |) are geometrically interpreted as principal vessel curvatures and provide a vesselness measure on the basis of this observation. They introduce two measures to describe structures in images:

R b = λ 1 λ 2
is the blob-like structure measure and S = H F = λ 2 1 + λ 2 2 , the Frobenius norm of the Hessian matrix, is the second-order structureness measure. These measures were combined in a vesselness function as given by Eq.5.

V σ (p) =      0 if λ 2 < 0 exp(- R 2 b 2β 2 )(1 -exp(-S 2 2c 2 )) otherwise.
(

) 5
Where parameters β and c are thresholds which control the filter's sensitivity to R b and S respectively. The two vesselness terms exp(-

R 2 b 2β 2 ) and (1 -exp(-S 2 2c 2
)) are combined using their product to ensure that the response of the filter is maximal only if both criteria are fulfilled.

According to scale space theory, V σ (p) will be the maximum only when the width of the vessel in pixel p matches a suitable scale factor σ . Thus, In the multi-scale vessel enhancement algorithm, V σ is computed for different scales then the largest one is taken as the final filter output.

V max (p) = max σ ∈Σ σ γ * V σ (p) (6)
The normalization factor σ γ [START_REF] Lindeberg | Principles for automatic scale selection[END_REF] is introduced to normalize the intensity values across the scale space. The parameter γ can be adjusted to show preference for particular scales. It is equal to one if no scale is preferred. Taking γ < 1 favours small scales while γ > 1 favours large ones.

Vessel Direction Information

The eigenvectors e 1 and e 2 of the Hessian matrix of Eq.4 characterize the local shape of the intensity. e 1 gives the direction of lower second derivative, i.e. the direction of a potential linear structure, while e 2 gives the direction of the stronger second derivative, i.e. the normal direction to a potential linear structure. Thus the proposed direction information D σ can be simply given by the following expression:

D σ (p) =      0 if λ 2 < 0 e1 otherwise. ( 7
)
Fig. 2 shows a motivating example of incorporating the direction information. The green circle in Fig. 2(a) shows a stenosis in a coronary artery. When the maximum Frangi vesselness is computed, the filter response significantly decreases in this location. However the first eigenvector ensure the continuity of concerned pixels and their neighbours (see Fig. 2(d)). One can note that neighbouring pixels in this location exhibit close direction.

Multiscale Segmentation Method

Leveraging from the above-mentioned multiscale geometrical information, we design a robust multiscale region growing process able to extract both wide and thin vessels while avoiding undesirable structure overemphasis and early segmentation stop.

A general overview of the proposed method is depicted in Fig. 3. First, seed points are selected by seeking for points with the highest Frangi vesselness response (section 3.1). Subsequently, starting from the highest scale σ max , the vesselness (V σ ) and direction (D σ ) features are computed for each scale. Then a region growing process based upon a specific criterion combining those two features is performed. This criterion will be described in section 3.2. If all the scales are explored i.e σ min is reached, the segmented image is obtained else the segmentation process restarts with the next scale σσ step . The details of the process evolution over the different scales are described in section 3.3.

Initial seed points selection

Candidate seed points are automatically generated by the Frangi vesselness map analysis. Pixels exhibiting the maximum of vesselness are considered as a good starting point. Since the vesselness is always maximum along the vessels' centerlines and close to zero outside tube-like regions, the initial seed points are generally located in the center of vascular regions. To avoid points in the peripheral positions which may have high response, the search is done inside a sub-area excluding the image borders.

The red points in Fig. 4(a), 4(b) and 4(c) show some automatically generated seeds on three different images.

These points are located in the center of wide vessels. One can see that some regions especially thin vessels with very low illumination may be unreachable. Thus, if the user is interested in thin vessel extraction, he may add some more additional seeds. The added seeds are depicted in green in Fig. 4. Sometimes, the area thought to be had to reach is segmented even without added markers.

Region Growing Condition

The idea behind our approach is to incorporate the vessel direction information with the vesselness response in order to improve the vessel segmentation result. This information have to be included in a way that favours neighbours exhibiting directions close to the direction of a seed point even if their vesselness responses are lower than a given threshold. This will help to ensure the continuity of the segmentation along the artery direction, and reduce the number of broken branches as well. Based on this idea, we propose a region growing algorithm based on the following convergence condition: for all pixel p, if there exists q ∈ N 4 (p) such that B(q) = 1 and Here B is the binary image corresponding to the segmentation, N 4 (p) is the 4-neighbourhood of pixel p, and Ω σ is a correlation index between vessel orientations at pixels p and q, defined as follows:

V σ (p) > η (1 -Ω σ (p, q)) , then B(p) = 1. ( 8
) (a) (b) (c)
Ω σ (p, q) = D σ (p).D σ (q) D σ (p) D σ (q) (9)
Ω σ (p, q) becomes one when the two local direction vectors are parallel and becomes zero if they are orthogonal. As a result, the growing process may aggregate pixels with close directions even if their vesselness values are low. The η parameter is a basic threshold of the vesselness value needed to growing the vessel region. If we assume that the vesselness is rescaled in the interval [0,1], η will be between 0 and 1. The smaller value of the parameter η can result in the overgrowth of the vessel region, while the bigger one will limit the extension of the seed point set.

MultiScale Region Growing Algorithm

A pseudo-code for the MSRG algorithm is given in Algorithm 1. The initial seeds (binary image S) are used to initialise an iterative region growing segmentation process, whose output is the binary image B. Then, at each scale starting from the largest scale σ max , the vesselness (V σ ) and direction (D σ ) features are computed from the original image I.

The root of the coronary tree is found in the coarsest scale as a candidate with the highest Frangi vesselness.

Most of the automatically generated seed points are located in large vessels, thus the coarse-to-fine scheme makes the growing process more fluid and faster. This scheme is also less sensitive to noise often encountered in thin vessels.

At each scale, the region growing process is implemented using a processing queue Q, which is initialised with the current segmentation image B. Then, the first pixel of the queue (q ← Q.first) is picked, its 4-neighbours meeting the growing condition of Eq.8 are added to the queue and to the segmented image, then q is pulled out of the queue, the new first pixel of the queue is picked, and so on until the queue is empty. If all the scales have been explored, i.e σ min is reached, then the segmented process is done else it restarts with the next scale σσ step .

Algorithm 1 MULTISCALEREGIONGROWING (I, S)

Require: S = 0

1: B ← S 2: σ ← σ max 3: while σ ≥ σ min do 4: Q ← {q; B(q) = 1} 5: while Q = / 0 do 6: q ← Q.first 7:
for all p ∈ N 4 (q) do 8:

if

B(p) = 0 and V σ (p) > η (1 -Ω σ (p, q)) then 9: B(p) ← 1 10: Q ← Q ∪ {p} 11: end if 12:
end for 13:

Q ← Q \ {q} 14:
end while

Materials

Real clinical X-ray angiograms were acquired during routine cardiac catheter examinations using a single plane Siemens Artis Zee angiography system available at the Cardiology Department of The University Hospital Fattouma Bourguiba, Monastir, Tunisia.

For each intervention, Dicom sequences of different projections were obtained by varying the left or the right anterior oblique (LAO or RAO) angle and the caudal or cranial (CAUD or CRAN) angle in order to visualize both left and right coronary arteries. The acquired series were then extracted and interesting angiograms from both healthy and diseased patients with different stenosis severity degrees.

A total of 50 images obtained from 25 patients were selected by the clinicians. Each image has a size of 512*512 pixels and 256 gray levels with a resolution of 0.3x0.3 mm. This database was equally split into 2 datasets. The first one (DS1) contains relatively easy cases in terms of vessel illumination and stenosis degrees, whereas the second dataset (DS2) presents difficult situations that may occur in X-ray coronary angiograms, such as low signal to noise ratio, severe stenosis and the overlapping of many vessels.

In order to evaluate the different segmentation methods on these images, manually gold standard were outlined by two trained operators and verified by clinicians resulting in two sets 1 and 2. The manually segmented images in set 1 and set 2 noted GT1 and GT2 respectively are used as ground truth. The disagreement between the two observers occurs especially on vessels borders and endpoints.

Performance measures

Three evaluation metrics, namely Dice coefficient [START_REF] Dice | Measures of the amount of ecologic association between species[END_REF], sensitivity and precision were used. The DICE similarity index is a measure of region overlap defined by:

Dice i = 2 * |GT i ∩ Seg| |GT i | + |Seg| for i ∈ {1, 2} (10)
where GT i is the ground truth region depicted by the observer i and Seg is the automatically segmented region. Dice is always between 0 and 1, with higher values indicating better match between automatic and manual segmentations.

Precision is the rate of true-positives to the total number of positives. It is expressed by the following formula:

Precision i = T P T P + FP for i ∈ {1, 2} (11)
TP is the number of true positive i.e. correctly classified vessel whereas FP is the number of false positive i.e. ground truth non-vessel pixels overlapping detected pixels.

Sensitivity measures the proportion of ground truth positives which are correctly identified as such. It is given by:

Sensitivity i = T P T P + FN for i ∈ {1, 2} (12)
where FN is the number of false negatives i.e. incorrectly rejected pixels.

Parameter setting

All the results presented in this work use β = 0.5 and half the value of the maximum Hessian norm for c (c = 1 2 max H(x) F ) as proposed in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF].

Scales should be set according to the approximate width of the smallest and largest vessel to be detected. Scales should better match the radii of coronary arteries. In our database vessel radii vary from less than 1 pixel to more than 10 pixels wide in some images. The performance of our method, as any other multiscale method, highly relies on how accurate the scales estimation is, especially on the σ max parameter. In fact, a lower σ max will suffer from noise disturbance and bar-like profiles. In contrast, an over-large σ max might blur the adjacent structures and over-estimate the vessel caliber. In our work, for each image, the largest vessel diameter d max was calculated by averaging the maxima of the distance transform computed on the two ground truth vessel trees and σ max was set to 0.5 * d max . σ min was set to 1 since vessels with smaller radii are not of interest. Four logarithmically spaced scales between σ min and σ max were then generated for each image.

The η parameter is calculated as a fixed quantile of the accumulated histogram of the vesselness image V σ . We experimentally optimize this parameter and find that 95 percent is an acceptable threshold for our dataset.

Methods

The evaluation of our method robustness and efficiency is done in two stages. First, we seek to see the effect of the enhancement filter chosen in this work. Then, we compare the proposed segmentation with other state of the art methods.

Enhancement filters

In order to justify the choice of the Frangi filter for vessel detection within our approach, we build our segmentation algorithm upon different enhancement filters which are Sato filter, VED filter and Frangi filter. which incorporate an automatic point seeding method and a multiscale edgeness term to detect vessel borders. We used the implementation available in the author website 3 and the recommended parameters.

Results

We divide the results in enhancement (section 5.1) and segmentation (section 5.2) evaluation.

Enhancement filter

Dice, sensitivity and precision of the three MSRG versions are presented in Fig. 7 for the first dataset and Fig. 8 for the second dataset. The first row compares the second observer and the three methods to the first observer and the second row presents an evaluation against the second observer. Regarding the two datasets, the Sato filter has the lowest performance in terms of Dice similarity, sensibility and precision. This confirms the sensibility of this filter to noise. The Frangi and VED filters have very similar performance measures with a slight superiority for the Frangi filter. This is due to the impact of the blurring effect of VED filter especially on small vessels. One can see in the red rectangles in Fig. 6(c) and Fig. 6(d) that some discontinuities occur in thin vessels. Thus, distal parts of those vessels are missed in the segmentation step.

Segmentation methods

Fig. 9 shows dice, sensitivity and precision for the four different methods (HT, FC, RG, AQCA and MSRG) as well as each observer against the other. It is worth to note that the plots in the first and second row have roughly similar appearance. A small inter-observer variability occurs in the specificity plot and shows that the first observer detects more vessels pixels than the second one. Regarding the dice metric, the FC method has the lowest dice, HT achieves better than FC but lower than the three other methods. RG and AQCA have slightly equal dice values and the MSRG achieves the best overlap values. The method HT has a higher sensitivity than FC, RG and AQCA but it has the disadvantages to produce many false positives as confirmed by the low precision values. Our method which can be seen as an enhanced region growing method achieves better sensitivity while keeping the same precision level as the region growing. It is also worth to note that the AQCA has the highest precision but it remains lower than the inter-observer variability.

Regarding the second dataset (Fig. 10), the methods keeps the same order in term of segmentation overlap (dice). The inter-observer variability has increased in term of specificity and precision. The precision of the HT, FC and RG methods have a large standard deviation which means that they may introduce a lot of false positives. Our method still achieve the best sensitivity and an acceptable rate of precision (80% ) and AQCA has the best precision. The outliers (red marks) seen in the different plots are related to bad seed or threshold automatic initialization.

Moreover, Fig. 11 shows some resulting outputs of the five different segmentation methods. From top to down, we find three different real angiograms; ground truth; hysteresis thresholding; fuzzy clustering; region growing; AQCA and our method outputs. The first and second images belong to the first dataset whereas the third images is characterised by a difficult background and many thin vessels overlap and belongs to second dataset. All the shown results are obtained without any added seeds.

The strength of our method is that detects correctly most of the relevant vessels are segmented, even small ones. From the third line, it's obvious that the HT method has difficulties to segment the small vessels. The red circle in the first row shows a missed small vessel by the HT segmentation method. The same branch is also omitted by the FC and RG. This branch is also detected by the AQCA but not in totality.

The green circles in images of the second row demonstrate that FC, AQCA and our approach are robust when dealing with a moderated stenosis. The three other methods fail to reach distal parts located after the stenosis. However, the FC method suffers from some discontinuities and is sensitive to noise which is also a limitation of the HT segmen-tation method as we can see in the figures of the last column where the yellow circles show non-vascular structures erroneously segmented by both HT and FC. Contrarily, the RG, MSRG and AQCA are able to avoid noise far enough from vessels since they start from correct seeds. However, the RG method produces many false positives around the thin vessels. The AQCA achieves a good segmentation of the principal arteries but gives more false negatives than the MSRG.

Discussion and Conclusion

In this paper, we have presented a novel method for coronary artery segmentation based on Hessian features analysis in a multiscale fashion. We have introduced a growing rule combining both vesselness and direction information at different scales. First, the Frangi vesselness filter was adopted for vessel enhancement since it nicely captures tubular structure and significantly reduces background artefacts. Then, instead of extracting vascular structure only on the basis of their vesselness response, we have also incorporated the direction information. This information was included in a way that favours neighbours exhibiting directions close to the direction of a seed point even if their vesselness responses are lower than a given threshold. Thereby, the growth runs fluidly along the artery direction, and the number of broken branches is reduced as well. The multiscale region growing is robust, especially when dealing with discontinuity factors such as noise, stenosis and low contrast. Quantitative evaluation have been performed on real angiography images at two stages. First, the choice of the Frangi filter for vessel-like structures detection was justified by comparing the segmentation results obtained from three different filters. The proposed segmentation method has been also compared to other segmentation methods and it has been shown that it generates a better performance than fuzzy clustering, hysteresis thresholding and region growing in terms of vessel overlap, sensitivity and precision. The AQCA method, however, by using a multiscale edgness term that is able to accurately detect vessel border, reaches a precision greater than our method. But its sensitivity is smaller.

The results obtained by our algorithm (an average of DICE coefficient of 80% ± 5% for DS1 and 70% ± 5% for DS2 ) seem to represent a significant misalignment. Fig. 12 gives details of misalignment between the MSRG algorithm and the first ground truth. Almost the false negatives occur in vessels boundaries. For the large vessels a difference of 1 to 2 pixels between each ground truth and the segmentation is noted. Contrarily, false positives frequently appear on thin vessels which are generally overemphasized.

Our observations throughout the two datasets showed that the discordance is due to two main causes which are: 1) detection of non-vascular structures like background structures which are adjacent to vessels; and 2) false positives on vessel borders.

The first type of error is normally related with the vessel enhancement method used. In fact, the Frangi filter enhances not only the coronary tree but also any tubular structure including diaphragm, spinal disks, catheters and other background artefacts. Such non-vascular structures are erroneously segmented by our method only if they intersect or are very close to a vascular branch, otherwise they are enhanced by the Frangi filter but omitted by our segmentation method since the growth process starts from correctly classified seeds. Such error occurs in just few images of our dataset and it can be resolved if these structures are beforehand detected and eliminated from the images using an a priori knowledge for the shape of the diaphragm border [START_REF] Petkov | Diaphragm border detection in coronary x-ray angiographies: New method and applications[END_REF] or a learning-based method to detect the catheter [START_REF] Hernandez-Vela | Accurate coronary centerline extraction, caliber estimation, and catheter detection in angiographies[END_REF].

The second error consists in inaccurate border detection especially for thin arteries. In fact, these vessels may be overemphasized when they are processed at a scale larger then their sizes. This issue could be addressed if a break constraint is added to our algorithm in order to stop the growing when the vessel boundaries are reached. Obviously, this needs a robust vessel border detector such as the mulitscale Canny edge detector used in the AQCA method [START_REF] Hernandez-Vela | Accurate coronary centerline extraction, caliber estimation, and catheter detection in angiographies[END_REF].

In the future work, we would like to overcome these limitations by automatic detection of non vascular structures and better vessel border delineation. Future investigation of our segmentation method can be used as a computerassisted quantitative analysis of blood vessels from angiograms or as an input for other processing tasks such as centerline extraction [START_REF] Hernandez-Vela | Accurate coronary centerline extraction, caliber estimation, and catheter detection in angiographies[END_REF], and 3D coronary tree reconstruction [START_REF] Liao | 3-d reconstruction of the coronary artery tree from multiple views of a rotational x-ray angiography[END_REF].

Figure 1 :

 1 Figure 1: Vessel enhancement with two different filters: (a) original image; and images enhanced by the (b) Sato and (c) Frangi filters.

Figure 2 :

 2 Figure 2: Motivating sample: (a) original image; (b) Frangi vesselness response; (c) vesselness in the green circle; (d) corresponding first eigenvector field.

Figure 3 :

 3 Figure 3: MultiScale Region Growing Overview.

Figure 4 :

 4 Figure 4: Seed points selection : automatically generated points using the vesselness map are plotted in red color and manually added points are plotted in green color.

  This section presents experiments carried out to assess the behaviour of the proposed segmentation method. The materials used in our experiments is described in section 4.1. Section 4.2 introduces the performance measures used to evaluate the algorithm. The different parameters used in the experiments are described in section 4.3. Finally section 4.4 introduces the state of the art methods used in the evaluation.

Fig. 5 Figure 5 :

 55 Fig.5 depicts two different ground truth of an image belonging to the first dataset. The common part is plotted in white color in Fig.5(b) while the red color illustrates pixels detected only by the first human observer. The blue color shows pixels detected only by the second human observer.

Figure 6 :

 6 Figure 6: Vessel Enhancing Diffusion vs. Hessian-based filter : (a) original image; (b) Frangi vesselness response; (c) output of the VED filter; (d): corresponding vesselness image.

Figure 7 :

 7 Figure 7: Enhancement filter evaluation for DS1 relatively to the first human observer (first row) and the second human observer (second row) : Dice (left), sensitivity (middle) and precision (right).

Figure 8 :

 8 Figure 8: Enhancement filter evaluation for DS2 relatively to the first human observer (first row) and the second human observer (second row) : Dice (left), sensitivity (middle) and precision (right).

Figure 9 :

 9 Figure 9: Segmentation evaluation for DS1 relatively to the first human observer (first row) and the second human observer (second row) : dice (left), sensitivity (middle) and precision (right).

Figure 10 :

 10 Figure 10: Segmentation evaluation for DS2 relatively to the first human observer (first row) and the second human observer (second row) : dice (left), sensitivity (middle) and precision (right).

Figure 11 :

 11 Figure 11: Segmentation results. Top to down: real angiograms ; ground truth ; hysteresis thresholding; fuzzy clustering; region growing; AQCA and our method .

Figure 12 :

 12 Figure 12: Comparison between MSRG and the first ground truth: The white color depicts true positives. False negative points are marked in blue and false positives in red. The shown results are obtained with manual added seeds in images 1 and 3.

http://www.cvc.uab.es/ cgatta/CoronaryCenterlineExtraction.html

The vesselness proposed by Frangi et al. in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] is often considered as the current gold-standard for vessel enhancement [START_REF] Merveille | Tubular structure filtering by ranking orientation responses of path operators[END_REF]. The Sato filter [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF] is also a well-cited method in this field and is integrated in many medical tools such as Mevislab 1 and Vmtk 2 . The VED filter was presented as a diffusion based version of the Frangi filter and is gaining in popularity especially for vessel enhancement in 3D images. The results obtained by the VED filter, similarly to other diffusion-based methods cannot be directly compared to Hessian-based methods since the VED output is a smoothed image and not a vesselness image. Readers may note this difference in Fig. 6 where Fig. 6(c) is the smoothed image obtained by applying the VED filter to the original image in Fig. 6(a), whereas Fig. 6(b) is the output of the Frangi filter. Thus, we apply the Frangi filter to the VED result (Fig. 6(c)). The final output image is given by Fig. 6(d).

All these methods were tested using the same scale range already discussed. The parameters α and β were set to 1 for the Sato filter whereas the VED filter was applied using the default parameters recommend by the authors in [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF] (w = 25, s = 5 and ξ = 10 -2 ).

Segmentation methods

The MSRG algorithm was compared against four state-of-art segmentation approaches, which are hysteresis thresholding (HT), fuzzy clustering (FC), multi-seed region growing(RG) and accurate quantitative coronary artery segmentation (AQCA).

Classic segmentation method such as hysteresis thresholding (HT), fuzzy clustering (FC) and multi-seed region growing(RG) are often applied to the final Frangi filter obtained by the maximum Frangi vesselness response. They are used generally in cases in which no dense abnormalities are present or when the vessel segmentation is auxiliary to the main task such as centerline extraction and 3D reconstruction of the coronary tree from different projections [START_REF] Liao | 3-d reconstruction of the coronary artery tree from multiple views of a rotational x-ray angiography[END_REF][START_REF] Andriotis | A new method of three-dimensional coronary artery reconstruction from x-ray angiography: Validation against a virtual phantom and multislice computed tomography[END_REF].

Their output was thought to be sufficient to launch other applications. An interesting relatively recent work (AQCA) tailored to coronary artery extraction was proposed by Hernandez-Vela et al. in [START_REF] Hernández-Vela | Accurate and robust fully-automatic qca: method and numerical validation[END_REF] and extended in [START_REF] Hernandez-Vela | Accurate coronary centerline extraction, caliber estimation, and catheter detection in angiographies[END_REF].

The hysteresis thresholding is applied on the maximum Frangi vesselness response. This method retains connected components with all points that have a vesselness value greater than a low threshold and having at least one pixel with a vesselness value greater than a given threshold. We used the thresholds suggested in [START_REF] Liao | 3-d reconstruction of the coronary artery tree from multiple views of a rotational x-ray angiography[END_REF] which are around 90 and 96 percent of the accumulated histogram of the vesselness image.

The second method segments the maximum vesselness image into two partitions using standard Fuzzy c-means algorithm (FCM).The FCM algorithm classifies the image by grouping similar pixels in the image into clusters. This clustering is achieved by iteratively minimizing an objective function with respect to fuzzy membership.

The third method is based on simple region growing (RG) on the maximum vesselness image with multiple seeds.

Seeds and thresholds are manually set. The intial seed points for RG are the same used for the MSRG algorithm. No manual seed were added to any method.

The last method (AQCA) is also based on the maximum Frangi vesselness map but it uses a graph cuts algorithm 1 http://www.mevislab.de/ 2 www.vmtk.org/