# Active Contour Driven by Local Region Statistics and Maximum A Posteriori Probability for Medical Image Segmentation

Authors: Xiaoliang Jiang; Bailin Li; Qiang Wang; Jiajia Liu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper presents a novel active contour model in a variational level set formulation for simultaneous segmentation and bias field estimation of medical images. An energy function is formulated based on improved Kullback-Leibler distance (KLD) with likelihood ratio. According to the additive model of images with intensity inhomogeneity, we characterize the statistics of image intensities belonging to each different object in local regions as Gaussian distributions with different means and variances. Then, we use the Gaussian distribution with bias field as a local region descriptor in level set formulation for segmentation and bias field correction of the images with inhomogeneous intensities. Therefore, image segmentation and bias field estimation are simultaneously achieved by minimizing the level set formulation. Experimental results demonstrate desirable performance of the proposed method for different medical images with weak boundaries and noise.


---

## Body

## 1. Introduction

Image segmentation is an important procedure in various image processing and computer vision applications. Many algorithms such as histogram thresholding, region growing, fuzzy c-means (FCM), and active contour models (ACM) have been proposed for image segmentation. However, due to the presence of noise, low contrast, and intensity inhomogeneity in medical images, the segmentation is still a challenging problem in the majority of applications.Since being introduced by Kass et al. [1], active contour ideas have been widely used in the computer vision field. There are mainly two kinds of active contour models: edge-based methods [2–4] and region-based methods [5–8]. Edge-based methods utilize image gradients to drive the curve evolution process, which is usually sensitive to image noise and weak edges. On the other hand, instead of employing gradient information, region-based methods utilize image region information to capture targets of interest. One of the early efforts towards region-based active contours was made by the M-S model [9], which approximates a given image by a piecewise smooth image. Compared to the edge-based model, the region-based models are less sensitive to initial contours and have better performance for images with weak object boundaries. The region-based model is insensitive to initialization of the level set function and can recognize boundaries of the object efficiently. Therefore, region-based models, especially the C-V model [5], have been widely applied in image segmentation.Recently, Li et al. [10] proposed a variational level set approach to segment images and estimate the bias field, which is motivated by the weighted K-means clustering. The main advantages of this model are robustness to initialization, accurate segmentation of images with intensity inhomogeneity, and estimate of the bias field. However, only using intensity means as a criterion of classification, this method may not be able to precisely extract the object boundaries, especially when the object texture distributions have identical means but different variances. Wang et al. [11] proposed an active contour model driven by local Gaussian distribution fitting (LGDF) energy with a level set function and local means and variances as variables for more accurate segmentation. The LGDF model can distinguish regions with similar intensity means but different variances by using the local Gaussian distribution fitting energy. However, the drawbacks of LGDF model are that it may introduce many local minima and cannot estimate the bias field so as to correct the intensity inhomogeneity in the original image.In this paper, a new region-based active contour model based on improved Kullback-Leibler distance and local neighborhood information is presented for medical image segmentation. In our model, the statistics of image intensities belonging to each different object in local regions are characterized by Gaussian distributions with different means and variances. According to maximum a posteriori probability and Bayes’ rule, we first derive a local objective function for image intensities in a neighborhood around each pixel. Then this local energy function is integrated with respect to the neighborhood center over the entire image domain to give a global criterion. In a level set formulation, this global criterion defines an energy in terms of level set functions that represent a partition of the image domain and a bias field that accounts for the intensity inhomogeneity of the image. Experimental results demonstrate desirable performance of the proposed method for different medical images with weak boundaries and noise.The rest of the paper is structured as follows. In Section2, a review on some well-known region-based models is given. Section 3 shows the proposed model in detail. In Section 4, experimental results are given. Finally, conclusions are drawn in Section 5.


## 2. Background

### 2.1. The Description of C-V Model

The Chan and Vese model, integrating the Mumford-Shah model and level set, utilizes the region information for segmentation [5]. The level set formulation of the C-V mode, regarding the time evolution of the level set function φ, can be described as
(1)∂φ∂t=δε(φ)×[μdiv⁡(∇φ|∇φ|)-ν-λ1(I(x)-c1)2+λ2(I(x)-c2)2],
where I(x) is the original image, μ, ν, λ1, λ2 are the positive constants of each item and the δε(φ) is the Dirac function. The first term keeps the level set function smooth, while the third and fourth terms are the internal and external forces, respectively, that drive the contour towards object boundaries. In C-V model, c1 and c2 are defined as follows:
(2)c1=∫ΩI(x)Hε(φ(x))dx∫ΩHε(φ(x))dx,c2=∫ΩI(x)(1-Hε(φ(x)))dx∫Ω(1-Hε(φ(x)))dx.Obviously,c1 and c2 are related to the global properties of the image intensities inside and outside the contour. The C-V model has good performance in image segmentation due to its ability of obtaining a larger convergence range and being less sensitive to the initialization. But because this model does not include local information, contour curve will deviate from the target boundary and it may cause error segmentation results when the energy function reaches to the minima.


### 2.2. Method of Li et al.

Li et al. proposed a variational level set approach to estimate the bias field and to segment the images with inhomogeneous intensity, which is motivated by a weighted K-means clustering. It is based on the following model to describe images in [10]:
(3)I=bJ+n,
where I is the measured image, J is the true image, b is the bias field, n is the noise. In [10], the assumptions on the true image I and the bias field b can be stated as follows.(I)
The bias fieldb is slowly varying in the entire image domain.(II)
The true image intensitiesJ are approximately a constant within each class of tissue; that is, J(x)≈ci for x∈Ωi, with {Ωi}i=1N being a partition of Ω.LetOx={y:∥x-y∥≤ρ} be the circular neighborhood with a small radius ρ centered at each point x in the image domain Ω. Then, in each small region, an approximation is given by
(4)b(y)J(y)≈b(x)ci,y∈Ox∩Ωi.Considering the task of classifying the intensityI(y) in the neighborhood Ox into N classes, a local intensity energy function is defined as
(5)Ex=∑i=1N∫Ox∩ΩiK(x-y)|I(y)-b(x)ci|2dy,
where b(x)ci are the cluster centers to be optimized and K(x-y) is a nonnegative weighting function, which is defined as follows:
(6)K(x-y)={1ae-|u|2/2σ2,for|x-y|≤ρ0,otherwise,
where a is a constant. Thus, (5) can be rewritten as
(7)Ex=∑i=1N∫ΩiK(x-y)|I(y)-b(x)ci|2dy.In the method of Li et al., the case ofN=2 is mainly considered. And the image domain is divided into two regions by the zero level contour of a function φ; that is, φ>0 and φ<0. Using the Heaviside function Hε, we can get the energy function as follows:
(8)E(φ,b,c)=∫(∑i=12∫K(x-y)|I(y)-b(x)ci|2Mi(φ(y))dy)dx,
where M1(φ(x))=Hε(φ(x)) and M2(φ(x))=1-Hε(φ(x)).Since local image intensity information is embedded into the energy function, the method can deal with some kinds of images with intensities inhomogeneity and estimate of the bias field. However, only using intensity means as a criterion of classification, this method is hard to correctly segment the image when the mean values of intensities of background and object are rather close to each other.


### 2.3. The LGDF Model

Wang et al. [11] proposed an active contour model driven by local Gaussian distribution fitting (LGDF) energy with a level set function and local means and variances as variables for more accurate segmentation. The energy functional for the LGDF model for each pixel using the level set function φ is defined as follows:
(9)ELGDF(φ,u1,u2,σ1·σ2)=-∫ω(x-y)log⁡p1,x(I(y))M1(φ(y))dy-∫ω(x-y)log⁡p2,x(I(y))M2(φ(y))dy+ν∫|∇H(φ(x))|dx+μ∫12(|∇φ(x)|-1)dx,
where ν≥0 and μ≥0 are fixed parameters, ω(x-y) is a nonnegative weighing function, and pi,x(I(y)) is
(10)pi,x(I(y))=12πσi(x)exp⁡(-(ui(x)-I(y))22σi2(x)),hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhi=1,2,
where ui(x) and σi(x) are local intensity means and standard deviation, respectively.The LGDF model can distinguish regions with similar intensity means but different variances by using the local Gaussian distribution fitting energy. Thus the LGDF model is able to deal with both noise and intensity inhomogeneity. However, it still easily gets stuck in local minima and cannot estimate the bias field so as to correct the intensity inhomogeneity in the original image.


## 2.1. The Description of C-V Model

The Chan and Vese model, integrating the Mumford-Shah model and level set, utilizes the region information for segmentation [5]. The level set formulation of the C-V mode, regarding the time evolution of the level set function φ, can be described as
(1)∂φ∂t=δε(φ)×[μdiv⁡(∇φ|∇φ|)-ν-λ1(I(x)-c1)2+λ2(I(x)-c2)2],
where I(x) is the original image, μ, ν, λ1, λ2 are the positive constants of each item and the δε(φ) is the Dirac function. The first term keeps the level set function smooth, while the third and fourth terms are the internal and external forces, respectively, that drive the contour towards object boundaries. In C-V model, c1 and c2 are defined as follows:
(2)c1=∫ΩI(x)Hε(φ(x))dx∫ΩHε(φ(x))dx,c2=∫ΩI(x)(1-Hε(φ(x)))dx∫Ω(1-Hε(φ(x)))dx.Obviously,c1 and c2 are related to the global properties of the image intensities inside and outside the contour. The C-V model has good performance in image segmentation due to its ability of obtaining a larger convergence range and being less sensitive to the initialization. But because this model does not include local information, contour curve will deviate from the target boundary and it may cause error segmentation results when the energy function reaches to the minima.


## 2.2. Method of Li et al.

Li et al. proposed a variational level set approach to estimate the bias field and to segment the images with inhomogeneous intensity, which is motivated by a weighted K-means clustering. It is based on the following model to describe images in [10]:
(3)I=bJ+n,
where I is the measured image, J is the true image, b is the bias field, n is the noise. In [10], the assumptions on the true image I and the bias field b can be stated as follows.(I)
The bias fieldb is slowly varying in the entire image domain.(II)
The true image intensitiesJ are approximately a constant within each class of tissue; that is, J(x)≈ci for x∈Ωi, with {Ωi}i=1N being a partition of Ω.LetOx={y:∥x-y∥≤ρ} be the circular neighborhood with a small radius ρ centered at each point x in the image domain Ω. Then, in each small region, an approximation is given by
(4)b(y)J(y)≈b(x)ci,y∈Ox∩Ωi.Considering the task of classifying the intensityI(y) in the neighborhood Ox into N classes, a local intensity energy function is defined as
(5)Ex=∑i=1N∫Ox∩ΩiK(x-y)|I(y)-b(x)ci|2dy,
where b(x)ci are the cluster centers to be optimized and K(x-y) is a nonnegative weighting function, which is defined as follows:
(6)K(x-y)={1ae-|u|2/2σ2,for|x-y|≤ρ0,otherwise,
where a is a constant. Thus, (5) can be rewritten as
(7)Ex=∑i=1N∫ΩiK(x-y)|I(y)-b(x)ci|2dy.In the method of Li et al., the case ofN=2 is mainly considered. And the image domain is divided into two regions by the zero level contour of a function φ; that is, φ>0 and φ<0. Using the Heaviside function Hε, we can get the energy function as follows:
(8)E(φ,b,c)=∫(∑i=12∫K(x-y)|I(y)-b(x)ci|2Mi(φ(y))dy)dx,
where M1(φ(x))=Hε(φ(x)) and M2(φ(x))=1-Hε(φ(x)).Since local image intensity information is embedded into the energy function, the method can deal with some kinds of images with intensities inhomogeneity and estimate of the bias field. However, only using intensity means as a criterion of classification, this method is hard to correctly segment the image when the mean values of intensities of background and object are rather close to each other.


## 2.3. The LGDF Model

Wang et al. [11] proposed an active contour model driven by local Gaussian distribution fitting (LGDF) energy with a level set function and local means and variances as variables for more accurate segmentation. The energy functional for the LGDF model for each pixel using the level set function φ is defined as follows:
(9)ELGDF(φ,u1,u2,σ1·σ2)=-∫ω(x-y)log⁡p1,x(I(y))M1(φ(y))dy-∫ω(x-y)log⁡p2,x(I(y))M2(φ(y))dy+ν∫|∇H(φ(x))|dx+μ∫12(|∇φ(x)|-1)dx,
where ν≥0 and μ≥0 are fixed parameters, ω(x-y) is a nonnegative weighing function, and pi,x(I(y)) is
(10)pi,x(I(y))=12πσi(x)exp⁡(-(ui(x)-I(y))22σi2(x)),hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhi=1,2,
where ui(x) and σi(x) are local intensity means and standard deviation, respectively.The LGDF model can distinguish regions with similar intensity means but different variances by using the local Gaussian distribution fitting energy. Thus the LGDF model is able to deal with both noise and intensity inhomogeneity. However, it still easily gets stuck in local minima and cannot estimate the bias field so as to correct the intensity inhomogeneity in the original image.


## 3. Proposed Method

### 3.1. Energy Functional Formulation

Recently, a kind of segmentation ordinations based on Kullback-Leibler distance, also named entropy, such as the minimum entropy method [12] and the maximum posterior entropy method [13], has been proposed. When comparing with other segmentation methods, these methods can achieve better segmentation accuracy, especially for image with the speckles and noises. KLD provides a mathematical formulation for measuring the difference between two probability density functions (PDFs). And the KLD-based image segmentation method significantly depends on the PDF. For discrimination between two probability density functions P1,x(I(x)) and P2,x(I(x)), the KLD can be defined as follows:
(11)DKL(P1,x∣P2,x)=∫ΩP1,x(I(x))log⁡P1,x(I(x))P2,x(I(x))dx.We take the logarithm forP1,x(I(x)) in (11). So, the improved KLD formulations can be rewritten as
(12)DKL(P1,x∣P2,x)=∫Ωlog⁡P1,x(I(x))log⁡P1,x(I(x))P2,x(I(x))dx.In this paper, we give the definition of the KLD in a spatially continuous domain. For a fixed regular regionΩx⊂Ω, it is a neighborhood centered at a point x∈Ω. The KLD of the point x is defined by
(13)DKL(P1,x∣P2,x)=∫Ωx∩Ω1log⁡P1,x(I(y)∣y∈Ωx∩Ω1)hhhhhhhhh×log⁡P1,x(I(y)∣y∈Ωx∩Ω1)P2,x(I(y)∣y∈Ωx∩Ω2)dy.In (13), a kernel function is introduced. We define the following two objective functions and a local energy function:
(14)DKL(P1,x∣P2,x)=∫Ωx∩Ω1K(x-y)log⁡P1,x(I(y))log⁡P1,x(I(y))P2,x(I(y))dy.According to the Bayesian formulation and maximum a posteriori probability method [14, 15], we define the following objective function for each point x in the image domain Ω:
(15)E=∫Ω[λ1DKL(P1,x∣P2,x)+λ2DKL(P2,x∣P1,x)]dx=∫Ω[(P1,xP2,x)∫ΩxK(x-y)(λ1log⁡P1,x-λ2log⁡P2,x)×log⁡(P1,xP2,x)dy∫ΩxK(x-y)(λ1log⁡P1,x-λ2log⁡P2,x)]dx,
where Pi,x(I(y)) is given by
(16)Pi,x(I(y))=12πσi(x)exp⁡(-(ui-I(y))22σi2(x)),hhhhhhhhhhhhhhhhhhhhhhhhhhhhi=1,2.From (5) in the above section, we can prove that ui(x) can be approximated to b(x)ci:
(17)ui(x)=∑y∈Ox∩OiI(y)|Ω|≈∑y∈Ox∩Oib(y)J(y)|Ω|≈∑y∈Ox∩Oib(x)ci|Ω|=b(x)ci,
where |Ω| is the number of pixels in Ox∩Oi.Then (16) can be rewritten as
(18)Pi,x(I(y))=12πσi(x)exp⁡(-(b(x)ci-I(y))22σi2(x)),hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhi=1,2.Remark. The proposed model is different from the method in [11] for the following reason. We introduced the Kullback-Leibler distance with likelihood ratios as measure of classification into this model. It can enhance the accuracy of segmentation. Moreover, our model can estimate the bias field so as to correct the intensity inhomogeneity in the original image. From (15), we can see that the energy function has the following properties:(I)
E=0 denotes that two region distributions are equal;(II)
ifE is minimization, two region distributions are not equal completely; that is, the reparability is better.


### 3.2. Level Set Formulation

Using the Heavside functionHε, then the energy functional (15) can written as follows:
(19)E=∫Ω[(P2,xP1,x)∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy+∫Ωxλ2K(x-y)log⁡P2,xhhhhhh×log⁡(P2,xP1,x)(1-Hε(φ(y)))dy∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy]dx.In order to obtain accurate evolution of the level set function, we need a signed distance function (SDF) to regularize the level set function. The SDF can be characterized as follows:(20)P(φ)=∫12(|∇φ(x)|2-1)dx.To regularize the level set contour, we need its length to derive a smooth contour during evolution:(21)L(φ)=∫|∇Hε(φ(x))|dx.The regularized Heaviside functionHε is defined by
(22)Hε(x)=12(1+2πarctan(xε)).Using the following regularized Dirac function,(23)δε(x)=1πεε2+x2.Thus, the entire energy functional is defined as(24)ε(φ)=α∫Ω[(P2,xP1,x)∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy+∫Ωxλ2K(x-y)log⁡P2,xhhhhhhhhhhh×log⁡(P2,xP1,x)(1-Hε(φ(y)))dy∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy]dx+μ∫12(|∇φ(x)|2-1)dx+ν∫|∇Hε(φ(x))|dx,
where α, λ1, λ2, μ, ν are positive constants.


### 3.3. Gradient Descent Flow

Minimization of the energy functionε(φ) in (24) with respect to φ can be achieved by solving the gradient descent flow equation:
(25)∂φ∂t=-αδε(φ)(e1-e2-e3)+μ(Δφ-div⁡(∇φ|∇φ|))+νδε(φ)div⁡(∇φ|∇φ|),
where
(26)e1=∫Ωxλ1K(x-y)[log⁡(2πσ1)+(I-b(x)c1)22σ12]2dy,e2=∫Ωxλ2K(x-y)[log⁡(2πσ2)+(I-b(x)c2)22σ22]2dy,e3=∫Ωx(λ1-λ2)K(x-y)×[log⁡(2πσ1)+(I-b(x)c1)22σ12]×[log⁡(2πσ2)+(I-b(x)c2)22σ22]dy.By fixingφ, σ, c, we find an optimal bias field that minimizes ε(φ):
(27)b=∑i=12(IciMi(φ(y)))*K/σi2∑i=12(ci2Mi(φ(y)))*K/σi2.Using the similar optimal process, the global meansci and local variance σi2 can be written as
(28)ci=∫(b*K)IMi(φ(y))dx∫(b2*K)Mi(φ(y))dx,σi2=∫K(y-x)(I(y)-b(x)ci)2Mi(φ(y))dy∫K(y-x)Mi(φ(y))dy.Equation (25) is the level set evolution equation to be solved in the proposed model. The term -αδε(φ)(e1-e2-e3), which is derived from the first term of the energy ε(φ), plays a key role in the proposed model since it is responsible for driving the active contour towards object boundaries.


## 3.1. Energy Functional Formulation

Recently, a kind of segmentation ordinations based on Kullback-Leibler distance, also named entropy, such as the minimum entropy method [12] and the maximum posterior entropy method [13], has been proposed. When comparing with other segmentation methods, these methods can achieve better segmentation accuracy, especially for image with the speckles and noises. KLD provides a mathematical formulation for measuring the difference between two probability density functions (PDFs). And the KLD-based image segmentation method significantly depends on the PDF. For discrimination between two probability density functions P1,x(I(x)) and P2,x(I(x)), the KLD can be defined as follows:
(11)DKL(P1,x∣P2,x)=∫ΩP1,x(I(x))log⁡P1,x(I(x))P2,x(I(x))dx.We take the logarithm forP1,x(I(x)) in (11). So, the improved KLD formulations can be rewritten as
(12)DKL(P1,x∣P2,x)=∫Ωlog⁡P1,x(I(x))log⁡P1,x(I(x))P2,x(I(x))dx.In this paper, we give the definition of the KLD in a spatially continuous domain. For a fixed regular regionΩx⊂Ω, it is a neighborhood centered at a point x∈Ω. The KLD of the point x is defined by
(13)DKL(P1,x∣P2,x)=∫Ωx∩Ω1log⁡P1,x(I(y)∣y∈Ωx∩Ω1)hhhhhhhhh×log⁡P1,x(I(y)∣y∈Ωx∩Ω1)P2,x(I(y)∣y∈Ωx∩Ω2)dy.In (13), a kernel function is introduced. We define the following two objective functions and a local energy function:
(14)DKL(P1,x∣P2,x)=∫Ωx∩Ω1K(x-y)log⁡P1,x(I(y))log⁡P1,x(I(y))P2,x(I(y))dy.According to the Bayesian formulation and maximum a posteriori probability method [14, 15], we define the following objective function for each point x in the image domain Ω:
(15)E=∫Ω[λ1DKL(P1,x∣P2,x)+λ2DKL(P2,x∣P1,x)]dx=∫Ω[(P1,xP2,x)∫ΩxK(x-y)(λ1log⁡P1,x-λ2log⁡P2,x)×log⁡(P1,xP2,x)dy∫ΩxK(x-y)(λ1log⁡P1,x-λ2log⁡P2,x)]dx,
where Pi,x(I(y)) is given by
(16)Pi,x(I(y))=12πσi(x)exp⁡(-(ui-I(y))22σi2(x)),hhhhhhhhhhhhhhhhhhhhhhhhhhhhi=1,2.From (5) in the above section, we can prove that ui(x) can be approximated to b(x)ci:
(17)ui(x)=∑y∈Ox∩OiI(y)|Ω|≈∑y∈Ox∩Oib(y)J(y)|Ω|≈∑y∈Ox∩Oib(x)ci|Ω|=b(x)ci,
where |Ω| is the number of pixels in Ox∩Oi.Then (16) can be rewritten as
(18)Pi,x(I(y))=12πσi(x)exp⁡(-(b(x)ci-I(y))22σi2(x)),hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhi=1,2.Remark. The proposed model is different from the method in [11] for the following reason. We introduced the Kullback-Leibler distance with likelihood ratios as measure of classification into this model. It can enhance the accuracy of segmentation. Moreover, our model can estimate the bias field so as to correct the intensity inhomogeneity in the original image. From (15), we can see that the energy function has the following properties:(I)
E=0 denotes that two region distributions are equal;(II)
ifE is minimization, two region distributions are not equal completely; that is, the reparability is better.


## 3.2. Level Set Formulation

Using the Heavside functionHε, then the energy functional (15) can written as follows:
(19)E=∫Ω[(P2,xP1,x)∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy+∫Ωxλ2K(x-y)log⁡P2,xhhhhhh×log⁡(P2,xP1,x)(1-Hε(φ(y)))dy∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy]dx.In order to obtain accurate evolution of the level set function, we need a signed distance function (SDF) to regularize the level set function. The SDF can be characterized as follows:(20)P(φ)=∫12(|∇φ(x)|2-1)dx.To regularize the level set contour, we need its length to derive a smooth contour during evolution:(21)L(φ)=∫|∇Hε(φ(x))|dx.The regularized Heaviside functionHε is defined by
(22)Hε(x)=12(1+2πarctan(xε)).Using the following regularized Dirac function,(23)δε(x)=1πεε2+x2.Thus, the entire energy functional is defined as(24)ε(φ)=α∫Ω[(P2,xP1,x)∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy+∫Ωxλ2K(x-y)log⁡P2,xhhhhhhhhhhh×log⁡(P2,xP1,x)(1-Hε(φ(y)))dy∫Ωxλ1K(x-y)log⁡P1,xlog⁡(P1,xP2,x)Hε(φ(y))dy]dx+μ∫12(|∇φ(x)|2-1)dx+ν∫|∇Hε(φ(x))|dx,
where α, λ1, λ2, μ, ν are positive constants.


## 3.3. Gradient Descent Flow

Minimization of the energy functionε(φ) in (24) with respect to φ can be achieved by solving the gradient descent flow equation:
(25)∂φ∂t=-αδε(φ)(e1-e2-e3)+μ(Δφ-div⁡(∇φ|∇φ|))+νδε(φ)div⁡(∇φ|∇φ|),
where
(26)e1=∫Ωxλ1K(x-y)[log⁡(2πσ1)+(I-b(x)c1)22σ12]2dy,e2=∫Ωxλ2K(x-y)[log⁡(2πσ2)+(I-b(x)c2)22σ22]2dy,e3=∫Ωx(λ1-λ2)K(x-y)×[log⁡(2πσ1)+(I-b(x)c1)22σ12]×[log⁡(2πσ2)+(I-b(x)c2)22σ22]dy.By fixingφ, σ, c, we find an optimal bias field that minimizes ε(φ):
(27)b=∑i=12(IciMi(φ(y)))*K/σi2∑i=12(ci2Mi(φ(y)))*K/σi2.Using the similar optimal process, the global meansci and local variance σi2 can be written as
(28)ci=∫(b*K)IMi(φ(y))dx∫(b2*K)Mi(φ(y))dx,σi2=∫K(y-x)(I(y)-b(x)ci)2Mi(φ(y))dy∫K(y-x)Mi(φ(y))dy.Equation (25) is the level set evolution equation to be solved in the proposed model. The term -αδε(φ)(e1-e2-e3), which is derived from the first term of the energy ε(φ), plays a key role in the proposed model since it is responsible for driving the active contour towards object boundaries.


## 4. Experimental Results

In this section, the proposed method has been tested on both synthetic and real medical images. The results are all carried out by Matlab 2011a in the PC with Pentium CPU 2.50 and 4 GB of RAM. In our experiments, the parameters are set as follows: iteration time stepΔt = 0.1, standard deviation of the Gaussian kernel σ=3, and weighting coefficients λ1=1.05, λ2=1.0, μ=1.0.


### 4.1. Application on Synthetic Images

We first apply our method to segment two synthetic images, which are displayed in the first column of Figure1. These two images are corrupted by strong noise and intensity inhomogeneity. And the initial contours and the final contours are plotted as red contours and blue contours, respectively. The final contours obtained by running our method are shown in the second column. The computed bias corrected images and bias fields are simultaneously obtained, as shown in the third and fourth columns. It can be seen that our method yields satisfactory segmentation results owing to considering and exploiting the image local region information, which can successfully extracts the object boundaries for these images.Application of our method on synthetic images. Column 1: original images and initial contours. Columns 2: final contours. Column 3: bias corrected images. Column 4: estimated bias fields.
(a)
(b)
(c)
(d)


### 4.2. Application on Real Medical Images

We also evaluate the performance of the proposed method on a set of real medical images in Figure2. The first row shows an MR image of a human bladder. The second row shows the results of our method on infrared images, which include some rather weak boundaries. The third row shows an X-ray vessel images with intensity inhomogeneity. In addition, parts of the vessel boundaries are quite weak. Such properties render it a nontrivial task to segment the vessels in the images. The original images and initial contours, the final results, the bias corrected images, and the bias fields are shown from the left column to the right column in Figure 2. Our method can get satisfactory segmentation results because it considers the statistical information in a transformed domain, where the intensity of the object and background are less overlapping than that in the original domain, making our method have a very strong discriminative capability for the object and background.Application of our method on real medical images. Column 1: original images and initial contours. Columns 2: final contours. Column 3: bias corrected images. Column 4: estimated bias fields.
(a)
(b)
(c)
(d)


### 4.3. Application on MR Image Segmentation

In order to further demonstrate the capability of our method in dealing with intensity inhomogeneity, we also use our method to segment images on two MR brain images with intensity inhomogeneity due to the bias field effect, which is usually seen in medical images. The segmentation results, the bias corrected images, and the bias fields are shown from the left column to the right column in Figure3. It can be clearly seen that the image quality is significantly improved by our method. Some regions whose intensity contrast is too low to be identified are able to be distinguished clearly after the bias correction.Application of our method on real medical images. Column 1: original images and initial contours. Columns 2: final contours. Column 3: bias corrected images. Column 4: estimated bias fields.
(a)
(b)
(c)
(d)


### 4.4. Comparison with Other Methods

Figure4 shows the comparison of the proposed method with C-V model [5], Li’s model [10], and LGDF [11] model in real medical images with intensity inhomogeneity. The first row shows the results for ultrasound images of left ventricle with severe intensity inhomogeneity and noise problems. The second row shows the left ventricle in a tagged MR image. It is evident that the image is corrupted by noise, severe intensity inhomogeneity, and weak boundaries. The third row shows the CT images of heart. It can be observed that the C-V model assumes that the image intensity is piecewise constant and uses the global intensity means to fit the image intensity. Therefore, they do not perform well in images with intensity inhomogeneity. Li’s model uses the local intensity mean to fit the measured image, and thus it can yield better segmentation results than C-V model on images with intensity inhomogeneity. However, if the intensity inhomogeneity is severe, using only the local mean information may fail to discriminate the intensity between an object and its background, leading to inaccurate segmentation. The LGDF model fully exploits intensity distribution information in local regions. However, it still easily gets stuck in local minima and cannot estimate the bias field. Our method can produce much better segmentation results because it can not only distinguish regions with similar intensity means but different variances, but also estimate the bias field so as to correct the intensity inhomogeneity in the original image, making our method have a very strong discriminative capability for the object and background.Comparison with other methods. Column 1: original images and initial contours. Column 2: C-V model. Column 3: Li’s model. Column 4: LGDF model. Column 5: our method.
(a)
(b)
(c)
(d)
(e)


### 4.5. Results for Synthetic Noisy Images

To further demonstrate the ability of the proposed model when dealing with noisy images, we create images by adding Gaussian noise to a synthetic image in Figure5. The two images are contaminated with multiplicative Gaussian noise with standard deviations of (0.1, 0.15, 0.2, and 0.25). It can be seen from the bottom row that the results appear similar visual quality for the original and noisy images. By quantitative comparison, we can see that our model do produce perfect results for these original and noisy images.Figure 5
Results for the images corrupted by additive Gaussian noise.


## 4.1. Application on Synthetic Images

We first apply our method to segment two synthetic images, which are displayed in the first column of Figure1. These two images are corrupted by strong noise and intensity inhomogeneity. And the initial contours and the final contours are plotted as red contours and blue contours, respectively. The final contours obtained by running our method are shown in the second column. The computed bias corrected images and bias fields are simultaneously obtained, as shown in the third and fourth columns. It can be seen that our method yields satisfactory segmentation results owing to considering and exploiting the image local region information, which can successfully extracts the object boundaries for these images.Application of our method on synthetic images. Column 1: original images and initial contours. Columns 2: final contours. Column 3: bias corrected images. Column 4: estimated bias fields.
(a)
(b)
(c)
(d)


## 4.2. Application on Real Medical Images

We also evaluate the performance of the proposed method on a set of real medical images in Figure2. The first row shows an MR image of a human bladder. The second row shows the results of our method on infrared images, which include some rather weak boundaries. The third row shows an X-ray vessel images with intensity inhomogeneity. In addition, parts of the vessel boundaries are quite weak. Such properties render it a nontrivial task to segment the vessels in the images. The original images and initial contours, the final results, the bias corrected images, and the bias fields are shown from the left column to the right column in Figure 2. Our method can get satisfactory segmentation results because it considers the statistical information in a transformed domain, where the intensity of the object and background are less overlapping than that in the original domain, making our method have a very strong discriminative capability for the object and background.Application of our method on real medical images. Column 1: original images and initial contours. Columns 2: final contours. Column 3: bias corrected images. Column 4: estimated bias fields.
(a)
(b)
(c)
(d)


## 4.3. Application on MR Image Segmentation

In order to further demonstrate the capability of our method in dealing with intensity inhomogeneity, we also use our method to segment images on two MR brain images with intensity inhomogeneity due to the bias field effect, which is usually seen in medical images. The segmentation results, the bias corrected images, and the bias fields are shown from the left column to the right column in Figure3. It can be clearly seen that the image quality is significantly improved by our method. Some regions whose intensity contrast is too low to be identified are able to be distinguished clearly after the bias correction.Application of our method on real medical images. Column 1: original images and initial contours. Columns 2: final contours. Column 3: bias corrected images. Column 4: estimated bias fields.
(a)
(b)
(c)
(d)


## 4.4. Comparison with Other Methods

Figure4 shows the comparison of the proposed method with C-V model [5], Li’s model [10], and LGDF [11] model in real medical images with intensity inhomogeneity. The first row shows the results for ultrasound images of left ventricle with severe intensity inhomogeneity and noise problems. The second row shows the left ventricle in a tagged MR image. It is evident that the image is corrupted by noise, severe intensity inhomogeneity, and weak boundaries. The third row shows the CT images of heart. It can be observed that the C-V model assumes that the image intensity is piecewise constant and uses the global intensity means to fit the image intensity. Therefore, they do not perform well in images with intensity inhomogeneity. Li’s model uses the local intensity mean to fit the measured image, and thus it can yield better segmentation results than C-V model on images with intensity inhomogeneity. However, if the intensity inhomogeneity is severe, using only the local mean information may fail to discriminate the intensity between an object and its background, leading to inaccurate segmentation. The LGDF model fully exploits intensity distribution information in local regions. However, it still easily gets stuck in local minima and cannot estimate the bias field. Our method can produce much better segmentation results because it can not only distinguish regions with similar intensity means but different variances, but also estimate the bias field so as to correct the intensity inhomogeneity in the original image, making our method have a very strong discriminative capability for the object and background.Comparison with other methods. Column 1: original images and initial contours. Column 2: C-V model. Column 3: Li’s model. Column 4: LGDF model. Column 5: our method.
(a)
(b)
(c)
(d)
(e)


## 4.5. Results for Synthetic Noisy Images

To further demonstrate the ability of the proposed model when dealing with noisy images, we create images by adding Gaussian noise to a synthetic image in Figure5. The two images are contaminated with multiplicative Gaussian noise with standard deviations of (0.1, 0.15, 0.2, and 0.25). It can be seen from the bottom row that the results appear similar visual quality for the original and noisy images. By quantitative comparison, we can see that our model do produce perfect results for these original and noisy images.Figure 5
Results for the images corrupted by additive Gaussian noise.


## 5. Conclusion

In this paper, we present a novel region-based active contour model in a variational level set formulation for simultaneous segmentation and bias field estimation of medical images. Based on objects of images with intensity inhomogeneity, we characterize the statistics of image intensities belonging to each different object in local regions as Gaussian distributions with different means and variances. Benefiting from the local variance, our method can obtain more accurate segmentation results. Experiments on synthetic and real medical images show that our model is desirable. In the future, our model will be extended from two-phase to multiphase level set formulation, which would further enhance its capability in processing more complex medical images with intensity inhomogeneity.



---

*Source: 182415-2014-07-08.xml*