Contrast-enhanced Bias-corrected Distance-regularized Level Set Method Applied to Hippocampus Segmentation

Abstract

Recently, the level set has become a popular method in many research fields. The main reason is that it can be modified into many variants. One such case is our proposed method. We describe a contrast-enhancement method to segment the hippocampal region from the background. However, the hippocampus region has quite similar intensities to the neighboring pixel intensities. In addition, to handle the inhomogeneous intensities of the hippocampus, we used a bias correction before hippocampal segmentation. Thus, we developed a contrast-enhanced bias-corrected distance-regularized level set (CBDLS) to segment the hippocampus in magnetic resonance imaging (MRI). It shows better performance than the distance-regularized level set evolution (DLS) and bias-corrected distance-regularized level set (BDLS) methods in 33 MRI images of one normal patient. Segmentation after contrast enhancement and bias correction can be done more accurately than segmentation while not using a bias-correction method and without contrast enhancement.

1. INTRODUCTION

Assessing the hippocampus is an important step for the early diagnosis of mild cognitive impairment (MCI), Alzheimer’s disease (AD), and other brain disorders. The best way is to detect such a disorder early, follow it, and if still possible, prevent disease progression. Hopefully, this will reduce the number of patients with brain disorders. In this respect, the accurate, reliable segmentation of the hippocampus could play a major role in overcoming such diseases, because an early symptom of these diseases is the deformation of the shape of the hippocampus, which can be detected.

The hippocampus belongs to the limbic system and is a major component of the brain in humans. It is the most electrically excitable part of the brain, where new neurons continue to be created throughout human life. Consolidation of information from short-term memory to long-term memory and spatial navigation occur in the hippocampus. Early symptoms of AD include damage to the hippocampus, memory loss, and disorientation. Difficulties in forming new memories are a result of severe damage to the hippocampus in both hemispheres. Moreover, if both the hippocampus and the parahippocampus are damaged, complete amnesia can result [1]. The shape of the hippocampus depends on age but can also be changed by traumatic stress. Not all elderly people have a smaller hippocampus and bad memories but those who do tend to perform less well on memory tasks and some aerobic tasks. In addition, younger people show more activation in the hippocampus than older people.

 

2. MATERIALS AND METHODS

2.1 Image Acquisition

In this work, the images used (256 × 256) were captured with a 1.5-T General Electric (GE) Medical System MRI device at Haeundae Paik Hospital, Korea. From this, we obtained 33 slides from one normal patient for hippocampal segmentation. The segmentation program used Matlab. The binary ground truth was obtained with manual segmentation.

2.2 Related Works

The first level set method was proposed by Osher and Sethian in 1987 [2]. They used it to capture dynamic interfaces and shapes by a contour as the zero level set, called a level set function (LSF). The LSF moves by evolution until it reaches its expected boundary. In the active contour model, there have been many applications to image segmentation [3-5]. These studies formulated the dynamic parametric contour C(s, t): [0, 1] × [0, ∞) → ℜ2, which has a curve evolution as follows:

with a spatial parameter s in [0, 1] that parameterizes points in the contour. F is the speed function and N is the inward normal vector to the curve C. N can be formulated as

The curve evolution in (1) can be converted into a level set by adding the dynamic contour C(s,t) as the zero level set of a time-dependent ϕ(x,y,t). This ϕ(x,y,t), called the LSF, has negative values inside the zero level contour and the opposite, positive values at the outside [6,7]. The formula in (2) is a conversion from the curve evolution to a partial differential equation (PDE):

It is also called an implicit active contour or a geometric active contour (GAC) model, which is given as an LSF.

The LSF also has its drawbacks, such as irregularities during evolution, and it introduces numerical error and can destabilize level set evolution [8,9]. Hence, reinitialization is introduced by performing a termination of the evolution periodically and reshaping the degraded LSF as a signed distance function [10-12].

One adverse effect of reinitialization is that the LSF can incorrectly move the zero level set away from the expected position; thus, reinitialization should be avoided as much as possible [13]. It is usually applied in an ad hoc manner because there is currently no consensus on when or how best to apply it [14]. Reinitialization can be expressed as follows:

The steady-state solution of this equation is a signed distance function, where sign (ϕ) is the sign function.

Gomes and Faugeras [14] used three PDEs to avoid reinitialization. The first PDE’s role is to restrict the LSF to a signed distance function, and the second and third PDEs describe the motion of the zero-level contour. However, they still had issues with destabilization of the level set evolution and destruction of the signed distance property, and had to introduce a separate reinitialization phase. Another study proposed a method for avoiding separate reinitialization of the GAC [15].

We chose the level set as the method for this segmentation for several reasons. First, numerical computation is performed on a fixed Cartesian grid without point parameterization on a contour, as in the parametric active contour model (ACM) [16-19]. As a result, level set methods have been applied widely in many scientific and engineering areas [20-24]. Second, the level set can handle topological changes that represent contours of a complex topology that cannot be supported by parametric ACM, and can handle shape corners and cusps in a propagating solution as well as three-dimensional effects [25-26]. Third, level set methods are widely known as versatile, robust, accurate, and efficient techniques for solving many problems. A previous study proposed a variant of the LSF to maintain signed distance properties that did not need reinitialization; however, the method still had unwanted effects on the numerical accuracy under several conditions [27].

2.3 Intensity Inhomogeneity and Bias Correction

Intensity inhomogeneity is a smooth change in intensity inside originally homogeneous regions. It should be removed or corrected because it will degrade the segmentation ability of an algorithm. Almost all brain images captured using MRI devices have intensity inhomogeneities because of several factors. There are now many intensity inhomogeneity correction methods [28].

The histogram-based method is a high-frequency method that does not use any knowledge about an image or uses only information present in an image, without making assumptions on spatial and intensity distribution; it is fully automatic and general in that it needs no initialization or prior information. In this work, we use intensity inhomogeneity as an attribute of a component of an image, formulated as follows:

where I is an observed image, J is the true image, b is the intensity inhomogeneity or bias field that has a low variance, and n is additive noise (i.e., zero mean Gaussian noise) [29].

The bias field b can be approximated by a constant in the neighborhood of each point in the image continuous domain I:Ω→R. The true image J takes N distinct constant values C1, …, CN in disjoint regions Ω1, …, ΩN, where {Ωi}Ni=1, { Ci}Ni=1 and bias field b. It is assumed that b slowly changes.

2.4 Contrast Enhancement

There are several ways to enhance the contrast in MR images, such as histogram equalization (HE), brightness biHE (BBHE), and contrast-limited adaptive HE (CAHE). We used CAHE in preprocessing to readily segment the hippocampus region.

CAHE prevents the tendency to over-amplify noise in homogeneous regions, as adaptive HE does, by limiting its contrast enhancement. Before computing the cumulative distribution function (CDF), CAHE clips the histogram at a predefined value to limit its amplification as well as the slope of the CDF and function of the transformation. The sizes of neighborhood regions and histogram normalization influence the clip limit. The resulting amplification limit is usually between 3 and 4. The excess over the clip limit is redistributed among all histogram bins. Some bins are redistributed over the clip limit repeatedly, recursively, to an effective clip limit bigger than the prescribed limit, but the exact value affected by the image and also the excess are negligible. The biggest disadvantage of CAHE is its high computational cost, because of the need to compute different neighborhood histograms for each pixel in an image. It uses interpolation to improve this, to minimize the cost without degrading result quality.

A modeled digital image histogram has intensity levels [0, L-1] and can be formulated as follows:

where rk is the kth intensity value and nk is the number of pixels that have intensity rk, and has a discrete function.

The normalized histogram can be formulated as:

where Pr(rk) is an estimated occurrence probability of intensity level rk. The new distribution of HE can be formulated as:

The first step is to separate the region into several non-overlapping, almost equal, regions and to calculate the histogram of each region. Then a clip limit for clipping the histograms is obtained from the limit for contrast enhancement to redistribute the histograms in such a way that each histogram height does not exceed the clip limit. The clip limit (β) can be calculated as:

where a is the clip factor, and the resulting contrast- limited histogram CDFs are obtained for mapping the grayscale. Furthermore, the results from the mapping of the four nearest regions can be combined linearly, as explained in a previous study [30].

2.5 Contrast-Enhanced Bias-Corrected Distance-Regularized Level Set

Chunming et al. [6] studied DLS and proposed a variant with a distance-regularized term and an external energy term that drives the motion of the curve towards the expected position. They also proposed a double-well potential for the distance regularization term, a potential function that forces the gradient magnitude of the LSF to its minimum points by maintaining the shape of the function’s curve. Then level set evolution not only maintains its regularity by the double-well potential as a penalty term but also its forward-and-backward (FAB) diffusion. This FAB diffusion is derived from the distance regularization term. In this way, reinitialization can be avoided. Next, the DLS is applied in a narrowband domain, using an edge-based active contour model to perform the image segmentation.

The difference between our proposed method and the existing DLS is that we use a Heaviside step function, which is 1 when Φ ≥ 0 and 0 when Φ < 0.

The energy formulation with distance regularization can be written as

where ϕ: Ω→ℜ is an LSF in a domain Ω, Rp(ϕ) is a level set regularization term, μ is a constant greater than 0, and εext(ϕ) is external energy.

Next, Rp(ϕ) can be formulated as:

Let p be a potential (energy density) function, and p: [0, ∞) → ℜ. When the zero level set LSF ϕ is at the expected position (i.e., the boundary object to be segmented), then the external energy reaches its minimum.

Let |∇ϕ | be zero if p(s) = s2. To maintain the computational accuracy of curve evolution, a level set regularization term should be imposed where p(s) has a minimum point, s = 1 (it may have more than s = 1 as minimum points), so that LSF will be smooth, ℜp(ϕ) is minimized, and the signed distance property ∇ϕ =1 is maintained.

Then potential p can be expressed as

The energy functional P(ϕ) is a penalty term that can maintain the signed distance property in the entire domain, which can be expressed as

To maintain the LSF profile, potential function p(s) should have minimum points at s = 1 and s = 0, the so-called double-well potential. Next, to achieve the minimum energy is to find the steadystate solution of this gradient flow equation, which is an evolution equation of time-dependent function ϕ(x, t) with a spatial variable x in domain Ω and temporal variable t ≥ 0, expressed as:

The evolution starts with an initial function The evolution of the time-dependent function ϕ (x, t) equals -∂F/∂ϕ, which is why it is called a gradient descent flow or steepest descent flow.

From (11), we can get the functional Rp(ϕ) as Gateaux derivatives:

where

Equation (17) is derived from (10) and the gradient flow of energy ε(ϕ), defined as:

Next, with equation (15), we can obtain:

Equation (19) is called the DLS because of its ability to maintain the signed distance property of the distance regularization of the LSF. DLS does not need reinitialization because of the distance regularization. If distance regularization in DLS is seen from the gradient flow of energy μRp(ϕ), then

Equation (20) can be defined in a diffusion equation as:

where D is a diffusion rate, defined as D=μdp(|∇ϕ |) is a diffusion rate whose value can be positive or negative. If it is positive, the diffusion goes forward, so |∇ϕ | will decrease. If it is negative, the diffusion goes backwards, so |∇ϕ | will increase. This diffusion is called a forward and backward mechanism. The |∇ϕ | decreases or increases to make it near its minimum points of the potential function p(s), while preserving the shape of the curve.

Equation (19) can also be expressed as:

The div(∇ϕ /|∇ϕ |) term computes the mean curvature of the LSF.

The sign of dp(s)=1-(1/s) indicates the FAB diffusion term. When |∇ϕ | > 1, while the diffusion rate μdp(|∇ϕ |) is positive, the diffusion equation in (19) is going forward and decreases |∇ϕ |. When |∇ϕ | < 1, while the diffusion rate μdp(|∇ϕ |) is negative, the diffusion equation is going backward and increases |∇ϕ |.

FAB diffusion with potential p = p1 keeps |∇ϕ | equals to 1 to preserve the signed distance property. This FAB diffusion has an unbounded diffusion rate, which goes to infinity as |∇ϕ | approaches 0. This disadvantage can be handled by a double-well potential, p = p2, where the μdp(|∇ϕ |) is bounded by a constant.

Several works have already proposed improvements to the DLS model. In one such study, an isotropic nonlinear diffusion filter was proposed by adopting a regularized P-M filter to remove image noise while preserving an edge information equation filter. It used the image gradient modulus while filtering the noise, and at the same time, also preserved the edge detection. It used a Gaussian filter for smoothing, then applied search filtering of the image’s gradient modulus and used a regularized P-M equation to obtain smoothed image results of isotropic nonlinear filtering. They also proposed a new adaptive function u’(x,y):

If the target object for segmentation is inside the evolution curve, Dint > Dext, where Dint is the average of the image’s gradient modulus of the inward region of the curve and Dext is the average image’s gradient modulus or the outward region of the curve, and the result is u’(x,y) = u(x,y). This makes the normal vector take a negative value and the evolution moves inside the curve. It also applies in the opposite direction if Dint < Dext, where u’(x,y) = -u(x,y) and the normal vector takes a positive value and the evolution moves outside the curve towards the expected position.

More DLS enhancement was performed by W. Weifeng, W. Yuan and H. Qian [31]. They improved the penalty term of DLS, which still had an unwanted effect, using three potential wells as a diffusion rate for the penalty term and also different narrowband steps to overcome it.

The procedure to segment the hippocampus is described in Fig. 1.

Fig. 1.Pipeline of the CBDLS procedure.

 

3. RESULTS AND DISCUSSION

3.1 Results of Segmentation and Visualization

Fig. 2 shows a brain MR images used in our experiments; the slice name is SUB40053.bmp. For simplicity, the MR images are cropped at the region of the hippocampus. After the bias field is estimated, it is corrected and its contrast is enhanced. The image slice is segmented and the segmentation final contour result is used to make a binary image of each slice.

Fig. 2.(a) Original MR image and (b)Cropped region of interest.

The slice is cropped at the hippocampus to ease the segmentation process. Then initial contours are placed on the body of the hippocampus, as shown in Fig. 3. The initial contour involves a step-by-step process of going outside each iteration to the edge of hippocampus area until it reaches a final contour in certain iterations. The final contour from the final iteration is shown in Fig. 3d–f. We can see that a better segmentation result is achieved with CBDLS. BDLS segmented the hippocampus better than DLS, and CBDLS segmented the hippocampus better than BDLS.

Fig. 3.(a) Initial contour of the level set, (b) Initial contour of the bias-corrected level set, (c) Initial contour of the bias-corrected and CAHE level set, (d) Final contour of the level set hippocampus segmentation, (e) Final contour of bias-corrected level set hippocampus segmentation, and (f) Final contour of bias-corrected and CAHE level set hippocampus segmentation.

After obtaining the segmentation result from the final contour of all three methods, we return the cropped images to the original size. Then we make the final contours as binary images for simplicity in Fig. 4.

Fig. 4.(a) Binary ground truth, (b) Binary image of level set hippocampus segmentation, (c) Binary image of bias-corrected level set hippocampus segmentation, and (d) Binary image of bias-corrected and CAHE level set hippocampus segmentation.

3.2 Validation and Discussion

In this experiment, we used the Jaccard similarity index to compare our proposed CBDLS procedure to DLS and BDLS. Jaccard similarity can be formulated as:

It is the size of the intersection divided by the size of the union of finite sample sets. In this experiment, A is a ground truth and B is the segmented image under consideration. Table 1 compares the Jaccard similarity index among the CBDLS, BDLS, and DLS.

Table 1.The Jaccard Similarity Index of DLS, BDLS, and CBDLS in hippocampus segmentation

If A and B are empty, then J(A, B) = 1

The Jaccard distance can be defined as dissimilarity: 1 minus Jaccard similarity or the difference in the sizes of the union and intersection divided by the union of two sets [32-33].

From Table 1 we can see that the sum and the mean of CBDLS are higher than for BDLS and DLS. Higher is better and ideally the value approaches 1. If the Jaccard similarity achieves 1, then the binary segmentation result is equal to the binary ground truth. Thus, we can conclude that CBDLS is more accurate than BDLS, and that BDLS is more accurate than DLS in the hippocampus segmentation of MR images.

From Table 1 and Fig. 5, we can see that our method is the best method for segmenting the hippocampus, with the highest accuracy. The addition of bias correction and contrast enhancement before segmentation using a modified level set can be a considerable advance to improve the image quality for ease of the segmentation process. If the segmentation process is easier, then costs can be lower.

Fig. 5.Comparison of the Jaccard similarity among methods.

The easier the segmentation, the more accurate the segmentation result. A drawback of our proposed method is that preprocessing the MR images until the final result as binary images is more time consuming than other methods. This is because two additional steps are needed before segmentation.

From Table 1 and Fig. 5, we can see that the sum and mean of CBDLS are higher than in BDLS and DLS. There are some segmentation results that are better than CBDLS in hippocampus segmentation (as can be seen in Fig. 5). This is because the image, after contrast enhancement and bias correction, makes the level set contour difficult to move inside the boundary of the hippocampus. The reason is that the intensity of each pixel inside the hippocampus has more contrast and more heterogeneity.

Consequently, there is a need to set the initial contour more in the beginning than with the other methods. Another drawback is that it takes more time and iterations to shape the hippocampus until reaching the ‘correct’ boundary. However, overall, the sum and mean of the Jaccard similarity index for 33 image slides of hippocampus segmentation show the best segmentation results with CBDLS, followed by BDLS and DLS. This is because in our method the contour is difficult to move to the background region, making it difficult to over-segment, and the opposite of the DLS method, which often over-segments the hippocampus area. BDLS is better than DLS because it removes the intensity inhomogeneity that makes segmentation of the hippocampus difficult to achieve. Thus, we can conclude that CBDLS is more accurate than BDLS, and BDLS is more accurate than DLS in segmenting hippocampus MR images.

 

4. CONCLUSION

We propose the use of contrast enhancement, bias correction, and a level set method to handle hippocampus segmentation, which suffers from intensity inhomogeneity and the nature of hippocampus where the edges have quite similar intensities to the neighboring pixels. Our experimental results indicate that our method outperforms both BDLS and DLS. The DLS method has a tendency to over-segment the hippocampus region, and BDLS has a tendency to under-segment it. Our proposed method meets the boundaries of the hippocampus almost exactly. This is because the bias correction makes the intensity inhomogeneity of the hippocampus more homogeneous, and the contrast enhancement provides the edges of the hippocampus region with more contrast versus the background. This helps prevent moving the contour beyond the real boundary of the hippocampus.

Although our method requires more time, it can be used as an alternative to existing methods to segment the hippocampus. In the future, we will seek to reduce the time needed to finish the segmentation of CBDLS, to make it more efficient.