Region-based Active Contour Model based on Markov Random Field to Segment Images with Intensity Non-Uniformity and Noise.

This paper represents a new region-based active contour model that can be used to segment images with intensity non-uniformity and high-level noise. The main idea of our proposed method is to use Gaussian distributions with different means and variances with incorporation of intensity non-uniformity model for image segmentation. In order to integrate the spatial information between neighboring pixels in our proposed method, we use Markov Random Field. Our experiments on synthetic images and cerebral magnetic resonance images show the advantages of the proposed method over state-of-art methods, i.e. local Gaussian distribution fitting.

INTRODUCTION

Image segmentation is an important step in many image analyses. The goal of image segmentation is to partition an image into a number of regions, each having similar characteristic. In the recent years, different methods have been developed for image segmentation. One of the most popular methods is active contour.

Active contours are based on the theory of surface evolution under a speed function determined by local, global and independent properties. They are able to provide smooth and closed contours to extract the boundaries of objects in image. They can be divided into two main groups: geometric and parametric. Parametric active contours are explicitly represented as parameterized contour/surface in a Lagrangian formulation.[1] One of main drawbacks of parametric active contours is the difficulty of naturally handling topological changes for splitting and merging of contours.

Osher and Sethian[2] proposed a level set frame work for geometric active contour to solve this problem. The main idea of the level set framework is to implicitly represent an evolving curve/surface C in the image plane (Ω ∈ ℜn) as zero level of a function φ defined in a higher dimension. According to defined energy functional, they can be broadly categorized into two classes: Edge-based models[34] and region-based models.[5-8]

Edge-based active contour models use image gradient as a constraint to stop the evolving curve on the desired object boundary. However, these models are not appropriate for segmenting images with high-level noise and weak edges.[8] Region-based models use global information such as region intensity homogeneity to evolve the curve. Therefore, they provide performance improvement over edge-based methods.

Active contour based on piecewise constant (PC) model originally proposed by Chan and Vese[5] is a popular region-based active contour model. This model assumes that each image region has statistically homogeneous intensity. In this model, Chan and Vese proposed a level set formulation for the approach proposed by Mumford-Shah to segment an input image domain Ω into two separated homogeneous regions, {Ω}, defined by the evolving contour C that minimizes the following functional:[9]

where I(y): Ω → ℜ is an image defined on Ω and the constants c and c approximate the image intensity inside and outside of the contour C. Obviously, because of intensity non-uniformity in real images, the constants c1 and c2 can be far from original image intensities in each region.

In order to overcome difficulty of segmentation of images with intensity non-uniformity, Li et al.[610] proposed the local binary fitting (LBF) model to use local intensity information for image segmentation. They locally approximate the intensities outside and inside the evolving contour using two fitting functions f1(x) and f2(x). The energy functional of the LBF model is defined as follows:

Where λ1, λ2, ν and ∝ are positive weighting constants and W(u) is a kernel function with localization property. Typically, W(u) is chosen as truncated Gaussian kernel, i.e. W(u) = 0, u > σ.[6]

In the presence of intensity non-uniformity and noise, image segmentation based on local intensity mean cannot provide accurate results. Thus, Wang et al.[78] proposed the local Gaussian distribution fitting (LGDF) model in which the local intensity is modeled by Gaussian distribution. The energy functional for LGDF model for each pixel is defined as follows:

where P(I(y)) is:

In this energy functional, u and σ2 are local mean and variance for intensity of the image, respectively.[811]

However, the region based-models are based on the assumption that pixels in each region are statistically independent of each other. This assumption ignores the correlations between pixels, which holds for most real images. To address this problem, Markov random fields can be used to model the correlations between neighboring pixels. Recently, Xu et al.[12] applied the Markov random field model into the region-based active contour proposed by Chan and Vese.[5] However, in this model, the intensity non-uniformity correction was not addressed in the proposed energy functional. Furthermore, the intensity means and standard deviations were simply calculated globally for each region without using a Gaussian kernel or any localized energy.

In this paper, we propose a novel region-based active contour model in a level set formulation that can be used to segment images with intensity non-uniformity and high-level noise. The main idea of our proposed method is to use Gaussian distributions with different means and variances while incorporating a model of intensity non-uniformity (bias) for simultaneous segmentation and intensity non-uniformity correction. Furthermore, by incorporating Markov random field to integrate the spatial information between neighboring pixels in the active contour model, the proposed generalized region-based model can be used to segment images in the presence of image artifacts.

The rest of the paper is organized as follows. The proposed method is introduced in following Section. Detailed simulation results are provided in Result Section and concluded remarks are given at the end.

PROPOSED METHOD

Let Ω ∈ ℜ2; be image domain and I(y): Ω → ℜ be the observed gray level intensity. Let h: Ω → {Ω1, Ω2} be a function that maps each point y in image domain into one of two regions {Ω1, Ω2} where Ω1∩Ω2 = Ω and ω ∈ Ω implies hω = Ω. The goal of image segmentation is to determine the function h that optimally partitions the image domain. In statistical framework, an optimal partitioning P(H = h) can be obtained by maximizing a posteriori probability P(H = h|I). This conditional probability can be expressed using the Bayes rule as follows:[13]

maxP(H = h|I)∝P(I|H = h)P(H = h)      (5)

where conditional probability P(I|H = h) is the a posteriori segmentation probability.

The P(H = h) is the prior probability about regions to be segmented and can be estimated using the Markov random field to model the image. A random field is Markov random field if all of its local conditional probability density functions satisfy the Markov property as follows: P(hω| h–ω) = P(hω| hη), where h–ω = {h(r): r ∈ Ω, r ≠ ω}, hη = {h(r): r ∈ ηω, r ≠ ω}, and the neighborhood ηω of ω is any bounded region such that s ∈ ηω ↔ ω ∈ ηs, ω ∉ ηω and ω, s ∈ Ω. Thus, the prior probability P(H = h) can be approximated by the normalized pseudo-likelihood as follows:[14]

where

and Z and Zω are normalization constants, β is an appropriately selected constant, dc represents the bin volume that guarantees the correct continuum limit,[11] A(ηω, Ω) is the area of the region defined by ηω∩Ω where ηω is selected to be a rectangular window that defines on initial contour.

To estimate the conditional probability P(I|H = h), we assume that the pixels within each region are independent and identically distributed realizations of the same random process. Thus, P(I|h) can be written as follows:

A common choice for conditional probability for a given value I with respect to h is to use Gaussian distributions with mean ∝ and variance σ2 as follows:

One of the main problems in image segmentation, e.g. medical images, is intensity non-uniformity or bias field that often arises from the imperfection of imaging devices. The bias field can be modeled as a smooth and spatially varying field b(x) multiplied by the constant true signal of the same tissue c. Thus, image intensity corresponding to the object domain Ωi (equation 9) can be approximated by Gaussian distribution with local mean bc and variance Ω2 as follows:[1516]

Maximization of the a posteriori probability is equivalent to minimization of its negative logarithm.

–log(P(H = h|I)) = –{logP(I|H = h) + log(P(H = h))}      (11)

After integrating over the entire image domains, the following energy functional is obtained:

where W(x – y) is a truncated Gaussian kernel with Ω > 0. With equations (7) and (10), energy functional in (12) can be re-written as follows:

Note that proposed energy functional has two differences from the functionals proposed by Wang et al.[78] and Xu et al.[12] First, in the proposed functional, a Gaussian kernel and a model for intensity non-uniformity are used to localize energy functional and simultaneously perform segmentation and bias correction. Second, integrating Markov random filed in the proposed energy functional improves the image segmentation in the presence of artifices. Thus, the proposed functional is of general form in comparison with other methods.

Level Set Formulation

The goal is that the image domain Ω can be partitioned into two regions {Ω1, Ω2} background and foreground, respectively, by using functional E(Ω1, Ω2, C). These two regions can be represented as the regions outside and inside of the zero level set φ. Thus, using M(φ(y)) as the characteristic function of region Ω, the energy functional[14] can be expressed using level set function φ(y) as follows:

where β, ν and μ are positive weighting constants. Using the Heaviside function H, the characteristic function of regions M(φ(y)) are defined as follows:

The Heaviside function and the Delta function, the derivation of Heaviside function, can be approximated as follows:

where ε is a constant.

For fix φ, by calculus of variations with gradient descent flow, it can be shown that the parameters c, σ2, b that minimize the energy functional in equation (14) can be obtained as follows:

Minimization of energy functiona E(.) with respect to φ can be achieved by solving the gradient decent flow equation as follows:

where

and * indicates the convolution operation.

The pseudocode for image segmentation is given in Table 1.

Table 1

The pseudocode for image segmentation

RESULTS

To validate the proposed model, an extensive validation was performed based on two imaging data types: (i) synthetic images and (ii) simulated and real medical images. In this section, the following parameters σ = 3, μ = 1, ε = 1, ν = 0.00001×255×255, β = 0.1, Δt = 0.01 were set.

Synthetic Images

Figure 1 shows the application of the proposed method for segmentation of synthetic images with noise and intensity non-uniformity. The segmentation results using the proposed method achieves satisfactory result.

Figure 1

Application of the proposed method for segmentation of synthetic image with noise and intensity non-uniformity. (a) Original image and initial contours; (b) curve evolution after 50 iterations; and (c) final zero-level contour using the proposed method

Medical Images

The proposed method is applied to simulated and real medical images with intensity non-uniformity and noise. The result of the vessel image segmentation is shown in Figure 2. As can be seen, the image has weak boundaries and is also corrupted by intensity non-uniformity. However, the method can segment the desired tissue appropriately.

Figure 2

Segmentation of vessel image using the proposed model. (a) Original image with initial contour; (b) curve evolution after 20 iterations; and (c) final zero-level contour

In order to study the robustness of the proposed method for brain tissue segmentation from magnetic resonance (MR) images with different noise (n) and intensity non-uniformity (RF), the white matter (WM) is segmented from simulated T1-weighted MR images with different noise and intensity non-uniformity (bias) level. They have a matrix size of 181×272×181 and voxel dimensions of 1 mm×1 mm×1 mm. Figure 3 shows the results of WM segmentation from simulated brain MR images using the Brainweb MR image simulator. The results obtained using the proposed method were compared with those obtained by the LGDF method. The first column in Figure 3 shows the selected axial slices with n=5%, RF=40; n=7%, RF=40 and n=9%, RF=40%, while the second and third columns illustrate the segmented WM using LGFD and the proposed method, respectively. As can be seen, the proposed method can properly extract WM in the presence of high-level intensity non-uniformity and noise in comparison with the LGDF method.

Figure 3

Comparison between the proposed method and local Gaussian distribution fitting (LGDF) method brain magnetic resonance images obtained from Brainweb. From top to down, the images are corrupted with noise = 5%, RF = 40; noise = 7%, RF = 40 and noise = 9%, RF = 40%. (a) Original image with initial contour; (b) final zero-level contour of LGDF method; (c) final zero-level contour of the proposed method

Furthermore, we applied our proposed model to segment WM from T1-weighted brain images of 20 normal subjects provided by the Center for Morphometric Analysis at Massachusetts General Hospital¥ (IBSR). IBSR provides T1-weighted coronal MR images where each slice has 256×256 pixels and the voxel dimensions are 1.17 mm× 1.17 mm×3.1 mm. The real MR brain data and their hand-guided expert segmentation results are available at these datasets. This database consists of clinical datasets that are a result of real MRI scanning. The images involve various problems related to segmentation of real MRI data, such as intensity non-uniformity and noise.

In order to quantitatively evaluate the proposed method, the Dice similarity metric is used,[417] which is defined as follows:

where L and L represent the obtained segmentation and ground truth, respectively, and n (.) indicates the number of voxels. This metric represents spatial overlap between two binary images and its value range between 0 (no overlap) and 1 (perfect agreement). On the other hand, to assess the success and error rate of the proposed method, the sensitivity (true positive fraction) and specificity (true negative fraction) were calculated. Sensitivity refers to the ability to correctly identify appropriate tissue in the segmented mask and specificity refers to the ability of the proposed segmentation method to correctly remove non-desired voxels. Sensitivity and specificity are the measures computed based on true positive (TP), false negative (FN) and false positive (FP).

Figure 4 illustrates the Dice similarity values (DSM) of the proposed method for WM segmentation from simulated MR images using the Brainweb MR image simulator in comparison with the LGDF method. It can be seen that the similarity values of our method are generally higher than those of LGDF, especially in high-level noise and intensity non-uniformity. In addition, the average results for WM segmentation from the eight simulated and 20 real datasets using the proposed method and LGDF are shown in Table 2. As can be seen, the method can segment the WM with mean similarity (Dice) of 99% and 98% in case of simulated and real images, respectively. If a segmentation result can be classed as high quality with an arbitrarily high threshold on the Dice coefficient, such as DSM >90%, then the proposed method demonstrates high performance for segmented WM from simulated and real IBSR datasets.

Figure 4

Dice similarity values for the segmented white matter from the simulated magnetic resonance images using the Brainweb simulator using the proposed method and local Gaussian distribution fitting method

Table 2

Quantitative evaluation of the proposed method for WM segmentation from simulated and real dataset IBSR

CONCLUSION

In this paper, we proposed a new method that can be used to segment images with intensity non-uniformity and noise. The segmentation framework is based on region-based active contour that uses Gaussian distributions with different means and variances. It incorporates a model of intensity non-uniformity to segment the image and correct the intensity non-uniformity. Furthermore, we have used the Markov random field to integrate the spatial information between neighboring pixels in our proposed method.

Recently, different region-based active contours have been proposed that use the intensity non-uniformity model and Markov random field separately to overcome bias and noise. However, in this research, we proposed a generalized region-based active contour model that integrates Markov random field and intensity non-uniformity model into energy functional. Furthermore, in the proposed model, we utilized local means and standard deviations by using Gaussian kernel. Thus, while the proposed model allows tissue segmentation and intensity non-uniformity correction to be performed simultaneously, it handles noise efficiently in real images such as cerebral MR images. Moreover, the proposed model is robust against initializations, allowing fully automatic applications.

We applied our method to the synthetic and cerebral MR images. Assessment of synthetic image segmentation was performed qualitatively and the brain tissue segmentation was performed quantitatively. As this method can be used to segment medical images, it is applied to segment brain tissues from simulated and real MR images. Experimental results show desirable performance of our proposed method to segment synthetic images and brain MR images in comparison with state-of-the-art methods such as LGDF even for images with high-level noise and bias.

BIOGRAPHIES

Zahra Shahvaran received the B.S. and M.S. degrees in Electrical engineering from Shiraz University of Technology, Shiraz, Iran, in 2009 and 2012, respectively.

Kamran Kazemi received the B.S. and M.S. degrees in Electrical engineering from Ferdowsi University of Mashhad and K. N. Toosi University of technology, Tehran in 2000 and 2002, respectively. He performed his Ph.D. degree in Electrical Engineering and Biomedical Engineering as cooperation between K. N. Toosi University of Technology, Tehran, Iran and University of Picardie Jules Verne (UPJV) Amiens, France. Actually, he is working as assistant professor in department of Electrical and Electronics Engineering, Shiraz University of Technology, Shiraz, Iran and also as invited researcher in the GRAMFC (Groupe de Recherche sur l’Analyse Multimodale de la Fonction Cérébrale) lab at the University of Picardie Jules Verne (UPJV), Amiens, France.

Mohammad Sadegh Helfroush received the B.S. and M.S. degrees in Electrical engineering from Shiraz University of Shiraz and Sharif University of technology, Tehran in 1993 and 1995, respectively. He performed his Ph.D. degree in Electrical Engineering from Tarbiat Modarres University, Tehran, Iran. Actually, he is working as associate professor in department of Electrical and Electronics Engineering, Shiraz University of Technology, Shiraz, Iran.

Nassim Jafarian received the B.S. and M.S. degrees in Electrical engineering from Babol University of Babol and Shiraz University of Technology, Shiraz, Iran, in 2008 and 2011, respectively.