Automatic segmentation of myocardium from black-blood MR images using entropy and local neighborhood information.

By using entropy and local neighborhood information, we present in this study a robust adaptive Gaussian regularizing Chan-Vese (CV) model to segment the myocardium from magnetic resonance images with intensity inhomogeneity. By utilizing the circular Hough transformation (CHT) our model is able to detect epicardial and endocardial contours of the left ventricle (LV) as circles automatically, and the circles are used as the initialization. In the cost functional of our model, the interior and exterior energies are weighted by the entropy to improve the robustness of the evolving curve. Local neighborhood information is used to evolve the level set function to reduce the impact of the heterogeneity inside the regions and to improve the segmentation accuracy. An adaptive window is utilized to reduce the sensitivity to initialization. The Gaussian kernel is used to regularize the level set function, which can not only ensure the smoothness and stability of the level set function, but also eliminate the traditional Euclidean length term and re-initialization. Extensive validation of the proposed method on patient data demonstrates its superior performance over other state-of-the-art methods.

Introduction

According to the World Health Organization, an estimated 17.5 million people died from cardiovascular diseases in 2005, representing 30% of all global deaths [1]. Being able to provide an early diagnosis and treatment will dramatically reduce this death toll. Recent advances in novel imaging and computing technology and their introduction into clinical routine have shown tremendous potential towards achieving such an ambitious goal. Over the diverse range of imaging modalities, cardiovascular magnetic resonance (CMR) imaging is a unique technique that is ionizing radiation free and can provide clear anatomy of the heart. Early iron detection using myocardial T2* values derived from CMR was developed and validated in large patient populations [2]. Extracting the anatomical information of the cardiac is the fundamental step for the development of clinical applications, and obtaining reproducible and unbiased quantitative measurement of the anatomy is indeed of central importance for the success of these applications.

Manual delineation in all the images of a subject is prohibitively a labor intensive, tedious and time-consuming task. Automatic or semi-automatic algorithms are therefore highly desired. Several approaches have been reported in the literatures for medical image segmentation. These approaches can be roughly categorized into three groups, namely, clustering algorithms [3,4,5], Markov random field (MRF) methods [6,7,8], and active contour-based methods [9,10,11,12,13,14]. Clustering methods iterate between segmenting the image and characterizing the properties of each class. Three commonly used clustering algorithms are the k-means [3], fuzzy c-means algorithm [4] and the expectation- maximization (EM) algorithm [5]. Without incorporating spatial information, these methods are sensitive to noise and image intensity inhomogeneity. MRF methods [6] incorporate spatial information using the Bayesian model, which is robust to noise. The segmentation is obtained using iterative methods such as iterated conditional modes [15] or simulated annealing [16]. However, a difficulty associated with MRF models is proper selection of the parameters controlling the strength of spatial interactions. Active contour models (ACMs) have been widely used for automatic image segmentation and object tracking [17,18]. The basic idea of the ACM is to evolve a curve to fit the desired object boundary according to predefined energy function [19,20]. According to the type of adopted image features, the existing active contour models can be broadly divided into edge-driven models [9,10,11] and region-driven models [12,13,14].

For edge-driven models, energy function is constructed using edge-based image features, i.e., local intensity gradients, outputs of edge indicator, and so on. This type of methods can be applied to images with intensity in homogeneities. However, these models are not only very sensitive to image noise, but also difficult to detect weak boundaries. Moreover, the segmentation result is highly dependent on the initial contour placement.

Alternatively, region-driven models depend on the region-based image features, i.e., some statistics of the primary image features in a region. These approaches have better performance in images with weak object boundaries due to the independence on image gradients. Chan–Vese (CV) model is a classical region-based active contour model [13]. Utilizing the global information of the image, the CV model can extract the contour that has an unobvious change of gradient in the image. Therefore, the CV model can obtain a good result even if the image has weak object boundaries. Furthermore, the CV model can reserve the topological structure in the evolutionary process automatically. However, image intensity inhomogeneity often occurs in CMR images, and the CV model cannot accurately detect the boundary of the object in such images.

To segment images with image intensity heterogeneity, Li et al. proposed the local binary fitting (LBF) model that is able to utilize image information in local regions [21]. A kernel function was introduced to define a local binary fitting energy in a variation formulation, and the local image intensity information can be embedded into the region-based active contour model. The LBF model has better performance than the CV model in segmentation accuracy. Nevertheless, the LBF model is susceptible to the initial contour placement and configuration of controlling parameters due to the limitation of the localized energy.

Shawn Lonkton and Allen Tannenbaum presented a novel framework based on the localizing region-based active contours to segment objects with heterogeneous feature profiles [22]. This method takes the local image information into account, and has significant improvement in accuracy for segmenting heterogeneous images. The limitations of this model are its sensitivity to the initialization parameters and the window size of the local region. Wang et al. proposed a new local Chan-Vese model for image segmentation [23], which is based on the techniques of curve evolution, local statistical function and level set method. The limitation of this model is its sensitivity to the window size of the local region.

The black-blood T2* has been shown to be effective by yielding high contrast images, providing superior myocardial border definition, and largely reducing blood signal contamination from the myocardium [24]. In recent years, a number of left ventricular (LV) segmentation methods from CMR cine images have been proposed [25,26,27]. However, there are limited reports on LV segmentation from the black-blood double inversion recovery images. We present a robust adaptive Gaussian regularizing CV model using the entropy and local neighborhood information for automatic LV segmentation from cardiac MR images, namely SMLV (LV segmentation method) for brevity. First, the circular Hough transformation (CHT) is used to locate the epicardial and endocardial contours of LV as the initialization. Second, an adaptive window is utilized to reduce the sensitivity to initialization rather than a fixed window using a fixed neighborhood window size in the previous works [28,29]. The interior and exterior energies are weighted by the entropies of interior and exterior regions, so that the homogeneity proportions of the interior and exterior regions are adjusted adaptively, reducing the effect of the optimal configuration of controlling parameters. Third, local region energies were computed over neighborhoods of points close to the curve. To improve the local neighborhood information (LNI), an adaptive window is used to reduce the sensitivity to initialization. Fourth, the Gaussian regularization is used to not only keep the level set function smooth and stability, but also remove the traditional Euclidean length term and re-initialization.

Qualitative and quantitative evaluations of the proposed method were carried out on both synthetic images and real medical images of clinical patients in terms of accuracy and robustness. The results show that the present method can segment the LV from cardiac images accurately and efficiently.

The remaining of the paper is organized as follows: “Material and methods” describes the image data, CV model, the classic CHT and the present method in details. The experiments setup and related evaluation are given in “Results and Discussion”. Finally, some concluding remarks are included in “Conclusions”.

Materials and Methods

1. Image Data

All CMR were performed at a 1.5T scanner (Siemens Sonata, Erlangen, Germany) with a four-element cardiac phased-array coil and standard electrocardiogram gating and were collected from 2006 to 2009. The final data set comprised of 111 × 8 images of 111 patients (57 male /54 female). The patients’ ages ranged from 11 to 51 years old. A single short-axis midventricular slice positioned halfway between the base and the apex of the LV was obtained by using a black-blood multi-echo gradient-echo sequence (flip angle 20°, sampling bandwidth 810 Hz/pixel, voxel size 1.56 × 1.56 × 10 mm3 and variable field of view (20∼30) × 40cm2 depending on patient size). The short-axis images were acquired at 8 echo times (TE) from 2.54ms to 17.90ms with an increment of approximately 2ms in a single breath-hold. The repetition time between radiofrequency pulses was 20ms. Black-blood sequence images were acquired in diastole after a double inversion recovery preparation pulse. For quantitative analysis of the algorithm, we compared the automatic segmentation results to the Gold Standard, which is defined by the medical experts using the computer-aided software. The study was approved by the local research ethics committee and the subjects gave their written informed consent.

2. Methodology

2.1 Circular Hough Transformation (CHT)

Proposed by Duda et al. [30], the CHT is one of the modified versions of the Hough transformation. The CHT aims to locate circular patterns in images. It is used to transform a set of feature points in the image space into a set of accumulated votes in a parameter space. Then, for each feature point, votes are accumulated in an accumulator array for all parameter combinations. The array elements that contain the highest number of votes indicate the presence of the shape. The radius and the center of the circle can also be computed by the CHT. The example of CHT is shown in . In this study, the detected circles are utilized as the initialization in order to reduce the time of the curve evolved to the object and to improve the robustness.

A CHT example.

(a) Original image. (b) The space of CHT. (c) Detected circles. (d) The circles imposed on the image.

2.2 Chan–Vese model (CV model)

Let I denotes a given image defined on domain Ω and let C be a closed curve dividing the image into object and background. The object denoted by Ω is the interior of the curve, and the background denoted by Ω is the exterior of the curve. The energy functional is given by Where μ is the coefficient of the length L(C) of curve C, ν denotes the coefficient of the area S(C) of Ω, and z is the spatial variable; c and c stand for the means inside and outside C, respectively; λ and λ represent the weights of the interior and exterior of the contour C; The boundary of the object is obviously the minimum of ε(C,c ,c ), which can be solved using the level set method [31].

In the level set method, C is considered as the zero level set of a signed distance function ϕ, i.e., C = {z|ϕ(z) = 0}, ϕ(Ω) > 0 and ϕ(Ω) < 0. The energy function expressed by level set function ϕ is where H(z) denotes the Heaviside function and δ(z) is the Dirac function.

F(ϕ(z)) is the fitting term of the energy functional as follows

Minimizing the above energy function by using the steepest descent method, we obtain the following variational formulation Where div(∙)is the divergence.

From the energy functional of CV model, the evolution of the curve is influenced by three terms. The first two terms regularizes the contour to keep the contour smooth during the evolution. The last term in equation (2) is called the data term, which has influence on the evolution of the contour. This model does not use the edge factor to stop the evolving curve on the boundary of the object, and is based on the assumption that the interior and exterior of the curve simultaneously are statistically homogeneous. The advantage of this model is robust to noise. As is well known to all, it is very suitable for segment images with two regions that have a distinct mean of pixel intensity. The performance of the segmentation is subject to optimal configuration of controlling parameters. Moreover, this model often fails to segment the objects with image intensity inhomogeneity and complex background [21]. The above deficiencies motive us to present a new method.

2.3 Entropy

Entropy is a measure of disorder, or more precisely unpredictability, which is an important concept of information theory. The entropy and its several modified versions [32,33] have been recognized as robust techniques for image segmentation. The more homogeneous of the region is, the smaller the entropy of the region. The definition of the entropy is where p(x) represents the probability density functions of the variable the x (e.g., the feature such as brightness, color, gradient). It is convenient to estimate the entropy of the interior region within the given image, E (C) is defined by

Similarly, the entropy of the exterior E is given as

For the level set formulation of the CV model, we replace the unknown variable C by the unknown variable ϕ. Using the Heaviside and Dirac functions, we express the equations (6) and (7) in the following way, respectively:

2.4 Entropy weighted CV model (ECV)

In the CV model, μ and ν are weight that represent the smoothness and tautness degrees of the curve, respectively. λ and λ are the weight of the homogeneities for the interior and exterior of the curve, respectively. The weight parameter is difficult to set for segmentation of heterogeneous images. CV model does not provide the strategy that can choose these parameters reasonably. The requirement of weight parameter is loose in the segmentation of images that have desired piecewise constant structure. However, the medical imaging and the property of organs will lead to heterogeneity. Moreover, the adjoining tissues have little gray level differences. Therefore, the traditional CV model cannot achieve accurate segmentation and the controlling parameters are difficult to set in medical images.

Consider that when E (C) < E (C), the interior homogeneity is greater than that of the exterior. Thus, enhancing the weight of the exterior homogeneity in the energy function is necessary to reduce the influence of lesser homogeneity to energy. When E (C) = E (C), enhancing the weight of the interior homogeneity is necessary. The energy function will be minimal when the homogeneity of the two regions is balanced.

To address the difficulties mentioned above, this paper utilizes entropy as a weight. The fitting energy is given by In the equation (10), when E (ϕ) < E (ϕ), the weight of the exterior homogeneity is enhanced. When E (ϕ) = E (ϕ), the energy function will acquire the minimum. As entropy changes along with the evolving curve, the weight of the homogeneity over two regions in the energy function will adjust automatically and adaptively. These basic remarks are illustrated in , the mauve arrow point to the dominant direction of the weight of interior and the exterior region.

The dominant direction of the weights of the interior and exterior.

(a) E (ϕ) > E (ϕ). (b) E (ϕ) < E (ϕ). (c) E (ϕ) = E (ϕ).

2.5 Adaptive window of local energy of points along the curve

In the CV model, the image intensities are assumed to be statistically homogeneous in the foreground and background. However, the assumption does not hold for some general images. Due to image intensity inhomogeneity and complex construction in medical images, the CV and ECV models often fail to achieve accurate segmentation results. illustrates the considerable challenge in medical image segmentation. It is obvious that the image intensity distributions of both the foreground and the background vary sharply; so it is difficult for the CV to find appropriate constants to fit the foreground and the background. Meanwhile, the foreground and background share the similar intensity to some degree. The heterogeneity of regions and the presence of nearby structures of similar image intensity affect the curve evolution, leading to the detection of false boundaries. CV and ECV models are unable to segment the entire object accurately due to the difficulty to find two optimal constants used as the means of the interior and exterior. These models do not include any local neighborhood information which is crucial for segmentation of medical images with image intensity inhomogeneity.

The considerable challenge in medical image segmentation.

(a) The background and the foreground in the cardiac MR image. (b) The histograms of the background and the foreground.

To solve the above problems, the information of the local regions inside and outside the curve in neighborhoods of points near the evolving curve [34] can be utilized. The object can be described in terms of smaller homogeneous local regions. In each of this small local region, the image intensities can be assumed to be statistically homogeneous. Let C be a closed contour in the image domain Ω, which separates each small region into two regions. As shown in , the irregular contour is the evolving curve and the square is the neighborhood of a point near the curve. Arrows point to the local regions inside and outside the contour in the neighborhood. The exterior and interior of the local region can be represented with local statistics.

Segmentation result is sensitive to the size of the local neighborhood. To avoid the effect of different localization sizes, we use an adaptive window which is estimated by the local entropy of the neighborhood. The evolving curve cannot move to minimize the energy computed in the heterogeneous region. If the entropy of the region is small, this region is homogeneous and we will increase the size of the window until the entropy of the region is more than the threshold. As shown in , the blue rectangle is the initial neighborhood window and the bigger green rectangle is the adaptive neighborhood window, which is determined by the entropy of the window. Instead of selecting a fixed size for all the neighborhood windows, we calculate a local window size for each point on the curve adaptively.

The adaptive neighborhood window.

The blue rectangle is the initial and the bigger green rectangle is the adaptive neighborhood window.

The mean of the local interior u and local exterior v , respectively, can be written as where Ω denotes the neighborhood domain of j-th point in an image I(x), and x is the spatial variable. The fitting energy of the neighborhood centered on the point j is shown as follow:

Given a center point j, the fitting energy F (ϕ) can be minimized when the contour C is exactly on the object boundary and the fitting values u and v optimally approximate the local image intensities on the two sides of C. To obtain the entire object and optimize the fitting energy of the image, each point is considered separately and moved to minimize the energy computed in its neighborhood. The fitting energy of the image is the sum of a family of energy computed in the neighborhood of points near the evolving curve and is accomplished by multiplication with the Dirac function δ(ϕ), as shown in the following equation

2.6 Gaussian Regularizing and Partial differential equation

As pointed out by Shi et al. [35], the evolution of a function according to its Laplacian is equivalent to Gaussian filtering the initial condition of the function. So the previous iteration result of the level set function can be viewed as the initial condition for the next iteration:

In order to enhance the smoothing capacity, we choose the ξ larger than the square root of the time-step Δt. As we use a Gaussian kernel to regularize the level set function the traditional regularized term div(∇ϕ/|∇ϕ|)δ(ϕ) can be re-moved [36].

The final energy function of the adaptive Gaussian regularizing CV model using entropy and LNI is: and the corresponding curve evolving equation is: where .

Results and Discussions

In this section, we evaluate and compare the proposed model (SMLV) with the CV model, the LBF model [21], and recently developed local neighborhood region-based CV models (NCV) [22] using both synthetic and in vivo medical images. The results of these experiments demonstrate the desirable properties of the new model. All the experiments are conducted in a computer with an Intel Core Duo processor, 3.10 GHz, and 4 GB RAM. In the NCV model, we choose a fixed and predetermined neighborhood size n × n of an image size M × N, where the optimal size is calculated as n = (M + N)/16 according to our previous studies.

1. Analyzing the entropy-weighted CV

The first experiment is designed to evaluate the entropy-weighted CV model, which is used to segment the synthetic image, a typical example of image with high contrast (). To show the effect of the entropy, we set μ = 0.01 × 2552, ν = 0, λ = 1 and λ = 1 in CV model, and μ = 0.01 × 2552, ν = 0 in entropy-weighted CV model. From the segmentation process shown in , we can see that the entropy-weighted CV model can segment the object with less iteration compared with CV model. To further study the behavior of entropy, the shows the entropy can adaptively adjust the proportion of the interior and exterior of the curve

The segmentations of CV model (Row 1) and entropy-weighted model (Row 2).

Column 1: the initializations. Column 2, 3, 5 and 6 shows the evolution process for iterations 50, 100, 150, 200, and 250 respectively. Column 7: the final segmentations of CV (iterations 320) and entropy-weighted model (iterations 270).

The weights of the interior and exterior.

(a) CV model. (b) The entropy-weighted CV model.

2. Segmentation of LV contours(the endocardium and epicardium of the left ventricle)

In this section, we validated the proposed model using a synthetic image and mass CMR images for the LV endocardium and epicardium segmentation. The myocardium and background area are inhomogeneous. LVs and adjoining organs overlap in some images presented. In addition, some areas in LV regions have similar gray level to surrounding tissues. These results partly explain why the CV model results in poor LV segmentation.

shows the results of automatic segmentation using the four models (LBF, CV, NCV and the present method), compared with the reference manual segmentation (as Ground Truth) by medical experts. The first column shows three cardiac MR images and their initializations, which are detected through the CHT, from the second to the fifth columns show results for the automatic segmentation using the four models, and the last column presents the manual segmentation. From , the boundaries of the LVs were automatically detected accurately by the present model proposed in this paper, which exhibited better results than the other two models. The NCV model got some wrong boundaries in the walls of the LV areas, and the high-intensity regions were considered as objects. On the other hand, the LBF and CV model may detect some wrong contours to some extent. Hence, shows the various advantages of the present model in medical image segmentation. Table 1 lists the controlling parameters of the four methods. It is obvious that the present method not only avoids to the optimal configuration of controlling parameters due to the action of entropy, but also takes shortest time to achieve accurate segmentations.

Table 1

The controlling parameters of the CV, LBF, NCV models and the present method.

Method	Patient	λ in	λ out	μ	ν
	1	1	1	0.5 × 2552	0.2
CV	2	1	1	0.5 × 2552	0.2
	3	1	1	0.5 × 2552	0.2
	1	1	1	0.03 × 2552	0.2
LBF	2	1	1	0.03 × 2552	0.2
	3	1	1	0.05 × 2552	0.2
	1	1	1	0.1	0.2
NCV	2	1	1	0.2	0.1
	3	1	1	0.2	0.15
	1	—	—	—	0.05
SMLV	2	—	—	—	0.05
	3	—	—	—	0.05

Comparing segmentations.

(a) Initial contour set by the CHT. (b) LBF. (c) CV. (d) NCV. (e) SMLV. (f) Ground truth.

For quantitative evaluation among the four automatic segmentation methods, we utilize the DICE [37], ME [38] and Hausdorff distance (HD) [39] to test the various algorithms, as follows where F and B are the foreground and background of the manual algorithm respectively, as the F and B are foreground and background of the automatic algorithm. The |F | denotes the number of pixels assigned to the foreground by the medical expert. The higher the value of DICE, the better the overall performance of the segmentation will be. DICE Index equal to 1 suggests a good match between manual and automatic segmentation. The lower the ME ratio, the fewer the misclassified pixels are. HD measures the maximum distance between two contours. illustrates the DICE, HD and ME values of the SMLV method. The DICE values are ranged from 0.9203 to 0.9607, the ME 0.0019 to 0.0052 and the HD from 1.5000 to 1.6623. To compare the automatic segmentation of the LV, we show in Fig. 10 the Dice values, scaled to a range of the same size (0.35). A quick look at the Dice similarity measures, ME and HD can be done also using a boxplot, see . Note that there is substantially less variability in the ratings of our method than others. According to , it can be observed that the myocardium extraction by the SMLV is most accurate with least misclassified pixels. The segmentation results demonstrate that the present approach outperforms the other traditional methods and has good agreement with the manual method.

Fig 10

Boxplot representations of Dice, ME and HD values for the 111 patients segmented by various methods.

3. Computation time

shows the comparing the time efficiency of the Gaussian regularizing with the re-initialization level set function. The Gaussian regularizing removes re-initialization and the traditional regularized term, so is capable of keeping the curve evolution stable and smooth while reducing the computation time of the algorithm.

In , we show more examples of accurate segmentations by the present method.

The other patients’ segmentations by the present method (10 iterations).

Row 1: the initial contours set by the CHT. Row 2: the segmentation endocardium and epicardium of the left ventricle.

Conclusions

We have presented an automatically method for myocardium segmentation from black-blood cardiac images. By utilizing the CHT, the present method can detect the LV endocardium and epicardium automatically. Traditional CV model is unsuitable for medical segmentation as it has limitation in parameter setting and suffers from heterogeneity intensities. To address the first issue, we utilize entropy to weight the energies of inside and outside the contour. For the second one, we incorporate the local neighborhood information to reduce the effect of inhomogeneity inside regions. Instead of selecting a fixed neighborhood window in NCV, we calculated the local window size for each point on the curve adaptively to reduce the model sensitivity to initialization. In addition, we use the Gaussian to regularize the level set function for removing traditional regularized term and the re-initialization, which not only ensure the smoothness and stability of the level set function, but also reduce the computation time. Encouraging experimental results on the real medical images demonstrated that the proposed algorithm is very robust, efficient and much less sensitive to the initial contour.

The limitation of this method is that the CHT may be failed in some images with the bright region of the fat and muscle. In the future, we will try to preproccess these images for improving the robust of the initialization by the CHT. MRI-T2* has been accepted as a clinical tool for monitoring iron overload in thalassemia patients. Our future studies will focus on the assessing the myocardial iron loading in CMR images. Due to the universality of bright-blood images, the present method will be improved in order to segment these bright-blood images for T2* assessment.