Quantification of multicontrast vascular MR images with NLSnake, an active contour model: In vitro validation and in vivo evaluation

The definitive version of the full paper is available from: http://onlinelibrary.wiley.com/doi/10.1002/mrm.10722/abstract

Quantification of multi-contrast vascular MR Images with the NLSnake, an active contour model: in vitro validation and in vivo evaluation Catherine Desbleds Mansard 1, Emmanuelle P. Canet Soulas 1, Alfred Anwander 1, Linda Chaabane 2, Bruno Neyran 1, Jean-Michel Serfaty 1, Isabelle E. Magnin 1, Philippe C. Douek 1, Maciej Orkisz 1

1 CREATIS CNRS Research Unit 5515, affiliated to INSERM, Lyon, France. 2 Laboratoire RMN CNRS Research Unit 5012, Villeurbanne, France.

Grant sponsors: This work has been supported by Rhone-Alpes Region project ADeMo and by Incentive Concerted Action CIVAREM. It is within the scope of scientific topics of GdR PRC ISIS.

Author for correspondence:

Maciej Orkisz CREATIS, INSA de Lyon, bât. Blaise Pascal 69621 Villeurbanne Cedex, France.

e-mail: {maciej.orkisz,catherine.mansard}@creatis.insa-lyon.fr Fax: (33) 472438526 Phone: (33) 472438782

Short Title: Quantification of MR Images by Active Contour

26/03/13 – page 1 of 31

Abstract Vessel wall measurements from multi-contrast MRI provide information on plaque structure and evolution. This requires an extraction of numerous contours. In this paper a contourextraction method is proposed, which uses an active contour model (NLSnake) adapted for a wide range of MR vascular images. Its originality resides in length normalization for the purpose of deformation computation. This implementation confers to the model several interesting properties: simplified parameter tuning, fast convergence and minimum user interaction. The model can be initialized far from the boundaries of the region to be segmented, even by only one pixel. Accuracy and reproducibility of NLSnake endoluminal contours was assessed on vascular phantom MR angiography and on high-resolution in vitro MR images of rabbit aorta. In vivo evaluation was performed on rabbit and clinical data for both internal and external vessel wall contours. In phantoms with 95% stenoses NLSnake gave 94.3  3.8% and the accuracy was even better for milder stenoses. In rabbit aorta images variability between NLSnake and experts was less than inter-observer variability, while maximum intra-variability of NLSnake was equal to 1.25%. In conclusion, NLSnake was successfully applied to quantify the vessel lumen in multi-contrast MR images using unchanged parameters.

KEY WORDS image processing, active contour model, MR vascular images, quantification.

26/03/13 – page 2 of 31

INTRODUCTION Lumen narrowing determined by an angiographic technique is still the reference measurement for the evaluation of atherosclerosis. However, measurements in the vessel wall would be more predictive. High-resolution MRI has emerged as a powerful noninvasive technique for the evaluation of vessel wall abnormalities. Its usefulness for the in vitro and in vivo study of plaque evolution has been demonstrated in humans and in animal models [1-7]. Multi-contrast MRI with measurements of lumen and outer vessel wall circumferences, areas, wall thickness, contrast and signal-to-noise ratio in the vessel wall provides information on plaque structure and evolution [8, 9]. Up-to-now, these measurements are generally based on a manual extraction of numerous contours, requiring at least two medical experts to preserve objectivity. Thus, with the increasing number of MRI exams, a fast, automated post-processing algorithm is necessary for multi-center longitudinal studies. However, there are several difficulties due to the characteristics of the MR data. Acquisition is performed in a multi-contrast mode, i.e. typically MR angiography (MRA), high-resolution T1, T2 and proton density with lumen blood appearing black or white. Signal-to-noise and contrast-to-noise with the surrounding tissues highly depend on the acquisition parameters, i.e. spatial resolution, fat saturation…, and on vessel-wall composition. Vessel-wall signal is often heterogeneous in pathological cases. Thus, automatic inner and outer contour detection is a challenging task. Moreover, a minimal user interaction is required and the process should be fast and reproducible. We propose in this article a method based on an active contour model adapted for a wide range of MR vascular images. Active-contour models, also known as snakes [10], are the basis of most tools currently used for a semi-automatic segmentation of medical images [11]. An active contour is a curve evolving from an initial shape towards a final solution, under external-forces action and internal-forces reaction. This curve is usually approximated either by a set of nodes 26/03/13 – page 3 of 31

that form the vertices of a polygon [12], or by a parametric curve described using a polynomial (spline) function [13] or by a harmonic model [14]. The success of these models is due to their great adaptability: active contours can conform to a large variety of shapes and they include mechanisms for interactive correction of errors. However, the best results are obtained when the initialization is close to the actual location of the boundaries. It is difficult to find a trade-off between quick convergence towards the final solution, accuracy of this solution and its independence from the initialization. Hence, operator interaction is usually required for initialization, for tuning of the model's parameters or for corrections. In the vascular cross-sectional image analysis the interactive initialization of the active contour is usually done either by clicking several points close to the desired boundary [12, 15] or by manually setting the center and the radius of a circle approximating this boundary [16]. Several attempts to the automation of the initialization process have been published. A multi-scale scanning of the entire image in order to find the most circular shape [17] leads to full automation but may fail when the cross-sectional shape is deformed by the pathology. A pre-segmentation of the image, based on the Markov Random Fields theory, has also been proposed for the initialization of an active contour [18]. The initial nodes of the active contour can also be sought radially, as the maximum contrast points, starting from an interactively indicated point within the vessel lumen [13, 19]. In all these developments a distinct algorithm, more or less complicated, carries out the initialization. We have developed an original implementation of closed active contour that starts from a single point, without any additional initialization algorithm [20]. From the methodological point of view, its originality resides in the length normalization of the snake when computing the internal forces. From the point of view of medical interest, the user interaction is reduced to a single-point initialization, while the parameters are pre-tuned and do not need to be modified for a large variety of MR images. Another possible strategy would be to use level-sets approach 26/03/13 – page 4 of 31

[21] that is also able to make evolve active contours towards the final solution starting from a distant initialization. However, this choice is usually done to cope with complex-topology objects in low-noise images, while vascular cross-sectional topology is relatively simple and the high-resolution MR images often are noisy. In this article, a validation of our model was performed on both phantoms and in vitro vascular MR images. NLSnake was also evaluated on various in vivo MR images from animal and human studies.

METHODS AND MATERIALS An active contour is a curve evolving from an initial shape towards a final solution, under external-forces action and internal-forces reaction. The main external force attracts the curve towards the boundaries in the image, while the internal forces tend to preserve its expected shape. In the absence of more detailed a priori knowledge of the expected shape, the internal forces control the continuity and the smoothness of the resulting curve. The final solution corresponds to equilibrium between these forces, which is attained by an iterative minimization of a weighted sum of associated potential energies. Mathematical description of snakes Let the active contour be expressed by the following parametric curve: s, t   xs, t , ys, t  , T

(1)

where t represents the time and s  [0,1] is the arc-length parameter. Energy E   is associated with  . This energy is a sum of an internal component Eint   and an external component

Eext   :

26/03/13 – page 5 of 31

E    Eint    Eext   .

(2)

The energy Eint   , corresponding to the internal forces, imposes constraints on the first and second derivatives of the curve: Eint    E elast    E flex    



1

0

 s, t  ds   s 2

1

 2 s, t 

0

s 2



2

ds ,

(3)

where   0,1 controls its elasticity, while   0,1 controls its flexibility. Setting a large value of  prevents excessive local elongation of the contour, while a large value of  prevents strong local curvature (it has a smoothing effect). The potential energy Eext () corresponds to external forces. The link between potentials P  and forces is given by: F    P  .

(4)

There are two kinds of external forces: image force and balloon force, respectively weighted by the coefficients wg and wb . The image force is associated with the intensity gradient by means of the potentials Pg. The balloon force [22], associated with the potentials Pb, inflates the contour outwards. To see the usefulness of the balloon force, let us consider an initialization within the vessel lumen, i.e. in a uniform region. The intensity gradient is then almost zero and so is the force supposed to attract the contour towards the boundary. Without the balloon force, the contour would not grow, due to its internal tension force. The external forces can be expressed as: Fext    wg Pg    wbPb   ,

(5)

and the corresponding external energy is: 1

Eext ()   wg Pg ()  wb Pb () ds . 0

26/03/13 – page 6 of 31

(6)

The contour moves until equilibrium is reached between the external and internal forces, i.e. until the associated potential energy E   is minimum. This equilibrium state is expressed by the Euler-Lagrange equation:      2   2      Fext () .    s  s  s 2  s 2 

(7)

The motion towards the energy minimum is described by the Lagrange equation:



 2       2   2      Fext () .      t s  s  s 2  s 2  t 2

(8)

The physical meaning of  is the mass density of the contour, while the damping parameter  represents the viscosity of the environment. In practice however, these are two adjustable parameters of the model, as well as , , wb and wg. To simplify the behavior of the model  is generally set equal to zero. For the sake of computations, the model has to be discretized. This leads to representing the contour by a set of N points (i, k ) (nodes), parameterized by an arc-length parameter i (i = 1, 2… N) and by a time parameter k (iteration index). To this purpose, the partial spatial derivatives in the equations (3), (7) and (8) are approximated by finite differences:   i  1, k   i  1, k / 2hk , s

(9)

where hk is the distance between nodes, called spatial discretization step. Spatial and temporal discretization of the equation (8) leads to an iterative shifting of the nodes, according to the following evolution equation:

(i, k )   I  A     i, k  1  Fext  i, k  1  . 1

(10)

This equation corresponds to  = 0. The matrix A is a discretized formulation of the internal forces. Note that the damping parameter  controls the magnitude of the deformations (shift per iteration). 26/03/13 – page 7 of 31

The discretization step hk within the matrix A depends on the number of snake points N and, initially, is equal to the distance between the snake points. The evolution equation (10) is only valid if the snake points are equally spaced and the step hk is unchanged. However, after each iteration the snake length grows due to the external forces, and hk does not correspond to the real distance between the snake points. The internal energy, especially the term Eelast   associated with the contour's tension, grows with the snake's length and can stop the snake before the boundary is reached. This is particularly apparent if the initialization is far from the final position. In this case, the model needs to be resampled with a new (larger) number of snake points N, or/and with a new step hk+1. The N × N matrices,  I  A and its inverse, have to be recomputed, which is a time-consuming task. The number of these computations can be significantly reduced when using a modified version of the evolution equation (10) with a careful implementation that avoids numerical instability [23] . Originality of the NLSnake The originality of our implementation resides in: 1) initializing the model with a single pixel and 2) in computing the deformations on a normalized-length copy ' i, k  1 of the actual contour (fig. 1): ' i, k  1   i, k  1 / hk 1 .

(11)

The nodes are initially distributed in equal distances on the perimeter of a small circle (diameter = 1 pixel) around the initialization point. Note that one can have many nodes within a single pixel, as their locations are real numbers, while the pixels' locations are integers. The deformations are then computed on the normalized-length copy of the contour, according to the following equation:

26/03/13 – page 8 of 31

(12)

Figure 1: One iteration of the NLSnake deformation. Step (I): normalization of the current contour (k-1) to '(k-1). Step (II): deformation of the normalized contour. Step (III): addition of the deformation (k) to obtain the updated contour (k).

(i, k )   I  A    ' i, k  1  Fext i, k  1   ' (i, k  1) . 1

These deformations are applied directly to the actual contour  i, k  1 :

i, k   i, k  1  i, k  .

(13)

After the deformation, the nodes are evenly redistributed along the contour in order to ensure a spatially constant discretization step hk. Note that the matrix A is computed with a normalized discretization step h = 1. The normalization process only changes the length of the contour for the computation of the internal energy, while the external forces remain unchanged. The number of nodes N has to be large enough to permit a good approximation of the largest

26/03/13 – page 9 of 31

contours in the image series. In our experimentation, N was set equal to 20 and never changed. Unlike the existing implementations, the NLSnake keeps the number of nodes constant from one iteration to another, despite the growth of the actual contour  i, k  . Therefore, the N  N matrices,  I  A and its inverse, are computed only once, thus making the entire process very fast. Indeed, as the discretization step of the normalized-length contour ' i, k  is spatially and temporally constant (h = 1), the evolution equation (12) remains valid and the re-discretization is not required. Application to vessel-wall contour extraction The above-described NLSnake model is first used to extract the endoluminal contour. The initialization is interactive: the user clicks a point within the vessel lumen. Once the NLSnake has converged to the inner vessel-wall boundary it is stored, and then it can be pushed outwards by one pixel so as to initialize the outer-contour search. To avoid collapsing towards the endoluminal contour, the inwards-directed displacements are inhibited and the balloon force coefficient is increased. Thus the model only grows until it reaches the outer boundary of the vessel wall. The model parameters were empirically pre-tuned, i.e. the values that had given the best results on a training set of images were kept unchanged during all the hereafter described experimentation on multi-contrast images:  = 0.058,  0.02,  0.03, wg 0.4, wb 0.020 for the endoluminal contour and wb 0.025 for the outer contour. The gradient intensity was calculated by simple finite difference. Phantom data Vascular phantoms were manufactured using Computer Assisted Design [24]. They represent the endoluminal shape of vascular segments. Each of them has a reference diameter of 6 mm and comprises two stenoses with different shapes (circular, elliptic, semi-lunar, concentric or eccentric) and severities (50%, 75% and 95% of area reduction) [25].

26/03/13 – page 10 of 31

During image acquisition the phantoms were filled with a dilution of Gadolinium [26] corresponding to the arterial peak after intra-vascular injection of 0.1 mmol/kg GdDOTA (DOTAREM, Guerbet, France). CE MRA images of phantoms were acquired on a clinical MRI system (1.5 T, Vision, Siemens, Erlangen Germany) using a body coil and 3D Gradient Echo FISP (Fast Imaging with Steady state Precession) sequence used with two different sets of acquisition parameters (TR/TE 5 ms/2.1 ms, resolution: 0.780.781 mm3) and (TR/TE: 4.4 ms/1.4 ms, resolution: 0.780.780.75 mm3). Animal Data A group of watanabe heritable hyperlipidaemic (WHHL) rabbits underwent a comparative study of atherosclerotic-plaque evolution with and without lipid-rich diet [27]. High-resolution MR images of the thoracic aorta were first acquired in vivo. Animals were then sacrificed and their entire aorta, heart and kidneys were removed after fixation under perfusion with paraformaldehyde. Before proceeding to histological preparation, in vitro high-resolution MRI was performed. This study complied with our institutional guidelines for the care and use of laboratory animals. In vivo high-resolution MR images of the thoracic aorta were acquired on a clinical MRI system (1.5 T, Vision, Siemens, Erlangen Germany). T1 weighted MR images were obtained with a 2D spin-echo sequence (TR/TE=855 ms/20 ms) and a spatial resolution of 0.20.22 mm3. T2 and proton density images were acquired with a fast 2D spin-echo sequence with TR/TE=2500 ms/54 ms and TR/TE=2500 ms/15 ms, respectively, with an in plane spatial resolution of 390 m and 310 µm, respectively, and with a slice thickness of 2 mm. All MR images were acquired with fat saturation and ECG triggering. In vitro high-resolution MR imaging of the aorta was performed on a research MRI system (2 T, Oxford magnet, MRRS console). T1 and T2 weighted images were acquired with a multislice 2D spin-echo sequence with TR/TE: 600ms/21ms and 1800 ms/50 ms, respectively. In 26/03/13 – page 11 of 31

both cases, the spatial resolution was from 78 to 97m and the slice thickness varied from 0.8 to 1 mm, depending on the arterial wall size. Human Data Human data were drawn from a multi-center (10 hospital centers in France) clinical trial CARMEDAS [28] (supported by a grant of the French Public Health Ministry), aiming at the definition of the best cost-effective diagnostic strategy for the assessment of symptomatic and asymptomatic carotid stenoses. The main inclusion criterion is a significant stenosis in Doppler-ultrasound (greater than 50 % for symptomatic patients and greater than 60 % for asymptomatic patients). In all patients CE MRA and high-resolution MRI are performed. Data analysis Data analysis was first carried out on cross-sectional endoluminal contour area and shape in MRA images of phantoms and on 46 T2 in vitro high-resolution MR images of rabbit aorta. Then both internal and external vessel wall contours were evaluated on in vivo data. On animal in vivo data we quantitatively evaluated the endoluminal area and vessel-wall thickness. Human data were assessed qualitatively, by visual inspection, as a preliminary step towards a thorough evaluation. Accuracy and reproducibility Phantom data Accuracy of the automatically extracted contours was assessed by comparing their actual endoluminal diameter with the diameters measured in the reference sections (fig. 2) and by comparing theoretical and measured stenosis degrees of each phantom. The reproducibility was evaluated by measuring the endoluminal diameter variations within the reference sections. As all the reference sections have the same diameter 6 mm), the accuracy and reproducibility of the diameter measurements were evaluated simultaneously (mean, standard deviation and maximum difference). 26/03/13 – page 12 of 31

Figure 2: Reference and minimum cross-sections for stenosis quantification.

The stenosis severity was computed as the following ratio, involving the minimum crosssectional area S min and the average of the measured reference areas S reference (fig. 2):

%S 

S reference  S min 100% . S reference

(14)

Animal data Accuracy of the automatically extracted contours on the animal data was assessed by comparing them with manually traced contours (as done by [16]). Two experts did the manual tracing in a blinded manner. To simultaneously capture the reproducibility the automatic extraction was run three times with varying initialization. Both experts carried out the manual tracing twice with at least one-week delay. Agreement of the endoluminal contour shape was assessed on in vitro images by measuring distances between automatic and manual contours. The distances were measured radially, starting from each point of the NLSnake. The vessel-wall thickness was evaluated on in vivo data. We compared the average thickness calculated in 8 sectors defined starting from the endoluminal contour gravity center and from the sagittal orientation (fig. 7). Statistical Analysis To study the agreement of the area difference between manual tracing and the NLSnake, the statistical method described by Bland and Altman [29] has been used. The results calculated by

26/03/13 – page 13 of 31

the Bland and Altman method are the mean difference, the standard deviation of the difference and the 95% confidence interval. The normal distribution of the data was verified. Operating and processing time The operating time was compared with that needed by the experts for manual tracing of the contours on the same computer and for the same data set. It included interaction time, clicking in the vessel lumen (NLSnake initialization) or contour tracing, but also (in both scenarios) image loading from disk. The processing time was evaluated by an empirical comparison, on the same computer and for the same data, between the computational time of NLSnake and of a conventional snake with re-discretization followed by a numerical inversion of the system matrix at each iteration.

RESULTS Phantom data An example of NLSnake segmentation and quantification is shown in figure 3. The diameter of reference sections (6 mm, i.e. reference area of 28.27 mm²), measured by NLSnake, was equal to 6.01  0.10 mm (n = 144 slices). The maximum difference between the true and measured diameters was equal to 0.22 mm, i.e. 3.7 % of the true diameter. This difference represents 0.28 pixel compared to the image resolution. The results of stenosis quantification are displayed in figure 4. The absolute difference between the true and estimated stenosis degrees was less than 2% on average and the maximum difference was equal to 6.4%. The largest errors occurred for the most severe (95%) stenoses where the diameter was less than one pixel.

26/03/13 – page 14 of 31

b)

a) c) Figure 3: Segmentation/quantification of an MRA image of a phantom with two 95% stenoses: translucent rendering of the phantom with rings corresponding to the extracted contours (a), example of a contour extracted in a reference slice with automatically detected maximum and minimum diameter (b), resulting quantification curves (c) with abscissa representing the arc-length and two ordinate scales: area (left) and stenosis percentage (right).

26/03/13 – page 15 of 31

Figure 4: Results (mean and standard deviation) of stenosis quantification in phantoms, represented separately for three different stenosis degrees (50, 75 and 95%). The central box represents the values from the lower to upper quartile (25 to 75 percentile). The middle line represents the median. The horizontal lines represent the minimum and the maximum values.

In vitro animal data Examples of NLSnake contours superimposed onto in vitro images are displayed in figure 5. The inter-variability of area measurements between the experts was equal to 5.56% (n = 46 slices), while the intra-variability was equal to 2.30% and to 6.41% for experts 1 and 2 respectively. The Bland and Altman test was performed using the measurements done by the expert having the smallest intra-variability (fig.6). These results show a good agreement between the NLSnake measurement and the manual tracing (mean difference: 0.065 mm², 95% confidence interval: from -0.25 to 0.12 mm², agreement interval: from -1.26 to 1.13 mm²). The reproducibility of NLSnake is expressed through its intra-variability: 1.23%. 26/03/13 – page 16 of 31

a) b)

d) c)

e) f) Figure 5: Examples of contours extracted with NLSnake in high resolution MR in vitro T2-weighted images of rabbit aorta for various signal-to-noise levels (a = 8.15, b = 8.05, c = 6.44, d = 6.00, e = 5.76, f = 3.71) and lumen-shape complexities 26/03/13 – page 17 of 31

Figure 6: Results of the Bland and Altman test between the Expert1 and NLSnake on endoluminal contours area (mm2) of 46 in vitro slices.

The average distance between NLSnake and the contours traced by the expert was equal to 0.67  0.64 pixels (n = 920 points). The largest distance was equal to 4.7 pixels, which corresponds to 272 m. In vivo data Examples of internal and external NLSnake contours superimposed onto in vivo images of rabbit aorta are shown in figure 7. The inter-variability and the intra-variability of the two criteria on T2-weighted images (n = 20 slices) are presented in table 1. The Bland and Altman test (fig.8) on endoluminal area measured by the expert with the smallest intra-variability shows a good agreement between the NLSnake measure and the manual tracing (mean difference: 0.187 mm², 95% confidence interval: from –0.62 to 0.99 mm², agreement interval: from –3.6 to 3.2 mm²). The reproducibility of NLSnake is expressed through its intra-variability: 1.25%.

26/03/13 – page 18 of 31

a)

c)

b)

d) Figure 7: NLSnake applied to internal contour (plain line) and external contour (dots) extraction in multi-contrast in vivo images of a rabbit aorta: a) T1, SNR = 3.01 b) T2, SNR = 3.63, c) proton-density, SNR = 2.83, d) sectors used for local measurements, superimposed onto a T2 image, SNR = 3.55.

Figure 9 shows typical results obtained on a human data set: CE MRA and high-resolution MRI. Operating Time For 15 in vitro images (slices) from a single vascular segment the operating time of the NLSnake (endoluminal contour) extraction was 30 seconds, compared to 20 minutes required by the manual tracing. Figure 10 illustrates an empirical comparison, on the same computer, between the computational time of the NLSnake and of a conventional snake.

26/03/13 – page 19 of 31

Figure 8: Results of the Bland and Altman test between the Expert1 and NLSnake on endoluminal contours area (mm2) of 20 in vivo high resolution MR images of rabbit aorta.

Table 1: Intra- and inter-variability on T2-weighted in vivo images (n = 20 slices). Criteria

Intra-

Intra-

Intra-

Inter-

Inter-

variability

variability

variability

variability varibility

Expert1

Expert 2

NLSnake

Experts

ExpertNLSnake

Endoluminal area

4.5%

Sectorial vessel wall 38.5%

9.26%

1,25%

11.3%

5.52%

41.44%

26.7%

47.2%

36.3%

thickness

26/03/13 – page 20 of 31

b)

a)

c) Figure 9: NLSnake applied to a human data set: a) – b) surface-rendered MRA image of a patient's carotid artery with rings corresponding to the extracted contours and the resulting stenosis quantification curves with abscissa scaled in mm and two ordinate scales: area in mm² (left) and stenosis percentage (right), c) high resolution MR T2 weighted in vivo image of a patient's carotid artery (internal and external boundaries extracted by NLSnake), SNR = 4.26.

DISCUSSION An active contour model (NLSnake) was implemented for the purpose of semi-automatic segmentation of vascular MR images. It was validated on endoluminal boundary segmentation in MRA images of phantoms and in high-resolution MR T2-weighted in vitro images of rabbit 26/03/13 – page 21 of 31

aorta. It was further evaluated on internal and external vessel-wall boundary extraction in multicontrast in vivo images of rabbit aorta and of human carotid arteries. An automated segmentation tool is useful if the results are reproducible, accurate and available more quickly than with an interactive tool. NLSnake demonstrated its reproducibility in two ways: 1) its intra-variability with different initializations was significantly lower than the experts’ intra-variability, 2) cross-sectional area measurements along the reference sections (constant diameter) of the phantoms were very stable. As for the accuracy, variations between areas measured with NLSnake and those measured by the experts were not larger than the variations between the experts. Differences between measurements with NLSnake and true dimensions of the phantoms were significantly smaller than the variations between NLSnake and the experts in high-resolution MR vascular images. This can be explained both by a better contrast-to-noise ratio in the MRA images and by the variability of the experts on highresolution MRI. In particular, detection of the outer vessel wall boundary on in vivo highresolution MRI is difficult due to local absence of contrast and to low signal-to-noise ratio. This partly explains the large variability of wall thickness measurements. Furthermore, this thickness typically corresponds to 3 or 4 pixels, i.e. placing the contour at one-pixel distance inside or beyond the actual boundary, all along its circumference, leads to 25% error, at least. There are two timesaving mechanisms, when using NLSnake: reduction of user interaction and fast convergence of the algorithm. Typical interaction scenario of active-contour initialization consists of several mouse-clicks that roughly place the snake's nodes along the boundary to be extracted. Obviously, a single click within the vessel lumen is faster than manual tracing, even approximate, of an entire contour. Further reduction of the interaction is achieved by using the current contour's center to initialize the contours in the neighboring slices and so on. In some existing general-purpose software tools the snakes in the neighboring slices are initialized by the final contour from the current slice. However, for oblique vessels, the contour location from 26/03/13 – page 22 of 31

one slice to another may significantly vary in the context of typical vascular high-resolution MRI, where the slice thickness is not negligible compared to the vessel diameter. Nevertheless, the projection of the current contour center usually remains within the lumen boundaries in the neighboring slices and NLSnake is able to correctly expand from this point to the actual boundary. Furthermore, the timesaving increases significantly with the length of the vessel (number of slices) to be segmented, e.g. for a sequence of 15 slices, interaction and computation time with NLSnake has been 40 times shorter (30 seconds) than with manual tracing (20 minutes).

Figure 10: Comparison of the computational time (on a PC Pentium III 600 MHz, Windows NT) between NLSnake deformation and a conventional implementation of active contour (with re-computation of the matrix  I  A and of its inverse at each rediscretization) 26/03/13 – page 23 of 31

Computational time is negligible, compared to the interaction time. Nevertheless, it is to be mentioned that NLSnake converges faster than the conventional implementations of snakes and the gain quickly increases with the initial number of nodes (fig.10). Indeed, without lengthnormalization, the snake has to be frequently re-discretized during its growth, as the evolution equation (10) only holds when the discretization step h is constant. The system matrices

 I  A and its inverse, have to be re-computed when the number of nodes is changed. Difficulty in using active contour models resides in the choice of weighting coefficient values

, , wb and wg of each part of the energy (3), (6). In conventional implementations of the snakes this choice is particularly difficult for the energy parts associated with two opposed forces: tension and balloon force, as these energies respectively depend on the snake length and number of nodes. In large cross-sections the contour growth may be stopped before reaching the boundary if the balloon coefficient is not large enough. Indeed, equilibrium between the balloon force and the internal tension is then reached too early. Conversely, in a small cross-section, the balloon force may blow the contour beyond the boundary, if its coefficient is too large. The damping factor  has also to be appropriately chosen according to the current discretization, as its value determines the growth speed of the contour (shift per iteration). It should make fast convergence possible but avoid stepping over the boundary (if the shift per iteration is too large). With NLSnake the number of nodes and the contour length that is used to compute the evolution are constant. Therefore, it is easier to set the coefficients. In particular, the damping factor can be chosen so that the shift magnitude is approximately one pixel per iteration. Consequently, the same set of experimentally determined parameters was successfully used for vessel wall boundaries segmentation in various dark blood and white blood images such as MRA and multi-contrast high-resolution MRI (T1, T2, proton density) in vitro and in vivo, in animals and in humans.

26/03/13 – page 24 of 31

Our objective is the extraction of both boundaries, internal and external, of the vascular wall in high-resolution MRI. Two interdependent NLSnakes are used. The internal boundary usually has better contrast, hence it is extracted at first. It is then used to initialize the external contour. The resulting segmentation (fig. 7) gives access to local thickness measurements and is a prerequisite for automated studies of plaque composition [8, 15, 16]. However, reliability of the external-boundary extraction is limited by local absence of contrast and by the low signal-tonoise ratio of the vessel wall in vivo. This is confirmed by the large variability of the vessel wall external boundary tracing done by the experts. Consequently, although the pre-tuned model parameters give satisfactory results for the endoluminal boundary in a wide range of images it may be useful to seek a specific set of model parameters adapted to the external boundary. Furthermore, in the cases where only a small part of the boundary is contrasted, an applicationspecific growth-stopping criterion is to be found in order to replace the gradient-based boundary detection criterion.

CONCLUSION We have developed a tool for the extraction of vessel wall boundaries from MR images, based on an original implementation of active contour model (NLSnake). One of the interesting aspects of this tool is minimum user-interaction. The contour-length normalization process used for the computation of the active-contour forces makes the model easier to use. NLSnake can therefore be used in a wide variety of images (MRA, multi-contrast high-resolution MRI with various spatial resolutions and SNR) with a fixed parameter set. Furthermore, initialization by a single pixel is highly timesaving in the context of blood-vessel segmentation and quantification, where large series of planar images along the vessels are involved. Although the model is initialized far from the final solution, the results are reproducible and accurate. However, the accuracy of the outer boundary extraction decreases in very noisy images of small arteries. The 26/03/13 – page 25 of 31

extraction of the vessel-wall boundaries is the first step towards an analysis of its composition from multi-contrast and dynamic contrast high-resolution MRI. Reliability of the outer boundary detection could be improved by finding a specific parameter set or a specific edgedetection criterion better suited to the external boundaries. Another strategy will also be considered, which simultaneously exploits all the available multi-contrast images of the same slice, provided that registration of these images can previously be carried out.

REFERENCES

[1]

C. Yuan, S.-X. Zhang, N. L. Polissar, D. Echelard, G. Ortiz, J. W. Davis, E. Ellington, M. S. Ferguson, and T. S. Hatsukami, “Identification of fibrous cap rupture with Magnetic Resonance Imaging is highly associated with recent transient ischemic attack or stroke,” Circulation, vol. 105, pp. 181-185, 2002.

[2]

R. Corti, V. Fuster, Z. A. Fayad, S. G. Worthley, G. Helft, D. Smith, J. Weinberger, J. Wentzel, G. Mizsei, M. Mercuri, and J. J. Badimon, “Lipid lowering by simvastatin induces regression of human atherosclerotic lesions: two years' follow-up by HighResolution noninvasive Magnetic Resonance Imaging,” Circulation, vol. 106, pp. 28842887, 2002.

[3]

W. Y. Kim, M. Stuber, P. Bornert, K. V. Kissinger, W. J. Manning, and R. M. Botnar, “Three-dimensional black-blood cardiac Magnetic Resonance coronary vessel wall imaging detects positive arterial remodeling in patients with nonsignificant coronary artery disease,” Circulation, vol. 106, pp. 296-299, 2002.

[4]

C. M. Kramer, “Magnetic resonance imaging to identify the high-risk plaque,” Am J Cardiol, vol. 90, pp. L15-L17, 2002. 26/03/13 – page 26 of 31

[5]

J. M. Serfaty, L. Chaabane, A. Tabib, J. M. Chevallier, A. Briguet, and P. C. Douek, “Atherosclerotic plaques: classification and characterization with T2-weighted highspatial-resolution MR Imaging - an in vitro study,” Radiology, vol. 219, pp. 403-410, 2001.

[6]

L. Chaabane, E. Canet Soulas, F. Contard, A. Salah, D. Guerrier, and P. C. Douek, “High Resolution Magnetic Resonance imaginig at 2T: potential for artherosclerotic lesions exploration in the apolipoprotein E-knockout mouse,” Invest Radiol, 2003.

[7]

B. D. Coombs, J. H. Rapp, P. C. Ursell, L. M. Reilly, and D. Saloner, “Structure of plaque at carotid bifurcation. High resolution MRI with histological correlation,” Stroke, vol. 32, pp. 2516-2521, 2001.

[8]

X. Kang, N. L. Polissar, C. Han, E. Lin, and C. Yuan, “Analysis of the measurement precision of arterial lumen and wall areas using high-resolution MRI,” Magn Reson Med, vol. 44, pp. 968-972, 2000.

[9]

D. Xu, J.-N. Hwang, and C. Yuan, “Atherosclerotic plaque segmentation at human carotid artery based on multiple contrast weighting MR images,” presented at Int. Conf. Image Processing, 2001.

[10]

M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active Contour Models,” Int J Comput Vision, vol. 1, pp. 321-331, 1987.

[11]

T. McInerney and D. Terzopoulos, “Deformable models in medical image analysis: a survey,” Med Image Analysis, vol. 1, pp. 91-108, 1996.

[12]

J. Thomas, B. Rutt, H. Ladak, and D. Steinman, “Effect of black blood MR image quality on vessel wall segmentation,” Magn Reson Med, vol. 46, pp. 299-304, 2001.

[13]

A. Sebbahi, A. Herment, A. de Cesare, and E. Mousseaux, “Multimodality cardiovascular image segmentation using a deformable contour model,” Comput Med Imaging Graph, vol. 21, pp. 79-89, 1997. 26/03/13 – page 27 of 31

[14]

L. H. Staib and J. S. Duncan, “Boundary finding with parametrically deformable models,” IEEE Pattern Analysis Machine Intell., vol. 14, pp. 1061-1075, 1992.

[15]

S. Zhang, T. S. Hatsukami, N. L. Polissar, C. Han, and C. Yuan, “Comparison of carotid vessel wall area measurements using three different contrast-weighted black blood MR imaging techniques,” Magn Reson Imaging, vol. 19, pp. 795-802, 2001.

[16]

C. Yuan, E. Lin, J. Millard, and J.-N. Hwang, “Closed contour edge detection of blood vessel lumen and outer wall boundaries in black-blood MR images,” Magn Reson Imaging, vol. 17, pp. 257-266, 1999.

[17]

D. Rueckert, P. Burger, S. M. Forbat, R. D. Mohiaddin, and G. Z. Yang, “Automatic tracking of the aorta in cardiovascular MR images using deformable models,” IEEE Trans. Med. Imaging, vol. 16, pp. 581-590, 1997.

[18]

W. Kerwin, C. Han, B. Chu, D. Xu, Y. Luo, J.-N. Hwang, T. Hatsukami, and C. Yuan, “A quantitative vascular analysis system for evaluation of atherosclerosis lesions by MRI,” presented at MICCAI, Utrecht (NL), 2001.

[19]

S. Kozerke, R. Botnar, S. Oyre, M. Scheidegger, E. Pedersen, and P. Boesiger, “Automatic vessel segmentation using active contours in cine phase contrast flow measurements,” J Magn Reson Imaging, vol. 10, pp. 41-51, 1999.

[20]

C. Desbleds-Mansard, A. Anwander, L. Chaabane, M. Orkisz, B. Neyran, P. C. Douek, and I. E. Magnin, “Size independent active contour model for blood vessel lumen quantification in High-Resolution Magnetic Resonance Images,” presented at MICCAI, Utrecht (NL), 2001.

[21]

K. C. Wang, R. W. Dutton, and C. A. Taylor, “Level sets for vascular model construction in computational hemodynamics,” IEEE Engineering in Med and Biol, vol. 18, pp. 33-39, 1999.

26/03/13 – page 28 of 31

[22]

L. D. Cohen, “On active contour models and balloons,” Comput Vision Graphics Image Process: Image Understand, vol. 53, pp. 211-218, 1991.

[23]

N. Ray, B. Chanda, and J. Das, “A fast and flexible multiresolution snake with a definite termination criterion,” Pattern Recogn, vol. 34, pp. 1483-1490, 2001.

[24]

C. P. Renaudin, B. Barbier, R. Roriz, D. Revel, and M. Amiel, “Coronary arteries: new design for three-dimensional arterial phantom,” Radiology, vol. 190, pp. 579-582, 1994.

[25]

M. Hernández-Hoyos, M. Orkisz, P. Puech, C. Mansard-Desbleds, P. C. Douek, and M. I.E.,

“Computer

assisted

analysis

of

three-dimensional

MR

angiograms,”

RadioGraphics, vol. 22, pp. 421-436, 2002. [26]

B. Marchand, P. C. Douek, P. Robert, C. Corot, J. P. Roux, P. Adeleine, M. HernándezHoyos, Y. Crémillieux, M. Orkisz, and E. Canet, “Standardized MR protocol for the evaluation of MRA sequences and/or contrast agents effects in high-degree arterial stenosis analysis,” MAGMA, vol. 14, pp. 259-267, 2002.

[27]

L. Chaabane, F. Contard, P. C. Douek, C. Corot, D. Guerrier, A. Briguet, and E. Canet, “MR contrast-enhancement of atherosclerotic plaque in watanabe rabbits: effects of a lipid-rich diet,” presented at ISMRM 10th Scientific Meeting, Honolulu, Hawaii, USA, 2002.

[28]

M. Nonent, J.-M. Serfaty, C. Pachai, J. Linard, A. Pasco-Papon, J.-F. Heautot, Y. Le Bras, I. Jars, V. Buthion, M. Lamure, and P. Douek, “Intermediate results of the CARMEDAS Multicenter study: Comparison of different imaging strategies for the diagnosis of carotid artery stenosis,” presented at XIV Annual International Workshop on MR Angiography, Essen, Germany, 2002.

[29]

J. Bland and D. Altman, “Statistical method for assessing agreement between two methods of clinical measurement,” The Lancet, vol. I, pp. 307-310, 1986.

26/03/13 – page 29 of 31

Figure captions Figure 11: One iteration of the NLSnake deformation. Step (I): normalization of the current contour (k-1) to '(k-1). Step (II): deformation of the normalized contour. Step (III): addition of the deformation (k) to obtain the updated contour (k). Figure 12: Reference and minimum cross-sections for stenosis quantification. Figure 13: Segmentation/quantification of an MRA image of a phantom with two 95% stenoses: translucent rendering of the phantom with rings corresponding to the extracted contours (a), example of a contour extracted in a reference slice with automatically detected maximum and minimum diameter (b), resulting quantification curves (c) with abscissa representing the arc-length and two ordinate scales: area (left) and stenosis percentage (right). Figure 14: Results (mean and standard deviation) of stenosis quantification in phantoms, represented separately for three different stenosis degrees (50, 75 and 95%). The central box represents the values from the lower to upper quartile (25 to 75 percentile). The middle line represents the median. The horizontal lines represent the minimum and the maximum values. Figure 15: Examples of contours extracted with NLSnake in high resolution MR in vitro T2weighted images of rabbit aorta for various signal-to-noise levels (a = 8.15, b = 8.05, c = 6.44, d = 6.00, e = 5.76, f = 3.71) and lumen-shape complexities Figure 16: Results of the Bland and Altman test between the Expert1 and NLSnake on endoluminal contours area (mm2) of 46 in vitro slices. Figure 17: NLSnake applied to internal contour (plain line) and external contour (dots) extraction in multi-contrast in vivo images of a rabbit aorta: a) T1, SNR = 3.01 b) T2, SNR = 3.63, c) proton-density, SNR = 2.83, d) sectors used for local measurements, superimposed onto a T2 image, SNR = 3.55. Figure 18: Results of the Bland and Altman test between the Expert1 and NLSnake on endoluminal contours area (mm2) of 20 in vivo high resolution MR images of rabbit aorta. 26/03/13 – page 30 of 31

Figure 19: NLSnake applied to a human data set: a) – b) surface-rendered MRA image of a patient's carotid artery with rings corresponding to the extracted contours and the resulting stenosis quantification curves with abscissa scaled in mm and two ordinate scales: area in mm² (left) and stenosis percentage (right), c) high resolution MR T2 weighted in vivo image of a patient's carotid artery (internal and external boundaries extracted by NLSnake), SNR = 4.26. Figure 20: Comparison of the computational time (on a PC Pentium III 600 MHz, Windows NT) between NLSnake deformation and a conventional implementation of active contour (with re-computation of the matrix  I  A and of its inverse at each re-discretization)

Table 1: Intra- and inter-variability on T2-weighted in vivo images (n = 20 slices).

26/03/13 – page 31 of 31