This article has been accepted for inclusion in a future issue of this journal.

Content is final as presented, with the exception of pagination.

IEEE TRANSACTIONS ON CYBERNETICS 1

A Level Set Approach to Image Segmentation

With Intensity Inhomogeneity
Kaihua Zhang, Lei Zhang, Senior Member, IEEE, Kin-Man Lam, Senior Member, IEEE,
and David Zhang, Fellow, IEEE

Abstract—It is often a difficult task to accurately segment well-known Mumford–Shah (MS) model [7], which assumes
images with intensity inhomogeneity, because most of repre- that the image intensity is piecewise smooth (PS), is suitable
sentative algorithms are region-based that depend on intensity to model images with intensity inhomogeneity. The MS model
homogeneity of the interested object. In this paper, we present
a novel level set method for image segmentation in the pres- uses a set of contours S to separate different regions. However,
ence of intensity inhomogeneity. The inhomogeneous objects are it is difficult to minimize its energy functional because the
modeled as Gaussian distributions of different means and vari- set S of low dimension is unknown and the problem is
ances in which a sliding window is used to map the original nonconvex [8]. Some simplified versions of the MS model
image into another domain, where the intensity distribution of have been proposed, such as the Chan–Vese (CV) model [1]
each object is still Gaussian but better separated. The means
of the Gaussian distributions in the transformed domain can be and the PS model [8], which represent contour S as the zero
adaptively estimated by multiplying a bias field with the original level of a function called level set function, and then segmen-
signal within the window. A maximum likelihood energy func- tation proceeds by evolving a level set equation. However, the
tional is then defined on the whole image region, which combines CV model is not applicable to images with intensity inho-
the bias field, the level set function, and the piecewise constant mogeneity because it models images as piecewise constant
function approximating the true image signal. The proposed level
set method can be directly applied to simultaneous segmentation functions. The PS model is able to achieve a desirable seg-
and bias correction for 3 and 7T magnetic resonance images. mentation result for an image with intensity inhomogeneity.
Extensive evaluation on synthetic and real-images demonstrate However, it needs to iterate two partial differential equations,
the superiority of the proposed method over other representative which is very time-consuming and thereby limits its practical
algorithms. application.
Index Terms—Active contour model, bias field correction, Recently, some local region-based (LRB) level set methods
intensity inhomogeneity, level set method, segmentation. have been proposed to deal with images with intensity inho-
mogeneity, such as LRB method [9], local binary fitting (LBF)
model [10], local intensity clustering (LIC) method [11], local
I. I NTRODUCTION
region model (LRM) [12], patch driven level set method [13]
NTENSITY inhomogeneity caused by imperfection of
I imaging devices or illumination variations often occurs
in real-world images, which leads to serious misclassifica-
based on sparse representation [14]–[16], and edge driven
level set method [17], etc. However, they have some draw-
backs: the LRB method has two drawbacks. First, it used Dirac
tion by intensity-based segmentation algorithms that assume functional is restricted to a neighborhood around zero level set,
a uniform intensity [1]–[6]. Statistically, misclassification is which makes the level set evolution act locally. Therefore, the
caused by the prolonged tails of intensity distribution of each evolution can be easily trapped into local minima [1]. Second,
object so that it is difficult to extract the desired objects accu- its region descriptor is only based on region mean informa-
rately based on their respective intensity distributions. The tion without considering region variance and thereby may lead
to inaccurate segmentation. This drawback also holds for the
Manuscript received May 29, 2013; revised December 14, 2013,
March 18, 2014, and February 10, 2015; accepted February 27, 2015. This LBF model, which uses a similar energy functional. The LIC
work was supported in part by the Natural Science Foundation of China under method can be considered as a locally weighted K-means
Grant 61402233, and in part by the Startup Foundation for Introducing clustering method [11], which does not consider the cluster-
Talent of Nanjing University of Information Science and Technology
under Grant S8113049001. This paper was recommended by Associate ing variance; similar drawback often exists in the K-means
Editor X. Li. clustering-based methods [18]. The LRM exploits local region
K. Zhang is with the Jiangsu Key Laboratory of Big Data Analysis statistics, i.e., local region means and variances, to interpret
Technology and Jiangsu Collaborative Innovation Center on Atmospheric
Environment and Equipment Technology, Nanjing University of the MS model. However, the local region means and variances
Information Science and Technology, Nanjing 210044, China (e-mail: are only defined empirically, but not derived from minimizing
[email protected]). the MS energy.
L. Zhang and D. Zhang are with the Department of Computing, Hong Kong
Polytechnic University, Hong Kong. In this paper, we present a level set method for image
K.-M. Lam is with the Department of Electronic and Information segmentation. By exploiting the local image region statis-
Engineering, Hong Kong Polytechnic University, Hong Kong. tics, we define a mapping from the original image domain to
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org. another domain in which intensity probability model is more
Digital Object Identifier 10.1109/TCYB.2015.2409119 robust to noise while suppressing the intensity overlapping to
2168-2267 c 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

2 IEEE TRANSACTIONS ON CYBERNETICS

some extent. We then devise a maximum likelihood energy the set of level set function. H represents the Heaviside func-
functional based on the distribution of each local region in the tion and δ denotes the Dirac function. B denotes bias field
transformed domain, which combines the bias field, the level function and its optimal solution is denoted as B̃. The symbol J
set function, and the piecewise constant function approximat- denotes true image signal. N represents Gaussian white noise.
ing the true image signal. Analysis of the proposed approach σ and σi are standard deviations of the Gaussian distribution,
shows that it is a soft classification method, which means that respectively. The symbol θ i is the set of estimated parameters
each pixel can be assigned to more than one class. In con- {ci , σi } where θ = {θ i }, and their optimal solutions are denoted
trast, the hard classification [1] assigns each pixel to only one as θ̃ , c̃i , and σ̃i . The symbol Ox denotes a local region cen-
class [19]. In addition, the proposed method can be applied tering at location x. The symbol Li is the number of pixels.
to simultaneous tissue segmentation and bias correction for N denotes the Gaussian distribution. P denotes the probability
magnetic resonance (MR) images. density function (PDF). D denotes the set of image intensity in
The main differences between this paper and its prelimi- the transformed domain. Kρ is an indicator function of a local
nary presented in [20] are summarized as follows: 1) more region and Gσ denotes a Gaussian function. Mi denotes the
related works and technique details about the principle of the ith region indicator function. l denotes iterations. t and t1
proposed algorithm are provided such as the detailed devia- denote time steps, respectively. K denotes a kernel. E denotes
tions about how to optimize the proposed objective function an energy functional.
are provided in Appendixes A and B; 2) a new two-phase level
set models are included; 3) the relationships among the pro- III. BACKGROUND
posed methods and other related works are discussed in detail;
Mumford and Shah [7] approximated an image with a
and 4) extensive evaluations on synthetic and real images
PS function U(x), such that U varies smoothly within each
are performed, including comparisons with the related works
sub-region, and abruptly across their boundaries. An energy
such as CV [1], global CV (GCV) [21], LRB [9], LBF [10],
functional is defined as [7]
LIC [11] models.  
The main contributions of this paper are summarized as MS
EU,S = (I − U)2 dx + μ |∇U|2 dx + v|S| (1)
follows.
1) It is difficult to use local region statistics to well  \S
segment images with severe intensity inhomogeneity where μ and v > 0 are two fixed parameters and |S| represents
because the regions must have sharp discontinuities in the length of contour S. Image segmentation can be performed
the statistics [22]. To handle this problem, we propose a by minimizing (1) with respect to U and S. However, it is
very simple method by transforming the pixel intensities MS in practice, due to the unknown set S
difficult to minimize EU,S
into another domain (via averaging the pixel intensity of lower dimension and nonconvexity of the functional [8].
in a local region), in which we have theoretically vali- Many methods have been proposed to simplify or modify the
dated that the intensities in the transformed domain have functional [1], [8], which will be reviewed afterward. Here,
less overlapping in the statistics, thereby achieving bet- we mainly discuss the two-phase case in which region  is
ter segmentation results than some representative level separated by a contour S where  = in(S) ∪ out(S) and
set methods [1], [9]–[11], [21] for images with severe in(S) ∩ out(S) = ∅.
intensity inhomogeneity. Chan and Vese [1] proposed an active contour model that
2) The proposed method can yield the closed-form solu- is a special case of (1) in which U(x) in (1) is replaced by
tions for the estimated parameters in the distribution, a piecewise-constant function. They proposed to minimize the
which significantly reduces computation effort. following energy functional:
3) We show that some representative level set methods  
(e.g., CV [1], LBF [10], and LIC [11]) are the special EcCV = (I − c ) 2
dx + (I − c2 )2 dx + v|S| (2)
1 ,c2 ,S 1
cases of the proposed method.
in(S) out(S)

where c1 and c2 are two constant functions which approx-

II. N OTATIONS imate the average intensities inside and outside the con-
In the following sections, bold italic variables (e.g., x) tour S, respectively. With the zero level set of the level
denote vectors; small letters (e.g., n) denote scalars; and capital set function  to represent the contour S, this model
letters (e.g., I) denote functions. achieves favorable results for images with piecewise constant
All notations are listed as follows: pixel coordinates x, y ∈ intensities.
 ⊂ R2 , where  represents image domain. The symbol R The PS model [8] aims to minimize the energy func-
denotes the set of real number and ∅ is an empty set. The tional (1) using the level set method. The contour
symbol I denotes an input image and I denotes an image in S(p) : R →  is implicitly represented by a level set
the transformed domain. The symbols U, U1 , and U2 are PS function (x) :  → R, i.e., S = {x ∈ |(x) = 0}.
functions, respectively. The symbols μ > 0, v > 0, and ρ are By approximating an input image with two smooth func-
fixed parameters. Constant c1 and c2 approximate image inten- tions U1 (x) and U2 (x) in the sub-regions 1 = {x ∈
sities inside and outside the contour S, respectively. The level  : (x) > 0} and 2 = {x ∈  : (x) < 0},
set functions are denoted as , 1 , and 2 . The symbol  is respectively, the energy functional of the PS model is
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: LEVEL SET APPROACH TO IMAGE SEGMENTATION WITH INTENSITY INHOMOGENEITY 3

Fig. 1. Illustration of the image model of intensity inhomogeneity. Left to

right: the true signal J(x), the bias field function B(x), the noise N(x), and
the observed image I(x).

Fig. 3. Distributions of adjacent regions in the original image intensity

domain (blue solid curves) and in the transformed domain (red dashed curves).

intensity I can be approximated by a Gaussian distribution

with mean BJ and variance σ 2 . However, using only one
Gaussian model is not accurate enough to describe the statisti-
Fig. 2. Graphical representation of the local region i ∩ Ox , i = 1, 2. The cal characteristics of image intensity. Often multiple Gaussian
white arrows point to the regions belonging to 2 ∩ Ox while the black ones
point to the regions belonging to 1 ∩ Ox . The blue arrows point to the
probability distributions are adopted, with each one modeling
regions of Ox which are inside the yellow circles centered at each location x the distribution of image intensity in each object domain. The
(denoted by yellow dots). distribution corresponding to the object domain i is [22]
 
1 (I(y) − Ui (x))2
P(I(y)|θ i , B, x) = √ exp − (5)
defined as follows: 2π σi 2σi2
 
where Ui (x) is the spatially varying mean that is estimated
EU1 ,U2 , = (I − U1 ) H()dx + μ |∇U1 |2 H()dx
PS 2
at the local region i ∩ Ox centered at each location x
  (see Fig. 2), and σi is the standard deviation. Since B(x)

+ (I − U2 )2 (1 − H())dx varies slowly, it can be assumed to be a constant in a small

window [11]. That is Ui (x) ≈ B(x)ci for all x ∈ i ∩ Ox .
  We symbolize θ i = {ci , σi } and θ = {θ i , i = 1, . . . , n} in the
+ μ |∇U2 | (1 − H(φ))dx + v |∇H()|dx (3)
2
following discussions.
  Using the local region statistics (5) is able to well describe
the gradually varying regions that commonly exist in images
where H(·) is the Heaviside functional. The optimal U1 and U2 with intensity inhomogeneity. However, the regions must have
must be obtained by solving two damped Poisson equations sharp discontinuities in the statistics for better separation [22].
before each iteration of the level set evolution and thereby the In the following discussions, we propose a simple method to
computational cost is very high. Moreover, U1 and U2 must map the image intensity to another domain in which all regions
be extended to the whole image domain, which is difficult to have sharp discontinuities in the statistics. The main idea of
implement in practice. our approach is to reduce the overlapping tails of the distri-
In the next section, we propose a statistical level set method bution in the transformed domain by compressing its profile
in which the CV and PS models can be viewed as two special with mean Ui fixed [2] (see the red dashed curves in Fig. 3),
cases based on some proper assumptions (see Section IV-E). thereby equipping it with better separation capability for all
regions in the statistics.
IV. M ETHODOLOGY
A. Statistical Model of Intensity Inhomogeneity B. Principle of the Proposed Method
Let B(x) :  → R be an unknown bias field, J(x) :  → R As shown by Fig. 2, we denote Ox a neighboring region
be the true signal to be restored, and N(x) :  → R be noise. centered at location x, i.e., Ox = {y||y − x| ≤ ρ}, where ρ
As illustrated by Fig. 1, we consider the following model of is the radius of the region Ox . The whole image domain 
intensity inhomogeneity [11], [23]: can be represented as  = ∪i=1,...,n i with i ∩ j = ∅,
for all i = j, where i is the ith object region. We define
I(x) = B(x)J(x) + N(x). (4) a mapping T : I(x|θ i , B) → I(x|θ i , B) from original image
intensity domain D(T ) to another domain R(T ) by averaging
Suppose there are n objects in the image domain , and i
is the domain of the ith object. The true signal J(x) is often image intensities

assumed to be piecewise constant within each object domain, 1
I(x|θ i , B) = I(y|θ i , B, x)dy (6)
i.e., J(x) = ci for x ∈ i , where ci is a constant. The bias field Li (x)
B(x) is often assumed to be smooth in the image domain . i ∩Ox

The noise N(x) is assumed to be Gaussian-distributed with where Li (x) = |i ∩ Ox | is the number of pixels in region
zero mean and variance σ 2 . Thus, the distribution of the image i ∩ Ox .
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

4 IEEE TRANSACTIONS ON CYBERNETICS

Remark: In this paper, (6) is not a preprocessing step of We have the following joint likelihood function:
a test image because each region i is unknown at first,

n 
n 
and our objective is to extract i by optimizing a statisti- P(D|θ, B) = P(D|θ i , B) = P(I(x|θ i , B))
cal model based on I, which is computed with the statistical i=1 i=1 x∈
information of I via their linkage in (6). This is signifi- (11)
n  
cantly different from that uses a bilateral filter [24] on the ∝ P(I(y|θ i , B, x)),
test image as a preprocessing step (see comparison results in i=1 x∈ y∈i ∩Ox
Section VI-C), Furthermore, the bilateral filter combines image
where θ = {θ i , i = 1, . . . , n}.
intensities based on both their geometric closeness and pho-
Our objective is to estimate the parameter set θ̃ and the bias
tometric similarity, but the proposed method just averages the
field function B̃ that maximize the joint likelihood P(D|θ, B)
image intensities in a local region without considering their
geometric or photometric similarity. {θ̃ , B̃} = arg max P(D|θ, B) (12)
The intensity of pixel x is assumed to be inde- θ,B
pendently distributed with a Gaussian distribution [2], which corresponds to the values of θ and B that best agree with
i.e., P(I(y|θ i , B, x)) = N (I|Ui (x), σi2 ), where Ui is the spa- the actually observed training samples {I(x), x ∈ }. From the
tial varying mean and σi is the standard deviation subject to definition of P(D|θ, B) in (11), it is intuitively observed that
the object in region i . Thus, for all I(x|θ i ) ∈ R(T ), the the solution {θ̃, B̃} in (12) can also maximize the joint like-
corresponding PDF is still Gaussian [18] as lihood Q(I(x)|θ , B) in (17) and the likelihood P(I(x)|θ i , B)
  in (9).
σi2 We define an energy functional EL (θ, B) that is the inverse
P(I(x|θ i , B)) = N I|Ui (x), . (7)
Li (x) log-likelihood function of P(D|θ, B) in (11)

Referring to the red dashed curves in Fig. 3, the overlapping EθL,B  − log P(D|θ , B)
n  
tails of the distributions are suppressed to some extent, thereby
the pixel intensities in the transformed domain are much easier = constant − log(P(I( y|θ i , B, x)))dydx.
to be classified than that in the original domain. i=1   ∩O
i x
In the following discussion, we first bridge the statistical (13)
linkage between the intensities I and I via (9) because the
PDF of I is assumed to be known [see (5)]. Then, we build Let Kρ (x, y) be the indicator function of region Ox (see
an objective function (11), which is the likelihood function by Fig. 2)

using the statistical information of I. Finally, we use the log- 1, |y − x| ≤ ρ
likelihood function [see (15)] with the statistical information Kρ (x, y) = (14)
0, else.
of I as the total objective function.
Since the intensity inhomogeneity manifests itself as a Putting (5) and (14) into (13), and eliminating the trivial
smooth intensity variation across an image (see Fig. 1), sim- constant term, EL (θ , B) can be rewritten as
ilar to [11], we can reasonably assume that I(y|θ i , B) ≈ n    
 (I( y) − B(x)ci )2
 i , B), for all y ∈ i ∩ Ox . L (x)
I(x|θ Thus, we have
EθL,B = Kρ (x, y) log(σi ) + dydx.
y∈i ∩Ox P(I(y|θ i , B)) ≈ P(I(x|θ i , B)) i . Because the 2σi2
i=1  i
product of Gaussian PDFs is still Gaussian. That is
P(I(x|θ i , B))Li (x) ∝ N (I|Ui (x), (σi2 /Li (x))). Therefore, we (15)
have
  C. Discussion
 σi2
P(I(y|θ i , B)) ∝ N I|Ui (x), . (8) In this section, we discuss some advantages of the proposed
Li (x) method, including achieving a soft classification, robustness
y∈i ∩Ox
against noise, and alleviating side effect of oversmoothing
Putting (7) into (8), we have object boundaries.
The joint likelihood function (11) can be rewritten as
 
P(I(x|θ i , B)) ∝ P(I( y|θ i , B, x)) (9) P(D|θ , B) = Q(I(x|θ, B)) (16)
y∈i ∩Ox x∈

where
which is the likelihood of θ i with respect to the samples
{I(y), y ∈ i ∩ Ox }. 
n

Let D = {I(x|θ i , B), x ∈ , i = 1, . . . , n}, we have the Q(I(x|θ, B)) = P(I(x|θ i , B)). (17)
following likelihood function for the ith object [18]: i=1

 The joint likelihood function P(D|θ, B) in (16) can be viewed

P(D|θ i , B) = P(I(x|θ i , B)). (10) as the likelihood function of parameter set θ associated with
x∈ the intensity I whose PDF is Q(I(x|θ, B)).
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: LEVEL SET APPROACH TO IMAGE SEGMENTATION WITH INTENSITY INHOMOGENEITY 5

Put (8) and (9) into (17), Q(I(x|θ, B)) can be rewritten as where Fi (y)   Kρ (x, y)(log(σi ) + ((I(y) − B(x)ci )2 /
a univariate Gaussian distribution 2σi2 ))dx, and n = 2 or 4, which constructs the force to drive
n the level set function to evolve. When the zero level set keeps
Q(I(x|θ , B)) = P(I(x|θ i , B)) ∝ N (I|U(x),
(x)) (18) stable, we obtain the final segmentation results in which each
i=1 region is indicated by its membership function Mi in (22).
where

n
Li (x)Ui (x) 
n
Li (x) E. Relation to Other Models
U(x) =
(x) ,
−1 (x) = . (19)
i=1
σi2 i=1
σi2 In this section, we explain the relationships of our method
with six representative active contour models, i.e., the CV
As shown by the PDF in (18) whose parameters are a model [1], the PS model [8], the LRM method [12], the LBF
combination of the parameters subject to all class intensities model [10], the LIC model [11], and the Local Gaussian
[see (19)], each pixel (voxel) intensity I(x|θ, B) is subject distribution (LGD) model [26] in detail.
to multiple classes. Therefore, our method achieves a soft We only discuss the two-phase level set method while simi-
classification, satisfying the condition of the partial volume lar discussions can be readily extended to the four-phase case.
effect [25] (i.e., the intensity of each volume voxel is mixed Let us revisit the proposed energy functional in (23). When
from multiple classes [25]). Moreover, as shown from (6), the the variance σi2 = 1, i = 1, 2, and the bias field B(x) = 1, (23)
intensity in the transformed domain is obtained by averaging can be rewritten as
the neighboring pixel intensities belonging to the same class 2  
that is a low-pass filter process. Thus, the classification result 1
Eθ , =
L
Kρ (x, y)(I(y) − ci )2 Mi ((y))dxdy
is less sensitive to noise. Finally, our method alleviates the side 2
i=1  
effect of oversmoothing object boundaries due to the follow-
2  
ing reasons. First, only the neighboring intensities belonging 1
to the same class contribute to each class of I in (6); sec- = (I(y) − ci )2 Mi ((y))dy Kρ (x, y)dx
2
i=1  
ond, the overlapping parts of the statistical distributions among
different classes of intensities are suppressed (see Fig. 3). 2 2 
πρ
= (I(y) − ci )2 Mi ((y))dy. (24)
2
D. Energy Functional Using the Level Set Method i=1 

One level set function  can only represent two regions, The above energy functional is similar to that of the
inside and outside the contour S, as 1 = in(S) = { > 0} well-known CV model [1] except for the trivial constant
and 2 = out(S) = { < 0}, respectively. This is called (πρ 2 /2) [see (2)].
the two-phase model. If there exist more than two different In (23), if we assume that Ui (x) = B(x)ci is a PS func-
regions, two level set functions, i.e., 1 and 2 , can be used tion and σi = 1, the objective functional EL (θ, B, ) can be
to represent all of the different regions based on the four-color rewritten as that in (24), where c1 and c2 are replaced by two
theorem [8] such that any two adjacent regions can be indi- PS functions, respectively, which is similar to the objective
cated by different colors. This is called the four-phase model. functional of the PS model [see (3)] [8].
We define the phase indicators as follows [8]: The LRM [12] handles inhomogeneous intensities well by
 using spatially varying means and variances to replace the con-
M1 () = H()
two phase: (20) stant means and variances in (2). However, the LRM method
M2 () = 1 − H() directly introduces a Gaussian kernel to compute the varying
⎧ means and variances, which is inconsistent with theory [26].
⎪ M () = H(1 )H(2 )
⎪ 1
⎪ The data-fitting term of the energy functional in the LBF
⎨ M () = H( )(1 − H( ))
2 1 2 model [10] is as follows:
four phase: (21)  
⎪ M3 () = (1 − H(1 ))H(2 )
⎪
⎪
⎩ LBF
E,U ,U = Gσ (y − x)(I(y) − U1 (x))2 H((y))dxdy
M4 () = (1 − H(1 ))(1 − H(2 )) 1 2
 
where H(·) is the Heaviside functional and  represents the set  
of level set functions, i.e.,  = {} for the two-phase model + Gσ (y − x)(I(y) − U2 (x))2 (1 − H((y)))dxdy
and  = {1 , 2 } for the four-phase model. Each phase indic-  
tor function Mi () in (20) or (21) is the membership function (25)
of the region i satisfying where Gσ (·) is a truncated
 Gaussian kernel with standard devi-
1, y ∈ i ation σ , which satisfies  Gσ (x)dx = 1, and U1 (·) and U2 (·)
Mi ((y)) = (22) are two smoothing functions that approximate the weighted
0, else.
average image intensities in a Gaussian window in the regions
Putting (22) into (15), the energy functional EL (θ, B) is 1 = {(x) > 0} and 2 = {(x) < 0}, respectively.
extended to the whole image domain It is worth noting that the principles of the proposed
n 
 model, the LBF model [10] and LIC model [11] are different.
EθL,B, = Fi (y)Mi ((y))dy (23) If we set σi2 = 1, B(x)ci = Ui (x) in (23), and use a trun-
i=1  cated Gaussian kernel Kρ with standard deviation ρ satisfying
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

6 IEEE TRANSACTIONS ON CYBERNETICS

 Kρ (x)dx = 1 to replace the constant kernel, then (23) can Note that B̃ is actually the normalized convolution [27], which
be written as naturally leads to a smooth approximation of the bias field B.
 
L 1 3) Minimization With Respect to σi : By fixing the other
EU 1 ,U2 ,
= Kρ (y − x)(I(y) − U1 (x))2 H((y))dxdy
2 variables in (23), we can obtain the minimizer of σi , denoted

  by σ̃i , as follows:
1 
+ Kρ (y − x)(I(y) − U2 (x))2 (1 − H((y)))dxdy. 
2  K (y, x)Mi ((y))(I(y) − B(x)ci )2 dydx
   ρ
(26) σ̃i =  . (29)
K (y, x)M ((y))dydx
ρ i
This LBF
is the same as E,U except for the trivial constant 
1 ,U2
1/2. If we set Ui = Bci in (26), it is the energy functional of
For an explanation of how to derive above solutions, please
the LIC model [11]. It should be noted that in the theoretical
refer to Appendixes A and B.
analysis of our method in Section IV-B, the constant kernel
The solutions θ̃ = {c̃i , σ̃i , i = 1, . . . , n, n = 2, or 4} and B̃
is reasonably used as a local region indicator. However, both
are then put into (23) in which the solution of Fi is
of the LBF and LIC models use a Gaussian kernel as the
  
locally spatially weighted function to relate the pixel x and its (I(y) − B̃(x)c̃i )2
neighboring pixel y. The closer y is to x, the larger weight is F̃i (y)  Kρ (x, y) log(σ̃i ) + dx. (30)
2σ̃i2
assigned, thereby representing the higher similarity between 
the intensities of pixels y and x.
The energy functional of the LGD model [26] is slightly dif- B. Two-Phase Level Set Evolution Formula
ferent from the LBF model which replaces the 2 -norm terms Minimizing the energy functional EL with respect to ,
θ̃ ,B̃,
in (26) with ((I(y) − Ui (x))2 /σi2 (x)), i = 1, 2, where Ui (x) we have the corresponding gradient descent formula as fol-
and σi (x) are spatially locally varying mean and variance of lows:
a Gaussian distribution. However, as pointed out by [22], the
∂ ∂E(θ̃ , B̃, )
spatially varying variance may be unstable due to its local =− = (F̃2 − F̃1 )δ() (31)
property. Our model is different from the LGD model due to ∂t ∂
the following reasons. First, the variances of the Gaussian dis- where δ() is the Dirac functional.
tributions in our model are piecewise constant in each region, In order to keep the numerical implementation stable, the
thereby making our method much more stable than the LGD level set function should be regularized during the iteration
model; second, we use a constant kernel to indicate a local of (31). Li et al. [28] proposed a signed distance regulariza-
region which is based on a solid theoretical foundation as tion formula. However, it produces some unnecessary valleys
explained in Section IV-B. However, the Gaussian kernel is and peaks, which makes level set evolution easy to fall into
used in LGD model whose physical meaning is similar to the some local minima. In this paper, we adopt a simple and sta-
LBF model; third, our model can be used for simultaneous ble method [29] to regularize the level set function during
segmentation and bias correction while the LGD model can iteration. After each iteration of the level set evolution, we
be only used for segmentation. diffuse the level set function using the following formula:

V. E NERGY M INIMIZATION AND L EVEL S ET l+1 = l + t · ∇ 2 l (32)

E VOLUTION F ORMULAS where l represents the level set function obtained from the
Minimizing with respect to each variable in θ =
EθL,B, lth iteration of (31), ∇ 2 represents the Laplacian operator, and
{ci , σi , i = 1, . . . , n} and B, we can be obtained their closed- t represents the diffusion strength. The l+1 in (32) can
form solutions, which are described below. also be approximated by l+1 = K ⊗ l , where K is either a
Gaussian kernel [30] or a constant kernel [20].
A. Closed-Form Solutions for Different Variables
1) Minimization With Respect to ci : By fixing the other C. Four-Phase Level Set Evolution Formula
variables in (23), we obtain the minimizer of ci , denoted by
c̃i , as follows: Minimizing the energy functional EL (θ̃ , B̃, ) with respect
to 1 and 2 , respectively, we have the corresponding
(Kρ ⊗ B)IMi ((y))dy gradient-descent formulas

c̃i = (27) ⎧  
(Kρ ⊗ B2 )Mi ((y))dy ⎪ ∂1

⎨ = − (F̃1 − F̃2 − F̃3 + F̃4 )H(2 ) + F̃2 − F̃4 δ(1 )
∂t
where ⊗ denotes the convolution operator. ⎩ ∂2 = −(F̃1 − F̃2 − F̃3 + F̃4 )H(1 ) + F̃3 − F̃4 δ(2 ).
⎪
2) Minimization With Respect to B: By fixing the ∂t
other variables in (23), we obtain the minimizer of B, (33)
denoted by B̃, as follows:
n ci Each level set function is regularized by the following
i=1 Kρ ⊗ (IMi ((x))) · σ 2 formula after each iteration of (33):
B̃(x) =  i
. (28)
n c2i
K
i=1 ρ ⊗ M i ((x)) · il+1 = li + t · ∇ 2 li , i = 1, 2. (34)
2 σi
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: LEVEL SET APPROACH TO IMAGE SEGMENTATION WITH INTENSITY INHOMOGENEITY 7

Algorithm 1 Locally Statistical Level Set Method

1: Initialization: B̃ = 1, σ̃i = i (i is taken as both variable
value and index), i = 1, . . . , n, n = 2, or 4, and the level
set function li = 0i , i = 1, . . . , n
2: Update ci to c̃i , i = 1, . . . , n, by (27)
3: Update B to B̃ by (28)
4: Update σi to σ̃i , i = 1, . . . , n, by (29)
5: Update Fi to F̃i , i = 1, . . . , n by (30)
6: Evolve the level set function according to (32) or (34)
once
7: Regularize the level set function according to (35)
8: If il+1 satisfies the stationary condition, stop; otherwise,
l = l + 1 and return to Step 2.

D. Numerical Implementation Fig. 4. Segmentation results on two synthetic images and real-vessel images
with intensity inhomogeneity. From top to bottom: segmentation results by
We only need to approximate the temporal derivatives the CV model [1], the GCV model [21], and our method. The blue circles
in (32) and (33) as a forward difference because there are represent the initial contours and the red lines represent the final segmentation
no partial derivatives. The Laplacian operator ∇ 2  is approx- contours. We set ρ = 12 for the left image and ρ = 6 for the other images.
imated as ∇ 2  ≈ K ⊗ , where ⊗⎛is a convolution⎞ operator,
0 1 0
and K is a kernel defined as K = ⎝ 1 − 4 1 ⎠. Therefore, are some representative level set methods for image
0 1 0 segmentation. The MATLAB source codes and some
the solution of the diffusion equation (32) can be discretized examples of the proposed method can be down-
as follows: loaded at http://www.comp.polyu.edu.hk/∼cslzhang/LSACM/
⎛ ⎞
0 t 0 LSACM.htm
l+1 = l ⊗ ⎝ t 1 − 4t t ⎠. (35) We initialize B̃(x) = 1, σ̃i = i, i = 1, . . . , n, and then the
0 t 0 initializations of c̃i , i = 1, . . . , n can be calculated by (27).
The standard Von Neumann analysis [29] can be used to ana- We have tested different values to initialize B̃(x) and c̃i , and
found that their corresponding results are similar, which can be
 step t. Putting
lyze the stability for the time
√
l (i, j) =
attributed to the convexity of our objective functional (23) with
rn eiξ1 +jξ2 into (35), where ξ12 + ξ22 = −1 denotes the
respect to these variables. The time step for level set evolution
imaginary unit, we obtain the amplification factor as follows:
  is set t1 = 1, the time step for regularization is set t = 0.1,
r = 1 + 2t · cos(ξ1 ) + cos(ξ2 ) − 2 . (36) and ε = 1 for all the experiments except for Figs. 11–13, in
which we set t = 0.01. Our method is stable for a wide range
Therefore, we have 1 − 8t ≤ r ≤ 1. By solving the
of ρ, e.g., 5 < ρ < 25. In most cases, we set ρ = 6. A small
inequality |1 − 8t| ≤ 1, we have
ρ makes computation in each iteration more efficient, but the
0 ≤ t ≤ 0.25. (37) convergence of the algorithm is slower. On the other hand,
The Heaviside functional H(z) is approximated by a smooth a large ρ increases computational burden in each iteration.
function Hε (z) as However, the convergence rate can be increased because infor-
   mation from larger regions is exploited. Therefore, the total
1 2 z
Hε (z) = 1 + arctan ,z ∈ R (38) computation burden is comparable for different ρ.
2 π ε
where ε > 0 is a constant. The Dirac functional δ(z) is
approximated by δε (z) as A. Comparisons With the CV [1] and GCV [21] Models
d(Hε (z)) 1 ε Both the CV [1] and GCV [21] models assume that image
δε (z) = = , z ∈ R. (39) intensities are piecewise constant and thereby use global inten-
dz π ε + z2
2
sity means to fit them. Therefore, they do not perform well
The procedures of the proposed algorithm are summarized
on images with intensity inhomogeneity. In this section, we
in Algorithm 1.
compare our method with the CV and GCV models on some
real and synthetic images with intensity inhomogeneity.
VI. E XPERIMENTAL R ESULTS
Fig. 4 demonstrates the segmentation results on two syn-
In this section, we compare our method with the CV thetic images and two real-vessel images with intensity inho-
model [1], the GCV model [21],1 the LRB model [9],2 mogeneity by using the CV model [1], the GCV model [21],
the LIC model [11], and the LBF model [10],3 which and our method, respectively. It can be seen that our method
1 http://www.math.ucla.edu/~xbresson/code.html achieves satisfying results because we exploit local image
2 http://www.shawnlankton.com/2007/05/active-contours/ region information, which can better separate each object from
3 http://www.engr.uconn.edu/~cmli/ its background. For the CV and GCV models, they use global
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

8 IEEE TRANSACTIONS ON CYBERNETICS

Fig. 5. Segmentation on a synthetic image (first row) and two real-MR

images (second and third rows) with severe intensity inhomogeneity. From
left to right: segmentation results using our method, the LBF model [10],
the LIC model [11] and the LRB model [9], and the histograms of image
intensity. The blue circles represent the initial contours, and the red lines
represent the final segmentation results. We set ρ = 6 for the first two images
and ρ = 20 for the last image.

Fig. 7. Comparison results on five images preprocessed by bilateral filtering.

Top to bottom: test images preprocessed by bilateral filtering; results of the CV
model [1], the GCV model [21], the LBF model [10], the LIC model [11],
and the LRB model [9], respectively.

method can achieve satisfying segmentation results on them,

Fig. 6. Segmentation results for natural images. The blue rectangles represent
the initial contours, and the red lines represent the final segmentation results.
which demonstrates its favorable generality.
We set ρ = 20 for all test images.

C. Results on Images Preprocessed by Bilateral Filtering

region information and thereby result in severe misclassifica- In Fig. 7, we evaluate the performance of the compet-
tion on these images because there exist severe overlapping ing methods on five synthetic and real-images selected from
intensities between each object and its background. Figs. 4 and 5 in which our method achieves favorable results.
The test images are first preprocessed by bilateral filtering [24]
(see the top row of Fig. 7). Compared with their original
B. Results on Synthetic and Real Images
images in Figs. 4 and 5, it is obvious that intensity contrast
We first use one synthetic image and two real-MR images between object and background is more obvious after prepro-
with severe intensity inhomogeneity to demonstrate the supe- cessing, which in turn results in an easier separation of the
rior performance of our method to the LBF model, the LIC object from background. Therefore, for the three real images,
model and the LRB model, all of which use local means the LBF, LIC and LRB methods achieve better results on the
to fit image intensities and thereby perform well on images preprocessed images than the original ones. However, some
with light intensity inhomogeneity. The LBF model [10] uses small noisy objects also have been extracted, thereby mak-
the local intensity mean to fit the intensities of the measured ing the results less accurate than our method. Furthermore,
images and thereby achieves much better segmentation results the intensities of the preprocessed images are still inhomo-
than the CV and GCV models on images with intensity inho- geneous and thereby the CV and GCV methods cannot work
mogeneity. However, if the intensity inhomogeneity is severe, well on them (see the second and third rows of Fig. 7) because
using only the local mean information may fail to discriminate both the CV and GCV methods are only suitable to segment
the intensities of an object from its background, thereby lead- images with piecewise constant intensities.
ing to inaccurate segmentation. Some experimental results are
shown in the second column of Fig. 5. Similar drawbacks exist
for the LIC model [11] (see the third column of Fig. 5) and D. Application to Simultaneous Segmentation
LRB model [9]. Moreover, as shown in the fourth column of and Bias Correction
Fig. 5, the localized Heaviside and Dirac functionals adopted In this section, we apply our level set method for simul-
by the LRB model make the level set evolution easily fall into taneous segmentation and bias correction, especially for
local minima [1]. Our method can produce much better seg- MR images. We compare our method with the LIC model [11]
mentation results because it takes into account the statistical because only the LIC model is applicable for simultaneous
information of intensities in a transformed domain in which segmentation and bias correction while the other competing
the intensities of object and background are less overlapping models can only be used for segmentation.
than those in the original domain. Then, we test our method on Fig. 8 shows the results by the two competing methods on
three natural images with complicate textures (see Fig. 6). Our synthetic images corrupted with additive Gaussian noises of
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: LEVEL SET APPROACH TO IMAGE SEGMENTATION WITH INTENSITY INHOMOGENEITY 9

Fig. 8. Experiments on a synthetic image with different additive Gaussian

noises. The test images are corrupted with Gaussian white noises of standard
deviations σ = 1 and σ = 5, respectively. From left to right: results by our
method, the LIC model [11], the ground-truth bias field, the estimated bias
field using our method, and the estimated bias field using the LIC model [11].
We set ρ = 10 for the experiments. Fig. 10. Bias-correction results on two 7T MRI images using the proposed
method. From left to right: original images, estimated bias fields, and the
bias-corrected images. We set ρ = 10 for both experiments.

whereas the segmentation result using our method is much

better. This is because our model takes into account the dif-
ferent probability distributions of various objects. Furthermore,
the estimated bias fields by our method are visually much more
similar to the ground truth than the LIC method (see the right
two columns of Fig. 8).
Fig. 9 shows the joint segmentation and bias correction
results on a 3T MRI image, which has four classes of tis-
sues: 1) whiter matter; 2) gray matter; 3) cerebrospinal fluid;
and 4) background. Because one level set function can only
represent two classes of tissues, we need to evolve two level
set functions 1 and 2 according to (33) to represent four
classes of tissues [8]. It can be seen that the tissue segmenta-
tion results by our method (the first and third rows in Fig. 9)
are much more accurate than those by the LIC method (the
second and fourth rows in Fig. 9). It should be noted that it is
very easy to initialize the level set functions in our method in
which the initial contours can be set inside, outside or across
the object boundary. The two initializations in the first and
third rows are very different, but the final segmentation results
are very similar, which demonstrates the robustness of our
method to different initializations.
Finally, Fig. 10 shows the bias-correction results using our
method on two 7T MRI images. This experiment aims to
obtain the bias corrected MRI images. The original images,
estimated bias fields, and the bias-corrected images are shown
in the left, middle, and right columns, respectively. It can be
Fig. 9. From left to right: the initializations of level set functions 1 and
clearly seen that the image quality is significantly improved
2 (the red lines represent the initial zero level sets of 1 , while the blue by our method. Some regions (inside the red circles) whose
lines represent the initial zero level sets of 2 ), estimated bias fields, tissue- intensity contrast is too low to be identified are able to be
classification results, and bias corrected images. The first and third rows are
the results by our method with two different initializations, while the sec-
distinguished clearly after the bias correction.
ond and fourth rows are the results by the LIC model [11] with the same
initializations. We set ρ = 10 for all the experiments.
E. Quantitative Evaluation
different levels. The evaluation criteria are the final segmen- In this section, we use Jaccard similarity (JS) [31], [32]
tation results and the similarity between the estimated bias to quantitatively evaluate the performance of the competing
field and the ground-truth bias. The top-left image in Fig. 8 methods on images whose intensity inhomogeneity has differ-
was added with Gaussian white noise of zero mean and unit ent strength. The JS index between two regions S1 and S2 is
standard deviation. Since the noise level is low, both models calculated as JS(S1 , S2 ) = (|S1 ∩ S2 |/|S1 ∪ S2 |), which is the
yield satisfying segmentation results, while our method out- ratio between the intersectional area of S1 and S2 and their
performs a little the LIC model. However, with increasing the union. Obviously, the closer the JS value is to 1, the more
noise level, the segmentation results by the LIC method are similar S1 is to S2 . In our experiments, S1 is the segmented
becoming very noisy (see bottom-left second figure in Fig. 8), object region, and S2 is the ground-truth.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

10 IEEE TRANSACTIONS ON CYBERNETICS

Fig. 13. Segmentation accuracy by our method with different level set ini-
tializations and region scale parameter ρ. Top: results with different level set
initializations by setting ρ = 10.5. The blue lines represent different initial
Fig. 11. Quantitative comparisons among our method, the CV model [1], contours, and the red lines represent the final segmentation contours. Bottom:
the GCV model [21], the LRB model [9], the LBF model [10], and the LIC the JS values for 20 different initial contours (the first six values are for the
model [11]. Top: test images, where the strength of intensity inhomogeneity six segmentations in the top row), while the right figure shows the JS values
is gradually increased from left to right. Bottom: the corresponding JS values for different region scale parameters ρ. The initial contour is the same as that
yielded by the competing methods on the five images. We set ρ = 10.5 for in the top right-most image.
all the experiments.

Fig. 12. From left to right: segmentation results on image with strong inten-
sity inhomogeneity by our method, the CV model [1], the GCV model [21],
the LRB model [9], the LBF model [10], the LIC model [11], and Canny edge
detector [33] (available from the MATLAB command edge). We set ρ = 10.5
for all the experiments.
Fig. 14. Segmentation results by our method for images with different
types of noises. Top: added Gaussian white noise with zero means and dif-
We test the competing methods on five synthetic images ferent variances (σ = 0.01, 0.02, 0.03, 0.04). Bottom: speckle noises that
are multiplicative noises with zero means and different variances (σ =
with different intensity inhomogeneity. Their corresponding JS 0.01, 0.02, 0.03, 0.04). We use the MATLAB command imnoise to add these
values are shown in the bottom image of Fig. 11. Obviously, noises to a clean image.
the JS values obtained by our method have little difference
for intensity inhomogeneity with different strength, which
demonstrates that our method is very robust to image inten-
sity inhomogeneity. For the CV and GCV models, when the
strength of intensity inhomogeneity is not strong (see images
1 and 2 in Fig. 11), both of them can yield a high JS value.
However, when the strength of inhomogeneity is not low (refer
to images 3–5), the performance of the CV and GCV models
degenerates rapidly. For the LRB models (i.e., the LRB model,
the LBF model and the LIC model), when the strength of
intensity inhomogeneity is strong (refer to images 4 and 5), the
performance of these methods also degenerates severely. This Fig. 15. Illustration of the JS values for images shown by Fig. 14.
is because only using local region means cannot discriminate
object from its background satisfactorily when the intensi- to measure segmentation accuracy. The top row in Fig. 13
ties between them overlaps severely. Our method consistently demonstrates the segmentation results using different initial
yields the highest JS values because it pursues segmentation contours (the blue lines). We can see that there are no obvious
in a transformed domain where the overlapping intensities are visual differences for these segmentation results (the JS val-
effectively suppressed. The final segmentation results for the ues of these results correspond to the first six values shown
top right-most image in Fig. 11 by the six competing meth- in the bottom-left figure of Fig. 12). We adopt 20 different
ods are shown in Fig. 12. The right-most segmentation results initial contours (the region scale parameter is set ρ = 10.5).
are generated by canny edge detector [33]. We can observe From the bottom-left figure of Fig. 13, we can observe that
that double contours near the boundary and some unnecessary the JS values change only from 0.97 to 0.98, which clearly
contours inside the object appear. demonstrate that our method can yield very high and stable
segmentation accuracies with different level set initializations.
F. Robustness to Initializations, Region Scale Parameters, The bottom-right figure of Fig. 13 shows the JS values com-
and Different Types of Noises puted by changing the region scale parameters ρ from 5.5 to
We use the test image in Fig. 12 to demonstrate the 22.5. The high and stable JS values again demonstrate that
robustness of our method to different level set initializations our method is very robust to the region scale parameters in a
and region scale parameters ρ. We again use the JS index wide range.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: LEVEL SET APPROACH TO IMAGE SEGMENTATION WITH INTENSITY INHOMOGENEITY 11

Fig. 14 illustrates the segmentation results for images with A PPENDIX B

different types of noises. The images on the top row of Fig. 14 In (23), we assume that the optimal smoothing function
are with different levels of Gaussian white noises, while the B(x) is B̃(x). Then, we add a variation function η(x) such
images on the bottom row are with different levels of multi- that B(x) = B̃(x) + η(x). Differentiating EL with respect to
plicative noises. We can observe that segmentation results are B and letting  → 0+ , we have
visually satisfying although the noises are very strong, thereby
validating the robust performance of our method to different δEL dEL
= lim
types of noises. Furthermore, as shown by Fig. 15, the JS val- δ B̃ →0+ d
 
 n 
ues on images with multiplicative noises are a little smaller I(y) − B̃(x)ci
than those with Gaussian white noises, which can be attributed = Kρ (x, y)Mi ((y)) ci dyη(x)dx.
i=1
σi2
to the assumption of Gaussian distribution (7) in our method.  
However, the differences (changes from 0.973 to 0.968) of JS (43)
values are so small that the visual effects of the segmentation
Finally, we obtain
results in Fig. 14 are not obvious. n ci
i=1 Kρ ⊗ (IMi ((x))) ·
σi2
B̃(x) =  (44)
VII. C ONCLUSION n c2i
i=1 Kρ ⊗ Mi ((x)) · σi2
In this paper, by exploiting local image region statistics,
we presented a level set method for segmenting images with which is (28).
intensity inhomogeneity. Our method combines the informa-
tion from the neighboring pixels belonging to the same class, R EFERENCES
which equips it with a strong capability to separate the desired [1] T. F. Chan and L. A. Vese, “Active contours without edges,” IEEE Trans.
object from its background. Moreover, the proposed method Image Process., vol. 10, no. 2, pp. 266–277, Feb. 2001.
[2] S. C. Zhu and A. Yuille, “Region competition: Unifying snakes,
yields a soft classification, which can, to some extent, satisfy region growing, and Bayes/MDL for multiband image segmentation,”
the condition of the partial volume effect. In addition, the seg- IEEE Trans. Pattern Anal. Mach. Intell., vol. 18, no. 9, pp. 884–900,
mentation results are insensitive to different initializations of Sep. 1996.
[3] K. Zhang, L. Zhang, H. Song, and W. Zhou, “Active contours with
the level set function, making it useful for automatic appli- selective local or global segmentation: A new formulation and level set
cations. Comparisons with several representative methods on method,” Image Vis. Comput., vol. 28, no. 4, pp. 668–676, 2010.
synthetic and real images have demonstrated the effectiveness [4] B. Wang, X. Gao, D. Tao, and X. Li, “A nonlinear adaptive level set for
image segmentation,” IEEE Trans. Cybern., vol. 44, no. 3, pp. 418–428,
of the proposed algorithm. Mar. 2014.
[5] B. Wang, X. Gao, D. Tao, and X. Li, “A unified tensor level set for image
segmentation,” IEEE Trans. Syst., Man, Cybern. B, Cybern., vol. 40,
A PPENDIX A no. 3, pp. 857–867, Jun. 2010.
[6] S. Balla-Arabé, X. Gao, and B. Wang, “A fast and robust level set method
In (23), we assume that the optimal solution of ci is c̃i , for image segmentation using fuzzy clustering and lattice Boltzmann
and then we add a variation ηi to the variable c̃i such that method,” IEEE Trans. Cybern., vol. 43, no. 3, pp. 910–920, Jun. 2013.
ci = c̃i + ηi . Keeping other variables except for ci fixed, [7] D. Mumford and J. Shah, “Optimal approximations by piecewise smooth
differentiating EL with respect to ci and letting  → 0+ , we
functions and associated variational problems,” Commun. Pure Appl.
Math., vol. 42, no. 5, pp. 577–685, 1989.
have [8] L. A. Vese and T. F. Chan, “A multiphase level set framework for image
segmentation using the Mumford and Shah model,” Int. J. Comput. Vis.,
δEL dEL vol. 50, no. 3, pp. 271–293, 2002.
= lim [9] S. Lankton and A. Tannenbaum, “Localizing region-based active
δc̃i →0+ d
  contours,” IEEE Trans. Image Process., vol. 17, no. 11, pp. 2029–2039,
(I(y)B(x) − B2 (x)c̃i ) Nov. 2008.
= −ηi Kρ (x, y) Mi ((y))dxdy [10] C. Li, C.-Y. Kao, J. C. Gore, and Z. Ding, “Implicit active contours
σi2 driven by local binary fitting energy,” in Proc. IEEE Conf. Comput. Vis.
 
Pattern Recognit., Minneapolis, MN, USA, 2007, pp. 1–7.
= 0. (40) [11] C. Li et al., “A level set method for image segmentation in the pres-
ence of intensity inhomogeneities with application to MRI,” IEEE Trans.
Therefore, we have Image Process., vol. 20, no. 7, pp. 2007–2016, Jul. 2011.
[12] T. Brox and D. Cremers, “On local region models and a statistical inter-
 
I(y)B(x) − B2 (x)c̃i pretation of the piecewise smooth Mumford–Shah functional,” Int. J.
Kρ (x, y) Mi ((y))dxdy = 0. (41) Comput. Vis., vol. 84, no. 2, pp. 184–193, 2009.
σi2 [13] L. Wang et al., “Patch-driven neonatal brain MRI segmentation with
  sparse representation and level sets,” in Proc. IEEE 10th Int. Symp.
Biomed. Imag. (ISBI), San Francisco, CA, USA, 2013, pp. 1090–1093.
From the above equation, we obtain [14] H. Song, B. Huang, Q. Liu, and K. Zhang, “Improving the spatial res-
olution of landsat TM/ETM+ through fusion with SPOT5 images via
(Kρ ⊗ B)IMi ((y))dy learning-based super-resolution,” IEEE Trans. Geosci. Remote Sens.,
 vol. 53, no. 3, pp. 1195–1204, Mar. 2015.
c̃i = (42) [15] H. Song, G. Wang, and K. Zhang, “Hyperspectral image denoising
(Kρ ⊗ B2 )Mi ((y))dy via low-rank matrix recovery,” Remote Sens. Lett., vol. 5, no. 10,
 pp. 872–881, 2014.
[16] H. Song, G. Wang, and K. Zhang, “Multiple change detection for
which is (27). multispectral remote sensing images via joint sparse representation,”
The σ̃i in (29) can be obtained using a similar method. Opt. Eng., vol. 53, no. 12, 2014, Art. ID 123103.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

12 IEEE TRANSACTIONS ON CYBERNETICS

[17] H. Song, “Active contours driven by regularized gradient flux flows for Lei Zhang (SM’14) received the B.Sc. degree from
image segmentation,” Electron. Lett., vol. 50, no. 14, pp. 992–994, 2014. the Shenyang Institute of Aeronautical Engineering,
[18] C. M. Bishop and N. M. Nasrabadi, Pattern Recognition and Machine Shenyang, China, in 1995, and the M.Sc. and
Learning, vol. 1. New York, NY, USA: Springer, 2006. Ph.D. degrees in control theory and engineering
[19] S. Yan, X. Xu, D. Xu, S. Lin, and X. Li, “Image classification with from Northwestern Polytechnical University, Xi’an,
densely sampled image windows and generalized adaptive multiple China, in 1998 and 2001, respectively.
kernel learning,” IEEE Trans. Cybern., vol. 45, no. 3, pp. 395–404, From 2001 to 2002, he was a Research Associate
Mar. 2015. at the Department of Computing, Hong Kong
[20] K. Zhang, L. Zhang, and S. Zhang, “A variational multiphase level set Polytechnic University, Hong Kong. From 2003
approach to simultaneous segmentation and bias correction,” in Proc. to 2006, he was a Post-Doctoral Fellow at the
IEEE Int. Conf. Image Process., Hong Kong, 2010, pp. 4105–4108. Department of Electrical and Computer Engineering,
[21] X. Bresson, S. Esedoglu, P. Vandergheynst, J.-P. Thiran, and S. Osher, McMaster University, Hamilton, ON, Canada. In 2006, he joined the
“Fast global minimization of the active contour/snake model,” J. Math. Department of Computing, Hong Kong Polytechnic University, as an Assistant
Imag. Vis., vol. 28, no. 2, pp. 151–167, 2007. Professor, where he has been an Associate Professor since 2010. His current
[22] T. Brox, “From pixels to regions: Partial differential equations in image research interests include image and video processing, computer vision, pat-
analysis,” Ph.D. dissertation, Dept. Comput. Sci., Saarland University, tern recognition, and biometrics. He has published about 200 papers in the
Saarbrücken, Germany, 2005. above areas.
[23] K. Zhang, Q. Liu, H. Song, and X. Li, “A variational approach to simul- Dr. Zhang was a recipient of the 2012–2013 Faculty Award in Research
taneous image segmentation and bias correction,” IEEE Trans. Cybern., and Scholarly Activities. He is currently an Associate Editor of the IEEE
Oct. 2014. T RANSACTIONS ON C IRCUITS AND S YSTEMS FOR V IDEO T ECHNOLOGY
[24] C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color and Image and Vision Computing.
images,” in Proc. 6th Int. Conf. Comput. Vis. Pattern Recognit., Bombay,
India, 1998, pp. 839–846.
[25] M. Styner, C. Brechbuhler, G. Szckely, and G. Gerig, “Parametric esti-
mate of intensity inhomogeneities applied to MRI,” IEEE Trans. Med.
Imag., vol. 19, no. 3, pp. 153–165, Mar. 2000.
[26] L. Wang, L. He, A. Mishra, and C. Li, “Active contours driven by local
Gaussian distribution fitting energy,” Signal Process., vol. 89, no. 12,
pp. 2435–2447, 2009.
[27] H. Knutsson and C.-F. Westin, “Normalized and differential convolu-
tion,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., Kin-Man Lam (SM’14) received the Associateship in electronic engineering
New York, NY, USA, 1993, pp. 515–523. (with distinction) from Hong Kong Polytechnic University (formerly called
[28] C. Li, C. Xu, C. Gui, and M. D. Fox, “Distance regularized level set Hong Kong Polytechnic), Hong Kong, the M.Sc. degree in communication
evolution and its application to image segmentation,” IEEE Trans. Image engineering from the Department of Electrical Engineering, Imperial College
Process., vol. 19, no. 12, pp. 3243–3254, Dec. 2010. of Science, Technology, and Medicine, London, U.K., and the Ph.D. degree
[29] K. Zhang, L. Zhang, H. Song, and D. Zhang, “Reinitialization-free from the Department of Electrical Engineering, University of Sydney, Sydney,
level set evolution via reaction diffusion,” IEEE Trans. Image Process., NSW, Australia, in 1986, 1987, and 1996, respectively.
vol. 22, no. 1, pp. 258–271, Jan. 2013. From 1990 to 1993, he was a Lecturer with the Department of Electronic
[30] K. Zhang, H. Song, and L. Zhang, “Active contours driven by local Engineering, Hong Kong Polytechnic University. He joined the Department of
image fitting energy,” Pattern Recognit., vol. 43, no. 4, pp. 1199–1206, Electronic and Information Engineering, Hong Kong Polytechnic University,
2010. as an Assistant Professor in 1996 and became an Associate Professor in 1999,
[31] K. Zhang, L. Zhang, and M.-H. Yang, “Fast compressive tracking,” where he has been a Professor since 2010. His current research interests
IEEE Trans. Pattern Anal. Mach. Intell., vol. 36, no. 10, pp. 2002–2015, include human face recognition, image and video processing, and computer
Oct. 2014. vision.
[32] H. Song, “Robust visual tracking via online informative feature selec- Dr. Lam is currently a General Co-Chair of the IEEE International
tion,” Electron. Lett., vol. 50, no. 25, pp. 1931–1933, 2014. Conference on Signal Processing, Communications, and Computing. He also
[33] J. Canny, “A computational approach to edge detection,” IEEE Trans. serves as an Associate Editor of the IEEE T RANSACTIONS ON I MAGE
Pattern Anal. Mach. Intell., vol. PAMI-8, no. 6, pp. 679–698, Nov. 1986. P ROCESSING, the Asia Pacific Signal and Information Processing Association
(APSIPA) Transactions on Signal and Information Processing, and the
EURASIP International Journal on Image and Video Processing. He was a
member of the organizing and program committees of several international
conferences. He is a BoG Member of the APSIPA and the Director-Student
Services of the IEEE Signal Processing Society.

Kaihua Zhang received the B.S. degree in tech- David Zhang (F’08) received the B.Sc. degree in computer science from
nology and science of electronic information from Peking University, Beijing, China, the M.Sc. degree in computer science from
the Ocean University of China, Qingdao, China, the the Harbin Institute of Technology, Harbin, China, in 1982, and the Ph.D.
M.S. degree in signal and information processing degrees in computer science and electrical and computer engineering from
from the University of Science and Technology of the Harbin Institute of Technology and the University of Waterloo, Waterloo,
China, Hefei, China, and the Ph.D. degree from the ON, Canada, in 1985 and 1994, respectively.
Department of Computing, Hong Kong Polytechnic From 1986 to 1988, he was a Post-Doctoral Fellow at Tsinghua University,
University, Hong Kong, in 2006, 2009, and 2013, Beijing, and then an Associate Professor at Academia Sinica, Beijing. He is
respectively. currently a Chair Professor with Hong Kong Polytechnic University, Hong
He is a Professor with the School of Information Kong, where he is the Founding Director of the Biometrics Technology Centre
and Control, Nanjing University of Information (UGC/CRC) supported by the Hong Kong Special Administrative Region
Science and Technology, Nanjing, China. From 2009 to 2010, he was a Government in 1998.
Research Assistant at the Department of Computing, Hong Kong Polytechnic Prof. Zhang is a Croucher Senior Research Fellow, the International
University. His current research interests include image segmentation, level Association of Pattern Recognition Fellow, and a Distinguished Speaker of
sets, and visual tracking. the IEEE Computer Society.