Adaptive Regularized Level Set Method for Weak Boundary Object Segmentation

Abstract

An adaptive regularized level set method for image segmentation is proposed. A weighted 𝑝(𝑥)-Dirichlet integral is presented as a geometric regularization on zero level curve, which is used to diminish the influence of image noise on level set evolution while ensuring the active contours not to pass through weak object boundaries. The idea behind the new energy integral is that the amount of regularization on the zero level curve can be adjusted automatically by the variable exponent 𝑝(𝑥) to fit the image data. This energy is then incorporated into a level set formulation with an external energy term that drives the motion of the zero level set toward the desired objects boundaries, and a level set function regularization term that is necessary for maintaining stable level set evolution. The proposed model has been applied to a wide range of both real and synthetic images with promising results.

1. Introduction

Image segmentation is a key initial step before performing high-level tasks such as objects recognition and tracking [1, 2] in most computer vision applications. Until now, image segmentation is yet a difficult task although it has been studied extensively in past decades. Typical difficulties result from facts that most natural images include noise, low intensity contrast with weak edges, and intensity inhomogeneity. A number of segmentation techniques are developed to overcome these difficulties, in which the level set methods [3–6] have been proved to be a class of efficient techniques.

The level set method, originally introduced by Osher and Sethian [7], is a general framework for the computation of evolving interfaces using implicit representations. In image processing and computer vision applications, geometric active contours (GAC) [3–5], that is, active contours implemented via level set methods [7], were first introduced independently by Caselles et al. [3] and Malladi et al. [4] for image segmentation. Afterwards, many works constitute very interesting applications of the level set method within the active contour framework [8–10]. The basic idea is to represent an active contour as a zero level set of an implicit function in a higher dimension, called level set function, and to deform the level set function according to an evolution partial differential equation (PDE). They are appealing for the ability to handle topological changes automatically, which is generally impossible in the traditional parametric active contours [11]. The evolution PDE for the level set function can be derived from the related evolution equation of a parameterized contour [4–6]. Alternatively, it can be directly derived from the minimization problem for an energy functional defined on the level set function [8–10, 12]. This type of methods is known as variational level set (VLS) methods. There are several desirable advantages of the VLS methods over classical image segmentation methods, such as thresholding and region grow. First, the VLS methods allow incorporation of various prior knowledge, such as shape and intensity distribution, under a principled energy minimization framework. Second, they can provide closed and continuous contours for further applications, such as image analysis and object tracking. Third, they are amenable to the introduction of constraints on level set function; smoothnesses of level set function are typical constraints. Generally, the constraints within a model can be categorized into two parts: level set function regularization and zero level curve regularization.

In conventional level set methods [3–6], the level set function may develop shocks during the evolution, which cause numerical errors and eventually destroies the stability of the level set evolution. To overcome this difficulty, level set function regularization, commonly known as reinitialization [13, 14], was introduced to restore the regularity of the level set function and maintain stable level set evolution. Although re-initialization as a numerical remedy is able to maintain the regularity of the level set function, there are serious theoretical and practical problems, as pointed out by Gomes and Faugeras [15]. Recently, Li et al. [16, 17] proposed a variational level set formulation with an intrinsic mechanism of maintaining the signed distance property of the level set function.

A problem for GAC models is the extraction of weak boundaries in noise and/or intensity inhomogeneity images. As we know, Gaussian kernel filter is a valid denoising method which smoothes away the isolated points. For example, Zhang et al. [18] proposed a local image fitting (LIF) energy based on Gaussian filtering for variational level set to regularize the level set function. Wang et al. [19] proposed a region-based tensor level set model for image segmentation, in which a three-order tensor involving the Gaussian filter bank was introduced to comprehensively depict features of pixels. These models are robust to noise. However, the level set evolution only depending on Gaussian kernel filter cannot achieve an accurate separation between weak objects and noise in complex images since the Gaussian kernel filter often eliminates weak boundaries or details. Some regularity must be imposed on the zero level curve to diminish the influence of image noise on level set evolution.

In this paper, we focus on the issue of zero level curve regularization with variational method. Length regularization [6, 8, 9, 18] is an initially popular choice of the geometric constraint on zero level curve, in the spirit of the Mumford-Shah functional [20]. But it is less robustness to noise. In [21], two smoother regularizations ∫_Ω|∇𝜙|²𝑑𝑥𝑑𝑦 and ‖|∇𝜙|‖_𝐿∞(Ω) were introduced. However, the smoother regularization may cause the active contours to pass through weak objects boundaries. Recently, there were many different choices of regularization in the literature, for example, p-Dirichlet integral regularization [22] and weighted p-Dirichlet integral regularization [23]. Different value of 𝑝≥1 results in a constrain which is somewhere between length-based and smoother regularization. However, the constant exponent 𝑝 cannot reflect the local property of image, thus the p-regularization do not adapt the exponent to fit the data automatically. This problem has limited their application.

This paper proposes an adaptive variational level set formulation with a weighted 𝑝(𝑥)-Dirichlet integral, an external energy, and a level set regularization term. The weighted 𝑝(𝑥)-Dirichlet integral integrating the gradient information is designed as a geometric regularization on zero level curve, which is used to diminish the influence of image noise on level set evolution while ensuring the active contours not to pass through weak object boundaries. To demonstrate the effectiveness of the weighted 𝑝(𝑥)-Dirichlet integral, we apply it to an edge-based GAC model for image segmentation. So an external energy based on Laplacian of Gaussian (LoG) filter is defined, and then it drives the level set function to deform in opposite direction (up or down) on either side of edge. The level set regularization term makes level set function behave approximately like a signed distance function, which ensures stable level set evolution. The resulting evolution of the level set function is the gradient flow that minimizes the overall energy functional. Due to the image data fitting in the weighted 𝑝(𝑥)-Dirichlet integral term, intensity information in local regions is extracted to guide the regularization of contours. Thereby our model can extract weak boundaries in noisy and/or intensity inhomogeneity images. An added benefit of the proposed model is that the level set function can be initialized to a constant function. The constant function is more easier to use in practice than the widely used signed distance function or binary step function.

The remainder of this paper is organized as follows. In Section 2, we simply review level set method and some typical regularizations on zero level set curve. The proposed model is introduced in Section 3. Numerical algorithm and experimental results are presented in Section 4. This paper is summarized in Section 5.

2. Background

2.1. Level Set Method and Level Set Function Regularization

In the level set formulation, a moving curve 𝐶(𝑡) is represented by the zero level set of a Lipschitz function 𝜙(𝑥,𝑦,𝑡) defined on the entire image domain. The evolution of the curve 𝐶(𝑡) along its normal direction with speed 𝐹 is described by the following evolution PDE [7]:𝜕𝜙𝜕𝑡+𝐹∇𝜙=0 (2.1) with the initial condition 𝜙(𝑥,𝑦,0)=𝜙₀(𝑥,𝑦). For image segmentation, the speed function 𝐹 depends on both image data and the level set function 𝜙.

The function 𝜙 in (2.1) may develop shocks during the evolution. As a result, some regularities must be imposed on 𝜙 in order to prevent 𝜙 to be too steep or too fat near the zero level curve. A common means is to initialize and periodically reinitialize the level set function to a signed distance function so as to keep steady level set evolution and ensure usable results. The reinitialization equation is𝜕𝜙𝜙,𝜕𝑡=sign1−∇𝜙 (2.2) where 𝜙 is the function to be re-initialized, and sign (⋅) is the sign function. Although re-initialization as a numerical remedy is able to maintain the regularity of the level set function, it still remains a fundamental problem as when and how to apply the re-initialization [14, 15].

In [16], Li et al. introduced an internal energy to keep the regularity of the level set function. It is formulated in a variational framework on the level set function 𝜙 as follows:𝐸𝜙=𝜇𝑃𝜙+𝐸_ext𝜙, (2.3) where 𝑃(𝜙) is the internal energy and 𝐸_ext(𝜙) is a certain external energy that would drive the motion of the zero level curve of 𝜙.

Let Ω be the image domain, the internal energy term is defined as𝑃𝜙=12_Ω∇𝜙−1²𝑑𝑥𝑑𝑦 (2.4) which makes the level set function 𝜙 behave approximately like a signed distance function, and so eliminates the re-initialization step during level set evolution. With this internal energy, level set evolution can be implemented by simple finite difference scheme and is initialized to a binary step function [8, 9, 16–18].

2.2. Some Typical Regularizations on Zero Level Curve

In order to control the smoothness of the zero level curve and further avoid the occurrence of small, isolated regions in the final segmentation, the regularization on zero level curve is very crucial in level set methods. The length regularization [6, 8, 9, 18] is to minimize the following energy functional:𝐿𝜙=_Ω∇𝐻𝜙𝑑𝑥𝑑𝑦=_Ω𝛿𝜙∇𝜙𝑑𝑥𝑑𝑦, (2.5) where 𝐻(⋅) is Heaviside function and 𝛿(⋅) is Dirac delta function. The energy functional 𝐿(𝜙) in (2.5) computes the length of the zero level curve of 𝜙 in the conformal metric 𝑑𝑠=|𝐶′(𝑝)|𝑑𝑝. The length regularization imposes a penalty on the length of the curve that smoothes the zero level curve and diminishes some false contours. But the smoothing is only along the tangent direction of each level line, so this regularization is less robustness to noise.

In [21],Chung and Veseproposed a smoother regularization:𝑅𝜙=_Ω∇𝜙²𝑑𝑥𝑑𝑦. (2.6) This regularization means isotropic smoothing at every point (𝑥,𝑦), so the level lines tend to maintain smoothing and further penalize the false contours in noise image. But the smoother regularization cannot stop the active contours to pass through some weak object boundary.

In [23],Zhou and Muproposed a weighted p-Dirichlet integral regularized level set evolution. The geometric regularization has the following form:𝐿_𝑝𝜙=_Ω𝛿𝜙∇𝜙^𝑝𝑑𝑥𝑑𝑦. (2.7) Different value of 𝑝≥1 results in a tradeoff between length regularization and smoother regularization. However, if the image intensities representing objects are nonuniform or if an image is highly degraded, this regularization may become sensitive to exponent 𝑝.

3. Weighted 𝑝(𝑥)-Dirichlet Integral Regularized Level Set Evolution

In this section, we propose a new variational level set formulation for image segmentation, in which a weighted 𝑝(𝑥)-Dirichlet integral is presented to regularize the zero level curve.

Let Ω⊂𝑅² be a image domain. For a given image 𝐼∶Ω→𝑅 and a level set function 𝜙(𝑥,𝑦)∶Ω→𝑅, we define an energy functional 𝐸(𝜙) by𝐸𝜙=𝐿_𝑝(⋅)𝜙+𝜈𝐸_ext𝜙+𝜇𝑃𝜙, (3.1) where 𝜈,𝜇>0 are constants, 𝐿_𝑝(⋅)(𝜙) is the zero level curve regularization term, 𝐸_ext(𝜙) is a certain external energy that would drive the motion of the zero level curve of 𝜙, and 𝑃(𝜙) is the level set function regularization term that controls the smoothness of the level set function during the level set evolution.

3.1. Weighted 𝑝(𝑥)-Dirichlet Integral

The zero level curve regularization term 𝐿_𝑝(⋅)(𝜙) is defined as follows:𝐿_𝑝(⋅)𝜙=_Ω1𝑝∇𝐺_𝜎𝛿𝜙∗𝐼∇𝜙^{𝑝(|∇𝐺𝜎∗𝐼|)}𝑑𝑥𝑑𝑦, (3.2) where ∇ is gradient operator and 𝐺_𝜎∗𝐼 is the convolution of the image 𝐼 with the Gaussian function 𝐺_𝜎 with standard deviation 𝜎. The exponent 𝑝(𝑠)∶[0,+∞)→[1,1.5] is a monotonically increasing function with lim_𝑠→0𝑝(𝑠)=1 and lim_𝑠→+∞𝑝(𝑠)=1.5. A simple example is𝑝(𝑠)=32−12+𝑠. (3.3)

The functional (3.2) is in fact a weighted 𝑝-Dirichlet integral with variable exponent 𝑝(|∇𝐺_𝜎∗𝐼|), thus it is called the weighted 𝑝(𝑥)-Dirichlet integral in this paper.

Remarks. (1) If 𝑝(|∇𝐺_𝜎∗𝐼|)=1, the functional 𝐿_𝑝(⋅)(𝜙) is simplified to 𝐿_𝑝(⋅)𝜙=_Ω𝛿𝜙∇𝜙𝑑𝑥𝑑𝑦. (3.4)It is well known that the energy functional computes the length of the zero level curve of 𝜙.(2)When 𝑝(|∇𝐺_𝜎∗𝐼|)=𝑝>1 is a constant, the functional 𝐿_𝑝(⋅)(𝜙) is the weighted 𝑝-Dirichlet integral. It is chosen as one of geometric regularization on zero level curve in [23].
However, the constant exponent 𝑝 is not an intelligent choice. An important reason is that the constant exponent 𝑝 does not reflect the local property of image, and so does not adapts the exponent to fit the image data automatically. We let 𝑝=𝑝(|∇𝐺_𝜎∗𝐼|) depend on local intensity information. This way the regularization ensures weak regularization (𝑝(|∇𝐺_𝜎∗𝐼|)≈1) in regions where the image almost has a constant intensity (i.e., where the intensity gradient almost is zero) to avoid the disappearance of weak boundaries, while it ensures strong regularization (𝑝(|∇𝐺_𝜎∗𝐼|)≈1.5) in other regions to force the false contours to vanish in final segmentation. We will show this by a simple experiment in Section 3.4. Therefore, the variable exponent 𝑝(𝑥) controls the tradeoff between weak regularization and strong regularization automatically.

3.2. External Energy Term 𝐸_ext(𝜙)

In image segmentation, an external energy depending on image information must be defined to move the zero level curve toward the objects boundaries. The weighted 𝑝(𝑥)-Dirichlet integral in (3.2) can be used in various applications with different definitions of the external energy 𝐸_ext(𝜙). In this subsection, we only provide an application of the weighted 𝑝(𝑥)-Dirichlet integral to an active contour model using edge-based information in the external energy, as a demonstration of the effectiveness of the weighted 𝑝(𝑥)-Dirichlet integral formulation.

Next, we define the external energy 𝐸_ext(𝜙) based on the Laplacian of a Gaussian (LoG) filter as follows, which drives the level set function to deform in opposite direction (up or down) on either side of edge:𝐸_ext𝜙=_ΩΔ𝐺_𝜎∗𝐼⋅𝐻−𝜙𝑑𝑥𝑑𝑦, (3.5) where Δ𝐺_𝜎 is the LoG filter. The LoG filter calculates the second derivative of an image, which is often used for zero crossing edge detectors. It is well known that at the inflection point the second derivative vanishes and changes sign. The LoG response is zero in areas where the image has a constant intensity. In the vicinity of a change in intensity, the LoG response is positive on the darker side and negative on the lighter side. By incorporating the edge-based information (the LoG filter) into the external energy term, the level set function can move up or down on either side of edge and cause the sign of 𝜙 to flip around edges. And the objects boundaries can be extracted at image locations where two opposite directions of flow encounter. In this formulation, the level set function can be initialized to a constant function. Such initialization scheme is quite simple and computationally efficient in practice application.

3.3. The Energy Formulation

In order to keep the regularity of level set function, the energy 𝑃(𝜙) in (2.4) is adopted in our model. This functional makes our level set formulation that has an intrinsic mechanism of maintaining the desired shape of 𝜙.

With the energy functionals (3.2),(3.5), and (2.4), the proposed energy formulation (3.1) can be rewritten as𝐸𝜙=_Ω1𝑝∇𝐺_𝜎𝛿𝜙∗𝐼∇𝜙^{𝑝(|∇𝐺𝜎∗𝐼|)}𝑑𝑥𝑑𝑦+𝑣_ΩΔ𝐺_𝜎1∗𝐼⋅𝐻−𝜙𝑑𝑥𝑑𝑦+2𝜇_Ω∇𝜙−1²𝑑𝑥𝑑𝑦. (3.6)

In order to compute the associated Euler-Lagrange equation for the unknown function 𝜙, we consider slightly regularized versions of the function 𝐻 and 𝛿(𝜙), denoted here by 𝐻_𝜀 and 𝛿_𝜀, as 𝜀→0. Let 𝐻_𝜀 be a 𝐶^∞(Ω) regularization of 𝐻, and 𝛿_𝜀=𝐻^′_𝜀. In a dynamical scheme via steepest descent, minimizing the energy functional (3.6) with respect to 𝜙, we obtain the evolution PDE:𝜕𝜙𝜕𝑡=𝛿_𝜀𝜙div∇𝜙^{𝑝(|∇𝐺𝜎∗𝐼|)−2}−1∇𝜙𝑝∇𝐺_𝜎𝛿∗𝐼−1^′_𝜀𝜙∇𝜙^{𝑝(|∇𝐺𝜎∗𝐼|)}+𝜈𝛿_𝜀𝜙Δ𝐺_𝜎∗𝐼+𝜇Δ𝜙−div∇𝜙∇𝜙 (3.7) with the initial condition 𝜙(𝑥,𝑦,0)=𝜙₀(𝑥,𝑦).

In all numerical experiments, we choose the following functions:𝛿_𝜀(𝑠)=1𝜋𝜀𝜀²+𝑠², (3.8) for our evolution PDE (3.7).

3.4. Weighted 𝑝(𝑋)-Dirichlet Integral Effect

In this subsection, we demonstrate the weighted 𝑝(𝑥)-Dirichlet integral (3.2) effect by a simple experiment, just as shown in Figure 1. We will see that the variable exponent 𝑝(|∇𝐺_𝜎∗𝐼|) in 𝐿_𝑝(⋅)(𝜙) effectively controls the tradeoff between the weak regularization and the strong regularization automatically. The test images are an aerial image (128 × 128) with cognitive boundaries (i.e., discontinuous boundaries) and a blood vessel image (103 × 131) with intensity inhomogeneity, as shown in Column 1. In this experiment, we take 𝜎=2.0, 𝜀=1.5, 𝜇=0.04. The level set function is initialized to 𝜙₀(𝑥,𝑦)=0.5 and the level set function 𝜙 evolves according to (3.7) with different exponents 𝑝(|∇𝐺_𝜎∗𝐼|). We observe from Figures 1(b) and 1(f) (where 𝑝(|∇𝐺_𝜎∗𝐼|)=1, length regularization) that some false contours appear in final segmentation (see the interior of object in Figure 1(b) and the upper left and right corners in Figure 1(f)). We can see from Figures 1(c) and 1(g) (where 𝑝(|∇𝐺_𝜎∗𝐼|)=1.5, weighted 𝑝-Dirichlet integral regularization) that some weak edges cannot be extracted effectively although the false contours are suppressed (see the middle of Figure 1(c) and the bottom of Figure 1(g)). With the variable exponent 𝑝(|∇𝐺_𝜎∗𝐼|)=1.5−(2+|∇𝐺_𝜎∗𝐼|)⁻¹(weighted 𝑝(𝑥)-Dirichlet integral regularization), the proposed model not only successfully extracted the weak boundaries, but also effectively suppressed the false contours, as shown in Figures 1(d) and 1(h).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 1 

Weighted 𝑝(𝑥)-Dirichlet integral effect on real images (top: 𝜈=0.25, bottom: 𝜈=0.5). Column 1: original images. Column 2: length regularization. Column 3: weighted 𝑝-Dirichlet integral regularization. Column 4: weighted 𝑝(𝑥)-Dirichlet integral regularization.

4. Numerical Algorithm and Experimental Results

4.1. Numerical Algorithm

In this subsection, we briefly present the numerical algorithm and procedure to solve the evolution (3.7). Equation (3.7) can be implemented via a simple explicit finite difference scheme as in [16, 17] rather than a complex upwind scheme as in the traditional level set methods [14]. We consider the 2D case with a time-dependent level set function 𝜙(𝑥,𝑦,𝑡). The spatial derivatives 𝜕𝜙/𝜕𝑥 and 𝜕𝜙/𝜕𝑦 in our model are approximated by the central difference, and the temporal partial derivative 𝜕𝜙/𝜕𝑡 is discretized as the forward difference.

We recall first the usual notations: let Δ𝑡 be the time step, ℎ be the space step, and (𝑥_𝑖,𝑦_𝑖)=(𝑖ℎ,𝑗ℎ) be the grid points. Let 𝜙^𝑛_𝑖,𝑗=𝜙(𝑥_𝑖,𝑦_𝑗,𝑛Δ𝑡) be an approximation of 𝜙(𝑥,𝑦,𝑡), with 𝑛≥0, 𝜙⁰=𝜙₀, the central differences areΔ^𝑥𝜙_𝑖,𝑗=𝜙_𝑖+1,j−𝜙_{𝑖−1,𝑗}2ℎ,Δ^𝑦𝜙_𝑖,𝑗=𝜙_𝑖,𝑗+1−𝜙_{𝑖,𝑗−1}2ℎ. (4.1) Set 𝑛=0, and start with initial level set function 𝜙⁰_𝑖,𝑗, we first compute 𝑝(|∇𝐺_𝜎∗𝐼|) according to (3.3). Then, the numerical approximation to (3.7) is given by the following discretization:𝜙^𝑛+1_𝑖,𝑗−𝜙^𝑛_𝑖,𝑗Δ𝑡=𝛿_𝜀𝜙^𝑛_𝑖,𝑗𝑙^𝑛_𝑖,𝑗−1𝑝∇𝐺_𝜎𝛿∗𝐼−1^′_𝜀𝜙^𝑛_𝑖,𝑗Δ^𝑥𝜙^𝑛_𝑖,𝑗²+Δ^𝑦𝜙^𝑛_𝑖,𝑗^{2𝑝(|∇𝐺𝜎∗𝐼|)/2}+𝜈𝛿_𝜀𝜙^𝑛_𝑖,𝑗Δ^𝑥Δ^𝑥𝐺_𝜎∗𝐼+Δ^𝑦Δ^𝑦𝐺_𝜎Δ∗𝐼+𝜇^𝑥Δ^𝑥𝜙^𝑛_𝑖,𝑗+Δ^𝑦Δ^𝑦𝜙^𝑛_𝑖,𝑗−𝑘^𝑛_𝑖,𝑗, (4.2) where𝑙^𝑛_𝑖,𝑗=Δ^𝑥Δ^𝑥𝜙^𝑛_𝑖,𝑗Δ^𝑥𝜙^𝑛_𝑖,𝑗²+Δ^𝑦𝜙^𝑛_𝑖,𝑗^{22−𝑝(|∇𝐼𝜎|)}+Δ^𝑦Δ^𝑦𝜙^𝑛_𝑖,𝑗Δ^𝑥𝜙^𝑛_𝑖,𝑗²+Δ^𝑦𝜙^𝑛_𝑖,𝑗^{22−𝑝(|∇𝐼𝜎|)}𝑘_{𝑛𝑖,𝑗},=Δ^𝑥Δ^𝑥𝜙^𝑛_𝑖,𝑗Δ^𝑥𝜙^𝑛_𝑖,𝑗²+Δ^𝑦𝜙^𝑛_𝑖,𝑗²+Δ^𝑦Δ^𝑦𝜙^𝑛_𝑖,𝑗Δ^𝑥𝜙^𝑛_𝑖,𝑗²+Δ^𝑦𝜙^𝑛_𝑖,𝑗². (4.3)(4.4) This system (4.2) is solved by an iterative method, and we use fixed space step ℎ=1 (implies pixel spacing) for computing spatial derivatives. The choice of the time step for this finite different scheme must satisfy the condition 𝜇Δ𝑡<0.25 for numerical stability, and for more details, we refer the reader to [17].

It is worth noting that the spatial derivatives in conventional level set formulation are discretized by upwind scheme [14]. Due to the added the level set regularization term in our model, Central difference scheme is valid for the PDE (3.7). all the spatial derivatives in (3.7) can be discretized by central difference scheme, and the corresponding numerical scheme is stable without the need for re-initialization. Moreover, the central difference scheme is more accurate and efficient than the first-order upwind scheme that is commonly used in conventional level set formulation, as pointed out by Li et al. in [17].

In summary, the main steps of the algorithm (3.7) are as follows.

(1) Initialize the level set function 𝜙₀= constant, set 𝑛=0.

(2) Compute 𝑝(|∇𝐺_𝜎∗𝐼|) according to (3.3).

(3) Solve the PDE in 𝜙 from (4.2), to obtain 𝜙^𝑛+1.

(4) Check whether the evolution has stationary. If not, 𝑛=𝑛+1 and repeat.

4.2. Experimental Results

This section shows the results of the proposed model for both synthetic and real images. The level set function 𝜙(𝑥,𝑦,𝑡) is simply initialized to 𝜙₀(𝑥,𝑦)=𝜌 (nonzero constant); we choose |𝜌|=0.5 for all experiments. Besides, we use the following default parameter setting for all experiments: 𝜎=2.0, 𝜀=1.5, 𝜇=0.04, Δ𝑡=5.0. The parameter 𝜈 is not the same in different experiments; we will give the exact value of 𝜈 in each experiment, together with the sign of 𝜌. For pixels on the borders of the image 𝐼, we take a mirror reflection in all experiments.

Figure 2 shows the segmentation process of the proposed model on a synthetic image (88 × 85), a T-shaped image (127 × 96), an X-ray image of vessels (176 × 167), and a hysterosalpingography (HSG) image (130 × 96) which are plotted in Column 1 of Figure 2. All of them are typical images with intensity inhomogeneity. In these images, parts of the objects boundaries are quite weak, which make it a nontrivial task to extract the objects from the background. Since the level set function is computed from the initial condition 𝜙₀(𝑥,𝑦)=𝜌 (top to bottom: 𝜌 = 0.5, −0.5, 0.5, and 0.5, resp.), there are no initial contours (zero-level curve) superimposed on the original images. The contours (zero-level curve) evolution processes are shown in Column 2 to Column 5. We see from Column 2 that the contours (zero-level curve) emerge automatically in a single iteration due to the introduction of the external energy term. And then the generated contours evolve toward objects boundaries, as shown in Column 3. Afterwards the evolving curve continues to propagate and the generated false small contours gradually disappear (see Column 4). This shrinking effect can be interpreted as the regularization of the weighted 𝑝(𝑥)-Dirichlet integral. Finally, the evolving curves convergence to the objects boundaries (see Column 5). These results show that our model, starting with a constant function, can detect weak boundary objects in intensity inhomogeneity images.

(a)

(b) 1iteration

(c) 15iteration

(d) 60iteration

(e) 130iteration

(f)

(g) 1iteration

(h) 15iteration

(i) 60iteration

(j) 160iteration

(k)

(l) 1iteration

(m) 5iteration

(n) 7iteration

(o) 15iteration

(p)

(q) 1iteration

(r) 2iteration

(s) 10iteration

(t) 70iteration

Figure 2 

Segmentation processes of our model (top to bottom: 𝜈 = 0.27, 0.46, 0.6, 0.46). Column 1: original images. Column 2 to 5: the curve evolution process from the automatically generated contours to the final contours.

Figure 3 shows that our model starting with a constant function can work well for images with multiple weak objects and is compared with the RSF model [8] and LIF model [18]. The RSF model and the LIF model represent the state of the art of variational level set methods which are able to handle intensity inhomogeneity efficiently and are less sensitive to noise. The codes of the RSF model and the LIF model are cited from http://www.engr.uconn.edu/~cmli/ and http://www4.comp.polyu.edu.hk/~cslzhang/papers.htm, respectively. Two test images, which are shown in Column 1, are a DNA channel image (229 × 168) and a bacteria image (173 × 173). These images intensities representing objects are typical nonuniform. To make fair comparison, we try our best to choose the optimal parameters and initial contours for RSF model. We can see from Column 2 and Column 3 that some unwanted contours were generated in the results of the RSF model and the LIF model. Our level set evolution starts with 𝜙₀(𝑥,𝑦)=𝜌 (top: 𝜌 = 0.5, bottom: 𝜌 = −0.5). Column 4 shows that our model successfully extracts all objects contours from background.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3 

Comparisons of the RSF model [8], the LIF model [18] and our model for images with multiple weak objects. Column 1: original images. Column 2: results of the RSF model. Column 3: results of the LIF model. Column 4: results of our model (top to bottom: 𝜈= 0.46, 0.15).

Figure 4 shows the robustness of our model to noise. For this purpose, we create four images by adding Gaussian noise to a synthetic dragon-like image (128 × 128) and a vascular biopsy image (94 × 123), as shown in Row 1. Then, the RSF model, the LIF model, and our model are applied to these original and noisy images. We observe that the three models are able to segment the objects in the original images (see Column 1 and Column 4). But for the noisy images, the contours of the RSF model collapse due to boundary leaks (see Row 2). Although some false contours are generated in the results of the LIF model, it still work by and large (see Row 3). We can see from Row 4 that our method successfully extracts all objects correctly in these noisy images (see Row 4).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

(v)

(w)

(x)

Figure 4 

Segmentation results on images with and without (Gaussian) noise. Row 1: original and noise images ((a) and (d) original images; (b) noisy image (zero mean, 0.04 variance); (c) noisy image (zero mean, 0.1 variance); (e) noisy images (zero mean, 0.01 variance); (f) noisy image (zero mean, 0.04 variance)). Row 2: results of the RSF model [8]. Row 3: results of the LIF model [18]. Row 4: results of our model (left to right: 𝜈 = 0.5, 0.28, 0.2, 0.4, 0.25, 0.21).

Figure 5 shows the comparison of the proposed model with the LIF model on bimodal images. The goal is to show the accuracy of our model. Five test images, which are shown in Row 1, are a synthetic image (84 × 84), wrench image (100 × 100), hand image (108 × 130), plane image (135 × 125), and rice image (128 × 128), respectively. The true objects can be extracted from the original images by a thresholding algorithm (see Row 2). The segmentation results obtained by the LIF model and our model are shown in Rows 3 and 4. It can be observed that the LIF model and our model have achieved similar final results for the first three images by visual comparison. For the plane image, the LIF model incorrectly identifies parts of the plane projection as the objects (see Figure 5(n)). For the rice images, the rice are very close to each other, and the LIF model fails to segment them (see Figure 5(o)). With the constant function initialization (𝜙₀ = 0.5), our method can achieve satisfactory results for the plane and rice images (see Figures 5(s) and 5(t)). We also demonstrate the accuracy of our model by quantitative comparison. The metric adopted in this paper is the dice similarity coefficient (DSC) [24] as follows:𝑆DSC=2𝑁₁∩𝑆₂𝑁𝑆₁+𝑁𝑆₂, (4.5) where 𝑁(⋅) indicates the number of pixels in the enclosed set, and 𝑆₁ and 𝑆₂ represent a given baseline foreground region (e.g., true object) and the foreground region found by the model, respectively. The closer the DSC value to 1, the better the segmentation. Table 1 shows the DSC values of the LIF model and our model. It can be clear that our model achieves more accurate results.

Table 1 

The DSC values of the LIF model [18] and our model for the images in Figure 5.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 5 

Segmentation results on bimodal images. Row 1: original images. Row 2: results of thresholding method (left to right: the threshold values are 140, 150, 180, 110, and 110, resp.). Row 3: results of the LIF model [18]. Row 4: results of our model (𝜈=0.4).

Figure 6 shows the segmentation results of our model for five real images with intensity inhomogeneity. Five medical images in different scenario are chosen to serve as the test images, which are two X-ray vessel images (111 × 110, 132 × 131), a bladder MR image (180 × 107), a brain MR images (120 × 160), and a wrist image (90 × 196), as shown in Row 1. Our level set evolution starts with 𝜙₀=𝜌 (left to right: 𝜌 = 0.5, 0.5, −0.5, 0.5 and 0.5, resp.). As can be seen in the second row of Figure 6, our model successfully extracted all objects contours from nonuniform background. This experiment also demonstrates that our model can handle more general images with intensity inhomogeneity.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 6 

Segmentation results of our model on images with intensity inhomogeneity (𝜈 = 0.50). Row 1: original images. Row 2: results of our model (left to right: 60, 100, 20, 20, and 160 iterations, resp.).

Figure 7 shows that our model can extract weak objects boundaries from images with complex background (250 × 180). In actual sky-mountain-water conflicts, the background of these real infrared images was quite complex (see Row 1). Moreover, parts of objects (boat, sun) are quite weak. Our level set evolution starts with 𝜙₀ = 0.5. As can be seen in the second row of Figure 7, our model successfully extracts the desired objects from the complex background.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7 

Segmentation results of our model on images with complex background (left to right: 𝜈 = 0.23, 0.23, 0.35, 0.28). Row 1: original images. Row 2: results of our model (left to right: 50, 50, 120, and 10 iterations).

Figure 8 shows that our model can handle real noisy image with different types of shapes. Five test images, which are shown in Row 1, are two breast cyst images (91 × 92, 157 × 110), an MR image of a human brain (159 × 122), a skin lesion image (180 × 180), and a ultrasonic image (100 × 100). These images are contaminated by texture tissue and clutter noise, thus it is very difficult to segment the objects in such images. Our level set evolution starts with 𝜙₀=𝜌 (left to right: 𝜌 = 0.5, 0.5, 0.5, −0.5 and −0.5, resp.). As can be seen in Row 2, our model successfully extracted the desired objects from the noisy images.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 8 

Segmentation results of our model on real medical images (left to right: 𝜈 = 0.22, 0.22, 0.22, 0.20, and 0.185, resp.). Row 1: original images;Row 2:results of our model (left to right: 120, 180, 300, 290, and 400 iterations).

5. Conclusion

We have proposed a novel variational level set formulation for image segmentation based on weighted 𝑝(𝑥)-Dirichlet integral and LoG filter. By incorporating local intensity information into the weighted 𝑝(𝑥)-Dirichlet integral, this regularization term preserves the good properties of both length regularization and 𝑝-Dirichlet integral regularization, as well as a combination of the two. The external energy based on the LoG filter drives the level set function up or down on either side of edge and locates objects edges. Due to the good properties of the weighted 𝑝(𝑥)-Dirichlet integral term, our model can extract weak edge in noise and/or intensity inhomogeneity images. The proposed model also allows the use of more general initialization of the level set function, that is, constant function. This implies that our model is free of manual initialization. Given its efficiency and accuracy, we expect that the proposed model will be useful for real-time applications. In our future work, we will extend the current model to motion tracking, where motion information can be included in the segmentation process.

Acknowledgments

The authors are very grateful to the anonymous reviewers for their careful reading and useful suggestions, which greatly improved the presentation of the paper. This work is supported by the Fundamental Research Funds for the Central Universities (CDJXS10100005) and Natural Science Foundation Project of CQ CSTC (cstcjjA40012) and CQ CSTC (cstc2011jjA40029).