A Multiphase Image Segmentation Based on Fuzzy Membership Functions and L1-norm Fidelity

A Multiphase Image Segmentation Based on Fuzzy Membership Functions and L1-norm Fidelity thanks:

Fang Li F. Li (fi) Department of Mathematics, East China Normal University, and Shanghai Key Laboratory of Pure Mathematics and Mathematical Practice, Shanghai, China.
22email: e-mail: fli@math.ecnu.edu.cn
S. Osher, J. Qin Department of Mathematics, University of California, Los Angeles, CA 90095, USA.
33email: sjo@math.ucla.edu,jxq@ucla.edu
and M. Yan Department of Computational Mathematics, Science and Engineering and Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA.
44email: e-mail:yanm@math.msu.edu
   Stanley Osher F. Li (fi) Department of Mathematics, East China Normal University, and Shanghai Key Laboratory of Pure Mathematics and Mathematical Practice, Shanghai, China.
22email: e-mail: fli@math.ecnu.edu.cn
S. Osher, J. Qin Department of Mathematics, University of California, Los Angeles, CA 90095, USA.
33email: sjo@math.ucla.edu,jxq@ucla.edu
and M. Yan Department of Computational Mathematics, Science and Engineering and Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA.
44email: e-mail:yanm@math.msu.edu
   Jing Qin F. Li (fi) Department of Mathematics, East China Normal University, and Shanghai Key Laboratory of Pure Mathematics and Mathematical Practice, Shanghai, China.
22email: e-mail: fli@math.ecnu.edu.cn
S. Osher, J. Qin Department of Mathematics, University of California, Los Angeles, CA 90095, USA.
33email: sjo@math.ucla.edu,jxq@ucla.edu
and M. Yan Department of Computational Mathematics, Science and Engineering and Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA.
44email: e-mail:yanm@math.msu.edu
   Ming Yan F. Li (fi) Department of Mathematics, East China Normal University, and Shanghai Key Laboratory of Pure Mathematics and Mathematical Practice, Shanghai, China.
22email: e-mail: fli@math.ecnu.edu.cn
S. Osher, J. Qin Department of Mathematics, University of California, Los Angeles, CA 90095, USA.
33email: sjo@math.ucla.edu,jxq@ucla.edu
and M. Yan Department of Computational Mathematics, Science and Engineering and Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA.
44email: e-mail:yanm@math.msu.edu
Received: date / Accepted: date
Abstract

In this paper, we propose a variational multiphase image segmentation model based on fuzzy membership functions and L1-norm 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 L1-norm based method is more robust to outliers such as impulse noise and keeps better contrast than its L2-norm counterpart. Theoretically, we prove the existence of the minimizer and analyze the convergence of the algorithm.

Keywords:
Image segmentation; fuzzy membership function; L1-norm; ADMM; segmentation accuracy.
journal: Journal of Scientific Computing

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 Mumford-Shah model MS (), a well-known variational segmentation model, penalizes the combination of the total length of the segmentation boundaries and the L2-norm error of approximating the observed image with an unknown piecewise smooth approximation. In other words, the Mumford-Shah model seeks an optimal piecewise smooth function with smooth boundaries to approximate the observed image.

However, the Mumford-Shah 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 Mumford-Shah 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) sum-to-one 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 two-phase 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 two-phase segmentation since more variables and constraints are involved in representing multiple regions and their boundaries effectively. The two-phase Chan-Vese 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 sum-to-one 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 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. Primal-dual 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 (), two-stage 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 L2-norm fidelity is adopted. However, when images are corrupted by non-Gaussian 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 L1-norm fidelity is used for both salt-and-pepper impulse noise and random-valued 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 L1-norm models is determined by the geometry such as area and length rather than the contrast in the L2-norm case chan2005aspects ().

Inspired by the fact that L1-norm 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 L1-norm 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 salt-and-pepper 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 closed-form 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 L1-norm fidelity is considered in the proposed method, while L2-norm 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 c-means 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 C-Means clustering method that solves

    where is a neighborhood of . FLICM is more robust to both Gaussian noise and impulse noise than FCM.

  • L2FS TVFCM () – L2-norm 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 () – L1-norm 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 four-phase model as follows:

    where the membership functions , , are represented by soft functions in the following way:

  • L2L0 storath2014fast () – L2-norm fidelity and L0-norm 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 L2-norm fidelity to measure the approximation error when the noise is the Gaussian white noise, we use the L1-norm fidelity. Same as in the Mumford-Shah 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 sum-to-one 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 semi-norm,

(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

Finally, since , a.e. and as , Fatou’s lemma yields

(7)

Combining inequalities (6) and (7) for all , on a suitable subsequence, we established

and hence must be a minimizer of the energy . This completes the proof.∎

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 one-to-one 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 Gauss-Seidel 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 right-hand 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 non-decreasing sequence , i.e.,

and a non-negative 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 non-increasing 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 non-increasing 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 closed-form 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).

Initialization: and are specified, . For , repeat until the stopping criterion is reached. Output: .
Algorithm 1 The proposed L1FS algorithm

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)

Let be a sequence generated by Algorithm 1 that satisfies the condition

Then any accumulation point of is a KKT point of problem (8). Consequently, any accumulation point of is a stationary point of problem (4).

Proof

From the updating formulas in Algorithm 1, we obtain the following inequalities related to the subsequent iteration:

(12a)
(12b)
(12c)
(12d)