Endocardial Border Detection in Cardiac Magnetic Resonance Images Using Level Set Method

Initialization with Binary Mask

The method that we propose in this paper is based on the initialization of the level set algorithm by using the binary image obtained from the thresholding operation. For the localization of the left ventricular cavity, we use a method based on the evaluation of a simple metric indicating the roundness of an object because we approximate that the left ventricular cavity has a round shape. To determine if an object is round, we estimate each object area and perimeter and this for all unconnected regions on the image. The metric roundness was defined by the following equation:

1

This metric is equal to one only for a circle and is less than one for any other shape.

On the short-axis images, the LV endocardial border is not always defined at locations with maximum intensity gradients. That is why we have proposed to smooth the input image with a Gaussian kernel filter to homogenate the LV endocardial border. In Fig. 1, we show an example that illustrates the difference between the left ventricular cavity detected from the original image and from the smoothing image.

We can clearly notice that the LV metric obtained before smoothing (0.6478) is less than the LV metric obtained after smoothing (0.8638). Figure 2 shows the metric of the LV for 18 images before and after smoothing. The proposed segmentation method in this paper starts by thresholding the input image using OTSU method. After this step, we eliminate the smallest regions (see Fig. 3).

Image Segmentation

In image segmentation, the methods based on active contours present dynamic curves that move toward object’s boundaries. To achieve this goal, these approaches explicitly define an external energy that can move the zero-level curve toward the object’s boundaries [12].

For our segmentation process, we propose to use a variational level set formulation of active contours without re-initialization. We describe here this method [12]. Let I be an image and g be the edge indicator function defined by:

2

where G_σ is the Gaussian kernel with standard deviation σ. We define an external energy for a function ϕ(x,y) as below:

3

where λ > 0 and ν are constants and the terms L_g(ϕ) and A_g(ϕ) are defined by

4

And

5

where δ is the univariate Dirac function and H is the Heaviside function. Now, we define the following total energy functional

6

7

The external energy drives the zero-level set toward the object boundaries, while the internal energy μP(ϕ) penalizes the deviation of ϕ from a signed distance function during its evolution. By variations calculating, the gâteaux [13] derivative (first variation) of the functional ε in Eq. 6 can be written as:

8

where Δ is the Laplacian operator. Therefore, the function ϕ that minimizes this functional satisfies the Euler–Lagrange equation . The steepest descent process for minimization of the functional ε is the following gradient flow:

9

This gradient flow is the evolution equation of the level set function in the proposed method. The second and the third term in the right-hand side of Eq. 9 correspond to the gradient flows of the energy functional λL_g(ϕ) and νA_g(ϕ), respectively, and are responsible of driving the zero-level curve toward the object boundaries. To explain the effect of the first term, which is associated to the internal energy μP(ϕ), we notice that the gradient flow Eq. 10 has the factor as diffusion rate.

10

If |Δϕ| >1, the diffusion rate is positive and the effect of this term is the usual diffusion, i.e., making φ more even and therefore reduce the gradient |Δϕ|. If |Δϕ| <1, the term has effect of reverse diffusion and therefore increases the gradient.

Image Acquisition

For our experiments, we used a free database available for research purposes [15]. The original images in this database consist of short-axis 2D sequences (between eight and 13 sequences of 25 images per patient), with breath-held, ECG-gated acquisitions from base to apex. The most basal slices included in the analysis were located just above mitral valve within LVC. To be included, the basal myocardium had to be visible in the entire circumference at end-systole. The most apical slice was chosen as the one with the smallest visible LVC at end-systole. The sequences were registered on the heart cycle, and thus, they can be stacked to construct 3D sequences. Each 3D sequence is made of 25 volumes and covers a complete heart cycle. A conventional 2D manual segmentation of the LV myocardium was performed by two independent and blinded expert cardiologists. They used a software package that is well established for the post-processing of medical images (Analyze R, Biomedical Imaging Resource, Mayo Clinic Foundation, Rochester, MN, USA). These two experts are called E1 and E2 in the sequel. The cine MR dataset was analyzed as a succession of 2D LV short-axis planes. For each slice location, the experts manually overlaid the endocardial and epicardial contours both at end-diastolic and end-systolic times. During manual tracing, papillary muscles and LV trabeculae were included within the LV myocardium. Then, the segmented slices were stacked to rebuild a 3D object for quantification. The time spent by the experts to manually segment one volume of a 3D + t sequence ranged from 15 to 20 min. For this reason, a manual segmentation is not available at every time-step.

Endocardial Border Detection

In this section, we present the results of endocardial border detection in cardiac magnetic resonance images using our segmentation approach. The proposed method produces good results, but sometimes some intensity inhomogeneities often occur in left ventricular MRI images and may cause considerable difficulties in image segmentation. As we show in Fig. 4a, heterogeneities in discussion can appear inside of the LVC or be attached with the myocardium muscle. In order to overcome the difficulties caused by intensity heterogeneities, we convoluted the original image with a Gaussian kernel filter. A typical example is given in Fig. 4 with (σ = 1). We notice here that it is necessary to examine the influence of the parameter sigma on the segmentation results.

For this purpose, we apply our method using different sigma parameters. The most accurate segmentation result is obtained with the parameter σ = 1.

Evaluation

The segmentation results were evaluated quantitatively using an error metric [16, 17].

11

where A is a binary images such all pixels inside the curves produced by a clinical expert and B is the set of all pixels inside the curves produced by a segmentation algorithm.

• ε = 3.5: the parameter in the definition of smoothed Dirac function

• μ = 0.1: coefficient of the internal (penalizing) energy term P(ϕ)

• λ = 1: coefficient of the weighted length term L_g(ϕ)

• ν = −0.5: coefficient of the weighted area term A_g(ϕ)

The Fig. 5 shows the results of the two expert segmentation (b and c) and illustrates the method segmentation results (d). The Fig. 6 shows the results of the segmentation method (d) and illustrates the initial and the final contour (the endocardial border) (b and c).

As shown in Fig. 6, the proposed method works efficiently. Others segmentation results are given in Fig. 7.

Table 1 gives the error metric of 18 patient selected from the data bases used in our study. The error metric was evaluated between the proposed method and the two experts (E1 for expert1 and E2 for expert2). As presented in Table 1, we obtained an overall average similarity area of 97.85% (expert1) and 97.89% (expert2).

We have also evaluated this error for the automatic algorithm proposed in [6], where the authors developed an automatic algorithm based on watershed segmentation of 4D cardiac MR images. The automatic segmentation results are available in the database used in this work. Table 2 gives the results of the error metric for the same patients selected in Table 1 at the end-diastole phase for this watershed segmentation algorithm.

We notice also that we obtained an overall average similarity area of 95.93% (expert1) for the automatic method proposed in [6] (Table 2) and 96.75% for expert2. The agreement between manual and the proposed method was assessed on LV area measures, as reported in Table 1. In order to gain a more meaningful data interpretation, the Bland–Altman [18] plots for the expert1 and the proposed method (Fig. 8a) and the expert2 and the proposed method (Fig. 8b) are shown in Fig. 8.

The Bland–Altman plot, or difference plot, is a graphical method used to compare two measurements techniques. In this graphical method, the differences (or alternatively the ratios) between the two techniques are plotted against the averages of the two techniques. Horizontal lines are drawn at the mean difference and at the limits of agreement, which are defined as the mean difference plus and minus 2 times the standard deviation of the differences [19].

We see clearly, from the tow plots, Fig. 8, the agreement between the two experts and the proposed method. This implies that the LV area measured with the new proposed method is in agreement with the LV area obtained with the two experts. The inter-observer variability was shown in Fig. 8c.

To show the robustness of our method, we illustrate the examples presented in Fig. 9. For all the examples, the inclusion of the left ventricular and the left myocardium and pathologies in the myocardium makes the separation of the two regions using any algorithms very difficult.

Generally pathologies in the myocardium cause a non-circular shape of the left ventricular. Meanwhile based on the metric used in this work, we can say that left ventricular cavity has the biggest value. For more explication, we illustrate the following example.

As shown in Fig. 10, the result of the thresholding operation shows three regions: in the middle the left (0.46) and the right (0.37) ventricular and the other region (0.13). So we see clearly that the left ventricular have not exactly a circular shape but they have the biggest value.

We have test our algorithm using MRI images with contrast enhancement and perfusion. Images were obtained with sequences with considerable less S/N and more inhomogeneities (Figs. 11 and 12). The results shown in Figs. 11 and 12 prove that our algorithm can successfully be implemented with this kind of image.