A Multiphase Image Segmentation Based on Fuzzy Membership Functions and L1-norm Fidelity
aa r X i v : . [ m a t h . O C ] F e b Journal of Scientific Computing manuscript No. (will be inserted by the editor)
A Multiphase Image Segmentation Based on FuzzyMembership Functions and L1-norm Fidelity
Fang Li · Stanley Osher · Jing Qin · Ming Yan
Received: date / Accepted: date
Abstract
In this paper, we propose a variational multiphase image segmen-tation model based on fuzzy membership functions and L1-norm fidelity. Thenwe apply the alternating direction method of multipliers to solve an equivalentproblem. All the subproblems can be solved efficiently. Specifically, we proposea fast method to calculate the fuzzy median. Experimental results and com-parisons show that the L1-norm based method is more robust to outliers suchas impulse noise and keeps better contrast than its L2-norm counterpart. The-oretically, we prove the existence of the minimizer and analyze the convergenceof the algorithm.
Keywords
Image segmentation; fuzzy membership function; L1-norm;ADMM; segmentation accuracy.
As a fundamental step in image processing, image segmentation partitions animage into several disjoint regions such that pixels in the same region sharesome uniform characteristics such as intensity, color, and texture. During thelast two decades, image segmentation methods using variational methods and
F. Li ( (cid:12) )Department of Mathematics, East China Normal University, and Shanghai Key Laboratoryof Pure Mathematics and Mathematical Practice, Shanghai, China.E-mail: e-mail: fl[email protected]. Osher, J. QinDepartment of Mathematics, University of California, Los Angeles, CA 90095, USA.E-mail: [email protected],[email protected] M. YanDepartment of Computational Mathematics, Science and Engineering and Department ofMathematics, Michigan State University, East Lansing, MI 48824, USA.E-mail: e-mail:[email protected] Fang Li et al. partial differential equations are very popular due to their flexibility in problemmodeling and algorithm design. There are two key ingredients of variationalsegmentation models. One is how to describe the regions or boundaries betweenthese regions, and the other is how to model the noise and describe the uniformcharacteristics of each region.The Mumford-Shah model [35], a well-known variational segmentationmodel, penalizes the combination of the total length of the segmentationboundaries and the L2-norm error of approximating the observed image withan unknown piecewise smooth approximation. In other words, the Mumford-Shah model seeks an optimal piecewise smooth function with smooth bound-aries to approximate the observed image.However, the Mumford-Shah model is hard to implement in practice be-cause the discretization of the unknown set of boundaries is very complex.Therefore, a parametric/explicit active contour method is used to represent thesegmentation boundaries [49]. In addition, the special Mumford-Shah modelwith a piecewise constant approximation is studied by Chan and Vese [15], andthe level set method [37] is applied to solve this problem. Using an implicitrepresentation of boundaries, the level set method has many advantages overmethods using explicit representations of boundaries. For instance, the levelset method handles the topological change of zero level set automatically [8,2,22]. Both the parametric/explicit active contour method and the level setmethod assume that each pixel belongs to a unique region. An alternativeway to represent various regions is to use fuzzy membership functions [14,6,34,25], which describe the levels of possible membership. Fuzzy membershipfunctions assume that each pixel can be in several regions simultaneously withprobability in [0,1]. These probabilities satisfy two constraints: (i) nonnegativ-ity constraint, i.e., the membership functions are nonnegative at all pixels; (ii)sum-to-one constraint, i.e., the sum of all the membership functions at eachpixel equals one. Then the length of boundaries can be approximated by theTotal Variation (TV) of fuzzy membership functions. The main advantages ofusing fuzzy membership functions over other methods include: i) the energyfunctional is convex with respect to fuzzy membership functions, guarantee-ing the convergence and stability of many numerical optimization methods. ii)fuzzy membership function has a larger feasible set, and it is possible to findbetter segmentation results.For two-phase segmentation, where there are only two regions, we only needone level set function or one fuzzy membership function. Multiphase segmenta-tion is more challenging than two-phase segmentation since more variables andconstraints are involved in representing multiple regions and their boundarieseffectively. The two-phase Chan-Vese model [15] has been generalized to multi-phase segmentation by using multiple level set functions to represent multipleregions [44]. Partitioning an image into N disjoint regions requires log N levelset functions. The advantage of using multiple level set functions is that it au-tomatically avoids the problems of vacuum and overlap of regions. However,the implementation is not easy, and special numerical schemes are needed toensure stability [32,33,42,27]. For fuzzy membership functions, the sum-to-one ultiphase Image Segmentation via Fuzzy and L1 3 constraint is not satisfied automatically in multiphase segmentation. However,this constraint is easy to deal with in many cases, e.g., Fuzzy C-Mean (FCM)and its variants have closed-form solutions for the membership functions andare widely used in data mining and medical image segmentation [4,1,38,16,29,31,30,24,7]. Variable splitting schemes are used in both [30] and [31] toget efficient numerical algorithms. The Alternating Direction Method of Mul-tipliers (ADMM) method is used in [24] to derive the algorithm with two setsof extra variables. Primal-dual type algorithm is derived in [10] to solve theTV regularized FCM segmentation model. Both [24] and [10] use projectionto simple to handle the constraints of membership functions. Other segmen-tation approaches include a convex approach [9], two-stage methods [11,43],one single level set function approach [18], et.al.Noise is unavoidable in images, and it is important to develop segmentationmethods that work on noisy images. Among many types of noise, the Gaus-sian white noise is frequently assumed, and the L2-norm fidelity is adopted.However, when images are corrupted by non-Gaussian noise, in order to obtaina faithful segmentation, one has to model the noise according to its specifictype [40,36,12,23,47,48,7]. Particularly, the L1-norm fidelity is used for bothsalt-and-pepper impulse noise and random-valued impulse noise in image pro-cessing [36,12,23]. In addition, it is robust to outliers and able to preservecontrast because the denoising process of L1-norm models is determined bythe geometry such as area and length rather than the contrast in the L2-normcase [13].Inspired by the fact that L1-norm is more robust to impulse noise andoutliers and can better preserve contrast, in this paper, we propose a vari-ational multiphase fuzzy segmentation model based on L1-norm fidelity and fuzzy membership functions . This model can also deal with missing data inimages. When there are missing pixels in an image, we randomly assign 0 or255 at these pixels by considering these pixels as corrupted by salt-and-pepperimpulse noise. ADMM [20,19], which was rediscovered as split Bregman [21]and found to be very useful for L1 and TV type optimization problems, isapplied to solve this nonconvex optimization problem. By introducing twosets of auxiliary variables, we derive an efficient algorithm with all the sub-problems having closed-form solutions. In the theoretical aspect, we prove theexistence of the minimizer and analyze the convergence of the algorithm. Wenote that the proposed method is closely related to the method in [24] sinceboth methods use TV regularization and ADMM algorithm. The difference isthat L1-norm fidelity is considered in the proposed method, while L2-normfidelity is used in [24].The outline of this paper is as follows. In Section 2, we review some closelyrelated existing works. In Section 3, the proposed model is described in detail,and the existence of the minimizer to the model is proved. In Section 4, anumerical algorithm is derived, and its convergence analysis is presented. InSection 5, experimental results and comparisons are presented to illustratethe effectiveness of the proposed method. Finally, the paper is concluded inSection 6.
Fang Li et al.
Let Ω ⊂ R be a bounded open subset with Lipschitz boundary, and let I : Ω → R s be the given clean or noisy image. Let s = 1 for grayscale imagesand s = 3 for color images. Our goal is to find an N -phase “optimal” partition { Ω i } Ni =1 such that Ω i T Ω j = ∅ for all i = j and S Ni =1 Ω i = Ω . Define the setof N -phase fuzzy membership functions as ∆ := ( ( u , ..., u N ) (cid:12)(cid:12)(cid:12) u i ∈ BV ( Ω ) , u i ( x ) ≥ , N X i =1 u i ( x ) = 1 , ∀ x ∈ Ω ) , where BV ( Ω ) is the space of functions with bounded variation [2]. The closelyrelated works are listed in the following and will be compared with our pro-posed method in Section 5. For the sake of simplicity, we use the notations U = ( u , · · · , u N ) and C = ( c , · · · , c N ) ∈ R sN , where c i ∈ R for grayscaleimages and c i ∈ R for color images. – FCM [4]– Fuzzy c-means clustering method that solvesmin ( U , C ) ∈ ∆ × R sN N X i =1 Z Ω | I ( x ) − c i | u pi ( x ) dx using the alternating minimization algorithm. Though p can be any numberno smaller than 1, it is commonly set to 2. – FCM S2 [16] – A variant of FCM that solvesmin ( U , C ) ∈ ∆ × R sN N X i =1 (cid:26)Z Ω | I ( x ) − c i | u pi ( x ) dx + α Z Ω (cid:12)(cid:12) ¯ I ( x ) − c i (cid:12)(cid:12) u pi ( x ) dx (cid:27) , where ¯ I is obtained by applying the median filter on I and α > – FLICM [29] – Fuzzy Local Information C-Means clustering method thatsolves min ( U , C ) ∈ ∆ × R sN N X i =1 (cid:26)Z Ω | I ( x ) − c i | u pi ( x ) dx + α Z Ω Z y ∈N ( x ) (1 − u i ( y )) p | I ( y ) − c i | u pi ( x ) dydx ) , where N ( x ) is a neighborhood of x . FLICM is more robust to both Gaus-sian noise and impulse noise than FCM. – L2FS [24] – L2-norm fidelity based Fuzzy Segmentation method that solvesmin ( U , C ) ∈ ∆ × R sN N X i =1 (cid:26)Z Ω k∇ u i ( x ) k dx + λ Z Ω | I ( x ) − c i | u i ( x ) dx (cid:27) , (1) ultiphase Image Segmentation via Fuzzy and L1 5 by ADMM. Here λ > R Ω k∇ u i ( x ) k dx denotes the TVof u i with k∇ u i ( x ) k := p ( ∇ x u i ( x )) + ( ∇ x u i ( x )) for x = ( x , x ). Forfixed C , [24] is related to the popular TV denoising method [39]. Note thatthe similar model is solved by other fast numerical methods in [30]. – L1SS [26] – L1-norm fidelity based Soft Segmentation method, in whichlog N soft functions are introduced to represent N phases. Since the modelfor multiphase segmentation is complicated for more than four phases, weshow the four-phase model as follows:min u ,u ∈ [0 , , C ∈ R sN X i =1 Z Ω k∇ u i ( x ) k dx + λ X j =1 Z Ω | I ( x ) − c j | M j ( x ) dx , where the membership functions M j , j = 1 , · · · ,
4, are represented by softfunctions u ( x ) , u ( x ) ∈ [0 ,
1] in the following way: M ( x ) = u ( x ) u ( x ) , M ( x ) = u ( x )(1 − u ( x )) ,M ( x ) = (1 − u ( x )) u ( x ) , M ( x ) = (1 − u ( x ))(1 − u ( x )) . – L2L0 [43] – L2-norm fidelity and L0-norm regularization based image par-tition model: min u k∇ u k + λ k u − I k . This model seeks a piecewise constant approximation of the original image I . Since this model can not specify the number of classes, a second stepis applied to combine some classes if this model returns more classes thanrequired. Here we apply FCM on the piecewise constant approximation u to obtain the final segmentation result. Remark:
There are two advantages of our proposed method over L1SS. Firstly,we use fuzzy membership functions to represent regions, where N fuzzy mem-bership functions are required for an N -phase segmentation. Hence, the solu-tion space is much larger than L1SS, which ensures the higher possibility toobtain optimal segmentation. Secondly, the proposed method can take use ofother commonly used segmentation methods such as FCM to gain good initial-ization of fuzzy membership functions. Multiphase segmentation is sensitiveto initialization, and a good initialization is very important for a successfulsegmentation. However, it is hard to use the existing segmentation methodsto get a good initialization for soft membership functions in L1SS. In this paper, we assume that the given image can be approximated by apiecewise constant function, i.e., I ( x ) = N X i =1 c i χ Ω i ( x ) + n ( x ) . Fang Li et al.
Here χ Ω i denotes the indicator function of region Ω i , i.e., χ Ω i ( x ) = (cid:26) , if x ∈ Ω i ;0 , otherwise ,c i is a constant that represents the given data in region Ω i , and n can beoutliers, impulse noise or other mixture types rather than Gaussian noise.Instead of using the L2-norm fidelity to measure the approximation errorwhen the noise is the Gaussian white noise, we use the L1-norm fidelity. Sameas in the Mumford-Shah model, we require the segmentation boundaries to besmooth. Then we have the following modelmin { Ω i } , C N X i =1 Z Ω k∇ χ Ω i ( x ) k dx + λ Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) I ( x ) − N X i =1 c i χ Ω i ( x ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dx, (2)where λ > χ Ω i in the first termis equal to the length of boundary ∂Ω i . An equivalent form of model (2) ismin { Ω i } , C N X i =1 (cid:26)Z Ω k∇ χ Ω i ( x ) k dx + λ Z Ω | I ( x ) − c i | χ Ω i ( x ) dx (cid:27) . (3)Because χ Ω i can take values 0 and 1 only and { χ Ω i } is a partition, ( χ Ω , . . . , χ Ω N )belongs to the set ∆ := ( ( u , ..., u N ) (cid:12)(cid:12) u i ∈ BV ( Ω ) , u i ( x ) ∈ { , } , N X i =1 u i ( x ) = 1 , ∀ x ∈ Ω ) . At any point x ∈ Ω , there is only one function having value 1, and all theother functions have value 0. Thus set ∆ is not continuous, which leads todifficulties and instabilities in numerical implementations. However, we canrelax binary indicator functions { χ Ω i } to fuzzy membership functions { u i } ,which satisfy the nonnegativity constraint and the sum-to-one constraint, i.e.,( u , ..., u N ) belongs to the set ∆ defined in (2). It is obvious that u i ( x ) ∈ [0 , ∆ is a simplex at any x ∈ Ω . So u i ( x ) can be understood as the probabilityof pixel x to be in the i th class. Then we can rewrite our model (3) asmin ( U , C ) ∈ ∆ × R sN E ( U , C ) = N X i =1 (cid:26)Z Ω k∇ u i ( x ) k dx + λ Z Ω | I ( x ) − c i | u i ( x ) dx (cid:27) . (4)Note that model (4) is convex with respect to U and C separately but notjointly. The difference between (4) and (1) is that the L2 fidelity term in (1)is replaced by the L1 fidelity term. The existence of a minimizer for E ( U , C )in (4) is stated in Theorem 1. Theorem 1 (Existence of a minimizer)
Given an image I ∈ L ∞ ( Ω ) , thereexists a minimizer of E ( U , C ) in ∆ × R sN for any parameter λ > . ultiphase Image Segmentation via Fuzzy and L1 7 Proof
Since E ( U , C ) is positive and proper, the infimum of E ( U , C ) must befinite. Let us pick a minimizing sequence ( U n , C n ) ∈ ∆ × R sN that satisfies E ( U n , C n ) → inf U , C E ( U , C ) as n → ∞ . Then there exists a constant M > E ( U n , C n ) = N X i =1 (cid:26)Z Ω k∇ u ni ( x ) k dx + λ Z Ω | I ( x ) − c ni | u ni ( x ) dx (cid:27) ≤ M. Then each term in E ( U n , C n ) is bounded, i.e., for each i = 1 , · · · , N , Z Ω k∇ u ni ( x ) k dx ≤ M, Z Ω | I ( x ) − c i | u ni ( x ) dx ≤ M. (5)Since u ni ( x ) ∈ [0 , R Ω u ni ( x ) dx ≤ | Ω | , where | Ω | is the area of Ω .Together with the first equality in (5), we have that u ni is uniformly boundedin BV ( Ω ) for all i = 1 , ..., N . By the compactness property of BV ( Ω ) and therelative compactness of BV ( Ω ) in L ( Ω ), up to a subsequence also denotedby { u i } n after relabeling, there exists a function u ∗ i ∈ BV ( Ω ) such that u ni → u ∗ i strongly in L ( Ω ) ,u ni → u ∗ i a.e. x ∈ Ω, ∇ u ni ⇀ ∇ u ∗ i in the sense of measure.Then by the lower semicontinuity of the TV semi-norm, Z Ω k∇ u ∗ i ( x ) k dx ≤ lim inf n →∞ Z Ω k∇ u ni ( x ) k dx. (6)Meanwhile since U n = ( u n , ..., u nN ) ∈ ∆ , we have U ∗ = ( u ∗ , ..., u ∗ N ) ∈ ∆ .It is easy to derive the first order optimality condition about c ni , which is0 ∈ Z Ω ∂ | I ( x ) − c ni | u ni ( x ) dx, where ∂ | · | is the subdifferential of | · | . Since u ni ( x ) ≥ R Ω u ni ( x ) dx > c ni satisfies | c ni | ≤ k I k ∞ . By the boundedness of sequence { c ni } , we can extract a subsequence also de-noted by { c ni } such that for some constant c ∗ i c ni → c ∗ i , as n → ∞ . Finally, since u ni ( x ) → u ∗ i ( x ), a.e. x ∈ Ω and c ni → c ∗ i as n → ∞ , Fatou’slemma yields Z Ω | I ( x ) − c ∗ i | u ∗ i ( x ) dx ≤ lim inf n →∞ Z Ω | I ( x ) − c ni | u ni ( x ) dx. (7) Fang Li et al.
Combining inequalities (6) and (7) for all i , on a suitable subsequence, weestablished E ( U ∗ , C ∗ ) ≤ lim inf n →∞ E ( U n , C n ) = inf U , C E ( U , C ) , and hence ( U ∗ , C ∗ ) must be a minimizer of the energy E . This completes theproof. ⊓⊔ The minimizer of E ( U , C ) is not unique due to the following hidden sym-metry property. Denote S N as the permutation group of { , ..., N } , i.e., eachpermutation σ ∈ S N is defined as a one-to-one map σ : { , ..., N } → { , ..., N } such that { σ (1) , ..., σ ( N ) } is a rearrangement of { , ..., N } . Denote U σ =( u σ (1) , ..., u σ ( N ) ), C σ = ( c σ (1) , ..., c σ ( N ) ). It is straightforward to show the fol-lowing theorem. Theorem 2 (Symmetry of minimizer)
For any ( U , C ) ∈ ∆ × R N andany σ ∈ S N , E ( U σ , C σ ) = E ( U , C ) . In particular, suppose that ( U ∗ , C ∗ ) is a minimizer of E ( U , C ) , i.e., ( U ∗ , C ∗ ) = arg min ( U , C ) ∈ ∆ × R N E ( U , C ) . Then, for any σ ∈ S N , we have ( U ∗ σ , C ∗ σ ) = arg min ( U , C ) ∈ ∆ × R N E ( U , C ) . In this section, we provide an efficient algorithm based on ADMM and discussits convergence.4.1 The algorithmADMM is applied, in this subsection, to solve the proposed fuzzy segmentationmodel (4). We introduce two sets of auxiliary variables D = ( d , ..., d N ) , W =( w , ..., w N ) such that ∇ u i = d i , u i = w i . Then the model (4) is equivalent tothe following minimization problem with equality constraints:min D , W , C , U N X i =1 (cid:26)Z Ω k d i ( x ) k dx + λ Z Ω | I ( x ) − c i | w i ( x ) dx (cid:27) + δ ∆ ( W ) , subject to ∇ u i = d i , u i = w i , ∀ i = 1 , . . . , N, (8)where δ ∆ is the indicator function of the set ∆ , i.e., δ ∆ ( W ) = (cid:26) , if W ∈ ∆, + ∞ , otherwise . ultiphase Image Segmentation via Fuzzy and L1 9 The augmented Lagrangian function for problem (8) is: L ( D , W , C , U ; Λ D , Λ W )= N X i =1 (cid:26)Z Ω k d i ( x ) k dx + λ Z Ω | I ( x ) − c i | w i ( x ) dx (cid:27) + δ ∆ ( W )+ N X i =1 (cid:26) h λ d i , ∇ u i − d i i + r Z Ω k∇ u i ( x ) − d i ( x ) k dx (cid:27) + N X i =1 (cid:26) h λ w i , u i − w i i + r Z Ω | u i ( x ) − w i ( x ) | dx (cid:27) , where Λ D = ( λ d , ..., λ d N ) , Λ W = ( λ w , ..., λ w N ) are the Lagrangian multipliersand r is a positive constant. Here h λ d i , ∇ u i − d i i = R Ω λ Td i ( x )( ∇ u i ( x ) − d i ( x )) dx ,and h λ w i , u i − w i i = R Ω λ w i ( x )( u i ( x ) − w i ( x )) dx .The ADMM solves the primal variables one by one in the Gauss-Seidelstyle and then updates the dual variables (Lagrangian multipliers). It can besummarized as follows: D k +1 = arg min D L ( D , W k , C k , U k ; Λ k D , Λ k W ) , W k +1 = arg min W L ( D k +1 , W , C k , U k ; Λ k D , Λ k W ) , C k +1 = arg min C L ( D k +1 , W k +1 , C , U k ; Λ k D , Λ k W ) , U k +1 = arg min U L ( D k +1 , W k +1 , C k +1 , U ; Λ k D , Λ k W ) ,λ k +1 d i = λ kd i + r (cid:0) ∇ u k +1 i − d k +1 i (cid:1) ,λ k +1 w i = λ kw i + r (cid:0) u k +1 i − w k +1 i (cid:1) . In the following, we show how to solve each subproblem and then give thealgorithm. D -subproblemThe subproblem for D is equivalent to D k +1 = arg min D N X i =1 Z Ω k d i ( x ) k dx + r Z Ω (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d i ( x ) − ∇ u ki ( x ) − λ kd i ( x ) r (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) dx . This is separable and the optimal solution of d k +1 i can be explicitly expressedusing shrinkage operators. We simply compute d k +1 i ( x ) = S ∇ u ki ( x ) + λ kd i ( x ) r , r ! , where S is the shrinkage operator defined as S ( v, τ ) := v k v k ∗ max ( k v k − τ, . For the sake of simplicity, we denote this step as D k +1 = S (cid:18) ∇ U k + Λ k D r , r (cid:19) . W -subproblemThe subproblem for W is equivalent tomin W N X i =1 r Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) w i ( x ) − u ki ( x ) − λ kw i ( x ) r + λ | I ( x ) − c ki | r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dx + δ ∆ ( W ) . Since ∆ is a convex simplex at any x ∈ Ω , the solution is given by W k +1 = Π ∆ " u ki + λ kw i r − λ | I − c ki | r Ni =1 , where Π ∆ is the projection onto the simplex ∆ , for which more details can befound in [17]. We denote the step as W k +1 = Π ∆ (cid:18) U k + Λ k W r − λ | I − C k | r (cid:19) . C -subproblemThe subproblem for C is C k +1 = arg min C N X i =1 Z Ω | I ( x ) − c i | w k +1 i ( x ) dx. It is separable, and c k +1 i can be solved independently. The optimality conditionfor each c k +1 i is 0 ∈ − Z Ω ∂ | I ( x ) − c k +1 i | w k +1 i ( x ) dx. (9)Because the right-hand side of (9) is maximal monotone [3], the bisectionmethod and ADMM are applied to solve it [28,26]. The next lemma providesa new way to find an optimal solution for discrete cases. ultiphase Image Segmentation via Fuzzy and L1 11 Lemma 1
Given a finite non-decreasing sequence { I [ j ] } nj =1 , i.e., I [1] ≤ I [2] ≤ ... ≤ I [ n ] , and a non-negative sequence { w [ j ] } nj =1 with A = P nj =1 w [ j ] > , the optimalsolution set for min c n X j =1 (cid:12)(cid:12) I [ j ] − c (cid:12)(cid:12) w [ j ] , (10) is (cid:2) I [ j ∗ ] , I [ j ∗ +1] ] (cid:3) , where j ∗ and j ∗ satisfy A − j ∗ X j =1 w [ j ] ≤ < A − j ∗ − X j =1 w [ j ] ,A − j ∗ +1 X j =1 w [ j ] < ≤ A − j ∗ X j =1 w [ j ] . The fuzzy median of { I j } nj =1 with respect to the weight { w [ j ] } nj =1 [28], whichis defined as ( I [ j ∗ ] + I [ j ∗ +1] ]) / , is an optimal solution. If, in addition, thereexists j ∗ such that A − j ∗ X j =1 w [ j ] < < A − j ∗ − X j =1 w [ j ] , then j ∗ = j ∗ − , and the optimal solution is unique.Proof The optimality condition of (10) is0 ∈ h ( c ) := n X j =1 ∂ | I [ j ] − c | w [ j ] . We can see that h ( c ) is non-increasing and it can be expressed as h ( c ) = A, c < I [1] , [ A − w [1] , A ] , c = I [1] ,A − w [1] , c ∈ ( I [1] , I [2] ) , · · · [ A − P sj =1 w [ j ] , A − P s − j =1 w [ j ] ] , c = I [ s ] ,A − P sj =1 w [ j ] , c i ∈ ( I [ s ] , I [ s +1] ) , · · ·− A, c > I [ n ] . The range of h is [ − A, A ], and h ( c ) can take multiple values when c = I [ j ] forany j with w [ j ] = 0. Therefore, we can always find j ∗ such that A − j ∗ X j =1 w [ j ] ≤ < A − j ∗ − X j =1 w [ j ] . Thus 0 ∈ h ( I [ j ∗ ] ) = h A − P j ∗ j =1 w [ j ] , A − P j ∗ − j =1 w [ j ] i , and I [ j ∗ ] is an opti-mal solution. In addition, we have that h ( c ) > c < I [ j ∗ ] . Similarly, wecan find j ∗ such that A − j ∗ +1 X j =1 w [ j ] < ≤ A − j ∗ X j =1 w [ j ] , and I [ j ∗ +1] is an optimal solution. In addition, we have that h ( c ) < c > I [ j ∗ +1] . Then h ( c ) being non-increasing gives us that the set of optimalsolutions for (10) is (cid:2) I [ j ∗ ] , I [ j ∗ +1] (cid:3) . When j ∗ = j ∗ + 1, the optimal solution isunique. ⊓⊔ Remark:
When there are missing pixels in images, we can put a mask onthe data fidelity term as in image inpainting problems [41] or assign a valueto each missing pixel. In [24], the missing pixels are assigned with zero, and C changes based on the percentage of missing pixels. While, this lemma tells usthat assigning 0 or 255 (for grayscale images) randomly to these missing pixelswill not change the value of c with a high probability because c i is a weightedmedian. Also, by assigning 0 or 255 randomly, this algorithm is able to dealwith cases where more than half of the pixels are missing. See the numericalexperiments in Section 5.Assume that I and w k +1 i are vectors in R n , where n is the total numberof pixels. We can reorder the components of I and w k +1 i such that the mono-tonicity of I in Lemma 1 is satisfied. If there are multiple optimal solutionsfor c i , we choose the smallest one. We denote this step as C k +1 = ψ ( W k +1 ) . U -subproblemThe subproblem for U is equivalent to U k +1 = arg min U N X i =1 Z Ω (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∇ u i ( x ) − d k +1 i ( x ) + λ kd i ( x ) r (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) dx + Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u i ( x ) − w k +1 i ( x ) + λ kw i ( x ) r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dx. It is separable for u k +1 i , and, from the following first order optimality conditionfor each u k +1 i : ∇ T ∇ u k +1 i ( x ) − d k +1 i ( x ) + λ kd i ( x ) r ! + u k +1 i ( x ) − w k +1 i ( x ) + λ kw i ( x ) r ! = 0 , ultiphase Image Segmentation via Fuzzy and L1 13 we can derive the closed-form solution of u k +1 i from the equation: (cid:0) ∇ T ∇ + I (cid:1) u k +1 i ( x ) = ∇ T d k +1 i ( x ) + w k +1 i ( x ) − ∇ T λ kd i ( x ) r − λ kw i ( x ) r . A diagonalization trick can be applied to find u k +1 i efficiently [45]. We denotethis step as U k +1 = (cid:0) ∇ T ∇ + I (cid:1) − (cid:18) ∇ T D k +1 + W k +1 − ∇ T Λ k D r − Λ k W r (cid:19) . Updating dual variables.
We denote the steps as Λ k +1 D = Λ k D + r (cid:0) ∇ U k +1 − D k +1 (cid:1) ,Λ k +1 W = Λ k W + r (cid:0) U k +1 − W k +1 (cid:1) . Combining all these steps together, we obtain the proposed L1 Fuzzy Seg-mentation algorithm (L1FS) in Algorithm 1 for solving (4).
Algorithm 1
The proposed L1FS algorithm – Initialization: U and C are specified, Λ D = , Λ W = . – For k = 0 , , , · · · , repeat until the stopping criterion is reached. D k +1 = S (cid:18) ∇ U k + Λ k D r , r (cid:19) , W k +1 = Π ∆ (cid:18) U k + Λ k W r − λ | I − C k | r (cid:19) , C k +1 = ψ ( W k +1 ) , U k +1 = (cid:0) ∇ T ∇ + I (cid:1) − (cid:18) ∇ T D k +1 + W k +1 − ∇ T Λ k D r − Λ k W r (cid:19) ,Λ k +1 D = Λ k D + r (cid:0) ∇ U k +1 − D k +1 (cid:1) ,Λ k +1 W = Λ k W + r (cid:0) U k +1 − W k +1 (cid:1) . – Output: C k +1 , U k +1 . Remark:
Though there are four variables D , W , C and U , they can begrouped into two variables ( D , W ) and ( C , U ) and the subproblems for thesetwo variables can be decoupled into the four subproblems. So it is essentiallya two block ADMM applied on a nonconvex optimization problem.Because problem (4) is nonconvex, a good initial guess of U and C ,which can be obtained from other methods using fuzzy membership functions,is helpful for the stability of L1FS. Thus we update D and W first becauseboth of them can use the initial guess of U . C is given, problem (4) is convex. Then, the algorithm is a standard ADMMby considering ( D , W ) together, and its convergence is guaranteed [5]. In thissection, we give a convergence result of Algorithm 1 for unknown C . To sim-plify notations, let us define the sextuples: Z := ( D , W , C , U , Λ D , Λ W ) . A point Z is a KKT point of (8) if it satisfies the following KKT conditions: ∂ k d ∗ i ( x ) k − λ ∗ d i ( x ) ∋ , (11a) λ | I − C | − Λ W + ∂δ ∆ ( W ) ∋ , (11b) Z Ω ∂ | I ( x ) − c ∗ i | w ∗ i ( x ) dx ∋ , (11c) ∇ T λ d i ( x ) + λ w i ( x ) = 0 , (11d) ∇ u i ( x ) − d i ( x ) = 0 , (11e) u i ( x ) − w i ( x ) = 0 . (11f)where ∂ k · k and ∂δ ∆ ( · ) are subdifferentials of k · k and δ ∆ ( · ), respectively. Theorem 3 (Convergence to stationary points of L1FS)
Let (cid:8) Z k (cid:9) ∞ k =1 be a sequence generated by Algorithm 1 that satisfies the condition lim k →∞ (cid:0) Z k +1 − Z k (cid:1) = . Then any accumulation point of (cid:8) Z k (cid:9) ∞ k =1 is a KKT point of problem (8).Consequently, any accumulation point of (cid:0) C k , U k (cid:1) is a stationary point ofproblem (4).Proof From the updating formulas in Algorithm 1, we obtain the followinginequalities related to the subsequent iteration: D k +1 − D k = S (cid:18) ∇ U k + Λ k D r , r (cid:19) − D k , (12a) W k +1 − W k = Π ∆ (cid:18) U k + Λ k W r − λ | I − C k | r (cid:19) − W k , (12b) C k +1 − C k = ψ ( W k +1 ) − C k , (12c) (cid:0) ∇ T ∇ + I (cid:1) (cid:0) U k +1 − U k (cid:1) = ∇ T ( D k +1 − ∇ U k ) + ( W k +1 − U k ) (12d) − ∇ T Λ k D r − Λ k W r , (12e) Λ k +1 D − Λ k D = r (cid:0) ∇ U k +1 − D k +1 (cid:1) , (12f) Λ k +1 W − Λ k W = r (cid:0) U k +1 − W k +1 (cid:1) . (12g) ultiphase Image Segmentation via Fuzzy and L1 15 By the assumption lim k →∞ (cid:0) Z k +1 − Z k (cid:1) = , the left-hand side and right-hand side of each equality in (12) go to zero as k → ∞ . Then we have S (cid:18) ∇ U k + Λ k D r , r (cid:19) − D k → , (13a) Π ∆ (cid:18) U k + Λ k W r − λ | I − C k | r (cid:19) − W k → , (13b) ψ ( W k +1 ) − C k → , (13c) ∇ T ( D k +1 − ∇ U k ) + ( W k +1 − U k ) − ∇ T Λ k D r − Λ k W r → , (13d) ∇ U k +1 − D k +1 → , (13e) U k +1 − W k +1 → . (13f)Assume Z ∗ = ( D ∗ , W ∗ , C ∗ , U ∗ , Λ ∗ D , Λ ∗ W ) is an accumulation point of Z k .(13) gives us S (cid:18) ∇ U ∗ + Λ ∗ D r , r (cid:19) − D ∗ = , (14a) Π ∆ (cid:18) U ∗ + Λ ∗ W r − λ | I − C ∗ | r (cid:19) − W ∗ = , (14b) ψ ( W ∗ ) − C ∗ = , (14c) ∇ T Λ ∗ D + Λ ∗ W = , (14d) ∇ U ∗ − D ∗ = , (14e) U ∗ − W ∗ = . (14f)By (14a), it follows that D ∗ is a solution of the minimization problemmin D N X i =1 (Z Ω k d i ( x ) k dx + r Z Ω (cid:13)(cid:13)(cid:13)(cid:13) d i ( x ) − ∇ u ∗ i ( x ) − λ ∗ d i ( x ) r (cid:13)(cid:13)(cid:13)(cid:13) dx ) . Thus D ∗ satisfies the first order optimality condition0 ∈ ∂ k d ∗ i ( x ) k + r (cid:18) d ∗ i ( x ) − ∇ u ∗ i ( x ) − λ ∗ d i ( x ) r (cid:19) . Using (14e), we simplify the above condition as0 ∈ ∂ k d ∗ i ( x ) k − λ ∗ d i ( x ) . (15)By (14b), W ∗ is a solution of the following minimization problem:min W N X i =1 ( r Z Ω (cid:12)(cid:12)(cid:12)(cid:12) w i ( x ) − u ∗ i ( x ) − λ ∗ w i ( x ) r + λ | I ( x ) − c ∗ i | r (cid:12)(cid:12)(cid:12)(cid:12) dx ) + δ ∆ ( W ) . Hence W ∗ satisfies the optimality condition ∈ ∂δ ∆ ( W ∗ ) + r ( W ∗ − U ∗ ) − Λ ∗ W + λ | I − C ∗ | . Together with the equality in (14f), the above equation can be simplified as ∈ ∂δ ∆ ( W ∗ ) − Λ ∗ W + λ | I − C ∗ | . (16)The equation (14c) implies that C ∗ is a solution of the minimization problemmin C N X i =1 Z Ω | I ( x ) − c i | w ∗ i ( x ) dx, and the optimal condition is0 ∈ Z Ω ∂ | I ( x ) − c ∗ i | w ∗ i ( x ) dx. (17)The equivalence of equations (14a)-(14c) with (15)-(17), together withequations (14d)-(14f) shows that the accumulation point Z ∗ satisfies the KKTcondition in equations (11a)-(11f), thus Z ∗ is a KKT point of problem (8).Since problem (8) and problem (4) are equivalent, the convergence of (cid:0) C k , U k (cid:1) in the statement follows directly. ⊓⊔ The convergence analysis is motivated by [46]. This theorem implies thatwhenever (cid:8) Z k (cid:9) ∞ k =1 converges, it converges to a KKT point of problem (8).However, since the minimization problem (8) is nonconvex, the KKT point isnot necessary to be an optimal solution. In order to demonstrate the effectiveness of the proposed method, we compareour method with some existing methods on both synthetic and real images.These methods (FCM, FCM S2, FLICM, L2FS, L2L0, and L1SS) are discussedin Section 2. Since the segmentation models in these methods are not jointlyconvex, and they are easy to get stuck in local minima, “good” initializationis crucial for all algorithms, especially when the given image is noisy. ForFCM, FCM S2, and FLICM, the initial U is uniformly distributed in [0 ,
1] andnormalized to satisfy the sum-to-one constraint. While for TV based methodsL2FS and L1FS, one can also use the results of FCM and FCM S2 as the initial U and C . Here we consider three ways for choosing the initial U and C : theresult of FCM, the result of FCM S2, and U as functions uniformly distributedin [0 ,
1] and C as the result of FCM. In all the experiments, we choose the onewith the highest performance among all the three initializations. For L1SS, weuse the initialization method as described in the original paper. The stoppingcriterion, which is the same for all the methods except L1SS, is defined as k U k +1 − U k k k U k k < ǫ ultiphase Image Segmentation via Fuzzy and L1 17 where ǫ is a very small number. For L1SS, this stopping criterion does not worksince L1SS leads to contrast loss due to the error in calculating class centers { c i } Ni =1 at early iterations. However, the contrast of L1SS will be enhancedgradually if the number of iterations increases. To gain satisfactory results, wechoose to stop L1SS by setting the maximum number of iterations to be 1000.The parameters of the methods being compared are set as follows. In FCM,we set p = 2. In FCM S2, we set p = 2 , α = 5, and the window size for themedian filter as 5 ×
5. In FLICM, we set p = 2 and the window size of theneighborhoods as 3 × ×
5, which depends on the noise level. However, theweight parameter λ for L2FS, L1SS, and L1FS are tuned for each experimentto achieve optimal results. In all experiments, the range of λ for L2FS is[0 . , . . , . , . ǫ = 10 − for the two-phase segmentation and ǫ = 10 − for themultiphase segmentation. We use the default parameters in L2L0 except thatwe set λ = 10.The clustering results of all methods are obtained by checking the max-imum value of their membership functions. Then we display the recoveredpiecewise constant image as the final result. To compare segmentation resultsquantitatively, we consider the Segmentation Accuracy (SA) defined by SA = . (a) (b) (c) Fig. 1
Synthetic piecewise constant test images. (a) A two-phase grayscale image withintensities 20 and 128, size 128 × × × The synthetic piecewise constant test images are displayed in Fig. 1. Thenoisy images are contaminated by three types of noise: Gaussian Noise (GN)with zero mean and standard deviation varying from 10 to 80, Salt-and-PepperImpulse Noise (SPIN) with 10% to 60% pixels contaminated, and Random-Valued Impulse Noise (RVIN) with 10% and 60% pixels contaminated.5.1 Test on Fig. 1aTab. 1 provides the SA comparison of these six algorithms, and Figs. 2-4 showsome of the corresponding results.
Table 1
SA performance of different methods applied on Fig. 1a contaminated by dif-ferent levels of GN, SPIN, and RVIN.
GN ( σ ) 10 20 30 40 50 60 70 80FCM L1SS
SPIN (%) 10% 20% 30% 40% 50% 60%FCM 0.9480 0.8983 0.8478 0.7979 0.7486 0.6980FLICM 0.9984 0.9921 0.9738 0.8002 0.7313 0.6554L2FS 0.9999 0.9998 0.9982 0.9983 - -L1SS 0.9998 0.9990 0.9977 0.9967 0.9956 0.9953L1FS
RVIN (%) 10% 20% 30% 40% 50% 60%FCM S2 0.9985 0.9972 0.9945 0.9862 0.9630 0.9042FLICM 0.9987 0.9985 0.9970 0.9958 0.9948 0.9919L2FS (a) GN σ =30 (b) 0.9642 (c) 0.9996 (d) 1 (e) 1 (f) 1(g) GN σ =60 (h) 0.8162 (i) 0.9968 (j) 1 (k) 0.9984 (l) 1 Fig. 2
Two-phase segmentation on the synthetic image Fig. 1a with Gaussian noise. Firstcolumn: images contaminated by Gaussian noise with standard deviations 30 and 60, re-spectively; Second column to last column: results of FCM, FLICM, L2FS, L1SS, and L1FS,respectively. The SA values are reported below each segmentation result.
For GN, all methods being tested give correct segmentation results forsmall standard deviations (e.g., σ = 10). As the standard deviation increases,the SA value of FCM decays very fast. All the other algorithms have verylarge SA values even when the standard deviation is 80. L1FS has the bestperformance for all cases, and it is able to give correct segmentation resultseven when σ ≤
70. L2FS is the second best algorithm, and it is able to givecorrect segmentation results when σ ≤ σ = 30 and 60 are displayed in Fig. 2. Theresults of FCM are relatively “noisy” (the second column). For FLICM (the ultiphase Image Segmentation via Fuzzy and L1 19(a) SPIN 20% (b) 0.8983 (c) 0.9921 (d) 0.9998 (e) 0.9990 (f) 1(g) SPIN 40% (h) 0.7979 (i) 0.8002 (j) 0.9983 (k) 0.9967 (l) 1 Fig. 3
Two-phase segmentation on the synthetic image Fig. 1a with SPIN. First column:images contaminated by 20% and 40% SPIN, respectively; Second column to last column:results of FCM, FLICM, L2FS, L1SS, and L1FS, respectively. The SA values are reportedbelow each segmentation result.(a) RVIN 20% (b) 0.9972 (c) 0.9985 (d) 1 (e) 0.9995 (f) 1(g) RVIN 40% (h) 0.9862 (i) 0.9958 (j) 0.9998 (k) 0.9979 (l) 1
Fig. 4
Two-phase segmentation on the synthetic image Fig. 1a with RVIN. First column:images contaminated by 20% and 40% RVIN, respectively; Second column to last column:results of FCM S2, FLICM, L2FS, L1SS, and L1FS, respectively. The SA values are reportedbelow each segmentation result. third column) and L1SS (the fifth column), the segmentation error occurs onthe middle edge.For SPIN, Tab. 1 shows that FLICM performs much better than FCM fornoise levels 10%-30%. However, if the noise level is higher than 30%, bothFCM and FLICM have very poor performance. L1SS achieves much betterperformance than FLICM, even when the noise level is higher than 30%. L2FSis slightly better than L1SS. Meanwhile, if the noise level is higher than 50%,L1SS fails to give a satisfactory result. L1FS has the highest SA among allmethods. It gives completely correct segmentation results for noise levels 10%-50%. Fig. 3 shows the results of all methods for noise levels 20% and 40%, respectively. For L2FS and L1SS, the segmentation errors occur along themiddle edge. We also observe that for high noise levels such as 40%, both L2FSand L1SS suffer from slight contrast loss, e.g., Fig. 3j and Fig. 3k. However,L1FS is still able to preserve contrast, e.g., Fig. 3f and Fig. 3l.For RVIN, FCM S2 is used to initialize U and C for TV based methodsL2FS and L1FS. Tab. 1 shows that FCM S2 has the worst performance amongall. L1SS is much better than FCM S2 especially for high noise levels. FLICMperforms slightly better than FCM S2. L1FS achieves the best performancewhich is slightly better than L2FS. L2FS gives correct segmentation results atnoise levels 10%-30%, while our method L1FS can give the correct segmenta-tion results at noise levels 10%-50%. Fig. 4 shows the results for noise levels20% and 40%. Again, we find that most of the errors occur along the middleedge for FLICM and L1SS. Moreover, the results of L2FS in Fig. 4j, L1SS inFig. 4e and Fig. 4k lose some contrast.Next we analyze the contrast problem for TV based methods. For L2FS,the estimated intensity c i in each segmented region roughly equals the meanvalue of the intensities in that region. In the Gaussian noise case, the noisehas zero mean and therefore it has almost no impact on the estimation of c i .However, for both impulse noise cases, the noise has a significant influenceon the estimation of c i by taking the average. More specifically, assumingthat the impulse noise follows the uniform distribution, its impact on theestimation of c i is like this. Given an image with intensity range [0 , Ω i with true intensity β ≤ β than those with intensity less than β , then c i ≥ β after taking the average and vice versa. Hence, the range of the image will beshrunk by applying L2FS even when the segmentation is correct, and therebythe recovered image will suffer from contrast loss. Note that the contrast lossproblem has also been reported for the TVL2 restoration model in impulsenoise removal [12]. For L1SS algorithm, one step of ADMM is used to solvethe C -subproblem approximately. Thus c i is not accurate for the first fewiterations. However, since L1SS uses the L1-norm fidelity, the loss of contrastbecomes more and more subtle as the number of iterations increases. In L1FS,we solve the C -subproblem exactly. Thus, L1FS can preserve contrast well inthe segmentation process.5.2 Test on Fig. 1bThe performance comparison for the multiphase synthetic piecewise constantgray image Fig. 1b is shown in Tab. 2 and Fig. 5.As shown in Tab. 2, FCM performs poorly for GN, while L2FS and L1FSperform relatively well with similar SA. For the noise level σ = 10, both L2FSand L1FS give a correct segmentation result. For SPIN, FCM also gives theworst performance in terms of SA, while L1FS achieves the best performance.L2FS fails to yield a correct segmentation result when noise levels σ ≥
50. ForRVIN, FCM S2 achieves high SA values since it can smooth out some noise ultiphase Image Segmentation via Fuzzy and L1 21
Table 2
SA comparison of different methods applied on Fig. 1b with different levels of GN,SPIN, and RVIN.
GN ( σ ) 10 20 30 40 50 60 70 80FCM 0.9987 0.8191 0.6634 0.5849 0.5233 0.4718 0.4319 0.4017L2FS L1FS - -RVIN (%) 10% 20% 30% 40% - - - -FCM 0.9922 0.9809 0.96672 0.9248 - - - -L2FS 0.9949 0.9923 0.9880 0.9731 - - - -L1FS - - - -
Table 3
SA comparison of different methods applied on Fig. 1c with different levels of GN,SPIN, and RVIN.
GN ( σ ) 10 20 30 40 50 60 70 80FCM L2L0 - -RVIN (%) 10% 20% 30% 40% 50% 60% - -FCM 0.8971 0.7992 0.6967 0.6085 0.5196 0.4294 - -L2FS 0.9974 0.9957 0.9899 0.9853 0.9841 0.8988 - -L2L0 0.9971 0.9955 0.9906 0.9856 0.9727 0.9488 - -L1FS - - in the segmentation process. L2FS performs much better than FCM S2, andL1FS performs the best among all methods.Fig. 5 shows some results corresponding to Tab.2. The first row is the noisyimages being tested. The second row shows the results of FCM (Fig. 5g-5j)and FCM S2 (Fig. 5k-5l). Most of them looks very “noisy” except Fig. 5k.L2FS and L1FS give very clean results in the third row and last row, respec-tively. Comparing the results by L2FS and L1FS for Gaussian noise, bothof them have high visual qualities. However, for SPIN and RVIN, it is obvi-ous that L1FS preserves the contrast much better than L2FS and has bettersegmentation results. σ =30 (b) GN σ =60 (c) SPIN 20% (d) SPIN 40% (e) RVIN 20% (f) RVIN 40%(g) 0.6634 (h) 0.4718 (i) 0.8431 (j) 0.6847 (k) 0.9809 (l) 0.9248(m) 0.9994 (n) 0.9950 (o) 0.9877 (p) 0.9673 (q) 0.9923 (r) 0.9731(s) 0.9993 (t) 0.9950 (u) 0.9948 (v) 0.9894 (w) 0.9957 (x) 0.9868 Fig. 5
Multiphase segmentation on the synthetic image Fig. 1b with different levels of GN,SPIN, and RVIN. First row: images contaminated different types of noise with differentlevels. Second row to last row: results of FCM, L2FS and L1FS, respectively. The SA valuesare reported below each segmentation result. σ ≤ σ ≤
40. When σ ≥
50, the performance of FCM decreases rapidly,while L2FS, L2L0, and L1FS still achieve very large SA values. Note that weinitialize U randomly for GN in this test. For SPIN and RVIN, FCM has theworst performance which is far lower than that of the other three methods. Itis also obvious that L1FS outperforms L2L0 and L2FS.Fig. 6 shows some results corresponding to Tab. 3. The first row showsthe tested noisy images. The results of FCM in the second row seems to be ultiphase Image Segmentation via Fuzzy and L1 23GN σ =30 GN σ =60 SPIN 20% SPIN 40% RVIN 20% RVIN 40%0.9998 0.7772 0.7294 0.5092 0.7992 0.60851 0.9992 0.9925 0.9822 0.9957 0.98531 0.9991 0.9880 0.9740 0.9955 0.98561 0.9994 0.9937 0.9854 0.9963 0.9906 Fig. 6
Multiphase segmentation on the synthetic color image Fig. 1c with different levels ofGN, SPIN, and RVIN. First row: images contaminated different types of noise with differentlevels. Second row to last row: results of FCM, L2FS, L2L0, and L1FS, respectively. The SAvalues are reported below each segmentation result. “noisy” in most cases. The results of L2FS (the third row), L2L0 (the fourthrow), and L1FS (the last row) are very clean. However, in terms of contrast, itis obvious that L1FS outperforms L2L0 and L2FS particularly for SPIN andRVIN. C , and thereforethe color is not blue any more. However, in Fig. 8i, since the background isclean, the correct color can be obtained. To sum up, Fig. 8 demonstrates thatthe proposed L1FS gives smoother segmentations than FCM. This paper presents a novel piecewise constant image segmentation modelbased on fuzzy membership functions and L1-norm fidelity. ADMM is appliedto derive an efficient numerical algorithm, in which each subproblem has aclosed-form solution. In particular, an efficient algorithm is proposed to findthe fuzzy median. The numerical results on both synthetic and real piecewiseconstant images demonstrate that the proposed method is superior to someexisting state-of-the-art methods since it is more robust to impulse noise andcan preserve better contrast. Even in the case of Gaussian noise, the proposedmethod can achieve similar results as its L2-fidelity counterpart. In this work,we assume that the image to be dealt with is piecewise constant, which workswell on images with homogeneous image features. The future work is to extendthis framework to piecewise smooth image segmentation. ultiphase Image Segmentation via Fuzzy and L1 25(a) N=2 (b) (c) (d) N=4 (e) (f)(g) N=2 (h) (i) (j) N=6 (k) (l)(m) N=8 (n) (o) (p) N=4 (q) (r)(s) N=3 (t) (u) (v) N=6 (w) (x)
Fig. 7
Segmentation on real images. First column and fourth column: real color images;Second column and fifth column: results of FCM; Third column and last column: results ofL1FS.(a) (b) (c) (d) (e) (f)(g) (h) (i) (j) (k) (l)
Fig. 8
Detailed comparison of FCM and L1FS on the woman image. First row: six seg-mented regions of FCM. Second row: six segmented regions of L1FS.6 Fang Li et al.
Acknowledgment
The research of F. Li was supported by the 973 Program 2011CB707104, theresearch of S. Osher and J. Qin was supported by ONR Grants N00014120838and N00014140444, NSF Grants DMS-1118971 and CCF-0926127, and theKeck Foundation, the research of M. Yan was supported by NSF Grants DMS-1317602. The work was done when the first author was visiting UCLA Depart-ment of Mathematics.
References
1. Mohamed N Ahmed, Sameh M Yamany, Nevin Mohamed, Aly A Farag, and ThomasMoriarty. A modified fuzzy c-means algorithm for bias field estimation and segmentationof MRI data.
Medical Imaging, IEEE Transactions on , 21(3):193–199, 2002.2. Gilles Aubert and Pierre Kornprobst.
Mathematical problems in image processing:partial differential equations and the calculus of variations , volume 147. Springer, 2006.3. H.H. Bauschke and P.L. Combettes.
Convex Analysis and Monotone Operator Theoryin Hilbert Spaces . CMS Books in Mathematics. Springer New York, 2011.4. James C Bezdek, LO Hall, and LP Clarke. Review of MR image segmentation techniquesusing pattern recognition.
Medical Physics , 20(4):1033–1048, 1992.5. Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Dis-tributed optimization and statistical learning via the alternating direction method ofmultipliers.
Foundations and Trends R (cid:13) in Machine Learning , 3(1):1–122, 2011.6. Xavier Bresson, Selim Esedo¯glu, Pierre Vandergheynst, Jean-Philippe Thiran, and Stan-ley Osher. Fast global minimization of the active contour/snake model. Journal ofMathematical Imaging and vision , 28(2):151–167, 2007.7. Xiaohao Cai. Variational image segmentation model coupled with image restorationachievements.
Pattern Recognition , 48(6):2029–2042, 2015.8. Vicent Caselles, Ron Kimmel, and Guillermo Sapiro. Geodesic active contours.
Inter-national Journal of Computer Vision , 22(1):61–79, 1997.9. Antonin Chambolle, Daniel Cremers, and Thomas Pock. A convex approach to minimalpartitions.
SIAM Journal on Imaging Sciences , 5(4):1113–1158, 2012.10. Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convexproblems with applications to imaging.
Journal of Mathematical Imaging and Vision ,40(1):120–145, 2011.11. Raymond Chan, Hongfei Yang, and Tieyong Zeng. A two-stage image segmentationmethod for blurry images with poisson or multiplicative gamma noise.
SIAM Journalon Imaging Sciences , 7(1):98–127, 2014.12. Raymond H Chan, Chung-Wa Ho, and Mila Nikolova. Salt-and-pepper noise removalby median-type noise detectors and detail-preserving regularization.
Image Processing,IEEE Transactions on , 14(10):1479–1485, 2005.13. Tony F Chan and Selim Esedoglu. Aspects of total variation regularized L1 functionapproximation.
SIAM Journal on Applied Mathematics , 65(5):1817–1837, 2005.14. Tony F Chan, Selim Esedoglu, and Mila Nikolova. Algorithms for finding global mini-mizers of image segmentation and denoising models.
SIAM Journal on Applied Math-ematics , 66(5):1632–1648, 2006.15. Tony F Chan and Luminita A Vese. Active contours without edges.
Image Processing,IEEE transactions on , 10(2):266–277, 2001.16. Songcan Chen and Daoqiang Zhang. Robust image segmentation using FCM withspatial constraints based on new kernel-induced distance measure.
Systems, Man, andCybernetics, Part B: Cybernetics, IEEE Transactions on , 34(4):1907–1916, 2004.17. Yunmei Chen and Xiaojing Ye. Projection onto a simplex. arXiv preprintarXiv:1101.6081 , 2011.ultiphase Image Segmentation via Fuzzy and L1 2718. Anastasia Dubrovina, Guy Rosman, and Ron Kimmel. Multi-region active contourswith a single level set function.
Pattern Analysis and Machine Intelligence, IEEETransactions on , 37(8):1585-1601, 2015.19. Daniel Gabay and Bertrand Mercier. A dual algorithm for the solution of nonlinearvariational problems via finite element approximation.
Computers & Mathematics withApplications , 2(1):17–40, 1976.20. Roland Glowinski and A Marroco. Sur l’approximation, par ´el´ements finis d’ordreun, et la r´esolution, par p´enalisation-dualit´e d’une classe de probl`emes de dirichletnon lin´eaires.
ESAIM: Mathematical Modelling and Numerical Analysis-Mod´elisationMath´ematique et Analyse Num´erique , 9(R2):41–76, 1975.21. Tom Goldstein and Stanley Osher. The split bregman method for L1-regularized prob-lems.
SIAM Journal on Imaging Sciences , 2(2):323–343, 2009.22. Weihong Guo, Jing Qin, and Sibel Tari. Automatic prior shape selection for imagesegmentation.
Research in Shape Modeling , 2015.23. Xiaoxia Guo, Fang Li, and Michael K Ng. A fast ℓ -TV algorithm for image restoration. SIAM Journal on Scientific Computing , 31(3):2322–2341, 2009.24. Yanyan He, M Yousuff Hussaini, Jianwei Ma, Behrang Shafei, and Gabriele Steidl. Anew fuzzy c-means method with total variation regularization for segmentation of imageswith noisy and incomplete data.
Pattern Recognition , 45(9):3463–3471, 2012.25. Nawal Houhou, J Thiran, and Xavier Bresson. Fast texture segmentation model basedon the shape operator and active contour. In
Computer Vision and Pattern Recognition,2008. CVPR 2008. IEEE Conference on , pages 1–8. IEEE, 2008.26. Miyoun Jung, Myeongmin Kang, and Myungjoo Kang. Variational image segmentationmodels involving non-smooth data-fidelity terms.
Journal of Scientific Computing ,59(2):277–308, 2014.27. Yoon Mo Jung, Sung Ha Kang, and Jianhong Shen. Multiphase image segmentation viaModica-Mortola phase transition.
SIAM Journal on Applied Mathematics , 67(5):1213–1232, 2007.28. P.R. Kersten. Fuzzy order statistics and their application to fuzzy clustering.
FuzzySystems, IEEE Transactions on , 7(6):708–712, 1999.29. Stelios Krinidis and Vassilios Chatzis. A robust fuzzy local information c-means clus-tering algorithm.
Image Processing, IEEE Transactions on , 19(5):1328–1337, 2010.30. Fang Li, Michael K Ng, Tie Yong Zeng, and Chunli Shen. A multiphase image segmen-tation method based on fuzzy region competition.
SIAM Journal on Imaging Sciences ,3(3):277–299, 2010.31. Fang Li, Chaomin Shen, and Chunming Li. Multiphase soft segmentation with totalvariation and H1 regularization.
Journal of Mathematical Imaging and Vision , 37(2):98–111, 2010.32. Johan Lie, Marius Lysaker, and Xue-Cheng Tai. A piecewise constant level set frame-work.
International Journal of Numerical Analysis and Modeling , 2(4):422–438, 2005.33. Johan Lie, Marius Lysaker, and Xue-Cheng Tai. A binary level set model and some ap-plications to Mumford-Shah image segmentation.
Image Processing, IEEE Transactionson , 15(5):1171–1181, 2006.34. Benoit Mory and Roberto Ardon. Fuzzy region competition: a convex two-phase seg-mentation framework. In
Scale Space and Variational Methods in Computer Vision ,pages 214–226. Springer, 2007.35. David Mumford and Jayant Shah. Optimal approximations by piecewise smooth func-tions and associated variational problems.
Communications on Pure and Applied Mthe-matics , 42(5):577–685, 1989.36. Mila Nikolova. A variational approach to remove outliers and impulse noise.
Journalof Mathematical Imaging and Vision , 20(1-2):99–120, 2004.37. Stanley Osher and James A Sethian. Fronts propagating with curvature-dependentspeed: algorithms based on Hamilton-Jacobi formulations.
Journal of ComputationalPhysics , 79(1):12–49, 1988.38. Dzung L Pham. Fuzzy clustering with spatial constraints. In
Image Processing. 2002.Proceedings. 2002 International Conference on , volume 2, pages II–65. IEEE, 2002.39. Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noiseremoval algorithms.
Physica D: Nonlinear Phenomena , 60(1):259–268, 1992.8 Fang Li et al.40. Alex Sawatzky, Daniel Tenbrinck, Xiaoyi Jiang, and Martin Burger. A variationalframework for region-based segmentation incorporating physical noise models.
Journalof Mathematical Imaging and Vision , 47(3):179–209, 2013.41. Jianhong Shen and Tony F Chan. Mathematical models for local nontexture inpaintings.
SIAM Journal on Applied Mathematics , 62(3):1019–1043, 2002.42. Jianhong Jackie Shen. A stochastic-variational model for soft Mumford-Shah segmen-tation.
International Journal of Biomedical Imaging , 2006:92329, 2006.43. Martin Storath and Andreas Weinmann. Fast partitioning of vector-valued images.
SIAM Journal on Imaging Sciences , 7(3):1826–1852, 2014.44. Luminita A Vese and Tony F Chan. A multiphase level set framework for image seg-mentation using the Mumford and Shah model.
International Journal of ComputerVision , 50(3):271–293, 2002.45. Yilun Wang, Junfeng Yang, Wotao Yin, and Yin Zhang. A new alternating minimizationalgorithm for total variation image reconstruction.
SIAM Journal on Imaging Sciences ,1(3):248–272, 2008.46. Yangyang Xu, Wotao Yin, Zaiwen Wen, and Yin Zhang. An alternating directionalgorithm for matrix completion with nonnegative factors.
Frontiers of Mathematics inChina , 7(2):365–384, 2012.47. Ming Yan. Restoration of images corrupted by impulse noise and mixed Gaussianimpulse noise using blind inpainting.
SIAM J. Imaging Sciences , 6(3):1227–1245, 2013.48. Ming Yan, Alex A. T. Bui, Jason Cong, and Luminita A. Vese. General convergent ex-pectation maximization (EM)-type algorithms for image reconstruction.
Inverse Prob-lems and Imaging , 7(3):1007–1029, 2013.49. Song Chun Zhu and Alan Yuille. Region competition: Unifying snakes, region growing,and Bayes/MDL for multiband image segmentation.