A Multiphase Image Segmentation Based on Fuzzy Membership Functions and L1norm Fidelity ^{†}^{†}thanks:
Abstract
In this paper, we propose a variational multiphase image segmentation model based on fuzzy membership functions and L1norm fidelity. Then we apply the alternating direction method of multipliers to solve an equivalent problem. All the subproblems can be solved efficiently. Specifically, we propose a fast method to calculate the fuzzy median. Experimental results and comparisons show that the L1norm based method is more robust to outliers such as impulse noise and keeps better contrast than its L2norm counterpart. Theoretically, we prove the existence of the minimizer and analyze the convergence of the algorithm.
Keywords:
Image segmentation; fuzzy membership function; L1norm; ADMM; segmentation accuracy.∎
1 Introduction
As a fundamental step in image processing, image segmentation partitions an image into several disjoint regions such that pixels in the same region share some uniform characteristics such as intensity, color, and texture. During the last two decades, image segmentation methods using variational methods and partial differential equations are very popular due to their flexibility in problem modeling and algorithm design. There are two key ingredients of variational segmentation models. One is how to describe the regions or boundaries between these regions, and the other is how to model the noise and describe the uniform characteristics of each region.
The MumfordShah model MS (), a wellknown variational segmentation model, penalizes the combination of the total length of the segmentation boundaries and the L2norm error of approximating the observed image with an unknown piecewise smooth approximation. In other words, the MumfordShah model seeks an optimal piecewise smooth function with smooth boundaries to approximate the observed image.
However, the MumfordShah model is hard to implement in practice because the discretization of the unknown set of boundaries is very complex. Therefore, a parametric/explicit active contour method is used to represent the segmentation boundaries zhu1996region (). In addition, the special MumfordShah model with a piecewise constant approximation is studied by Chan and Vese CV (), and the level set method osher1988fronts () is applied to solve this problem. Using an implicit representation of boundaries, the level set method has many advantages over methods using explicit representations of boundaries. For instance, the level set method handles the topological change of zero level set automatically caselles1997geodesic (); aubert2006mathematical (); Qin2014 (). Both the parametric/explicit active contour method and the level set method assume that each pixel belongs to a unique region. An alternative way to represent various regions is to use fuzzy membership functions chan2006algorithms (); bresson2007fast (); mory2007fuzzy (); houhou2008fast (), which describe the levels of possible membership. Fuzzy membership functions assume that each pixel can be in several regions simultaneously with probability in [0,1]. These probabilities satisfy two constraints: (i) nonnegativity constraint, i.e., the membership functions are nonnegative at all pixels; (ii) sumtoone constraint, i.e., the sum of all the membership functions at each pixel equals one. Then the length of boundaries can be approximated by the Total Variation (TV) of fuzzy membership functions. The main advantages of using fuzzy membership functions over other methods include: i) the energy functional is convex with respect to fuzzy membership functions, guaranteeing the convergence and stability of many numerical optimization methods. ii) fuzzy membership function has a larger feasible set, and it is possible to find better segmentation results.
For twophase segmentation, where there are only two regions, we only need one level set function or one fuzzy membership function. Multiphase segmentation is more challenging than twophase segmentation since more variables and constraints are involved in representing multiple regions and their boundaries effectively. The twophase ChanVese model CV () has been generalized to multiphase segmentation by using multiple level set functions to represent multiple regions vese2002multiphase (). Partitioning an image into disjoint regions requires level set functions. The advantage of using multiple level set functions is that it automatically avoids the problems of vacuum and overlap of regions. However, the implementation is not easy, and special numerical schemes are needed to ensure stability PCLSM (); lie2006binary (); softMS (); jung2007multiphase (). For fuzzy membership functions, the sumtoone constraint is not satisfied automatically in multiphase segmentation. However, this constraint is easy to deal with in many cases, e.g., Fuzzy CMean (FCM) and its variants have closedform solutions for the membership functions and are widely used in data mining and medical image segmentation FCM (); ahmed2002modified (); pham2002fuzzy (); FCMS2 (); FLICM (); li2010softseg (); li2010competition (); TVFCM (); cai2015variational (). Variable splitting schemes are used in both li2010competition () and li2010softseg () to get efficient numerical algorithms. The Alternating Direction Method of Multipliers (ADMM) method is used in TVFCM () to derive the algorithm with two sets of extra variables. Primaldual type algorithm is derived in CP () to solve the TV regularized FCM segmentation model. Both TVFCM () and CP () use projection to simple to handle the constraints of membership functions. Other segmentation approaches include a convex approach chambolle2012convex (), twostage methods chan2014two (); storath2014fast (), one single level set function approach dubrovina2015multi (), et.al.
Noise is unavoidable in images, and it is important to develop segmentation methods that work on noisy images. Among many types of noise, the Gaussian white noise is frequently assumed, and the L2norm fidelity is adopted. However, when images are corrupted by nonGaussian noise, in order to obtain a faithful segmentation, one has to model the noise according to its specific type sawatzky2013variational (); nikolova2004variational (); chan2005salt (); guo2009fast (); Yan13 (); Yan13b (); cai2015variational (). Particularly, the L1norm fidelity is used for both saltandpepper impulse noise and randomvalued impulse noise in image processing nikolova2004variational (); chan2005salt (); guo2009fast (). In addition, it is robust to outliers and able to preserve contrast because the denoising process of L1norm models is determined by the geometry such as area and length rather than the contrast in the L2norm case chan2005aspects ().
Inspired by the fact that L1norm is more robust to impulse noise and outliers and can better preserve contrast, in this paper, we propose a variational multiphase fuzzy segmentation model based on L1norm fidelity and fuzzy membership functions. This model can also deal with missing data in images. When there are missing pixels in an image, we randomly assign 0 or 255 at these pixels by considering these pixels as corrupted by saltandpepper impulse noise. ADMM ADMM1 (); ADMM2 (), which was rediscovered as split Bregman SB () and found to be very useful for L1 and TV type optimization problems, is applied to solve this nonconvex optimization problem. By introducing two sets of auxiliary variables, we derive an efficient algorithm with all the subproblems having closedform solutions. In the theoretical aspect, we prove the existence of the minimizer and analyze the convergence of the algorithm. We note that the proposed method is closely related to the method in TVFCM () since both methods use TV regularization and ADMM algorithm. The difference is that L1norm fidelity is considered in the proposed method, while L2norm fidelity is used in TVFCM ().
The outline of this paper is as follows. In Section 2, we review some closely related 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, a numerical algorithm is derived, and its convergence analysis is presented. In Section 5, experimental results and comparisons are presented to illustrate the effectiveness of the proposed method. Finally, the paper is concluded in Section 6.
2 Related works
Let be a bounded open subset with Lipschitz boundary, and let be the given clean or noisy image. Let for grayscale images and for color images. Our goal is to find an phase “optimal” partition such that for all and . Define the set of phase fuzzy membership functions as
where is the space of functions with bounded variation aubert2006mathematical (). The closely related works are listed in the following and will be compared with our proposed method in Section 5. For the sake of simplicity, we use the notations and , where for grayscale images and for color images.

FCM FCM ()– Fuzzy cmeans clustering method that solves
using the alternating minimization algorithm. Though can be any number no smaller than , it is commonly set to 2.

FCM_S2 FCMS2 () – A variant of FCM that solves
where is obtained by applying the median filter on and is a weight parameter. It can also be solved by the alternating minimization algorithm, and it is more robust to impulse noise than FCM.

FLICM FLICM () – Fuzzy Local Information CMeans clustering method that solves
where is a neighborhood of . FLICM is more robust to both Gaussian noise and impulse noise than FCM.

L2FS TVFCM () – L2norm fidelity based Fuzzy Segmentation method that solves
(1) by ADMM. Here is a parameter and denotes the TV of with for . For fixed , TVFCM () is related to the popular TV denoising method ROF (). Note that the similar model is solved by other fast numerical methods in li2010competition ().

L1SS TVSEGL1 () – L1norm fidelity based Soft Segmentation method, in which soft functions are introduced to represent phases. Since the model for multiphase segmentation is complicated for more than four phases, we show the fourphase model as follows:
where the membership functions , , are represented by soft functions in the following way:

L2L0 storath2014fast () – L2norm fidelity and L0norm regularization based image partition model:
This model seeks a piecewise constant approximation of the original image . Since this model can not specify the number of classes, a second step is applied to combine some classes if this model returns more classes than required. Here we apply FCM on the piecewise constant approximation 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 fuzzy membership functions are required for an phase segmentation. Hence, the solution space is much larger than L1SS, which ensures the higher possibility to obtain optimal segmentation. Secondly, the proposed method can take use of other commonly used segmentation methods such as FCM to gain good initialization of fuzzy membership functions. Multiphase segmentation is sensitive to initialization, and a good initialization is very important for a successful segmentation. However, it is hard to use the existing segmentation methods to get a good initialization for soft membership functions in L1SS.
3 The proposed model
In this paper, we assume that the given image can be approximated by a piecewise constant function, i.e.,
Here denotes the indicator function of region , i.e.,
is a constant that represents the given data in region , and can be outliers, impulse noise or other mixture types rather than Gaussian noise.
Instead of using the L2norm fidelity to measure the approximation error when the noise is the Gaussian white noise, we use the L1norm fidelity. Same as in the MumfordShah model, we require the segmentation boundaries to be smooth. Then we have the following model
(2) 
where is a tuning parameter. Note that the TV of in the first term is equal to the length of boundary . An equivalent form of model (2) is
(3) 
Because can take values 0 and 1 only and is a partition, belongs to the set
At any point , there is only one function having value 1, and all the other functions have value 0. Thus set is not continuous, which leads to difficulties and instabilities in numerical implementations. However, we can relax binary indicator functions to fuzzy membership functions , which satisfy the nonnegativity constraint and the sumtoone constraint, i.e., belongs to the set defined in (2). It is obvious that and is a simplex at any . So can be understood as the probability of pixel to be in the th class. Then we can rewrite our model (3) as
(4) 
Note that model (4) is convex with respect to and separately but not jointly. 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 in (4) is stated in Theorem 3.1.
Theorem 3.1 (Existence of a minimizer)
Given an image , there exists a minimizer of in for any parameter .
Proof
Since is positive and proper, the infimum of must be finite. Let us pick a minimizing sequence that satisfies as . Then there exists a constant such that
Then each term in is bounded, i.e., for each ,
(5) 
Since , we have , where is the area of . Together with the first equality in (5), we have that is uniformly bounded in for all . By the compactness property of and the relative compactness of in , up to a subsequence also denoted by after relabeling, there exists a function such that
Then by the lower semicontinuity of the TV seminorm,
(6) 
Meanwhile since , we have .
It is easy to derive the first order optimality condition about , which is
where is the subdifferential of . Since and , the above equation implies that the constant satisfies
By the boundedness of sequence , we can extract a subsequence also denoted by such that for some constant
The minimizer of is not unique due to the following hidden symmetry property. Denote as the permutation group of , i.e., each permutation is defined as a onetoone map such that is a rearrangement of . Denote , . It is straightforward to show the following theorem.
Theorem 3.2 (Symmetry of minimizer)
For any and any ,
In particular, suppose that is a minimizer of , i.e.,
Then, for any , we have
4 The numerical algorithm and its convergence analysis
In this section, we provide an efficient algorithm based on ADMM and discuss its convergence.
4.1 The algorithm
ADMM is applied, in this subsection, to solve the proposed fuzzy segmentation model (4). We introduce two sets of auxiliary variables such that . Then the model (4) is equivalent to the following minimization problem with equality constraints:
(8) 
where is the indicator function of the set , i.e.,
The augmented Lagrangian function for problem (8) is:
where are the Lagrangian multipliers and is a positive constant. Here , and .
The ADMM solves the primal variables one by one in the GaussSeidel style and then updates the dual variables (Lagrangian multipliers). It can be summarized as follows:
In the following, we show how to solve each subproblem and then give the algorithm.
subproblem
The subproblem for is equivalent to
This is separable and the optimal solution of can be explicitly expressed using shrinkage operators. We simply compute
where is the shrinkage operator defined as
For the sake of simplicity, we denote this step as
subproblem
The subproblem for is equivalent to
Since is a convex simplex at any , the solution is given by
where is the projection onto the simplex , for which more details can be found in chen2011projection (). We denote the step as
subproblem
The subproblem for is
It is separable, and can be solved independently. The optimality condition for each is
(9) 
Because the righthand side of (9) is maximal monotone bauschke2011convex (), the bisection method and ADMM are applied to solve it fuzzy_median (); TVSEGL1 (). The next lemma provides a new way to find an optimal solution for discrete cases.
Lemma 1
Given a finite nondecreasing sequence , i.e.,
and a nonnegative sequence with , the optimal solution set for
(10) 
is , where and satisfy
The fuzzy median of with respect to the weight fuzzy_median (), which is defined as , is an optimal solution. If, in addition, there exists such that
then , and the optimal solution is unique.
Proof
The optimality condition of (10) is
We can see that is nonincreasing and it can be expressed as
The range of is , and can take multiple values when for any with . Therefore, we can always find such that
Thus , and is an optimal solution. In addition, we have that when . Similarly, we can find such that
and is an optimal solution. In addition, we have that when . Then being nonincreasing gives us that the set of optimal solutions for (10) is . When , the optimal solution is unique.∎
Remark: When there are missing pixels in images, we can put a mask on the data fidelity term as in image inpainting problems shen2002mathematical () or assign a value to each missing pixel. In TVFCM (), the missing pixels are assigned with zero, and changes based on the percentage of missing pixels. While, this lemma tells us that assigning 0 or 255 (for grayscale images) randomly to these missing pixels will not change the value of with a high probability because is a weighted median. Also, by assigning 0 or 255 randomly, this algorithm is able to deal with cases where more than half of the pixels are missing. See the numerical experiments in Section 5.
Assume that and are vectors in , where is the total number of pixels. We can reorder the components of and such that the monotonicity of in Lemma 1 is satisfied. If there are multiple optimal solutions for , we choose the smallest one. We denote this step as
subproblem
The subproblem for is equivalent to
It is separable for , and, from the following first order optimality condition for each :
we can derive the closedform solution of from the equation:
A diagonalization trick can be applied to find efficiently wang2008new (). We denote this step as
Updating dual variables. We denote the steps as
Combining all these steps together, we obtain the proposed L1 Fuzzy Segmentation algorithm (L1FS) in Algorithm 1 for solving (4).
Remark: Though there are four variables , , and , they can be grouped into two variables and and the subproblems for these two variables can be decoupled into the four subproblems. So it is essentially a two block ADMM applied on a nonconvex optimization problem.
Because problem (4) is nonconvex, a good initial guess of and , which can be obtained from other methods using fuzzy membership functions, is helpful for the stability of L1FS. Thus we update and first because both of them can use the initial guess of .
4.2 Convergence analysis
If is given, problem (4) is convex. Then, the algorithm is a standard ADMM by considering together, and its convergence is guaranteed ADMM (). In this section, we give a convergence result of Algorithm 1 for unknown . To simplify notations, let us define the sextuples:
A point is a KKT point of (8) if it satisfies the following KKT conditions:
(11a)  
(11b)  
(11c)  
(11d)  
(11e)  
(11f) 
where and are subdifferentials of and , respectively.
Theorem 4.1 (Convergence to stationary points of L1FS)
Proof
From the updating formulas in Algorithm 1, we obtain the following inequalities related to the subsequent iteration:
(12a)  
(12b)  
(12c)  
(12d)  