Adaptive Regularized Level Set Method for Weak Boundary Object Segmentation

An adaptive regularized level set method for image segmentation is proposed. A weighted p x 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 p x 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 proposedmodel has been applied to awide range of both real and synthetic images with promising results.


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 Ω |∇φ| 2 dxdy and |∇φ| L ∞ Ω 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 p ≥ 1 results in a constrain which is somewhere between length-based and smoother regularization.However, the constant exponent p 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 p x -Dirichlet integral, an external energy, and a level set regularization term.The weighted p x -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 p x -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 p x -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.

Level Set Method and Level Set Function Regularization
In the level set formulation, a moving curve C t is represented by the zero level set of a Lipschitz function φ x, y, t defined on the entire image domain.The evolution of the curve C t along its normal direction with speed F is described by the following evolution PDE 7 : with the initial condition φ x, y, 0 φ 0 x, y .For image segmentation, the speed function F 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 where φ is the function to be re-initialized, and sign • is the sign function.Although reinitialization 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: where P φ is the internal energy and E 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 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 .

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: where H • is Heaviside function and δ • is Dirac delta function.The energy functional L φ in 2.5 computes the length of the zero level curve of φ in the conformal metric ds |C p |dp.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 Vese proposed a smoother regularization: R φ Ω ∇φ 2 dx dy.

2.6
This regularization means isotropic smoothing at every point x, y , 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 Mu proposed a weighted p-Dirichlet integral regularized level set evolution.The geometric regularization has the following form: Different value of p ≥ 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 p.

Weighted p x -Dirichlet Integral Regularized Level Set Evolution
In this section, we propose a new variational level set formulation for image segmentation, in which a weighted p x -Dirichlet integral is presented to regularize the zero level curve.
Let Ω ⊂ R 2 be a image domain.For a given image I : Ω → R and a level set function φ x, y : Ω → R, we define an energy functional E φ by where ν, μ > 0 are constants, L p • φ is the zero level curve regularization term, E ext φ is a certain external energy that would drive the motion of the zero level curve of φ, and P φ is the level set function regularization term that controls the smoothness of the level set function during the level set evolution.

Weighted p x -Dirichlet Integral
The zero level curve regularization term L p • φ is defined as follows: where ∇ is gradient operator and G σ * I is the convolution of the image I with the Gaussian function G σ with standard deviation σ.The exponent p s : 0, ∞ → 1, 1.5 is a monotonically increasing function with lim s → 0 p s 1 and lim s → ∞ p s 1.5.A simple example is The functional 3.2 is in fact a weighted p-Dirichlet integral with variable exponent p |∇G σ * I| , thus it is called the weighted p x -Dirichlet integral in this paper.
It is well known that the energy functional computes the length of the zero level curve of φ.
2 When p |∇G σ * I| p > 1 is a constant, the functional L p • φ is the weighted p-Dirichlet integral.It is chosen as one of geometric regularization on zero level curve in 23 .
However, the constant exponent p is not an intelligent choice.An important reason is that the constant exponent p does not reflect the local property of image, and so does not adapts the exponent to fit the image data automatically.We let p p |∇G σ * I| depend on local intensity information.This way the regularization ensures weak regularization p |∇G σ * I| ≈ 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 p |∇G σ * I| ≈ 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 p x controls the tradeoff between weak regularization and strong regularization automatically.

External Energy Term E 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 p x -Dirichlet integral in 3.2 can be used in various applications with different definitions of the external energy E ext φ .In this subsection, we only provide an application of the weighted p x -Dirichlet integral to an active contour model using edge-based information in the external energy, as a demonstration of the effectiveness of the weighted p x -Dirichlet integral formulation.
Next, we define the external energy E 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: where ΔG σ 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.

The Energy Formulation
In order to keep the regularity of level set function, the energy P φ 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

3.6
In order to compute the associated Euler-Lagrange equation for the unknown function φ, we consider slightly regularized versions of the function H and δ φ , denoted here by H ε and δ ε , as ε → 0. Let H ε be a C ∞ Ω regularization of H, and δ ε H ε .In a dynamical scheme via steepest descent, minimizing the energy functional 3.6 with respect to φ, we obtain the evolution PDE: with the initial condition φ x, y, 0 φ 0 x, y .In all numerical experiments, we choose the following functions: 3.8 for our evolution PDE 3.7 .

Weighted p X -Dirichlet Integral Effect
In this subsection, we demonstrate the weighted p x -Dirichlet integral 3.2 effect by a simple experiment, just as shown in Figure 1.We will see that the variable exponent

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 φ x, y, t .The spatial derivatives ∂φ/∂x and ∂φ/∂y in our model are approximated by the central difference, and the temporal partial derivative ∂φ/∂t is discretized as the forward difference.
We recall first the usual notations: let Δt be the time step, h be the space step, and x i , y i ih, jh be the grid points.Let φ n i,j φ x i , y j , nΔt be an approximation of φ x, y, t , with n ≥ 0, φ 0 φ 0 , the central differences are Set n 0, and start with initial level set function φ 0 i,j , we first compute p |∇G σ * I| according to 3.3 .Then, the numerical approximation to 3.7 is given by the following discretization: where This system 4.2 is solved by an iterative method, and we use fixed space step h 1 implies pixel spacing for computing spatial derivatives.The choice of the time step for this finite different scheme must satisfy the condition μΔt < 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.

Experimental Results
This section shows the results of the proposed model for both synthetic and real images.The level set function φ x, y, t is simply initialized to φ 0 x, y ρ 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, Δt 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 I, 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 φ 0 x, y ρ top to bottom: ρ 0.5, −0.5, 0.5, and 0.5, resp., there are no initial contours zerolevel 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 zerolevel 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 generat-ed false small contours gradually disappear see Column 4 .This shrinking effect can be interpreted as the regularization of the weighted p x -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.
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 φ 0 x, y ρ top: ρ 0.5, bottom: ρ −0.5 .Column 4 shows that our model successfully extracts all objects contours from background.
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 .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 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: where N • indicates the number of pixels in the enclosed set, and S 1 and S 2 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  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 φ 0 ρ 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.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 0.5.As can be seen in the second row of Figure 7, our model successfully extracts the desired objects from the complex background.

Conclusion
We have proposed a novel variational level set formulation for image segmentation based on weighted p x -Dirichlet integral and LoG filter.By incorporating local intensity information into the weighted p x -Dirichlet integral, this regularization term preserves the good properties of both length regularization and p-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 p x -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.

Table 1 :
The DSC values of the LIF model 18 and our model for the images in Figure5the LIF model and our model.It can be clear that our model achieves more accurate results.

Figure 8
Figure7shows 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 0.5.As can be seen in the second row of Figure7, our model successfully extracts the desired objects from the complex background.Figure8shows 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 φ 0 ρ 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.
1 Initialize the level set function φ 0 constant, set n 0. 2 Compute p |∇G σ * I| according to 3.3 .3 Solve the PDE in φ from 4.2 , to obtain φ n 1 .4 Check whether the evolution has stationary.If not, n n 1 and repeat.