Linkage between piecewise constant Mumford-Shah model and ROF model and its virtue in image segmentation
Xiaohao Cai Raymond Chan Carola-Bibiane Schonlieb Gabriele Steidl Tieyong Zeng


Thermonuclear explosions of rapidly differentially rotating
  white dwarfs: Candidates for superluminous Type Ia supernovae?
    M. Fink1,4
  M. Kromer2,3
  W. Hillebrandt4
  F. K. Röpke2,3
  R. Pakmor3
  I. R. Seitenzahl5,6
  S. A. Sim7
    Received 22 May 2018 / Accepted 17 July 2018
============================================================================================================================




The piecewise constant Mumford-Shah (PCMS) model and the Rudin-Osher-Fatemi (ROF) model are two important variational models
in image segmentation and image restoration, respectively.
In this paper, we explore a linkage between these models. 
We prove that for the two-phase segmentation problem a partial minimizer of the PCMS model 
can be obtained by thresholding the minimizer of the ROF model.
A similar linkage is still valid for multiphase segmentation under specific assumptions.
Thus it opens a new segmentation paradigm: image segmentation can be done via image restoration plus thresholding.
This new paradigm, which circumvents the innate non-convex property of the PCMS model, therefore improves
the segmentation performance in both efficiency (much faster than state-of-the-art methods based on PCMS model,
particularly when the phase number is high)
and effectiveness (producing segmentation results with better quality) due to the flexibility of the ROF model in tackling degraded images,
such as noisy images, blurry images or images with information loss. As a by-product
of the new paradigm, we derive a novel segmentation method,
called thresholded-ROF (T-ROF) method, to illustrate the virtue of managing image segmentation through image restoration techniques.
The convergence of the T-ROF method is proved, and elaborate experimental results and 
comparisons are presented.



Image segmentation, image restoration, Mumford-Shah model, piecewise constant Mumford-Shah model, Chan-Vese model, total variation ROF model, thresholding.







§ INTRODUCTION


Image segmentation aims to group parts of a given image with similar characteristics together, while
image restoration intends to remove image degradations such as noise, blur or occlusions.



The piecewise constant Mumford-Shah (PCMS) model (nonconvex, a special case of the Mumford-Shah model <cit.>)
and the Rudin-Osher-Fatemi (ROF) model (convex, <cit.>) are two of the most famous variational models in the
research areas of image segmentation and restoration, respectively.
Following earlier works e.g. <cit.>, in this paper, we show a linkage between the PCMS and ROF models,
which gives rise to a new image segmentation paradigm:
manipulating image segmentation through image restoration plus thresholding.

Let us first recall the PCMS and ROF models. 
Throughout this paper, let Ω⊂ℝ^2 be a bounded, open set with Lipschitz boundary,
and f: Ω→ [0,1] be a given (degraded) image.  
In 1989 Mumford and Shah  <cit.> proposed solving segmentation problems by minimizing over
Γ⊂Ω and u ∈ H^1(Ω\Γ) the energy functional

    E_ MS(u, Γ; Ω)= H^1(Γ)+ λ^'∫_Ω\Γ|∇ u|^2 dx+ λ∫_Ω(u-f)^2 dx,   λ^', λ > 0,

where H^1 denotes the one-dimensional Hausdorff measure in ℝ^2. 
The functional E_ MS contains three terms:
the penalty term on the length of Γ, the H^1 semi-norm that enforces  the smoothness of u in Ω\Γ,
and the data fidelity term controlling the distance of u to the given image f.
Related approaches in a spatially discrete setting were proposed in <cit.>.
An early attempt to solve the challenging task of finding a minimizer of the non-convex and non-smooth Mumford-Shah 
functional (<ref>) was done by approximating it using a sequence of simpler elliptic problems, see <cit.>
for the so-called Ambrosio-Tortorelli approximation.
Many approaches to simplify model (<ref>) were meanwhile proposed in the literature, for example, in <cit.>, a
convex relaxation of the model was suggested. 
Another important simplification is to restrict its solution to 
be piecewise constant, which leads to the so-called PCMS model. 

PCMS model.
This model is based on the restriction ∇ u = 0 on Ω\Γ,
which results in 

    E_ PCMS(u, Γ; Ω)= H^1(Γ) + λ∫_Ω(u-f)^2 dx.

Assuming that Ω=⋃_i=0^K-1Ω_i with pairwise disjoint sets Ω_i and constant functions u(x) ≡ m_i 
on Ω_i, i = 0,… K-1, model (<ref>) can be rewritten as

    E_ PCMS(Ω, m)=
    1/2∑_i=0^K-1 Per(Ω_i; Ω)
    + λ∑_i=0^K-1∫_Ω_i(m_i-f)^2dx,

where Ω := {Ω_i}_i=0^K-1, m := {m_i}_i=0^K-1, and 
Per(Ω_i;Ω) denotes the perimeter of Ω_i in Ω.
If the number of phases is two, i.e. K=2, the PCMS model is the model of the active contours without edges
(Chan-Vese model) <cit.>, 

    E_ CV(Ω_1,m_0,m_1) =  Per(Ω_1;Ω) +
    λ(∫_Ω_1
    (m_1-f)^2   dx + ∫_Ω\Ω_1 (m_0-f)^2   dx).

In <cit.> the authors proposed to solve (<ref>), where
it can easily get stuck in local minima. To overcome this drawback, a convex relaxation
approach was proposed in <cit.>.
More precisely, it was shown that a global minimizer
of E_ CV(·, m_0,m_1) for fixed m_0, m_1 can be found by solving

    u̅ = u∈ BV(Ω) argmin{TV(u) + λ∫_Ω( (m_0-f)^2 - (m_1-f)^2  ) u  dx },

and setting Ω_1 := {x ∈Ω: u̅(x) > ρ} for any choice of ρ∈ [0, 1),
see also <cit.>. Note that the first term of (<ref>) is known as the total variation (TV) 
and the space BV is the space of functions of bounded variation, see Section 2 for the definition.
In other words, (<ref>) is a tight relaxation of the Chan-Vese model with fixed m_0 and m_1.
For the convex formulation of the full model (<ref>), see <cit.>.

There are many other approaches for two-phase image segmentation based on the Chan-Vese model and its convex version,
see e.g. <cit.>. 
In particular, a  hybrid level set method was proposed
in <cit.>, which replaces the first term of (<ref>) by a boundary feature map and the data fidelity terms
in (<ref>) by the difference between the given image f and a fixed threshold chosen by a
user or a specialist. Method <cit.> was used in medical image segmentation.
However, since every time it needs the user to choose a proper threshold for its model, 
it is not automatic and thus its applications are restricted. 
In <cit.>, the TV term
of (<ref>) was replaced by a weighted TV term which
helps the new model to capture much more important geometric properties.
In <cit.>, the TV term of (<ref>) was replaced by
a wavelet frame decomposition operator which, similar to the model in
<cit.>,  can also capture important geometric
properties. 
Nevertheless, for its solution u, no similar conclusions as the ones in <cit.> can be 
addressed; that is, there is no theory to support that its segmentation result Ω_1  = {x: u(x) > ρ} for ρ∈ [0, 1) 
is a solution as to some kind of objective functional. In <cit.>, the Chan-Vese model 
was extended for 3D biopores segmentation in tomographic images. 

In <cit.>, Chan and Vese proposed a multiphase segmentation model based on the
PCMS model using level sets. However, this method can also  get stuck easily in local minima.
Convex (non-tight) relaxation approaches for the PCMS model were proposed, which are
basically focusing on solving

    min_m_i, u_i∈ [0,1]{∑_i=0^K-1∫_Ω |∇ u_i|dx + λ∑_i=0^K-1∫_Ω(m_i-f)^2 u_idx },
       s.t.  ∑_i=0^K-1 u_i =1.

For more details along this line, refer e.g. to <cit.> and the references therein.
We are interested in a relation between the PCMS model and the  ROF model for 
 image restoration which we introduce next.

ROF model.
In 1992, Rudin, Osher and Fatemi <cit.> proposed the  variational model

    min_u∈ BV(Ω){TV(u) + μ/2∫_Ω(u-f)^2dx },   μ >0.

which has been studied extensively
in the literature, see e.g. <cit.> and references therein.

Actually, a subtle connection between image segmentation and image restoration has been raised in <cit.>.
In detail,  a two-stage image segmentation method is proposed – smoothing and thresholding (SaT) method –
which finds the solution of a convex variant of the Mumford-Shah model  in the first stage followed by a thresholding step in the second one.
The convex minimization functional in the first stage (the smoothing stage)
is the ROF functional (<ref>) plus an additional smoothing term ∫_Ω | ∇ u|^2   dx.
The SaT method is very efficient and flexible: it performs excellently for degraded images
(e.g. noisy and blurry images and images with information loss);
the minimizer from the first stage is unique; and one can change the number of phases K without solving
the minimization functional again. The success of the SaT method indicates a new methodology 
for image segmentation: first smoothing and then thresholding. This approach was extended in <cit.> 
for images corrupted by Poisson and Gamma noises, and in <cit.> to degraded color images.

Our contribution.
In this paper we highlight a relationship between the PCMS model (<ref>) and the ROF model (<ref>).
We prove that thresholding the minimizer of the ROF model leads to a partial minimizer (cf. definition (<ref>)) 
of the PCMS model when K=2
(Chan-Vese model (<ref>)), which remains true under specific assumptions when K>2.
This linkage between the PCMS model and the ROF model validates the effectiveness of our proposed SaT method
in <cit.> for image segmentation. 
Due to the significance of the PCMS model and ROF model, respectively in image
segmentation and image restoration, this linkage  bridges  to some extent these two research areas
and  might serve as a motivation to improve and design better methods.
A direct  benefit is a new efficient segmentation method
– thresholded-ROF (T-ROF) method – proposed in this paper. 
The T-ROF method exactly follows the new paradigm to perform image segmentation through image restoration 
plus iterative thresholding, where these thresholds are selected automatically following certain rules.
This appears to be more sophisticated
than the SaT method <cit.> which is based on K-means.
We emphasize that we just need to
solve the ROF model once, and our method gives optimal segmentation results akin to the PCMS model.
We prove the convergence of our T-ROF method regarding thresholds automatic selection.

On the one hand, the T-ROF method can be regarded as a special case of our proposed SaT method. However, it is directly obtained
from the linkage between the PCMS model and the ROF model discovered in this paper and thus is more theoretically justified.
Moreover, the strategy of choosing the thresholds automatically and optimally in the T-ROF method is not covered in
the SaT method in <cit.>. The strategy makes our T-ROF method more effective particularly for degraded images whose phases
have close intensities. On the other hand, the T-ROF method inherits the advantages of the SaT method – fast speed
and computational cost independent of the required number of phases K. In contrast, methods solving the 
PCMS model become computational demanding as the required number of phases increases.

Numerical experiments and detailed comparisons to the start-of-the-art methods
are presented to demonstrate the great performance of the proposed T-ROF method.

Partial results of this paper have been  presented in the conference paper <cit.>. 

The paper is organized as follows. In Section <ref>, we introduce the basic notation.
The linkage between the PCMS and ROF models is presented in Section <ref>. 
In Section <ref>,
we present our T-ROF model and provide an algorithm to solve it together with its
convergence analysis. 
In Section <ref>,
we demonstrate the performance of our T-ROF method on various synthetic and real-world images 
and compare it with 
the representative related segmentation methods <cit.>.
Conclusions are given in Section <ref>.




§ BASIC NOTATION


We briefly introduce the basic notation which will be used in the followings, see <cit.> for more details.
By BV(Ω) we denote the space of functions of bounded variation defined on Ω,  i.e.,
the Banach space of functions u:Ω→ [0,1]
with finite norm u_BV := u_L^1(Ω) + TV(u), where

    TV(u) := sup{∫_Ω u(x)   divφ  dx: 
    φ∈ C_c^1(Ω,ℝ^2), φ_∞≤ 1 }.

The distributional first order derivative Du of u
is a vector-valued Radon measure having total mass 
|Du| (Ω) = ∫_Ω |Du| dx = TV(u). In particular, we have for u ∈ W^1,1(Ω) that
Du = ∇ u ∈ L^1(Ω) so that in this case
TV(u) = ∫_Ω |∇ u|   dx.

 A “set" is understood as a Lebesgue measurable
set in ℝ^2, where we mainly
consider equivalence classes of sets which are equal up to Lebesgue measure zero.
By |A| we denote the Lebesgue measure of a set A.
For a Lebesgue measurable set A ⊂Ω, the perimeter of A in Ω is
defined by

    Per( A;Ω) := TV(χ_A),

where χ_A is the characteristic function of A. Hence A is of finite perimeter
if its characteristic function has finite bounded total variation. If A has a C^1 boundary,
then Per( A;Ω) coincides with H^1(∂ A ∩Ω).
For A,B ⊆Ω the relation

    Per(A ∪ B;Ω) +  Per(A ∩ B;Ω) ≤ Per(A;Ω) +  Per(B;Ω)

holds true. 
We will also use the notation Per( A; Ω̂) for non open set Ω̂, 
in the sense
Per( A; Ω̂) =  Per( A;  int(Ω̂)),
where int(Ω̂) denotes the interior of Ω̂. 

We define the mean of f on A ⊂ℝ^2 by

    mean_f (A) :=
    {[ 1/|A|∫_A f   dx,     if  |A| > 0,;               0,       otherwise. ].





§ LINKAGES


We first propose our  T-ROF model, and then use it to
derive a linkage between the PCMS and ROF models using the T-ROF model.




 §.§ Thresholded-ROF (T-ROF) Model


To motivate our T-ROF model, we start by considering 
for fixed τ∈ (0,1)[Note that E(∅,τ) = 0
and E(Ω,τ) = μ∫_Ω (τ - f)   dx.
Since f maps into [0,1], the global minimizer of E(·, τ)
for fixed τ≤ 0 is Ω and for fixed τ≥ 1 is ∅.
Therefore we restrict our attention to τ∈ (0,1).], 
the minimization of

    E(Σ, τ) :=  Per(Σ;Ω) + μ∫_Σ (τ - f)   dx.

The following proposition gives a way to solve it.


 
For any fixed τ∈ (0,1), the minimizer Σ_τ of E(·,τ) in (<ref>)
can be found by solving the convex minimization problem

    u̅ = u∈ BV(Ω) argmin{ TV(u) + μ∫_Ω (τ-f)   u   dx }

and then setting
Σ_τ = { x ∈Ω: u̅(x) > ρ} for any ρ∈ [0,1).



For a proof of the proposition, we refer to Proposition 2.1 in the review paper <cit.>, where
the proof uses the same ideas as in <cit.>.
The functional (<ref>) is convex and it is well known
that there exists a global minimizer. 
Hence the proposition ensures the existence
of a global minimizer of (<ref>). 
Moreover, based on 
the following Lemma from <cit.> and
a smoothness argument, an explanation that the minimizing
set Σ_τ is unique was given in <cit.>.


 
For fixed 0 < τ_1 < τ_2 < 1,
let Σ_i = Σ_τ_i 
be minimizers of (<ref>), i=1,2.
Then
|Σ_2\Σ_1| = 0
is fulfilled, i.e.,
Σ_1⊇Σ_2 up to a measure zero set.



The following proposition gives another way of solving (<ref>) via the ROF function.
For the proof see <cit.>.


 
The set Σ_τ := {x ∈Ω : u(x) > τ} solves E(·;τ) in (<ref>) for every τ∈ (0,1)
if and only if
the function u ∈ BV(Ω) solves the ROF model (<ref>), i.e.,

    min_u∈ BV(Ω){TV(u) + μ/2∫_Ω(u-f)^2dx },   μ >0.




Using Lemma <ref>, after minimizing (<ref>)
for 0 < τ_1 < τ_2 < ⋯ < τ_K-1< 1, we have for the corresponding sets

    Ω⊇Σ_1⊇Σ_2⊇⋯⊇Σ_K-1⊇∅.


Setting Σ_0 := Ω and Σ_K := ∅,
we see that the sets

    Ω_i:= Σ_i\Σ_i+1,    i=0,…,K-1

are pairwise disjoint and fulfill ⋃_i=0^K-1Ω_i = Ω.

Let Σ := {Σ_i }_i =1^K-1, τ := {τ_i }_i =1^K-1 (τ_i < τ_j, i < j), and

    E(Σ, τ) := ∑_i=1^K-1(  Per(Σ_i;Ω)
    + μ∫_Σ_i (τ_i - f)   dx ),   μ > 0.

Our  T-ROF (thresholded-ROF) model aims to find a pair (Σ^*, τ^*) fulfilling the two conditions

    E(Σ^*, τ^*) ≤ E(Σ, τ^*)   for  all  Σ⊂Ω^K-1,

and

    τ_i^* = 1/2(m_i-1^* + m_i^*),    i=1,…,K-1,

where

    m_i^* :=  mean_f (Ω_i^*),   Ω_i^*:=Σ_i^* \Σ_i+1^*,  Σ_0^* := Ω,  Σ_K^* := ∅.

Since E(·,·) in (<ref>) is separable in each i and, for each i, it is precisely of the form of  (<ref>), a minimizer of
E(·,τ^*) in (<ref>) for some fixed τ^* can be found easily by componentwise minimization for each i 
using Proposition <ref>, i.e., thresholding
the minimizer of the ROF model (<ref>) with τ^*.
However, finding a pair of (Σ^*, τ^*) satisfying (<ref>) as well as the condition (<ref>)
is not straightforward. In Section <ref>, we will propose an efficient way to address the solution of the T-ROF
model (<ref>) that satisfies the condition (<ref>).
Once we have obtained (Σ^*, τ^*), the desired segmentation for the given image f is then given by Ω_i^*, i=0,…,K-1.
Finally, note that finding a pair (Σ, τ^*) which solves the T-ROF model differs from solving

    min_Σ, τ E(Σ, τ),     τ_i = 1/2(m_i-1 + m_i),    i=1,…,K-1,

since we do not minimize over all feasible τ. For an example see Remark 1 in <cit.>.




 §.§ Linkage of PCMS and ROF Models
 

Recall that 
(Σ^*, m^*) is a partial minimizer of some objective function E if

    E(Σ^*, m^*) ≤ E(Σ^*, m),     m, 
    
    E(Σ^*, m^*) ≤ E(Σ, m^*),    Σ.

We see immediately that a partial minimizer 
Ω^*= {Ω_i^*}_i=0^K-1, 
m^* = {m_i^*}_i=0^K-1
of the PCMS model (<ref>) has to fulfill

    m_i^* =  mean_f( Ω_i^*),    i=0,…,K-1.

We note that if E is differentiable on its domain, then every partial minimizer contained
in the interior of the domain is stationary, see e.g. <cit.>.
Note also that a partial minimizer must not be a local minimizer and conversely, see Fig. <ref>.






Case K=2.
We first give the relationship between our T-ROF model  and the
PCMS model for K=2, and then use it to derive
the relationship between the ROF and PCMS models.


 
(Relation between T-ROF and PCMS models for K=2)
For K=2, let (Σ_1^*,τ_1^*) with  0 < |Σ_1^*| < |Ω| 
be a solution of the T-ROF model (<ref>)–(<ref>). Then
(Σ_1^*,m_0^*,m_1^* ) is a partial minimizer of the PCMS
model (<ref>) with the parameter λ := μ/2(m_1^* - m_0^*).




Since E(Σ_1^*, τ_1^*)≤ E(∅, τ_1^*) = 0, we conclude
∫_Σ_1^* (τ_1^* - f)  dx < 0 which implies that 
τ_1^* <  mean_f(Σ_1^*) = m_1^*.
Similarly, since E(Σ_1^*, τ_1^*)≤ E(Ω, τ_1^*), we see that

    0 <  Per(Σ_1^*;Ω)  ≤μ∫_Ω\Σ_1^* (τ_1^* - f)   dx

and consequently m_0^* =  mean_f(Ω\Σ_1^*)≤τ_1^*. Therefore m_0^* < m_1^*.

Clearly, the set Σ_1^* is also a minimizer of E(·, τ_1^*) + C 
with the constant C := λ∫_Ω (m_0^* -f)^2   dx.
Setting τ_1^* := m_1^*+m_0^*/2, we obtain

    E(Σ, τ_1^*)  + C
    	   =
    	 Per(Σ;Ω) +μ∫_Σ (τ_1^* - f)   dx + C 
       =
    	  Per(Σ;Ω) + μ/2(m_1^* - m_0^*)∫_Σ[(m_1^* -f)^2  -  (m_0^* - f)^2 ]  dx  + C 
       =
    	 Per(Σ;Ω) + λ( ∫_Σ (m_1^* -f)^2 dx + ∫_Ω\Σ (m_0^* -f)^2 dx ).

By the definition of m_i, i=0,1 and (<ref>)
we get the assertion.




Since f∈ [0, 1], we have 0< m_1^* - m_0^* ≤1. The parameter λ = μ/2(m_1^* - m_0^*) in the
Chan-Vese model (<ref>) is larger than μ and increases dramatically if
(m_1^* - m_0^*) becomes smaller. Hence, this λ is adapted to
the difference between m_1^* and m_0^* and penalizes the data term more in the Chan-Vese model if this difference becomes smaller.
We now know that when (m_1^* - m_0^*) is very small, λ used in the Chan-Vese model should be large; however,
in practice, to solve the Chan-Vese model, λ is given beforehand with no knowledge about this kind of information. 
It is therefore very hard if not impossible for the Chan-Vese model to be given a good value of λ to obtain a high quality 
segmentation result.
In contrast, it is much easier for our  T-ROF model to get good results by just tuning the threshold τ_1^* 
(automatically, see Section <ref>), 
no matter how large or small the difference between m_1^* and m_0^* is.



Next we give the relationship between the ROF model  and the
PCMS model for K=2.


  (Relation between ROF and PCMS models for K=2)
Let K=2 and u^* ∈ BV(Ω) solve the ROF model (<ref>).
For given 0<m_0 < m_1 ≤ 1, let 
Σ̃ := {x ∈Ω : u^*(x) > m_1+m_0/2} fulfill 
0 < |Σ̃| < |Ω|.
 Then Σ̃ is a  minimizer of the PCMS model (<ref>) for
 λ := μ/2(m_1 - m_0) and fixed m_0, m_1.
 In particular, (Σ̃, m_0, m_1) is a partial minimizer of (<ref>)
 if m_0 =  mean_f(Ω\Σ̃) and m_1 =  mean_f(Σ̃).




Following (<ref>), we have for all m_0 < m_1, λ = μ/2(m_1 - m_0) and Σ⊂Ω that

    E(Σ, m_1+m_0/2 ) + λ∫_Ω (m_0 -f)^2   dx =
    E_ CV(Σ,m_0,m_1).

By Proposition <ref>, Σ̃ minimizes
E(·, m_1+m_0/2 ) and hence also E_ CV(·,m_0,m_1).
If m̃_0 =  mean_f(Ω\Σ̃) and m̃_1 =  mean_f(Σ̃), then
0 < m̃_0 < m̃_1 (cf. the first part of the proof of Theorem <ref>);
and by (<ref>) they minimize E_ CV(Σ̃, ·, ·).




Case K>2.
Now we consider K>2. For Σ in (<ref>) and {Ω_i}_i=0^K-1
in (<ref>), we know Ω_i∪⋯∪Ω_K-1 = Σ_i so that

    ∑_i=0^K-1 Per(Ω_i∪⋯∪Ω_K-1; Ω) =
    ∑_i=0^K-1 Per(Σ_i; Ω).

If int(Σ_i)⊃Σ̅_i+1, 
that is the closure of set Σ_i+1 is inside the interior of Σ_i and therefore 
∂Σ_i∩∂Σ_i+1 = ∅, then

    ∑_i=0^K-1 Per(Ω_i∪⋯∪Ω_K-1; Ω) =
    ∑_i=0^K-1 Per(Σ_i; Ω) =
    1/2∑_i=0^K-1 Per(Ω_i; Ω).

We consider the following variant (extension) of the PCMS model (<ref>)

    min_Ω_i, m_i{∑_i=0^K-1 Per(Ω_i∪⋯∪Ω_K-1; Ω)
    + ∑_i=0^K-1μ̃_i∫_Ω_i (m_i-f)^2dx },

where μ̃_i > 0 are regularization parameters.
Let us call model (<ref>) the PCMS-V model.
The relationship between the T-ROF model (<ref>)–(<ref>) and the PCMS-V model (<ref>)
is given in the following Theorem <ref>.


  (Relation between the T-ROF and PCMS-V models for K>2)
For K>2, let ( {Σ^*_i }_i=1^K-1, {τ_i^*}_i=1^K-1) 
be a solution of the T-ROF model
(<ref>)–(<ref>)  with m^*_i < m^*_i+1, i=0,…,K-1.
Then ({Ω_i^*}_i=0^K-1, {m_i^*}_i=0^K-1) is a partial minimizer of the
PCMS-V model (<ref>) with regularization parameters defined as

    μ̃_i =
    μ/2(m_1^*-m_0^*),     i=0, 
    μ/2(m_i^*-m_i-1^*)+μ/2(m_i+1^*-m_i^*),     i=1, …, K-2,
    μ/2(m_K-1^*-m_K-2^*),     i=K-1,

where {Ω_i^*}_i=0^K-1
is obtained by (<ref>) with {Σ^*_i }_i=1^K-1.



When minimizing E(Σ, τ^*) (<ref>) with respect to
Σ, i.e.,

    min_Σ_i{ Per(Σ_i, Ω)+ μ∫_Σ_i (τ_i^*-f)dx },
         i=1, …, K-1,

it is clear that finding Σ_i^* is independent with finding Σ_j^* for j≠ i.
Hence, Σ_1^*, …, Σ_i-1^*, Σ_i+1^*, …, Σ_K-1^*
can be regarded as fixed when finding Σ_i^*.
Note that, we also have
Σ_1^*⊇Σ_2^* ⊇⋯⊇Σ_K-1^* using Lemma <ref>.
From Theorem <ref>, we know the minimizer of
(<ref>) for each i is a partial minimizer of

    min_Σ_i,m_i-1, m_i{ Per(Σ_i; Ω) +
    μ/2(m_i^*-m_i-1^*)(∫_Σ_i
    (m_i-f)^2dx + ∫_Ω\Σ_i (m_i-1-f)^2dx) },

which is equivalent to

    min_Ω_i,m_i-1, m_i   { Per(Ω_i∪Ω_i+1^*∪⋯∪Ω_K-1^*; Ω) 
        + μ/2(m_i^*-m_i-1^*)(∫_Ω_i (m_i-f)^2dx
     + ∫_Ω_i-1 (m_i-1-f)^2dx) }.

The proof is completed by summing up the above objective functions for all i = 1, …, K-1.




Note that for the standard PCMS model in (<ref>), the regularization parameter μ is fixed. In contrast,
its variation PCMS-V model in (<ref>) derived from our T-ROF model (<ref>) has more flexible regularization
parameters. This kind of resetting for regularization parameters will avail the PCMS-V model (and our T-ROF method which is given in the
next section) for multiphase segmentation particularly for images containing phases with close intensities.
We demonstrate this fact in experimental results later.


In summary, we conclude that, when K=2, the T-ROF model gives a segmentation result of the PCMS model
for fixed λ and m_i, i=0, 1 using the way of defining λ in Theorem <ref>, and
vice versa. When K>2, Theorem <ref> tells us that, if ∂Σ_i∩∂Σ_i+1 = ∅, i=1, …, K-1,
then the T-ROF model gives a segmentation result of the PCMS model with μ redefined as in (<ref>);
otherwise, only an approximation to the PCMS model can be achieved by the T-ROF model. 
It is worth mentioning that the case of ∂Σ_i∩∂Σ_i+1≠∅ is due to jumps in the ROF solution,
which are a subset of the jumps of the original image f, see <cit.> for more details. 

Even though the equivalence between the PCMS model and the T-ROF model cannot hold anymore 
when K>2 and ∂Σ_i∩∂Σ_i+1≠∅,  
the way of quantifying, reducing and/or finally overcoming the gap between them will be of great interest for future research. 
Note, importantly, that the lack of equivalence for K>2  does not mean that the T-ROF model performs poorer than the PCMS model. 
In the experimental results, we will show that 
in some cases the T-ROF method is actually much better than the state-of-the-art methods based on the PCMS model.




§ T-ROF ALGORITHM AND ITS CONVERGENCE
 

Proposition <ref> implies that we can
obtain a minimizer Σ of E(·,τ) in (<ref>)
by minimizing the ROF functional and subsequently thresholding the minimizing function by τ_i, i=1,…,K-1.
This method is particularly efficient since the minimizer of the ROF functional remains the same and thus just need
to be solved once when we apply various thresholds τ^ (k).
Here, at iteration k, when we have (Σ^(k), τ^(k)), we use the following rule to
obtain τ^(k+1):

    m_i^(k) =  mean_f (Σ_i^(k)\Σ_i+1^(k)), τ_i^(k+1) = 1/2(m_i-1^(k) + m_i^(k)),   i=1,…,K-1,

where we assume the ordering
0 < τ_1^(k) <  ⋯ < τ_K-1^(k) < 1 and 

    |Ω_i^(k)| = |Σ_i^(k)\Σ_i+1^(k)| > 0,     i = 0, …, K-1.

We abbreviate the rule in (<ref>) by 

    τ^(k+1)= Φ(Σ^(k), τ^(k)).



If the above criterion (<ref>)  is not satisfied at step k, 
for example  |Ω_i^(k)| = 0, i∈{0, …, K-1}, then
we remove τ_i^(k) from the 
threshold sequence {τ_i^(k)}_i =1^K-1, which will not affect much of the segmentation result practically when not considering 
measure zero set. A new threshold sequence {τ_i^(k)}_i =1^K̂-1, K̂<K can be formed 
which is considered in (<ref>) now.

For a set  Ω̂⊂Ω and Ω̂⊃Σ≠∅, where Σ is the minimizer of E_Ω(·,τ)
in (<ref>), we want that segmenting  
Ω̂ into Σ and Ω̂\Σ 
does not produce a worse segmentation result than the
naive segmentation of Ω̂ into the set itself and an empty set. In other words, we require 

    E_Ω̂(Σ,τ) ≤ E_Ω̂(Ω̂,τ) =  μ(τ -  mean_f (Ω̂) ) |Ω̂|,
    	  
    	E_Ω̂(Σ,τ) ≤ E_Ω̂(∅,τ) = 0.

If the requirement (<ref>) is not fulfilled, this might imply that the difference between Σ and
Ω̂\Σ is too subtle to be distinguished and  we keep Ω̂ rather than segmenting it into 
Σ and Ω̂\Σ.


 
For i = 1, …, K-2, let Σ_i≠Ω and Σ_i+1≠∅ be the minimizers of E(·,τ_i) and
E(·,τ_i+1) for  0< τ_i < τ_i+1 < 1 appearing in the T-ROF model, where |Σ_i\Σ_i+1| > 0.
If the inequalities in (<ref>) hold, then 

    τ_i ≤ mean_f (Σ_i\Σ_i+1) ≤τ_i+1.



 By definition we have Σ_i⊃Σ_i+1.
Let Ω̂ = Σ_i and Σ = Σ_i+1. 
Then using the first inequality in (<ref>), we have

    Per(Σ_i+1;Σ_i) + μ∫_Σ_i+1 (τ_i+1 - f)   dx 
       ≤μ∫_Σ_i (τ_i+1 -f)   dx,

i.e.,

    0 ≤ Per(Σ_i+1;Σ_i)  ≤μ∫_Σ_i\Σ_i+1 (τ_i+1 - f)   dx

which implies 
mean_f ( Σ_i\Σ_i+1) ≤τ_i+1.

Let Ω̂ = Ω\Σ_i+1 and 
Σ = Σ_i\Σ_i+1.  
Using the second inequality in (<ref>), we have

    Per(Σ_i\Σ_i+1; Ω\Σ_i+1) 
    	+ μ∫_Σ_i\Σ_i+1 (τ_i - f)   dx ≤ 0,

which implies with |Σ_i\Σ_i+1| > 0 that

    μ∫_Σ_i\Σ_i+1 (τ_i - f)   dx ≤ 0,

i.e.,
τ_i ≤ mean_f (Σ_i\Σ_i+1).
This completes the proof.


In practice, for Σ_i⊃Σ_i+1 (0< τ_i < τ_i+1 < 1), if 
Σ_i\Σ_i+1 contains different texture from Σ_i+1, 
then it is highly likely that segmenting Σ_i into Σ_i\Σ_i+1
and Σ_i+1 is much better than the naive segmentation of Σ_i.
This gives us a useful criterion to design our T-ROF algorithm: at some step, if the naive segmentation of Σ_i is better
than segmenting it into Σ_i\Σ_i+1
and Σ_i+1, then 
removing Σ_i+1 (τ_i+1) from the set (threshold) sequence is more meaningful and 
will not influence the segmentation result much.

In sum, after obtaining (Σ^(k), τ^(k)) at iteration k, we apply 
the criterions in (<ref>) and (<ref>) to remove measure zero phases and ignore
unnecessary segmentation, respectively. To achieve this, we can particularly keep checking the following steps i)–iii) until they are all satisfied: 
i) if (Σ_i^(k), τ_i^(k)) and (Σ_i+1^(k), τ_i+1^(k)), 
0≤ i ≤ K-1, do not fulfil
the criterion  (<ref>), then remove (Σ_i+1^(k), τ_i+1^(k)) from 
(Σ^(k), τ^(k)), and form a new (Σ^(k), τ^(k)) by reordering the indices, and update K; 
ii) if (Σ_i^(k), τ_i^(k)), 1≤ i ≤ K-1, and one of its neighbors in (Σ^(k-1), τ^(k-1)) 
say (Σ_j^(k-1), τ_j^(k-1)), j ∈{i-1, i, i+1}, do not fulfil
the criterion  (<ref>), then replace the values of (Σ_i^(k), τ_i^(k)) by (Σ_j^(k-1), τ_j^(k-1));
and iii) if (Σ_i^(k), τ_i^(k)) and (Σ_i+1^(k), τ_i+1^(k)), 0≤ i ≤ K-1, do not fulfil
the criterion  (<ref>), then remove (Σ_i+1^(k), τ_i+1^(k)) from 
(Σ^(k), τ^(k)), and form a new (Σ^(k), τ^(k)) by reordering the indices, and update K.
We abbreviate the above steps i)–iii) conducting the criterions (<ref>) and (<ref>) by 

    (Σ^(k), τ^(k))  :=  C (Σ^(k), τ^(k)).


By combining the criterion in (<ref>) and  the rule of updating τ^(k) by (<ref>),
we obtain our T-ROF segmentation method in Algorithm <ref>, to find a solution of the T-ROF model (<ref>)–(<ref>).





From  Algorithm <ref>, we see that  the T-ROF method exactly follows the new paradigm which is: performing image segmentation
through image restoration plus a thresholding. The T-ROF method can therefore be regarded as a special case of our proposed
SaT method in <cit.>. Particularly, in addition to the theoretical bare bones in the T-ROF method and the fast speed it inherits from
the standard SaT method, the iterative way of selecting optimal thresholds – the strategy in (<ref>) – in the T-ROF 
method makes it more effective in multiphase segmentation when compared with the standard SaT method which uses the K-means to 
select thresholds. This advantage is more notable when segmenting images with phases which have close intensities.

There are many efficient methods to solve the ROF model, for example
the primal-dual algorithm <cit.>, alternating direction method with multipliers (ADMM) <cit.>,
or the split-Bregman algorithm <cit.>. In this paper,
we use the ADMM to solve the ROF model (<ref>).
The convergence property of Algorithm <ref> is discussed in the next section.




 §.§ Convergence Proof of the T-ROF Algorithm
  

We now prove the convergence of our T-ROF algorithm – Algorithm <ref>.


 
Let A,B,C,D ⊂ℝ^2 be bounded measurable sets with A ⊇ B ⊇ C ⊇ D
and let f: ℝ^2 →ℝ be a Lebesgue-integrable function.
Then the following implications hold true:

 
  i) if mean_f(A \ B) ≤ mean_f(B \ C),
then: 

mean_f(A \ C) ≤ mean_f(B \ C)
and
mean_f(A \ B) ≤ mean_f(A \ C);
 
  ii) if mean_f(A \ B) ≤ mean_f(B \ C) ≤ mean_f(C \ D),
then: 

 mean_f(A \ C) ≤ mean_f(B \ D).




i) We first prove the first assertion in i).
If |A \ C| = 0, then this assertion is clearly true.
Let |A \ C| > 0.
Since
A\ C = A\ B ∪ B \ C and A\ B ∩ B \ C = ∅,
we conclude |A\ C| = |A\ B| + |B \ C| and in particular
|A\ B| >0 if |B \ C| = 0.
Assume that mean_f(A \ C) >  mean_f(B \ C). Then

    ∫_A \ C f   dx / |A\ C|
    =
    ∫_A \ B f   dx + ∫_ B \ C f   dx / |A\ B|
    + |B \ C| 
    >  mean_f(B \ C) =
    {[                      0       if  |B \ C| = 0,; ∫_B \ C f   dx/|B \ C|       if  |B \ C| > 0. ].

Both cases yield a contradiction to the assumption mean_f(A \ B) ≤ mean_f(B \ C).

Concerning the second assertion in i) we are done if |A \ B| = 0.
If |A \ B| > 0, then by the above considerations |A \ C| > 0
and assuming
mean_f(A \ B) >  mean_f(A \ C)
we obtain

    ∫_A \ C f   dx / |A\ C|
    =
    ∫_A \ B f   dx + ∫_ B \ C f   dx / |A\ B|
    + |B \ C| 
    <
    ∫_A \ B f   dx / |A\ B|.

This yields again a contradiction to the assumption.

ii) Applying the first/second implication in i) with respect to the first/second inequality in the assumption of ii) we obtain

    mean_f(A \ C) ≤ mean_f(B \ C)
       and   mean_f(B \ C) ≤ mean_f(B \ D),

which lead to
mean_f(A \ C) ≤ mean_f(B \ D)
and we are done.



Recall that τ_i^(k+1) = 1/2(m_i-1^(k) + m_i^(k)),
where
m_i^(k) =  mean_f(Σ_i^(k)\Σ_i+1^(k)).
Using the criterion derived in (<ref>) and Lemma <ref>, we can prove the following lemma.

 
Our T-ROF Algorithm <ref> produces sequences
( τ^(k) )_k and ( m^(k))_k
with the following properties:


  i)
0 ≤ m_0^(k)≤τ_1^(k)≤ m_1^(k)≤⋯≤ m_K-2^(k)≤τ_K-1^(k)≤ m_K-1^(k).

  ii) Set
τ_0^(k) := 0 and τ_K^(k) := 1.
If τ_i^(k)≥τ_i^(k-1) and
τ_i+1^(k)≥τ_i+1^(k-1),
then
m_i^(k)≥ m_i^(k-1), i=0,…,K-1
and this also holds true if ≥ is replaced everywhere by ≤.



 i) At step k, let τ_i = τ_i^(k).
For fixed 0 ≤τ_1 < τ_2 < ⋯ < τ_K-1 < 1, let Σ_i be a minimizer of E(·,τ_i).
With Lemma <ref>, we only need to prove m_0 ≤τ_1 and τ_K-1≤ m_K-1.
Assume that
Ω_0 := Ω\Σ_1 and
Ω_K-1 := Σ_K-1 have positive measure.

Since |Ω_0|, |Ω_K-1|>0, neither Ω nor ∅ is minimizer of E(·,τ_i),
we verify that

    Per(Σ_1;Ω) + μ∫_Σ_1 (τ_1 - f)   dx    ≤μ∫_Ω (τ_1 -f)   dx,
    
     Per(Σ_K-1;Ω) + μ∫_Σ_K-1 (τ_K-1 - f)   dx    ≤ 0.

From (<ref>), we have

    Per(Σ_1; Ω) + μ∫_Σ_1 (τ_1 - f)   dx ≤μ∫_Ω_0 (τ_1 - f)   dx + μ∫_Σ_1 (τ_1 - f)   dx.

Note that m_i :=  mean_f (Ω_i).
Thus, 0 ≤μ∫_Ω_0 (τ_1 - f)   dx which gives m_0 ≤τ_1.
From (<ref>), we have
0 ≥∫_Ω_K-1 (τ_K-1 - f)   dx,
which means τ_K-1≤ m_K-1.

ii) We only prove for sign ≥. The proof for sign ≤ follows the same lines.

If τ_1^(k)≥τ_1^(k-1), let τ_1 := τ_1^(k-1), τ_2 := τ_1^(k), using the conclusion from
Lemma <ref> i), we get

    0 ≤ mean_f (Ω\Σ_1^(k-1)) ≤τ_1^(k-1)≤ mean_f (Σ_1^(k-1)\Σ_1^(k))
    ≤τ_1^(k).

Hence, from the second implication of Lemma <ref> i), we have m_0^(k)≥ m_0^(k-1).

If τ_K-1^(k)≥τ_K-1^(k-1), let τ_1 := τ_K-1^(k-1), τ_2 := τ_K-1^(k),
using the conclusion from Lemma <ref> i), we get

    0 ≤τ_K-1^(k-1)≤ mean_f (Σ_K-1^(k-1)\Σ_K-1^(k))
    ≤τ_K-1^(k)≤ mean_f (Σ_K-1^(k)).

Hence, from the first implication of Lemma <ref> i), we have m_K-1^(k)≥ m_K-1^(k-1).

For i ∈{1,…,K-2 } and if τ_i^(k)≥τ_i^(k-1)
and τ_i+1^(k)≥τ_i+1^(k-1), from Lemma <ref> i), we have

    τ_i^(k-1)≤τ_i+1^(k-1)   and  τ_i^(k)≤τ_i+1^(k).

Hence, we can only have one of the following orderings:


  a) τ_i^(k-1)≤τ_i+1^(k-1)≤τ_i^(k)≤τ_i+1^(k);

  b) τ_i^(k-1)≤τ_i^(k)≤τ_i+1^(k-1)≤τ_i+1^(k).

In case a) we obtain by Lemma <ref> i) that

    τ_i^(k-1)≤ m_i ^(k-1)≤τ_i+1^(k-1)≤τ_i^(k)≤ m_i^(k).

In case b) we conclude by Lemma <ref> i) with the settings
τ_1 := τ_i^(k-1), τ_2 := τ_i^(k)
and
τ_1 := τ_i^(k), τ_2 := τ_i+1^(k-1)
and
τ_1 := τ_i+1^(k-1), τ_2 := τ_i+1^(k), respectively,

    τ_i^(k-1)≤    mean_f (Σ_i^(k-1)\Σ_i^(k))
    ≤τ_i^(k)≤ mean_f (Σ_i^(k)\Σ_i+1^(k-1))
    ≤τ_i+1^(k-1)≤
    ≤    mean_f (Σ_i+1^(k-1)\Σ_i+1^(k))
    ≤τ_i+1^(k).

By Lemma <ref> ii) this implies m_i^(k-1)≤ m_i^(k), which completes the proof.


It is straightforward to verify that the generated sequences satisfying (<ref>) will 
fulfil the above Lemma <ref>.



To prove the convergence of the sequence (τ^(k))_k, we
define a sign sequence
ζ^(k) = (ζ_i^(k))_i=1^K-1 as follows:
If τ_i^(k)≠τ_i^(k-1),

    ζ_i^(k):=
    
    +1,    if τ_i^(k) > τ_i^(k-1), 
    
    -1,    if τ_i^(k) < τ_i^(k-1),

and otherwise

    ζ_i^(k):=
    ζ_j^(k),    if i=1,
    ζ_i-1^(k),    if i≠ 1,

where
j := min{ l  |  τ_l^(k)≠τ_l^(k-1)}.
Additionally, since τ_0^(k) and τ_K^(k) do not change with k, we set 

    ζ_0^(k) := ζ_1^(k),   ζ_K^(k) := ζ_K-1^(k).

By s_k we denote the number of sign changes in ζ^(k).
Writing for simplicity ± instead of ± 1, we obtain for example that

    ζ^(k) = {ζ_i^(k)}_i=1^K-1 = 
     (+  + | - |+ | -|   +  + + |-  -|  +| -| +|  - - |+ |-   |++| —| +  +| - |+ |-)

has s_k = 17.


 
i) The number of sign changes s_k is monotone decreasing in k. 

ii) If ζ_1^(k+1)≠ζ_1^(k), then we have the strict decrease s_k+1 < s_k.



 i) Let s_k = N -1
and consider the sequence

    ζ_0^(k), ⋯, ζ_l_1^(k)|_v_1^(k) ⋯  ,
    |ζ_i_j^(k), ⋯, ζ_l_j^(k)|_v_j^(k) ⋯  ,
    |ζ_i_N^(k), ⋯, ζ_K^(k)_v_N^(k),

where v_j^(k), named a sign block, contains those successive components with the same sign, and ζ_0^(k) and ζ_K^(k)
(defined in (<ref>)) are used for the boundary elements of the whole sequence
which by definition  do not belong to ζ^(k) (see (<ref>) and (<ref>)). 
Note that #v_1^(k)≥ 2 and #v_N^(k)≥ 2.



1. For j≥ 2, if #v_j^(k)≥ 3, we consider ζ_i^*^(k+1) with i_j≤ i^*-1 < i^* < i^*+1 ≤ l_j, i.e.,
ζ_i^*-1^(k) = ζ_i^*^(k) = ζ_i^*+1^(k).
WLOG let ζ_i^*^(k) = -1.
Then we obtain by Lemma <ref> ii) that
m_i^*-1^(k)≤ m_i^*-1^(k-1) and m_i^*^(k)≤ m_i^*^(k-1).
Therefore

    τ_i^*^(k+1) = m_i^*-1^(k) + m_i^*^(k)/2≤m_i^*-1^(k-1) + m_i^*^(k-1)/2
    = τ_i^*^(k)

and consequently ζ_i^*^(k+1) = -1,
 or
ζ_i_j^(k+1) = ζ_i_j+1^(k+1) = ⋯ = ζ_i^*^(k+1) =1 (which is obtained by definition
(<ref>) when τ_i^*^(k+1) = τ_i^*^(k)).

Specifically, for j = 1, if #v_1^(k)≥ 3, we consider ζ_i^*^(k+1) with 0≤ i^*-1 < i^* < i^*+1 ≤ l_1.
Following the same lines above, we have, if ζ_i^*^(k) = -1,
then ζ_i^*^(k+1) = -1,
or
ζ_0^(k+1) = ζ_1^(k+1) = ⋯ = ζ_i^*^(k+1) =1. Moreover, in the latter case, we can also prove 
ζ_i^*^(k+1) = ζ_i^*+1^(k+1) = ⋯ = ζ_l_1^(k+1)  =1,
using the definition in (<ref>) for ζ_1^(k+1) when τ_1^(k+1) = τ_1^(k)
(which is just opposite to the way of defining ζ_i^(k+1) when  τ_i^(k+1) = τ_i^(k) for i>1).

The case ζ_i^*^(k) = 1 can be handled in the same way.


Therefore we obtain for a sequence v_j^(k) with #v_j^(k)≥ 3 the following properties from step k to step k+1:


  - no sign changes if the boundary signs are kept;

  - maximal one additional sign change if the left or right boundary sign is kept;

  - maximal two additional sign changes if both boundary signs are not kept.

Here is an example:

    -|++++|-   →      |++++|,  |+++-|, 
      |—+|, |–++|, 
        |-++-|, |–+-|, |—-|.



2. Assume that there is a sequence  in (<ref>) with no j such that #v_j^(k) = 1.
Consider
ζ_l_j^(k) and ζ_i_j+1^(k) which are different by definition.
WLOG, let
ζ_l_j^(k) = -1 and ζ_i_j+1^(k) = +1.
Then ζ_l_j -1^(k) = -1 and ζ_i_j+1 +1^(k) = +1. From Lemma <ref> ii), we have

    m_l_j-1^(k)≤ m_l_j-1^(k-1),    m_i_j+1^(k)≥ m_i_j+1^(k-1).

If
m_l_j^(k) < m_l_j^(k-1) (or m_l_j^(k)≥ m_l_j^(k-1)),
then
τ_l_j^(k+1) < τ_l_j^(k) (or τ_i_j+1^(k+1)≥τ_i_j+1^(k)).
This means that
ζ_l_j^(k+1)≠ζ_l_j^(k)
and
ζ_i_j+1^(k+1)≠ζ_i_j+1^(k) are not possible at the same time.

Let us call a subsequence containing consecutive sign blocks #v_j^(k)≥ 2 for all j within the whole sequence
a non-oscillating subsequence, short nosc-sequence.
Then we can use the  above arguments to show that the number of sign changes  for an inner nosc-sequence from step k to step k+1 
behaves as follows:


  - no additional sign changes if the boundary signs are kept;

  - maximal one additional sign change if the left or right boundary sign is kept;

  - maximal two additional sign changes  if both boundary signs are not kept.

Here is an example of the sign changes of a nosc-sequence containing two sign blocks for the first two boundary conditions:

    |+++—|    →      
    |+++ – -|,  | +++ +- -|, |+ ++ ++  -|,
     |++- – -|,
       
    |+++ - - +|,    |+++ +- +|,  | +++ +++| ,  |++-  – +|,
       
    |- ++ - - -|, | - -+ - - -|, |- - + ++ -|, | - -+ + - -|, 
       |- +++ - -|,|- - ++ + -| , |- + - - - -| , |- - - - - -|.


In particular, accompanying the proof in part 1, we obtain that, if the whole sequence is a nosc-sequence, the sign changes of the whole sequence 
will not increase.

3. Now, we consider the case # v_j^(k) = 1 for all j_1 ≤ j ≤ j_2, where
#v_j_1 - 1^(k) > 1 and #v_j_2 + 1^(k) > 1. 
We prove that from step k to step k+1 the signs of

    ζ_i_j_1-1^(·),  ζ_i_j_1^(·),  …,  ζ_i_j_2^(·),  ζ_i_j_2+1^(·)

can not change at the same time. 
We call the above sequence (<ref>) oscillatory subsequence, short osc-sequence.
Here is an example, where the osc-sequence is underbraced and the inner pattern overbraced:

    … + +|- + - +|- - …,   or  … + +|- + - + -  |+ + ….

WLOG assume that
ζ_i_j_1^(k)=-1 so that ζ_i_j^(k)=(-1)^j-j_1+1 for
j_1≤ j ≤ j_2, and ζ_i_j_1-1^(k) = ζ_i_j_1-2^(k)=+1 and
ζ_i_j_2+1^(k) = ζ_i_j_2+2^(k)=(-1)^j_2-j_1.
From Lemma <ref> ii), we know that

    m_i_j_1-2^(k)≥ m_i_j_1-2^(k-1),    and  
    m_i_j_2+1^(k)≥ m_i_j_2+1^(k-1)   if j_2-j_1  is even, 
    
    m_i_j_2+1^(k)≤ m_i_j_2+1^(k-1)   if j_2-j_1  is odd.

Assume that ζ_i_j_1-1^(·), ζ_i_j_1^(·), ⋯, ζ_i_j_2^(·),
ζ_i_j_2+1^(·) all change signs at the same time from step k to step k+1, 
e.g. for the left example in (<ref>), we have

    … + + | - + - + | - - …  →  …  - | + - + - |+  … .

Then we deduce 

    τ_ i_j_1-1^(k+1)   =  m_i_j_1-2^(k) + m_i_j_1-1^(k)/2 ≤τ_i_j_1-1^(k) = m_i_j_1-2^(k-1) + m_i_j_1-1^(k-1)/2

and since 
m_i_j_1-2^(k)≥ m_i_j_1-2^(k-1)
this implies

    m_j_1 - 1^(k)≤ m_j_1 - 1^(k) .

Taking into account that τ_ i_j_1^(k+1) = τ_ i_j_1^(k) would require the same sign ζ_ i_j_1^(k+1)
as in ζ_ i_j_1-1^(k+1) (see definition (<ref>)), we get further

    τ_ i_j_1^(k+1)   =  m_i_j_1-1^(k) + m_i_j_1^(k)/2  
    >
    τ_i_j_1^(k) = m_i_j_1-1^(k-1) + m_i_j_1^(k-1)/2

and since 
m_i_j_1-1^(k)≤ m_i_j_1-1^(k-1)
this implies

    m_i_j_1^(k) > m_i_j_1^(k-1).

Continuing in this way, we get finally
the contradiction

    m_i_j_2+1^(k) < m_i_j_2+1^(k-1),    if j_2-j_1  is even, 
    
    m_i_j_2+1^(k) > m_i_j_2+1^(k-1),    if j_2-j_1  is odd.


Then the above arguments show that the number of sign changes  for an osc-sequence from step k to step k+1 
behaves as follows:


  - the number of sign changes cannot increase;

  - the number of sign changes either not change or at least decreases by two if both boundary signs are kept;

  - the number of sign changes at least decreases by two if both boundary signs are not kept;

  - the number of sign changes at least decreases by one if one boundary sign is kept.


4. Now we can decompose the sequence ζ^(k) into subsequences with osc-pattern  and nosc-pattern 
(note that the case of a sequence ζ^(k) which is a nosc-sequence has been proved in part 2), 
where these patterns overlap in the boundary signs. 
This is illustrated by an example, where the osc-patterns are underbraced and the nosc-pattern are overbraced:



    + + 
    +  - +  -  +   + + - 
     -   + - +  -
     - + -   +   + — + 
      +  - + -  -.




    + + +  - +  -  ++ + --   + - +  - - + -   ++ - - - +  +  - +  - -



If nothing changes in the boundary of the osc-pattern,
then the number of sign changes cannot increase for ζ^(·) from step k to step k-1.
If both signs in the boundary of an osc-sequence changes, there are at least two sign changes less within the osc-sequence
which can be at most added in the nosc-sequences, but we will not obtain more sign changes for the whole sequence ζ^(·).
If one sign at the osc-boundary is kept and the other is changed we can use counting arguments
to show that the number of sign changes of the whole sequence cannot increase.
This completes the proof.

ii)
If ζ_1^(k+1)≠ζ_1^(k), we first prove all the components
in v_1^(k) must change signs at step k+1. 
If #v_1^(k) = 2, we get the proof directly since ζ_0^(k+1) follows the 
sign of ζ_1^(k+1) according to the definition in (<ref>);
if #v_1^(k)≥ 3, the proof can be found in the proof of part 1.

By parts 1–4 of the proof, we have block v_1^(k) will then be merged with 
all or part of v_2^(k) and have the same sign, where the whole sequence ζ^(·) will 
have one sign change less, since the rest components in v_2^(k),
if they exist, will be merged with all or part of v_3^(k), keeping this process until 
merging with the last block v_N^(k). Therefore we have s_k+1 < s_k.
This completes the proof.


Now we can prove the convergence of our T-ROF algorithm.



 
The sequence (τ^(k))_k∈ℕ
produced by the T-ROF Algorithm <ref> converges to a vector τ^*.
Then (Σ^*,τ^*) is a solution of the T-ROF model (<ref>).




We prove the assertion by induction on the number of sign changes s_k at some iteration step k.

i) Assume that s_k = 0. WLOG let
ζ_i^(k) = +1, i= 1, …, K-1,
i.e.,
τ_i^(k)≥τ_i^(k-1).
From Lemma <ref> ii), we obtain
m_i^(k)≥ m_i^(k-1) and consequently τ_i^(k+1)≥τ_i^(k), i= 1, …, K-1.
Therefore s_k+1 = 0 and
ζ_i^(k+1) = +1, i= 1, …, K-1.
This means that each sequence (τ_i^(k))_k is monotone increasing.
Since the sequences are moreover bounded in [0, 1], we conclude that (τ^(k))_k
converges. Note that s_k=0 when K=2.

ii) Assume that (τ^(k))_k converges if
s_k ≤ N-1 for some k ∈ℕ.

iii) We prove that (τ^(k))_k converges in case of s_k = N.
If there exists a k̂ such that
ζ_1^(k̂)≠ζ_1^(k̂ -1),
we get s_k̂≤ N-1 directly from Lemma <ref> ii); therefore
(τ^(k))_k converges by ii).
If ζ_1^(k+1) = ζ_1^(k) for all k>k̂, then
(τ_1^(k))_k>k̂ is monotone and bounded and converges consequently
to some threshold τ_1^*. By the definition of Σ_1^(k),
we have Σ_1^(k) converges, and therefore
m_0^(k) also converges. Since τ_1^(k) = (m_0^(k) + m_1^(k))/2,
we have m_1^(k) converges.
Now we prove τ_2^(k) also converges (here we assume there does not exist a k̂
such that ζ_2^(k+1) = ζ_2^(k) for all k>k̂, since otherwise the convergence of
τ_2^(k) is obtained immediately).
If not, τ_2^(k) must contain at least two convergent subsequences converge
to say τ_2^* and τ_2^'. WLOG, let τ_2^* < τ_2^'.
Then we have Σ_2^* ⊇Σ_2^' from Lemma <ref>.
In what follows, we show by contradiction that it is impossible.


  - If Σ_2^* ⊃Σ_2^', then
from Lemmas <ref> and <ref>, we have
mean_f (Σ_1 \Σ_2^*) <  mean_f (Σ_1\Σ_2^'), which implies m_1^(k) diverges.


  - If Σ_2^* = Σ_2^',  then there exists a k̂ such that
Σ_2^(k) = Σ_2^' for k > k̂. Considering thresholds
τ_3^(k), ⋯, τ_K-1^(k) in set Σ_2^*, the convergence of
them is obtained immediately from ii) by their sign changes which are ≤ N-1. Therefore
m_2^(k) converges, with the convergence of m_1^(k), so τ_2^(k) also converges. 

Analogously, we can have τ_3^(k) converges from the convergence of
τ_2^(k) and m_2^(k).
Repeating this procedure up to the final index K-1, we obtain the assertion.






§ NUMERICAL RESULTS


In this section we test our proposed T-ROF method on many kinds of images. More precisely, we use
the T-ROF Algorithm <ref> with a discrete ROF model (see e.g. <cit.>)
whose minimizer is computed numerically by an ADMM algorithm with its inner parameter fixed to 2.
Speedups by using more sophisticated methods will be considered in future work.
The stopping criteria in the T-ROF algorithm for u and τ are

    u^(i)-u^(i-1)_2 / u^(i)_2 ≤ϵ_u 	   and  τ^(k) - τ^(k-1)_2 ≤ϵ_τ,

where ϵ_u and ϵ_τ are fixed to 10^-4 and 10^-5,
respectively. The initialization of {τ_i^(0)}_i=1^K-1 was computed by
the fuzzy C-means method <cit.> with 100 iterations.

We compare our method with the recently proposed multiphase segmentation methods
<cit.>. Note that the methods <cit.> work with
the fixed fuzzy C-means codebook {m_i}_i=0^K-1 which are not updated (all the codebooks used in the
following examples are given in Appendix). Such update is however involved in <cit.>.
The default stopping criterion used in <cit.> is the maximum iteration steps,
which will be shown in each example;
the default stopping criterion used in <cit.> is the relative error with tolerance set to 10^-4;
and the default stopping criterion used in <cit.> is the same as the one used in <cit.>,
together with maximum 300 iteration steps.
We choose the regularization parameter,  λ/μ (note that λ is used in the PCMS model and μ is
used in the SaT method <cit.> and our T-ROF model), in front of the fidelity term for all the methods
by judging the segmentation accuracy (SA) defined as

    SA := #correctly classified pixels/#all pixels

unless otherwise stated. 
We show the results for two two-phase and five multiphase images. 
Moreover, an extra example – one three-phase segmentation problem based on retina manual segmentation –
is provided to further demonstrate the superior performance of our proposed method.   
All computations were run on
a MacBook with 2.4 GHz processor and 4GB RAM.



 §.§ Two-phase Image Segmentation


Example 1. Cartoon image with some missing pixel values.
Fig. <ref> (a) is the clean two-phase image, with constant value in each phase.
Fig. <ref> (b) is the test corrupted image generated by
removing some pixel values (in this example, 80% pixels are removed) randomly from Fig. <ref> (a).
Fig. <ref> (c)–(g) are the results of the methods <cit.>, respectively.
Fig. <ref> (h)–(j) give the results of our T-ROF method respectively at iterations 1, 2 and 6 (final result)
with μ=1, which clearly shows the
effectiveness of the updating strategy of our T-ROF method on τ given in (<ref>).
For the T-ROF method, the iteration steps to find u and τ are 418 and 6 (simply represented as
418 (6) in Table <ref>), respectively.
From those results, we see that only
methods <cit.> and our T-ROF method give good results.
The quantitative comparison in terms of μ, iteration steps, computation time, and SA for each method is given
in Table <ref>. We can see that the T-ROF method gives the highest SA, which
shows the effectiveness and necessity of updating the threshold τ, compared with the slightly poor results of
the SaT method <cit.> which just uses the fixed thresholds selected by K-means.






Example 2. Two-phase image with close intensities.
Fig. <ref> (a) is an image generated by adding Gaussian noise with mean 0 and
variance 10^-8 onto a constant image with constant value 0.5. Fig. <ref> (b) is a
two-phase mask separating the whole domain into two parts (the black and the white color parts).
The testing noisy image Fig. <ref> (c) is generated from Fig. <ref> (a)
by keeping the pixel values belonging to the white part in the mask and reducing the pixel values
belonging to the black part by a factor of 2×10^-4.
This way of generating test images has the following two main features:
i) the noise pattern (e.g. Gaussian) presented in the noisy image is locally changed slightly, which 
will make the test more challenging so as to better evaluate the performance of each method and to 
do comparison; ii) the way of lower or increase the intensities of some specified areas can help to generate 
test images easily which contain phases with close intensities, where these generated images  
are good candidates to test the performance of different methods in classifying the close phases.  
Fig. <ref> (d)–(h) are the results of methods
<cit.>, respectively.
Fig. <ref> (i)–(j) are the results of our T-ROF method with μ=8 at iterations 1 and 6 (final result),
respectively, which again clearly show the effectiveness of the
updating strategy on τ given in (<ref>).
Obviously, except method <cit.>, all other methods can get good results.
The quantitative results given in Table <ref> show that our method is
the fastest and gives the highest SA, which again validates
 the necessity of updating the threshold τ given in (<ref>)
 against the way of obtaining thresholds by K-means used in the SaT method <cit.>.









 §.§ Multiphase Image Segmentation


Example 3. Five-phase noisy image segmentation.
Fig. <ref> (a) and (b) are the clean image and Gaussian noisy image with mean 0 and
variance 10^-2. Fig. <ref> (c)–(h) are the results of methods
<cit.> and our T-ROF method (with μ=8),
respectively. From there results, we see that all the results are very good except the result of method <cit.>.
From Table <ref>, we can see that the T-ROF method is the fastest.






Example 4. Four-phase gray and white matter segmentation for a brain MRI image.
In this example, we test the four-phase brain MRI image used in <cit.>,
see Fig. <ref> (a). The gray and white matter segmentation for this kind of image is
very important in medical imaging. Fig. <ref> (b)–(g) are the results of methods
<cit.> and our T-ROF method (with μ=40), respectively.
We can see that all the methods work well for this kind of image. In particular, the T-ROF method with 11 τ-value
updates is faster than other methods, e.g., three times faster than the algorithm of Pock et al. <cit.>
with assigned parameters. Note that the SaT method <cit.> is also very fast, due to the fact that it is akin to the
T-ROF method.









Example 5. Stripe image.
In this example, we test methods on segmenting the noisy stripe image in Fig. <ref> (b),
which is generated by imposing Gaussion noise with mean 0 and variance 10^-3 on
the clean image Fig. <ref> (a) with 30 stripes.
The results are generated by segmenting Fig. <ref> (b) into five, ten, and fifteen phases, respectively.
The rows two to seven of Fig. <ref> are the results of methods
<cit.> and our T-ROF method (with μ=8),
respectively. The quantitative comparison is shown in Table <ref>,
from which we can see that methods <cit.> and our T-ROF method give much better
results in terms of SA; note, importantly, that our method is always the fastest compared with methods <cit.>.
Moreover, Table <ref> shows clearly the great advantage of the T-ROF method and SaT method <cit.>
in computation time: their computation time is independent to
the required number of phases K, whereas this is not the case for other methods
(as the number of phases goes larger, their computation time significantly increases inevitably).
From Table <ref>, we also see that the T-ROF method gives much better results than the SaT method,
which again shows the excellent performance and necessity of updating the threshold τ
benefited from the rule proposed in (<ref>). 








Example 6. Three-phase image containing phases with close intensities.
In this example, we test a three-phase image, where two phases of it have very close intensities.
The test image in Fig. <ref> (c) is generated using the same way as that in Example 2 with
Gaussian noise of mean 0 and variance 10^-2, and the factors used in the black and
white color parts in the mask are 0.1 and 0.6, respectively.
For the results of methods <cit.>, in order to reveal the performance of these methods
clearly, we give two representative results for each method compared using different regularization
parameters, see Fig. <ref> (d)–(m).
From Fig. <ref>, we can see that our T-ROF method gives the best result, see 
Fig. <ref> (o) (with μ=8). In detail, we notice that method <cit.> gives very poor results, see Fig. <ref> (d)–(e).
For methods <cit.>, no matter how their parameters are tuned, they all cannot
achieve good results as the result of our T-ROF method given in Fig. <ref> (o). 
More specifically, methods <cit.> either give results
separating different phases unclearly, or give results not removing the noise successfully.
The main reason is that, the regularization parameters in methods <cit.> are constant, see 
the parameter λ used in the PCMS model (<ref>); however, constant parameter used to penalize all the phases equally is
obviously not appropriate for this case.
The results of methods <cit.>, Fig. <ref> (h) and (l), are better than the results of methods
<cit.>. Again, after comparing
Fig. <ref> (h), (l) and Fig. <ref> (o) from visual-validation and the quantitative results given in
Table <ref>, we see that our T-ROF method gives the best result
in terms of segmentation quality and computation time.
	
	

	
Example 7. Four-phase image containing phases with close intensities.	
In order to show the powerful of our T-ROF method in handing images containing phases with close intensities, 
in this example, we test these methods on a four-phase image where each two phases with close intensities.
Fig. <ref> (a) and (b) are respectively the clean image and the noisy image
generated by adding Gaussian noise with mean 0 and variance 3× 10^-2.
Similar to Example 6, this example also provides two results for each method compared using different
representative parameters, see Fig. <ref> (c)–(l).
From Fig. <ref> and Table <ref> (which gives quantitative results), we see that
methods <cit.> all give poor results compared with
our result in terms of segmentation quality and computation time.
In particular, our method gives much better results than that of its akin SaT method <cit.>,
see Fig. <ref> (k) and (n), which further verifies the excellent performance of updating the threshold τ
using the rule proposed in (<ref>).


	










 §.§ Discussion About the Threshold τ 


We report, in Fig. <ref>, the convergence history of τ of our T-ROF algorithm <ref> corresponding to iteration steps
for Examples 1–7. From Fig. <ref>, we can see that τ always
not only converges, but converges very quickly – just a few steps are enough (generally within ten iterations for the examples
shown in this paper).

It is worth mentioning that, for all the Examples 1–7, the rule (<ref>) and the criterion derived in (<ref>)
are always satisfied at each iteration. This means that these rules are indeed very easy to fulfil, if one uses meaningful number of phases K
(here K is our prior knowledge for each segmentation example) and 
initializations (here we use K-means to obtain initializations). The chances that these criterions are likely to be invoked based on our tests are the cases, 
for example, meaningless number of phases K and/or initializations. For example, given a test image which has two phases, using K = 3 would be highly likely to 
obtain an empty phase at some iteration and one would finally obtain a segmentation result with two phases. An empty phase that is also likely to be obtained
is when two thresholds are too close in the initialization. This performance is actually a kind of
automatic error correction, which is an innate ability of our T-ROF algorithm. In contrast, the methods solving the PCMS model 
(like the ones we compared in this paper) cannot have. 

  





 §.§ Real-World Example


In the previous examples, we mainly investigated the performance of the aforementioned segmentation methods
using synthetic images and quantified their accuracy using SA (segmentation accuracy) defined in (<ref>). 
Before closing this section, to complement the test, we hereby test those methods using a manual segmentation of a real-world image, 
the retina image, which is from the DRIVE data-set[http://www.isi.uu.nl/Research/Databases/DRIVE/]. 
Fig. <ref> (a) and (b) are the clean manual segmentation image and the noisy image
generated by adding Gaussian noise with mean 0 and variance 0.1. Note that in Fig. <ref> (a), 
we changed the original binary manual segmentation image to three phases by lowering the intensity of those vessels on the right hand side
from 1 to 0.3; the intensities of the background and the vessels on the left hand side are respectively 0 and 1. 
Obviously, segmenting the noisy three-phase image in Fig. <ref> (b) is extremely challenging 
due to those thin blood vessels which have a big chance of being smoothed out. 

To quantify the segmentation accuracy, in addition to the measure SA in (<ref>) which
quantifies the whole segmentation accuracy, we here also use the DICE score which is able
to quantify the accuracy for individually segmented phases, i.e., 

    DICE(Ω_i, Ω_i^') := 2 |Ω_i ∩Ω_i^'|/|Ω_i| + |Ω_i^'|,

where Ω_i and Ω_i^' are respectively the segmentation result and the ground truth 
corresponding to phase i, and |Ω_i| is the cardinality of set Ω_i in discrete setting.
Note that the DICE score not only involves the correctly segmented areas like the measure SA in (<ref>), but the incorrectly areas. 

Fig. <ref> shows the segmentation results of the methods compared, and Table <ref> gives the quantitative results
including the SA and the DICE score which quantifies the three individually segmented phases, i.e.,
the areas of the background, vessels on the right hand side with lower intensity and vessels on the left hand side.
We clearly see that all the methods successfully segmented the vessels on the left hand side, but
 failed to segment the vessels on the right hand side (please refer to
Section <ref> for the theoretical discussion about the drawback of the PCMS model for case K > 2), except for method <cit.> and 
our proposed T-ROF method (also with much faster speed). Particularly, to segment the vessels on the right hand side which have lower intensity, 
the proposed T-ROF method achieved higher DICE score compared to the SaT method <cit.> – 0.7749 versus 0.5673 – 
which further verifies the excellent performance of updating the threshold τ
using the rule proposed in (<ref>).

	








§ SUMMARY AND CONCLUSIONS


In this paper, we unveiled a linkage between the PCMS model and the ROF model,
which is important to build the connection between image segmentation and image
restoration problems. The built linkage theoretically validates our proposed novel segmentation
methodology – pursuing image segmentation by applying image restoration
plus thresholding.
This new segmentation methodology can circumvent the innate non-convex property of
the PCMS model, and thus improves the segmentation performance in both efficiency and effectiveness.
In particular, as a direct by-product of the built linkage, we proposed a segmentation method named
T-ROF method. The convergence of this method has also been proved.
Elaborate experimental results were presented which all support the excellent performance of
the proposed T-ROF method in terms of segmentation quality and computation time.
For the future work, considering using other thresholding rules, e.g. using the median instead of the mean in
the T-ROF method, may be worthwhile. Moreover, similar to the linkage built in this paper between the
PCMC and ROF models, investigating the relationship
between other models, such as those variants of the PCMS and ROF models
respectively in image segmentation and image restoration, will also be of interest. 






§ APPENDIX


The codebooks (m := {m_i}_i=0^K-1) computed by <cit.> and used for methods <cit.> and our T-ROF method
(as initializations) in Examples 1-7 are listed in Table <ref>.






10

ACC05
F. Alter, V. Caselles, and A. Chambolle, A characterization of
 convex calibrable sets in R^N, Mathematische Annalen, 332 (2005),
 pp. 329–366.

AFP00
L. Ambrosio, N. Fusco, and D. Pallara, Functions of Bounded
 Variation and Free Discontinuity Problems, Oxford University Press, Oxford,
 2000.

AT90
L. Ambrosio and V. Tortorelli, Approximation of functions depending
 on jumps by elliptic functionals via t-convergence, Communications in Pure
 and Applied Mathematics, 43 (1990), pp. 999–1036.

ABM06
H. Attouch, G. Buttazzo, and G. Michaille, Variational analysis in
 Sobolev and BV Spaces, SIAM, Philadelphia, 2006.

BCCJKM11
L. Bar, T. Chan, G. Chung, M. Jung, N. Kiryati, R. Mohieddine, N. Sochen,
 and L. Vese, Mumford and shah model and its applications to image
 segmentation and image restoration, in Handbook of Mathematical Imaging,
 Springer, 2011, pp. 1095–1157.
 
 
 BCPSS17
B. Bauer, X. Cai, S. Peth, K. Schladitz, and G. Steidl. 
Variational-based segmentation of biopores in tomographic images,
Computers & Geosciences, 98 (2017), pp. 1–8.

BPV91
G. Bellettini, M. Paolini, and C. Verdi, Convex approximations of
 functionals with curvature, 
 Natur. Rend. Lincei (9) Mat. Appl., 2 (1991), pp. 297–306.

BEF84
J. C. Bezdek, R. Ehrlich, and W. Full, FCM: The fuzzy c-means
 clustering algorithm, Computers & Geosciences, 10 (1984), pp. 191–203.

BZ87
A. Blake and A. Zisserman, Visual Reconstruction, MIT Press,
 Cambridge, MA, 1987.

BPCPE10
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, Distributed
 optimization and statistical learning via the alternating direction method of
 multipliers, Foundations and Trends in Machine Learning, 3 (2010),
 pp. 1–122.

BEVTO07
X. Bresson, S. Esedoglu, P. Vandergheynst, J. Thiran, and S. Osher, Fast global minimization of the active contour/snake model, Journal of
 Mathematical Imaging and Vision, 28 (2007), pp. 151–167.

BCB12
E. Brown, T. Chan, and X. Bresson, Completely convex formulation of
 the Chan-Vese image segmentation model, International Journal of Computer
 Vision, 98 (2012), pp. 103–121.

C15
X. Cai, Variational image segmentation model coupled with image
 restoration achievements, Pattern Recognition, 48 (2015), pp. 2029–2042.

CCNZ15
X. Cai, R. Chan, M. Nikolova, and T. Zeng, A three-stage approach
 for segmenting degraded color images: smoothing, lifting and thresholding
 (slat), Journal of Scientific Computing, 72 (2017), pp. 1313–1332.

CCZ13
X. Cai, R. Chan, and T. Zeng, A two-stage image segmentation method
 using a convex variant of the mumford-shah model and thresholding, SIAM
 Journal on Imaging Sciences, 6 (2013), pp. 368–390.

CFNSS15
X. Cai, J. Fitschen, M. Nikolova, G. Steidl, and M. Storath, 
Disparity and optical flow partitioning using extended Potts priors, 
Information and Inference: A Journal of the IMA, 4 (2015), pp. 43–62.

CS13
X. Cai and G. Steidl, Multiclass segmentation by iterated rof
 thresholding, in EMMCVPR, LNCS, Springer, 2013, pp. 237–250.

Ch04
A. Chambolle, An algorithm for total variation minimization and
 applications, Journal of Mathematical Imaging and Vision, 2004,
 pp. 89–97.
 
 Ch05
A. Chambolle, Total variation minimization and a class of binary MRF models, 
in EMMCVPR, LNCS, Springer, 2005, pp. 136–152.

 CNCP10
A. Chambolle, V. Caselles, M. Novaga, D. Cremers, and T. Pock,
An introduction to total variation for image analysis,
in Theoretical Foundations and Numerical Methods for Sparse Recovery, 
Radon Ser. Comput. Appl. Math., 9 (2010), pp. 263–340.

CP11
A. Chambolle and T. Pock, A first-order primal-dual algorithm for
 convex problems with applications to imaging, Journal of Mathematical
 Imaging and Vision, 40 (2011), pp. 120–145.

CYZ13
R. Chan, H. Yang, and T. Zeng, A two-stage image segmentation method
 for blurry images with poisson or multiplicative gamma noise, SIAM Journal
 on Imaging Sciences, 7 (2014), pp. 98–127.

CEN06
T. F. Chan, S. Esedoglu, and M. Nikolova, Algorithms for finding
 global minimizers of image segmentation and denoising models, SIAM Journal
 on Applied Mathematics, 66 (2006), pp. 1632–1648.

CEPY06
T. F. Chan, S. Esedoglu, F. Park, and A. Yip, Total variation image restoration:
Overview and recent developments, in Handbook of Mathematical Models in Computer Vision, 
edited by N. Paragios, Y. Chen, and O. Faugeras, Springer-Verlag, New York, 2006, 
pp. 17–31.

CV01-2
T. F. Chan and L. A. Vese, Active contours without edges, IEEE
 Transactions on Image Processing, 10 (2001), pp. 266–277.

DCS10
B. Dong, A. Chien, and Z. Shen, Frame based segmentation for medical
 images, Commun. Math. Sci., 32 (2010), pp. 1724–1739.

GG84
S. Geman and D. Geman, Stochastic relaxation, Gibbs distributions,
 and the Bayesian restoration of images, IEEE Transactions on Pattern
 Analysis and Machine Intelligence, 6 (1984), pp. 721–741.

GO09
T. Goldstein and S. Osher, The split Bregman method for
 l1-regularized problems, SIAM Journal on Imaging Sciences, 2 (2009),
 pp. 323–343.

GFK07
J. Gorski, F. Pfeiffer, and K. Klamroth, Biconvex sets and
 optimization with biconvex functions - a survey and extensions, Mathematical
 Methods of Optimization Research, 66 (2007), pp. 373–407.

HHMSS12
Y. He, M. Y. Hussaini, J. Ma, B. Shafei, and G. Steidl, A new fuzzy
 c-means method with total variation regularization for image segmentation of
 images with noisy and incomplete data, Pattern Recognition, 45 (2012),
 pp. 3463–3471.

LS12
J. Lellmann and C. Schnörr, Continuous multiclass labeling
 approaches and algorithms, SIAM Journal on Imaging Sciences, 44 (2011),
 pp. 1049–1096.

LNZS10
F. Li, M. Ng, T. Zeng, and C. Shen, A multiphase image segmentation
 method based on fuzzy region competition, SIAM Journal on Imaging Sciences,
 3 (2010), pp. 277–299.

MS89
D. Mumford and J. Shah, Optimal approximation by piecewise smooth
 functions and associated variational problems, Communications on Pure and
 Applied Mathematics, XLII (1989), pp. 577–685.

NEC06
M. Nikolova, S. Esedoglu, and T. F. Chan, Algorithms for finding
 global minimizers of image segmentation and denoising models, SIAM Journal
 on Applied Mathematics, 66 (2006), pp. 1632–1648.

PCCB09
T. Pock, A. Chambolle, D. Cremers, and H. Bischof, A convex
 relaxation approach for computing minimal partitions, IEEE Conference on
 Computer Vision and Pattern Recognition, 2009, pp. 810–817.

PCBC09a
T. Pock, D. Cremers, H. Bischof, and A. Chambolle, An algorithm for
 minimizing the piecewise smooth mumford-shah functional, in ICCV, 2009.

ROF92
L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation
 based noise removal algorithms, Physica D, 60 (1992), pp. 259–268.

VC02
L. Vese and T. Chan, A multiphase level set framework for image
 segmentation using the mumford and shah model, International Journal of
 Computer Vision, 50 (2002), pp. 271–293.

YBTB10
J. Yuan, E. Bae, X.-C. Tai, and Y. Boykov, A continuous max-flow
 approach to potts model, in European Conference on Computer Vision, 2010,
 pp. 379–392.

ZGFN08
C. Zach, D. Gallup, J.-M.Frahm, and M. Niethammer, Fast global
 labeling for real-time stereo using multiple plane sweeps, Vision, Modeling,
 and Visualization Workshop, 2008.

ZMSM08
Y. Zhang, B. Matuszewski, L. Shark, and C. Moore, Medical image
 segmentation using new hybrid level-set method, in 2008 Fifth International
 Conference BioMedical Visualization: Information Visualization in Medical and
 Biomedical Informatics, 2008, pp. 71–76.







