ImprovedMumford-Shah Functional for Coupled Edge-Preserving Regularization and Image Segmentation

An improved Mumford-Shah functional for coupled edge-preserving regularization and image segmentation is presented. A nonlinear smooth constraint function is introduced that can induce edge-preserving regularization thus also facilitate the coupled image segmentation. The formulation of the functional is considered from the level set perspective, so that explicit boundary contours and edge-preserving regularization are both addressed naturally. To reduce computational cost, a modified additive operator splitting (AOS) algorithm is developed to address diffusion equations defined on irregular domains and multi-initial scheme is used to speed up the convergence rate. Experimental results by our approach are provided and compared with that of MumfordShah functional and other edge-preserving approach, and the results show the effectiveness of the proposed method.


INTRODUCTION
Mumford-Shah (MS) functional is an important variational model in image analysis.It minimizes a functional involving a piecewise smooth representation of an image and penalizing the Hausdorff measure of the set of discontinuities, resulting in simultaneous linear restoration and segmentation [1,2].
However, the MS functional is based on Bayesian linear restoration, so the resultant linear diffusion not only smoothes all structures in an equal amount but also dislocates them more and more with the increasing scale that may blur true boundaries [2].The situation becomes worse for poor-quality images with artifacts and low contrast, making the coupled segmentation unreliable.To address this problem, many improvements on MS model from nonlinear diffusion perspective are developed.However, due to the unknown discontinuities set of lower dimension, most approaches solve the weak formulation of the improved functional.In [3], the smooth constraint and the data fidelity are defined by the norm functions instead of quadratic functions in the MS functional.The resultant diffusion is modulated by the magnitude of the gradient that can deblur the edges.In [4], an edge-preserving regularization model based on the half-quadratic theorem is proposed, where the diffusion is nonlinear both in intensity and edges.But these approaches solving weak formulation concentrate rather on image restoration than image segmentation.Therefore, exact boundary locations cannot be explicitly yielded.
To solve the MS functional in such a way that image segmentation can be explicitly yielded, level set and curve evolution formulations of the MS functional have been developed in recent years.By viewing an active contour as the set of discontinuities, active contours without edges model is proposed [5,6].It introduces Heaviside function and embeds the level set function into MS model for piecewise constant and piecewise smooth optimal approximations, leading to the coupled image restoration and level set evolution.Similar work can be found in [7], where numerical implementation of the MS functional from the level set perspective is derived in detail.And also in [8], the MS functional is formulated from curve evolution perspective.
In this paper, inspired by the nonlinear diffusion theory and level set method, an improved Mumford-Shah functional is presented from both the theoretical and numerical aspects.The main contribution of this paper is as follows.First, a nonlinear smooth constraint function is proposed and introduced into the functional that can induce edge-preserving regularization.Second, different from existing edge-preserving approaches that solve the weak formulation of the problem, we formulate the proposed functional from the level set perspective so that nonlinear edgepreserving regularization and explicit boundary contours are both addressed naturally.Third, to reduce the computational cost, a modified additive operator splitting (AOS) algorithm is developed to address the diffusion equations defined on irregular domains.Furthermore, multi-initial scheme is used to speed up the convergence rate.
The remainder of the paper is organized as follows.In Section 2, mathematical background is sketched; in Section 3, an improved Mumford-Shah functional is proposed and level set formulation of the functional is derived.In Section 4, numerical implementation of the improved functional is described in detail, where a modified AOS algorithm is proposed.In Section 5, experimental results are provided and comparisons are discussed.Finally in Section 6, conclusions are reported.

Mumford-Shah functional
Let Ω ⊂ R m be open and bounded image domain, and let f : Ω → R be the original image data, the linear restoration of the ideal image data u is formulated as where n ∼ N(0, σ 2 ) is the white noise.By introducing Markov random field (MRF) line process, the above restoration is expressed as the minimization problem [2]: where ω λ,μ (x) = min(λx 2 , μ), λ and μ are two positive weights.The continuous form of (2) is the MS functional [1]: where C ⊂ Ω is the set of discontinuities in the image and α, β, γ are weights that control the competition of the various terms.The first term comes from the piecewise smooth constraint.The second term is the data fidelity term that makes the restoration more like its original.And the third term stands for the (m − 1)-dimensional Haussdorf measure of the set of discontinuities.Due to the unknown discontinuities set of lower dimension, it is not easy to minimize MS functional in practice.Some approaches solve the weak formulation of the MS functional [9][10][11].But the weak formulation cannot explicitly yield the boundary contours.However, from curve evolution perspective [8], minimizing (3) with respect to u + , u − , C, we can obtain the coupled diffusion and curve evolution equations: where u + and u − represent the value inside and outside the current curve C, respectively.N is the normal vector of the curve and κ is the curvature of the curve C.

Nonlinear edge-preserving regularization
As was discussed above, MS functional is derived from linear restoration problem.However, in most cases, image restoration is a complex and nonlinear process where the linear relationship between f and u in (1) does not come into existence.Consider the nonlinear regularization such that the restored image can preserve edge while remove noise, then the restoration energy can be expressed as follows [12]: From variational method, the corresponding Euler-Lagrange equation minimizing (6) is the stable state of the following diffusion partial difference equation: where n denotes the normal to the image domain boundary ∂Ω, and ψ (|∇u|)/2|∇u| is the diffusion coefficient that controls the diffusion process.To assure edge preserving and the stability, Discussions on the choice of ψ(•) can be found in [4,12].
When the diffusion coefficient is equal to 1, (7) degrades to homogenous linear diffusion.In this view, ( 4) is homogenous linear diffusion inside and outside the current curve C, respectively.Although the diffusion is not across C, it smoothes the regions inside and outside C in an equal amount.Therefore, before the curve arrives at the true boundary, the true boundaries may become obscure more and more with the increasing scale that can dislocate the boundary position and also mislead the coupled curve evolution.

IMPROVED MUMFORD-SHAH FUNCTIONAL
In order to perform edge-preserving regularization and segmentation simultaneously, we present an improved Mumford-Shah functional.Inspired by the idea of nonlinear diffusion theory, we propose and introduce a nonlinear smooth constraint function ψ(|∇u|) to replace |∇u| 2 in (3).The improved Mumford-Shah functional is given by where ψ(|∇u|) is a nonlinear increasing function of |∇u| 2 for piecewise smooth constraint and should also satisfy the conditions (1), (2), and (3) for edge-preserving regularization.
Inspired by the level set method [13], we minimize (8) from the level set evolution perspective.Embed the current curve C as zero level set of a higher-dimensional level set function φ and introduce the Heaviside function and Dirac function δ 0 (φ) = (d/dφ)H(φ), then the length is approximated by L(φ) = Ω |∇H(φ)|dx dy.We can formulate the functional (8) as follows: Let μ = α/β, minimizing E IMS (u + , u − , φ) with respect to u + , u − , φ, we can obtain the coupled equations: where δ ε (φ) is a slightly regularized version of δ 0 (φ).The detailed derivation is provided in appendix.
Choosing proper ψ(|∇u|) such that ψ (|∇u|)/2|∇u| is close to 1 in the homogenous regions and close to zero on the edges, then (11) and ( 12) can perform edge-preserving diffusion inside and outside the current curve, which can remove noise while preserve the edges even with the increasing scale.Therefore, the coupled curve evolution can be precisely guided towards the true boundaries.Note that MS functional could be considered as the special case when ψ(x) = x 2 .In this case the diffusion coefficient ψ (x)/2x ≡ 1, thereby leading to homogenous linear diffusion.

Modified AOS algorithm for diffusion equations
In general, standard explicit difference scheme can be used to iteratively update u + and u − [7].However, the convergence rate is very slow.In [14], an efficient AOS scheme, which uses semi-implicit scheme discretization and uses an additive operator splitting to decompose complex process into multiple linear process, is proposed to solve diffusion equations defined on the rectangular image domain as expressed in (7).Its separability allows parallel implementations, making the scheme ten times faster than usual explicit one.However, diffusion equations ( 11) and ( 12) are defined on irregular domains divided by the current curve, making the question more complex and difficult.Aiming at this point, a modified AOS scheme is developed.
Let Ψ(X) = ψ (|∇u + + |.The diffusion equation ( 11) is equivalent to Use center difference scheme and let In that case, ( 14) can be discretized as Now we only consider pixels {i | φ i > 0} and store u i as a vector u + .Since H(φ i ) = 1 for φ i > 0, we can obtain where N l (i) is the set of the two neighbors of pixel i along l direction, and In vector-matrix notation and using semi-implicit discretization, (17) can be written as where A l (u n,+ ) is a matrix along l direction whose element is a n,+ i jl .The solution is given by Directly solving (20) leads to great computation effort.Using AOS scheme, we can obtain [14] The key of AOS algorithm is that A l (u n,+ ) is tridiagonal and diagonally dominant along the l direction.Therefore (21) can be converted to solve tridiagonal system of linear equations for m parallel processes that can be rapidly implemented by Thomas algorithm.Similarly, solving (12) we can obtain where A l (u n,− ) is obtained by replacing H(φ) by 1 − H(φ) in (18).

Extend u + on {φ < 0} and u − on {φ > 0}
For solving the level set equation ( 13), extending u + on {φ < 0} and u − {φ > 0} is necessary for calculating the jumps.This can be made by the modified five-point averaging scheme as follows [15]:

Let
Discretization of the level set equation ( 13) yields Consequently we have

Algorithm description for the coupled PDEs
The coupled PDEs (11), (12), and ( 13) should be alternatively iterated until convergence.The diffusion equations ( 11) and ( 12) are solved by the modified AOS algorithm.The level set evolution ( 13) is by standard explicit iterate scheme, where multi-initial scheme is used so that it can speed up the convergence rate and also it has the tendency to converge to a global minimizer [6].The complete algorithm is described as follows.
Step 1. Use multi-initial scheme to plant seed curves and set the initial front curve to be pixels on all seed curves.Initialize φ 0 as a signed distance function to the initial front curve.

RESULTS AND DISCUSSIONS
To evaluate the effectiveness of the proposed model, a representative CT pulmonary vessel image is chosen as an example.The vessels and their branches, which exhibit much variability with artifacts and low contrast, make segmentation and restoration very difficult.In the experiment, a regularization H ε of the Heaviside function is used.We choose 2 for our model and γ = 0.002 * 255 2 for MS model.We set , where λ is the contrast parameter separating the forward diffusion from backward diffusion, which is chosen to be the 50 percent quantile of which is close to 1 as x → 0 and 0 as x → ∞, leading to edgepreserving diffusion.
The first experiment compares the performance of the improved Mumford-Shah functional and original MS functional.Figures 1(A) and 1(B) show the coupled diffusion and curve evolution by the improved and original MS functional, respectively.From Figure 1(A), it can been seen that our model has succeeded in preserving the locality of vessel boundaries while smoothing the interior of the vessel, so that the restored image is edge-preserving regularization as shown in Figure 1(A)(h).Therefore, the final segmentation results are promising in that vessels, even their thin branches, could be located precisely, as shown in Figure 1(A)(g).Whereas in Figure 1(B), we can see that almost all structures are blurred with the increasing scale, thus the restored image is not satisfactory in that important vessel branches become more and more obscure and the boundaries are dislocated to some extent as shown in Figure 1(B)(h).Therefore, the resultant segmentation results as shown in Figure 1(B)(g) are not reliable and some small vessel branches cannot be extracted precisely as well.
Furthermore, the segmentation results are quantitatively evaluated by the criteria of intraregion uniformity and interregion discrepancy.The intraregion uniformity can be measured by the variance within each region and the interregion discrepancy by the difference of the mean between adjacent regions.The smaller the variance is, the better the intraregion uniformity is.And the larger the difference is, the better the interregion discrepancy is.The comparison results are illustrated in Table 1.We can conclude that our model achieves better segmentation results than MS functional in that the variance of each region is smaller and the difference between regions is larger by our approach.It also indicates that the regularization by our approach is more reliable because it can induce better segmentation result.
The second experiment compares our model with the other edge-preserving approach.Figure 2(A) shows the regularization and segmentation results by our approach, and Figure 2(B) by the approach proposed in [4].Both approaches can perform edge-preserving regularization as shown in Figures 2(A)(a) and 2(B)(a).However, the segmentation results are different.Our approach formulates the improved Mumford-Shah functional from the level set perspec-tive, thereby explicit boundary contours can be obtained naturally as shown in Figure 2(A)(b).The edge-preserving approach in [4] formulates its functional from the weak formulation, resulting in two coupled diffusion equations.One is for image intensity and the other is for edge function.Observe in Figure 2(B)(b) that only the image of the edge function is obtained and that cannot explicitly give the boundary location.The comparison shows that our approach is very effective both in regularization and segmentation results.

CONCLUSION
In this paper, an improved Mumford-Shah functional is presented that can perform edge-preserving regularization and image segmentation simultaneously.Different from existing edge-preserving approaches, the formulation of the proposed functional is considered from the level set evolution perspective thus explicitly yielding boundary contours and image restoration.Both are addressed naturally.A modified AOS algorithm is developed to address the diffusion equations defined on irregular domains.It is ten times faster than usual explicit scheme.The experimental results are evaluated by the criteria of intraregion uniformity and interregion discrepancy, which show that our model outperforms MS functional both in segmentation and regularization.Comparison with other edge-preserving approach shows that our approach is very promising.It can be applied to a variety of image segmentations and restorations.

APPENDIX
As ∇H(φ) = H (φ) • ∇φ = δ 0 (φ) • ∇φ, (10) can be rewritten as From variational method, fixing u the Euler-Langrange equation for φ is Coupled edge-preserving regularization and curve evolution by the improved Mumford-Shah functional Coupled linear diffusion and curve evolution by the MS functional

Figure 1 :Figure 2 :+ 2 −(A. 4 ) 9 )
Figure 1: Coupled diffusion and curve evolution for CT pulmonary vessel segmentation and restoration (A) by the improved Mumford-Shah functional and (B) by the MS functional.(a) Original CT pulmonary vessel.(b)-(f) Coupled curve evolution and diffusion.(g) and (h) show the final segmentation and restoration results, respectively.

Table 1 :
Comparison of segmentation results by improved Mumford-Shah functional and by MS functional.IMS: improved Mumford-Shah functional; MS: MS functional.