# Fast Mumford-Shah Two-Phase Image Segmentation Using Proximal Splitting Scheme

Authors: Nilima Shah; Dhanesh Patel; Pasi Fränti
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

The Mumford-Shah model is extensively used in image segmentation. Its energy functional causes the content of the segments to remain homogeneous and the segment boundaries to become short. However, the problem is that optimization of the functional can be very slow. To attack this problem, we propose a reduced two-phase Mumford-Shah model to segment images having one prominent object. First, initial segmentation is obtained by the k-means clustering technique, further minimizing the Mumford-Shah functional by the Douglas-Rachford algorithm. Evaluation of segmentations with various error metrics shows that 70 percent of the segmentations keep the error values below 50%. Compared to the level set method to solve the Chan-Vese model, our algorithm is significantly faster. At the same time, it gives almost the same or better segmentation results. When compared to the recent k-means variant, it also gives much better segmentation with convex boundaries. The proposed algorithm balances well between time and quality of the segmentation. A crucial step in the design of machine vision systems is the extraction of discriminant features from the images, which is based on low-level segmentation which can be obtained by our approach.


---

## Body

## 1. Introduction

Image is an entity having different objects. Image segmentation is one of the techniques to recognize the regions that show different objects in the image.Let the functiongx,y represent the intensity of light at any point x,y of domain Ω∈R2 when an image is captured. So, image is defined as g=gx,y∣x,y∈Ω⊂R2. The RGB color representation of an image is mathematically represented by gx,y∈R3. The aim is to find regions Ω1,Ω2,Ω3,⋯⋯Ωn for n∈N, of a domainΩ⊂R2, corresponding to different objects or their parts or shadows in the image. These regions Ωi,i=1,2,⋯⋯n, of the image will have boundaries. In summary, the segmentation problem has the following definition: suppose Ω be the domain of g. Compute and find Ω1,Ω2,⋯Ωn such that
(1)Ω=⋃i=1nΩi,and g varies smoothly over the regions Ωi,i=1,2,⋯⋯nas well as g varies discontinuously over the boundaries between Ωi,i=1,2,⋯n. The same problem can be represented in terms of approximation theory: find an approximation f of g in such a way that each fi represents the entire region Ωi via a piecewise smooth function.There are many approaches to image segmentation [1] including thresholding, clustering, or classification methods [2–4], region-based methods [5], and edge-based active contour methods [6–8]. Segmentation techniques based on deep convolutional neural networks have been developed for various medical images MRI, CT, and X-rays, which show promising results in comparison to conventional segmentation methods [9–12]. The segmentation problem has also been considered as an energy minimization problem. But till date, obtaining a generalized solution that works without prepriori facts about the object elements, its characteristics like shape, color, texture, appearance of shadows, and overlapping of objects are an open problem.Segmentation is essentially a clustering problem with additional connectivity constraint. k-means is widely used, its properties are well known, and it has also been applied to image segmentation with modest modifications. It is fast and simple to implement [13], which makes it suitable for large datasets like images. However, applying k-means merely based on the pixel values usually results in very fragmented segments. Regularized k-means (reg-KM) [14] attacks this by combining the intensity and the spatial information using two terms. The first term calculates the ratio of the segment perimeter and its area while the second term measures the intracluster dissimilarity as in the standard k-means.In 1989, Mumford and Shah introduced the Mumford-Shah model [15]. Since then, it has applications in the areas such as image denoising, image restoration, and image segmentation [16]. The Mumford-Shah model can be implemented using the Douglas-Rachford algorithm, ADI implicit method, and other numerical optimizations [17–22]. The improvement of local minima was achieved by graph-based methods, but these methods do not allow for open boundaries [23].A convex relaxation approach was presented in [24], using functional lifting which guarantees a globally optimal solution, is independent of initialization, and takes care of open boundaries. There are suitable algorithms to solve the simplified Mumford-Shah model which include convex boundaries for any kind of image segmentation [18, 22, 25–27]. The problems of the previous approaches are that those based on Mumford-Shah are rather slow, while many region and contour-based models provide poor segmentation result if the parameters are not carefully tuned.In this paper, we introduce a two-phase image segmentation algorithm. Instead of using region or block-based heuristics to solve the segmentation problem, we define it as a global optimization problem. We input the number of segments, select the initial centroids randomly, and apply the k-means clustering technique to obtain the initial segmented image, which is further improved by minimizing the Mumford-Shah model. k-means is known to be sensitive to the initialization. A typical solution is to try the algorithm several times with different initialization or use some more robust initialization technique, see [13]. Our main contribution is to show how the proximal splitting technique can be used to solve the two-phase Mumford-Shah model.We compare our segmentation with ground truth data for 100 color images from the Weizmann image dataset [28]. Statistically evaluating using error metrics defined in [29], we find 70% of our segmentation have an error value of 0.5 or less. The method works especially well with segments that are difficult to determine, there are multiple objects in the image, and objects have stripes, or complex objects posing so that they have artificial boundaries. For these troublesome images, we find our algorithm gives slightly better segmentation when compared to the classical Chan-Vese optimization algorithm [8], and much better than the k-means variant (reg-KM) [14]. Our implementation of the Douglas-Rachford (DR) algorithm is also a lot faster than the Chan-Vese algorithm. The Chan-Vese model jointly uses the reduced Mumford-Shah model and level set method which also uses active contours for image segmentation. This is time consuming because it is solved using reinitialization of a level set function [30, 31].The modern machine vision methods have higher-level goals like extracting an anatomical object of interest for diagnosis of diseases, scene classification, detection of human activity in visual surveillance, assigning a name to a human face in an image, and classifying handwritten characters [32]. A crucial step in the design of such machine vision systems is the extraction of discriminant features from the images [33, 34]. The application of the proposed approach is to obtain a low-level segmentation based on color features and spatial connectivity, which can be further fed to the above-mentioned higher-level semantic segmentation.The rest of the paper is organized as follows: Section2 introduces the Mumford-Shah image segmentation model and shows its discrete formulation. Section 3 explains the proximal splitting technique used for optimizing the Mumford-Shah functional and also contains pseudocode. In Section 4, we explain the evaluation criteria used and discuss the results. Conclusions are drawn in Section 5.


## 2. The Mumford-Shah Image Segmentation Model

The Mumford-Shah model is based on the functionalEf;Γ which is to be minimized and it depends on two parameters μ and λ. The functional Ef;Γ is given by
(2)Ef;Γ=μ2∫∫Ωf−g2dxdy+∫∫Ω−Γ∇f2dxdy+λLΓ,where g is an image, function f approximates g which is differentiable over the regions Ωi, and
(3)Ω=⋃i=1nΩi∪Γ,where Γ is made up of smooth arcs and it creates boundaries for Ωi. The overall length of all the smooth arcs of Γ is expressed by LΓ. The parameters μ and λ control the approximation quality and segment coarseness, respectively.The smaller value ofE gives good segmentation f,Γ and better representation of the imageg. The first term in equation (2) measures how well f approximates the image g. The second term measures how much f varies on each Ωi. The third term in equation (2) measures the length of the boundaries and ensures that the boundaries Γ are as small as possible. Mumford and Shah mentioned that [15], dropping any of the terms in equation (2), we have infE=0. If the first term is not in equation (2), we have f=0,Γ=∅. If the second term is not there in equation (2), we take f=g,Γ=∅. If the third term is not there in equation (2), we take Γ as a fine grid of N horizontal and vertical lines and hence Ωi=N2 number of small squares. Also, f is average of g on each Ωi.The presence of all the three terms leads to nontrivial optimal solutionEf,Γ. The optimal value f,Γ can be visualised as a cartoon of the actual image g. Both Mumford and Shah were uncertain for the well-posedness of minimization of E. However, they conjectured that for the continuous function g, E has a minimum in the set of all pairs f,Γ where f is differentiable on each Ωi and Γ a finite set of singular points joined by a finite set of C1 arcs.The energy functionalEf,Γ in equation (2) whose exact minimization is a nontrivial task. When f is a piecewise constant function on each Ωi, obtained as the average value of g on each Ωi, the second term in equation (2), ∫∫Ω−Γ∇f2dxdy=0. Therefore, equation (2) reduces to the simplified Mumford-Shah model [5, 21, 35] given as
(4)Ef,Γ=μ2∫∫Ωf−g2dxdy+λLΓ.This is a special case of equation (2). This is also called the minimal partition problem [6]. If λ is large, the emphasis is on the length of the boundaries Γ and the minimization of E leads to fewer boundaries and coarser segmentation. If λ is small, the emphasis is on the first term of equation (4) which leads to an approximation f that follows g much more closely and allows for more boundaries. Thus, a small λ leads to a finer segmentation. If λ=0, the minimum of E forces f=g leading to finest possible segmentation, which is a trivial segmentation, where every region consists of a single point (pixel).


### 2.1. Image Segmentation Model

Equation (4) is the reduced two-phase segmentation model. In the two-phase model, the image is divided into two segments, the object and the background. Chan-Vese used the level set method to implement the two-phase segmentation model [6, 8]. In our approach, we solve the reduced two-phase model by rewriting the problem using the Douglas-Rachford proximal splitting algorithm. We use a convexification technique similar to [18, 22].In this work, we have considered RGB images. The goal is to divide the given image into two regions, the background and the foreground object. We define the foreground object asΩ0, background as Ω1, and the entire image as Ω. So, image g=Ω0∪Ω1∪Γ, here Γ represent the arcs of the segment Ω0. The image g contains n×m pixels.Letc0 and c1 denote the color values of randomly selected pixels as foreground pixel and background pixel, respectively. From simplified equation (4), we have
(5)Ec0,c1,Γ=μ2∫∫Ω0ω0dxdy+μ2∫∫Ω1ω1dxdy+λLΓ,where ω0=gx,y−c02 and ω1=gx,y−c12, x=1⋯n; y=1⋯m. Let Ω0=x,y∣ω0x,y<ω1x,y and Ω1=complementof Ω0.Now, our aim is to find the regionsΩ0 and Ω1 that minimize functional Ec0;c1;Γ which is given in equation (5). Defining ω=ω0−ω1, we represent the domain Ω with binary indicator function f so that
(6)fx,y=χΩx,y1;x,y∈Ω0,i.e.,ω<0,0,otherwise.This functionf is the expected approximation of image g. Using f, the minimization of functional Ec0;c1;Γ from equation (5) can be seen as follows:
(7)minfx,y∈0,1f,ω+λfTV.The first term in equation (7) is the inner product of f and ω and is defined as f,ω=∫fω. The total variation pseudonorm is fTV, which equals to LΓ for binary indicator f=χΩ by coarea formula established by Fleming and Rishel [20].


### 2.2. Discrete Formulation of Segmentation Model

To discretize variational problem (7) on a grid of n×m pixels, we assume that N=n×m. Let vi∈R2 be a vector field for every i and suppose
(8)v→=vi→i=1N∈RN×2,then norm of v is given by
(9)v→1=∑i=1Nvi→.The total variation pseudonorm, forf=fiN/i=1∈RN, is defined as
(10)fTV=∇f1=∑i=1N∇fi,using equation (9). Thus, we obtain
(11)minfi∈0,1f,ω+λ∇f1.We assume that the pixels are indexed on a two-dimensional grid whereδ1=1,0 and δ2=0,1, the finite difference gradient operator is
(12)∇fii=1N=fi+δ1−fi,fi+δ2−fii=1N∈RN×2.We use periodic boundary conditions for simplicity. The inner productRN is
(13)f,ω=∑i=1Nfiωi.Equation (11) is a nonsmooth and nonconvex and usually direct minimization that leads to poor local minima. By replacing the binary constraint fi∈0,1 with a box constraint fi∈0,1, we obtain the convex formulation of the energy functional of equation (11). This defines the following finite dimensional convex problem:
(14)minfi∈0,1Nf,ω+λ∇f1.One can prove that this relaxation is exact, meaning that the minimizerf, when it is unique, is binary, f∈0,1. It means that Ω such that f=χΩ actually solves the original segmentation problem.


## 2.1. Image Segmentation Model

Equation (4) is the reduced two-phase segmentation model. In the two-phase model, the image is divided into two segments, the object and the background. Chan-Vese used the level set method to implement the two-phase segmentation model [6, 8]. In our approach, we solve the reduced two-phase model by rewriting the problem using the Douglas-Rachford proximal splitting algorithm. We use a convexification technique similar to [18, 22].In this work, we have considered RGB images. The goal is to divide the given image into two regions, the background and the foreground object. We define the foreground object asΩ0, background as Ω1, and the entire image as Ω. So, image g=Ω0∪Ω1∪Γ, here Γ represent the arcs of the segment Ω0. The image g contains n×m pixels.Letc0 and c1 denote the color values of randomly selected pixels as foreground pixel and background pixel, respectively. From simplified equation (4), we have
(5)Ec0,c1,Γ=μ2∫∫Ω0ω0dxdy+μ2∫∫Ω1ω1dxdy+λLΓ,where ω0=gx,y−c02 and ω1=gx,y−c12, x=1⋯n; y=1⋯m. Let Ω0=x,y∣ω0x,y<ω1x,y and Ω1=complementof Ω0.Now, our aim is to find the regionsΩ0 and Ω1 that minimize functional Ec0;c1;Γ which is given in equation (5). Defining ω=ω0−ω1, we represent the domain Ω with binary indicator function f so that
(6)fx,y=χΩx,y1;x,y∈Ω0,i.e.,ω<0,0,otherwise.This functionf is the expected approximation of image g. Using f, the minimization of functional Ec0;c1;Γ from equation (5) can be seen as follows:
(7)minfx,y∈0,1f,ω+λfTV.The first term in equation (7) is the inner product of f and ω and is defined as f,ω=∫fω. The total variation pseudonorm is fTV, which equals to LΓ for binary indicator f=χΩ by coarea formula established by Fleming and Rishel [20].


## 2.2. Discrete Formulation of Segmentation Model

To discretize variational problem (7) on a grid of n×m pixels, we assume that N=n×m. Let vi∈R2 be a vector field for every i and suppose
(8)v→=vi→i=1N∈RN×2,then norm of v is given by
(9)v→1=∑i=1Nvi→.The total variation pseudonorm, forf=fiN/i=1∈RN, is defined as
(10)fTV=∇f1=∑i=1N∇fi,using equation (9). Thus, we obtain
(11)minfi∈0,1f,ω+λ∇f1.We assume that the pixels are indexed on a two-dimensional grid whereδ1=1,0 and δ2=0,1, the finite difference gradient operator is
(12)∇fii=1N=fi+δ1−fi,fi+δ2−fii=1N∈RN×2.We use periodic boundary conditions for simplicity. The inner productRN is
(13)f,ω=∑i=1Nfiωi.Equation (11) is a nonsmooth and nonconvex and usually direct minimization that leads to poor local minima. By replacing the binary constraint fi∈0,1 with a box constraint fi∈0,1, we obtain the convex formulation of the energy functional of equation (11). This defines the following finite dimensional convex problem:
(14)minfi∈0,1Nf,ω+λ∇f1.One can prove that this relaxation is exact, meaning that the minimizerf, when it is unique, is binary, f∈0,1. It means that Ω such that f=χΩ actually solves the original segmentation problem.


## 3. Proximal Splitting Technique

Proximal algorithms are used for solving large-scale, constrained, or distributed optimization problems [36, 37]. Many image processing problems can be formulated as
(15)minx∈RNf1x+⋯+fnx,where f1,⋯,fn:RN→−∞;+∞. These algorithms use splitting technique such that the functions f1,⋯,fn in equation (15) can be used to get an algorithm which is easily implementable. Usually, these algorithms are called proximal because each function in equation (15) is involved via its proximity operator. We propose to use proximal splitting scheme to solve our optimization problem, i.e., simplified Mumford-Shah model of equation (14).


### 3.1. Splitting

We introduce the variableu=∇f (http://www.numericaltours.com). The optimization problem (14) can be seen as follows:
(16)minz=f,u∈Z=RN×RN×2Fz+Hz,where Fz and Hz are defined as follows:
(17)Fz=Ff,u=f,ω+i0,1Nf+λu,Hz=Hf,u=lDf,u,D=z=f,u:u=∇f.Here, the constraints are defined using indicator function as follows:(18)l0,1Ns=0ifs∈0,1N,+∞otherwise,lDf,u=0iff,u∈D=z=f,u:u=∇f,+∞otherwise.


### 3.2. Proximal Mapping

The proximal mapping,proxγT:RN→RN, for any convex function T with γparameter and t,s∈RN(19)proxγTt=argmins12t−s2+γTs.Here,proxγHt=s∗ is the minimizer of function T. Proximal operators can be viewed as generalized projections (see Appendix A). Similarly, reflection is defined as
(20)rproxγTt=2proxγTt−t.


### 3.3. Douglas-Rachford Algorithm

The Douglas-Rachford (DR) algorithm minimizes functionals iteratively as in equation (16):
(21)minzFnz+Hnz.It is given thatFz and Hz are convex functions. Hence, we can compute their proximal mappings proxγF and proxγH, respectively. It is not required that the functions Fz and Hz are smooth. A DR iteration is (see Appendix C):
(22)z~k+1=1−μ2z~k+μ2rproxγnHrproxγFz~k,zk+1=rproxγnFz~k,where
(23)rproxγnFz=2proxγnFz−z.We can show that for any value ofγ>0,0<μ<2, and z~0,zk⟶z∗is a minimizer of Fn+Hn.


#### 3.3.1. Proximal Operator ofH

The proximal mapping ofH is the orthogonal projection on the convex set D:
(24)fn~,u~n=proxγnHf,u=ProjDfn,un.The following linear system of equations givesf~ and u~ because
(25)u~n=∇f~n,where
(26)fn~=IdN−∆−1fn−divun(see AppendixB).Here, by convention,∆=div∘∇ and div=−∇∗. We have
(27)fn~^m=h^nmKm,∀m=m1,m2,where
(28)Km=1+4sinπm1n2+4sinπm2n2.As we have discretized variational problem (7) on a grid of n×m pixels, m1 takes values 1⋯n and m2 takes the values 1⋯m.Also,hn=fn−divun and where the two-dimensional discrete Fourier transform of an image h is hn.


#### 3.3.2. Proximal Operator ofF

The functionFf,u can be written as
(29)Ffn,u=F0fn+λu1,where
(30)F0fn=fn,ω+l0,1Nfn.The proximal mapping ofFn is
(31)proxγnFfn,u=proxγ0Ff;proxγλ.u.We define the value ofλ>0. The proximal operator of F0 is
(32)proxγ0Ff;Proj0,1Nfn−γω,where
(33)Proj0,1Nhn=max0,min1,hn.The proximal operator of thel1−l2 norm .1 is a soft thresholding of the amplitude of the vector field:
(34)proxγ.1ui=max0,γuiui.


### 3.4. Implementation

The implementation is sketched in the following pseudocode. We use the following notations:g: input image of size X×Ysegments
: number of segmentsiterations: number of iterationslambda
λ: control parameter in Mumford-Shah functionalr
c
0
,
r
c
1: color values of randomly selected pixel0 and pixel1c
0
,
c
1: color values of pixel0 and pixel1 after iterations of k-meansω
,
ω
0
,
ω
1
:indicate regions Ω, Ω0, and Ω1f: indicator functionu: gradient of fE: value of functional E at each iterationThe following subprograms are also needed:f
f
t: fast Fourier transformgrad: gradienti
f
f
t
: inverse fourier transformdiv: divergenceAlgorithm 1: Minimization of MS model using DR algorithm.
MS (g, segments, max iterations, lambda) → f
rc0⟵rgb values of randomly selected pixel0
rc1⟵rgb values of randomly selected pixel1
c0,c1⟵ kmeansrc0,rc1
mu⟵1, gamma⟵1,
ω0⟵diffc0,g
ω1⟵diffc1,g
ω⟵diffω0,ω1
f⟵indicatorω
u⟵grad(f)
REPEAT
term1⟵1−mu/2∗f
term2⟵mu/2∗2∗proxf02∗projc1f,u1−f,gamma,ω
term3⟵2∗projc1f,u1−f
f1⟵term1+term2−term3
term1⟵1−mu/2∗u
term2⟵mu/2∗2∗proxL12∗projc2u−u,gamma∗lambda
term3⟵2∗projc2u−u
u1⟵term1+term2−term3
f⟵projc1f1,u1
u⟵projc2u1
Efn⟵sumω∗f+lambda∗sumsqrtsumu2
UNTILn≥2 AND (∆Eiter<0.1) OR (iter=maxiter)
Procedures
diff
a
,
b → c
FOR all i∈1,X DO
FOR allj∈1,Y DO
ci,j=ai,j−bi,j
indicator
w
→
f
FOR all i∈1,X DO
FOR allj∈1,Y DO
IFwi,j>0 THEN
fi,j=1
ELSE
fi,j=0
Auxiliary functions
projc1
t
f
1
,
t
u
1
→
i
f
f
t
K⟵gridk()
ifft=ifftffttf1−divtu1/K
gridk() → K
FOR all i∈1,X DO
FOR allj∈1,Y DO
Ki,j=1+4∗sinX∗π/i2+4∗sinY∗π/j2
projc2
f → gf
gf=gradf(f)
proxf0
f
,
g
a
m
m
a
,
w
FOR all i∈1,X DO
FOR allj∈1,Y DO
mi,j=max0,min1,fi,j−gamma∗wi,j


## 3.1. Splitting

We introduce the variableu=∇f (http://www.numericaltours.com). The optimization problem (14) can be seen as follows:
(16)minz=f,u∈Z=RN×RN×2Fz+Hz,where Fz and Hz are defined as follows:
(17)Fz=Ff,u=f,ω+i0,1Nf+λu,Hz=Hf,u=lDf,u,D=z=f,u:u=∇f.Here, the constraints are defined using indicator function as follows:(18)l0,1Ns=0ifs∈0,1N,+∞otherwise,lDf,u=0iff,u∈D=z=f,u:u=∇f,+∞otherwise.


## 3.2. Proximal Mapping

The proximal mapping,proxγT:RN→RN, for any convex function T with γparameter and t,s∈RN(19)proxγTt=argmins12t−s2+γTs.Here,proxγHt=s∗ is the minimizer of function T. Proximal operators can be viewed as generalized projections (see Appendix A). Similarly, reflection is defined as
(20)rproxγTt=2proxγTt−t.


## 3.3. Douglas-Rachford Algorithm

The Douglas-Rachford (DR) algorithm minimizes functionals iteratively as in equation (16):
(21)minzFnz+Hnz.It is given thatFz and Hz are convex functions. Hence, we can compute their proximal mappings proxγF and proxγH, respectively. It is not required that the functions Fz and Hz are smooth. A DR iteration is (see Appendix C):
(22)z~k+1=1−μ2z~k+μ2rproxγnHrproxγFz~k,zk+1=rproxγnFz~k,where
(23)rproxγnFz=2proxγnFz−z.We can show that for any value ofγ>0,0<μ<2, and z~0,zk⟶z∗is a minimizer of Fn+Hn.


### 3.3.1. Proximal Operator ofH

The proximal mapping ofH is the orthogonal projection on the convex set D:
(24)fn~,u~n=proxγnHf,u=ProjDfn,un.The following linear system of equations givesf~ and u~ because
(25)u~n=∇f~n,where
(26)fn~=IdN−∆−1fn−divun(see AppendixB).Here, by convention,∆=div∘∇ and div=−∇∗. We have
(27)fn~^m=h^nmKm,∀m=m1,m2,where
(28)Km=1+4sinπm1n2+4sinπm2n2.As we have discretized variational problem (7) on a grid of n×m pixels, m1 takes values 1⋯n and m2 takes the values 1⋯m.Also,hn=fn−divun and where the two-dimensional discrete Fourier transform of an image h is hn.


### 3.3.2. Proximal Operator ofF

The functionFf,u can be written as
(29)Ffn,u=F0fn+λu1,where
(30)F0fn=fn,ω+l0,1Nfn.The proximal mapping ofFn is
(31)proxγnFfn,u=proxγ0Ff;proxγλ.u.We define the value ofλ>0. The proximal operator of F0 is
(32)proxγ0Ff;Proj0,1Nfn−γω,where
(33)Proj0,1Nhn=max0,min1,hn.The proximal operator of thel1−l2 norm .1 is a soft thresholding of the amplitude of the vector field:
(34)proxγ.1ui=max0,γuiui.


## 3.3.1. Proximal Operator ofH

The proximal mapping ofH is the orthogonal projection on the convex set D:
(24)fn~,u~n=proxγnHf,u=ProjDfn,un.The following linear system of equations givesf~ and u~ because
(25)u~n=∇f~n,where
(26)fn~=IdN−∆−1fn−divun(see AppendixB).Here, by convention,∆=div∘∇ and div=−∇∗. We have
(27)fn~^m=h^nmKm,∀m=m1,m2,where
(28)Km=1+4sinπm1n2+4sinπm2n2.As we have discretized variational problem (7) on a grid of n×m pixels, m1 takes values 1⋯n and m2 takes the values 1⋯m.Also,hn=fn−divun and where the two-dimensional discrete Fourier transform of an image h is hn.


## 3.3.2. Proximal Operator ofF

The functionFf,u can be written as
(29)Ffn,u=F0fn+λu1,where
(30)F0fn=fn,ω+l0,1Nfn.The proximal mapping ofFn is
(31)proxγnFfn,u=proxγ0Ff;proxγλ.u.We define the value ofλ>0. The proximal operator of F0 is
(32)proxγ0Ff;Proj0,1Nfn−γω,where
(33)Proj0,1Nhn=max0,min1,hn.The proximal operator of thel1−l2 norm .1 is a soft thresholding of the amplitude of the vector field:
(34)proxγ.1ui=max0,γuiui.


## 3.4. Implementation

The implementation is sketched in the following pseudocode. We use the following notations:g: input image of size X×Ysegments
: number of segmentsiterations: number of iterationslambda
λ: control parameter in Mumford-Shah functionalr
c
0
,
r
c
1: color values of randomly selected pixel0 and pixel1c
0
,
c
1: color values of pixel0 and pixel1 after iterations of k-meansω
,
ω
0
,
ω
1
:indicate regions Ω, Ω0, and Ω1f: indicator functionu: gradient of fE: value of functional E at each iterationThe following subprograms are also needed:f
f
t: fast Fourier transformgrad: gradienti
f
f
t
: inverse fourier transformdiv: divergenceAlgorithm 1: Minimization of MS model using DR algorithm.
MS (g, segments, max iterations, lambda) → f
rc0⟵rgb values of randomly selected pixel0
rc1⟵rgb values of randomly selected pixel1
c0,c1⟵ kmeansrc0,rc1
mu⟵1, gamma⟵1,
ω0⟵diffc0,g
ω1⟵diffc1,g
ω⟵diffω0,ω1
f⟵indicatorω
u⟵grad(f)
REPEAT
term1⟵1−mu/2∗f
term2⟵mu/2∗2∗proxf02∗projc1f,u1−f,gamma,ω
term3⟵2∗projc1f,u1−f
f1⟵term1+term2−term3
term1⟵1−mu/2∗u
term2⟵mu/2∗2∗proxL12∗projc2u−u,gamma∗lambda
term3⟵2∗projc2u−u
u1⟵term1+term2−term3
f⟵projc1f1,u1
u⟵projc2u1
Efn⟵sumω∗f+lambda∗sumsqrtsumu2
UNTILn≥2 AND (∆Eiter<0.1) OR (iter=maxiter)
Procedures
diff
a
,
b → c
FOR all i∈1,X DO
FOR allj∈1,Y DO
ci,j=ai,j−bi,j
indicator
w
→
f
FOR all i∈1,X DO
FOR allj∈1,Y DO
IFwi,j>0 THEN
fi,j=1
ELSE
fi,j=0
Auxiliary functions
projc1
t
f
1
,
t
u
1
→
i
f
f
t
K⟵gridk()
ifft=ifftffttf1−divtu1/K
gridk() → K
FOR all i∈1,X DO
FOR allj∈1,Y DO
Ki,j=1+4∗sinX∗π/i2+4∗sinY∗π/j2
projc2
f → gf
gf=gradf(f)
proxf0
f
,
g
a
m
m
a
,
w
FOR all i∈1,X DO
FOR allj∈1,Y DO
mi,j=max0,min1,fi,j−gamma∗wi,j


## 4. Experimental Results

### 4.1. Image Dataset

The performance of our algorithm is evaluated using the partial Weizmann dataset [28] images depicting one object in the foreground. Currently, the image database contains a hundred grey level and hundred color images along with their ground truth segmentations, out of which we have used only color images. The dataset comprises of images that clearly have an object that is different from the background. Each image has three corresponding ground truth segmentations as shown in Figure 1.Figure 1
Segmentation results for selected two images using the Douglas-Rachford algorithm.


### 4.2. Evaluation Criterion

The performance of our algorithm is based on region-based error metrics defined in [29]. Considering S1 as segmentation obtained by an algorithm, S2 as one of the ground-truth, we calculate
(35)ES1,S2,pi=RS1,pi∩RS2,pi¯RS1,pi,a value in the range0,1. For a given pixel pi, we consider segment S1 in the solution and segment S2 in the ground truth that contain the pixel. If one segment is a proper subset of the other, then the pixel lies in an area of refinement, and the local error should be zero. Otherwise, the two regions overlap in an inconsistent manner and we should calculate the corresponding error. Here, RS;pi is the set of pixels corresponding to the region in segmentation S that contains pixel pi, and the local refinement error is defined as ES1,S2,pi.This local error measure is not symmetric, and it gives refinement in one direction only. It is to be noted thatES1,S2,pi=0 if S1 is a refinement of S2 at pixel pi, but not the opposite. Considering this local refinement error in each direction at each pixel, there are two methods to combine the values into an error measure for the entire image: global consistency error GCE that forces all local refinements to be in the same direction and local consistency error LCE that allows refinement in different directions in different parts of the image [38]. For a given n as the number of pixels, we have
(36)GCES1,S2=1nmin∑inES1,S2,pi,∑inES2,S1,pi,LCES1,S2=1n∑inminES1,S2,pi,ES2,S1,pi.The region-based consistency measure LCE is tolerant to refinement in either direction at each image pixel. If we simply replace the pixelwise minimum with a maximum, we get a measure that does not tolerate refinement at all and penalizes dissimilarity between segmentations proportional to the degree of region overlap. Applying to segmentationsS1 and S2, the bidirectional consistency error (BCE) is defined as
(37)BCES1,S2=1n∑inmaxES1,S2,pi,ES2,S1,pi.In addition, we can ask if segmentationS1 obtained by our algorithm is consistent with the collection of all ground truth segmentations Sk of that image. The BCE∗ measure is obtained by computing the minimum error at each pixel over each ground truth segmentation Sk:
(38)BCE∗S1=1n∑inminkmaxES1,Sk,pi,ESk,S1,pi.Structural similarity index [39] measures the visual quality between two images. It is calculated as a weighted combination of three factors luminance, contrast, and structural similarity. It can be expressed as
(39)SSIMx,y=lx,yα.cx,yβ.sx,yγ,where l is the luminance, c is the contrast, s is the structural similarity, and α, β, and γ are positive constants. The three factors are calculated separately as follows:
(40)lx,y=2μxμy+C1μx2+μy2+C1,cx,y=2σxσy+C2σx2+σy2+C2,sx,y=σxy+C3σxσy+C3,where μx and μy are local means, σx and σy are the standard deviations, and σxy is the cross-covariance between the two images x and y, respectively. If α=β=γ=1, then SSIM simplifies to
(41)SSIM=2μxμy+C1+2σxσy+C2μx2+μy2+C1+σx2+σy2+C2.We use SSIM to calculate the similarity between the image generated by the segmentation result, and the image generated from the ground truth segmentation. Here, we have implemented the proposed method for two-phase segmentation, which gives us foreground (object) and background. We have used SSIM to evaluate the foreground.


### 4.3. Results

We have executed the DR algorithm over 100 images from the Weizmann dataset [29], for various values of λ as 0.10, 0.25, 0.50, 0.75, and 1.00.λ is a parameter which is positive, and unfortunately, there is no fitness function to adjust it automatically. This is a drawback of the method. reg-KM and Chan-Vese methods also have similar limitation.For each of the 100 images, the segmentation results are compared with their three ground truths to obtain the error measures GCE andBCE∗. The average value of the three error values is taken as the final error. Out of these λ values, we found that the results using λ=0.75 and λ=1.00 are the best in terms of accuracy in most of the images. We elaborate these in Figure 2. It is very important that the error values are more concentrated towards zero. Figure 2(a) illustrates that for λ=0.75 and λ=1.00, around 69% of the images have GCE value of 0.5 or less. Similarly, Figure 2(b) illustrates that for λ=0.75 and λ=1.00, around 70% of the images have BCE∗ value of 0.5 or less.Figure 2
Error metrics for different values ofλ for 100 images.
(a)
(b)Figure1 shows the original image, the three human segmentations comprising the ground truth, the initial segmentation obtained by the k-means clustering technique, and the result of the Douglas-Rachford algorithm for two selected images from the dataset. It is clear that when λ is large, the minimization of E leads to fewer boundaries and coarser segmentation and when λ is small, we obtain a finer segmentation. Minimization of E is achieved when the difference of E between two iterations is less than 0.1 tolerance. We vary λ from 0.1 to 1.00 to have segmentation from finer to coarser, as established in [15]. Figure 3 compares the performance and efficiency of the iterations with different λ values. It shows that the value of the minimization of energy functional keeps decreasing when iterated long enough. We observe that for λ values of 0.75 and 1.00, minimization of energy functional reaches the lowest values. At the same time for λ value of 0.75, it converges faster.Figure 3
Behaviour ofE for different values of λ with tolerance=0.1.
(a)
(b)Further, for comparison with other methods, we select three images from the dataset such that their segmentation is not an easy task. Figure4 shows these three selected images and their different ground truths. For the church image, it is difficult to determine which is the object, the church, or the signboard. This confusion is also seen in the ground truths. The same is observed in the case of a lion, where either the lion or the ball, or both, is considered as the object. The segmentation of caterpillar is the most difficult because it has stripes which can act as an artificial boundary.Figure 4
Images and their ground truths.Figure5 shows segmentations obtained by the reg-KM algorithm [14], Chan-Vese algorithm [8], and the proposed DR algorithm. The results are optimized as follows. Results of reg-KM are with the regularization parameter 0.5. The Chan-Vese algorithm has been used with various initial contours like a small circle, large circle, and user-defined box. We tested all these and chose to use small circles (see Figure 6) as it provided the best segmentation for these images. In our algorithm, we have used the k-means clustering method for obtaining initial segmentation and then Mumford-Shah functional minimized by Douglas-Rachford using the parameter λ=0.5.Figure 5
Segmentation results using reg-KM, Chan-Vese, and Douglas-Rachford (DR) algorithms.Figure 6
Result of the Chan-Vese algorithm.Table1 shows the numerical comparison of reg-KM, Chan-Vese, and the proposed DR algorithm. The results show that the reg-KM algorithm has worse accuracy in all three cases which we obtained through the error evaluation criteria mentioned in Section 4.2. The results of the Chan-Vese method and the proposed DR method have similar accuracy; Chan-Vese is better with caterpillar while DR is better with lion.Table 1
Comparison of reg-KM, Chan-Vese algorithm, and proposed DR algorithm in terms of BCE∗ error.
Image
reg-KM
Chan-Vese
Proposed DR
Church
0.2732
0.1563
0.1552
Lion
0.2682
0.1893
0.1435
Caterpillar
0.3919
0.3327
0.3723Table2 summarizes the SSIM values for the three methods: reg-KM, Chan-Vese, and the proposed DR. In case of all images, the proposed DR algorithm achieves significantly higher SSIM values which indicates its superiority over the compared methods.Table 2
Comparison of reg-KM, Chan-Vese algorithm, and proposed DR algorithm in terms of SSIM.
Image
reg-KM
Chan-Vese
Proposed DR
Church
0.6864
0.6878
0.7334
Lion
0.7026
0.6880
0.7948
Caterpillar
0.1238
0.1187
0.2059Table3 summarizes the processing times of these algorithms. We can see that the proposed DR algorithm is an order of magnitude faster than Chan-Vase. The three images are segmented in about 1.5 seconds, on average, while Chan-Vese takes 40 seconds per image, on average. The reg-KM technique does not give very good segmentation but is definitely fast in terms of execution time. However, the DR algorithm uses a minimization technique and therefore it takes more time but gives much better results than reg-KM and nearly the same results when compared to the popular Chan-Vese method.Table 3
Comparison of reg-KM, Chan-Vese algorithm, and proposed DR algorithm in terms of execution time (seconds).
Image
reg-KM
Chan-Vese
Proposed DR
Church
<1
60.00
2.02
Caterpillar
48.42
1.42
Lion
10.63
1.38


## 4.1. Image Dataset

The performance of our algorithm is evaluated using the partial Weizmann dataset [28] images depicting one object in the foreground. Currently, the image database contains a hundred grey level and hundred color images along with their ground truth segmentations, out of which we have used only color images. The dataset comprises of images that clearly have an object that is different from the background. Each image has three corresponding ground truth segmentations as shown in Figure 1.Figure 1
Segmentation results for selected two images using the Douglas-Rachford algorithm.


## 4.2. Evaluation Criterion

The performance of our algorithm is based on region-based error metrics defined in [29]. Considering S1 as segmentation obtained by an algorithm, S2 as one of the ground-truth, we calculate
(35)ES1,S2,pi=RS1,pi∩RS2,pi¯RS1,pi,a value in the range0,1. For a given pixel pi, we consider segment S1 in the solution and segment S2 in the ground truth that contain the pixel. If one segment is a proper subset of the other, then the pixel lies in an area of refinement, and the local error should be zero. Otherwise, the two regions overlap in an inconsistent manner and we should calculate the corresponding error. Here, RS;pi is the set of pixels corresponding to the region in segmentation S that contains pixel pi, and the local refinement error is defined as ES1,S2,pi.This local error measure is not symmetric, and it gives refinement in one direction only. It is to be noted thatES1,S2,pi=0 if S1 is a refinement of S2 at pixel pi, but not the opposite. Considering this local refinement error in each direction at each pixel, there are two methods to combine the values into an error measure for the entire image: global consistency error GCE that forces all local refinements to be in the same direction and local consistency error LCE that allows refinement in different directions in different parts of the image [38]. For a given n as the number of pixels, we have
(36)GCES1,S2=1nmin∑inES1,S2,pi,∑inES2,S1,pi,LCES1,S2=1n∑inminES1,S2,pi,ES2,S1,pi.The region-based consistency measure LCE is tolerant to refinement in either direction at each image pixel. If we simply replace the pixelwise minimum with a maximum, we get a measure that does not tolerate refinement at all and penalizes dissimilarity between segmentations proportional to the degree of region overlap. Applying to segmentationsS1 and S2, the bidirectional consistency error (BCE) is defined as
(37)BCES1,S2=1n∑inmaxES1,S2,pi,ES2,S1,pi.In addition, we can ask if segmentationS1 obtained by our algorithm is consistent with the collection of all ground truth segmentations Sk of that image. The BCE∗ measure is obtained by computing the minimum error at each pixel over each ground truth segmentation Sk:
(38)BCE∗S1=1n∑inminkmaxES1,Sk,pi,ESk,S1,pi.Structural similarity index [39] measures the visual quality between two images. It is calculated as a weighted combination of three factors luminance, contrast, and structural similarity. It can be expressed as
(39)SSIMx,y=lx,yα.cx,yβ.sx,yγ,where l is the luminance, c is the contrast, s is the structural similarity, and α, β, and γ are positive constants. The three factors are calculated separately as follows:
(40)lx,y=2μxμy+C1μx2+μy2+C1,cx,y=2σxσy+C2σx2+σy2+C2,sx,y=σxy+C3σxσy+C3,where μx and μy are local means, σx and σy are the standard deviations, and σxy is the cross-covariance between the two images x and y, respectively. If α=β=γ=1, then SSIM simplifies to
(41)SSIM=2μxμy+C1+2σxσy+C2μx2+μy2+C1+σx2+σy2+C2.We use SSIM to calculate the similarity between the image generated by the segmentation result, and the image generated from the ground truth segmentation. Here, we have implemented the proposed method for two-phase segmentation, which gives us foreground (object) and background. We have used SSIM to evaluate the foreground.


## 4.3. Results

We have executed the DR algorithm over 100 images from the Weizmann dataset [29], for various values of λ as 0.10, 0.25, 0.50, 0.75, and 1.00.λ is a parameter which is positive, and unfortunately, there is no fitness function to adjust it automatically. This is a drawback of the method. reg-KM and Chan-Vese methods also have similar limitation.For each of the 100 images, the segmentation results are compared with their three ground truths to obtain the error measures GCE andBCE∗. The average value of the three error values is taken as the final error. Out of these λ values, we found that the results using λ=0.75 and λ=1.00 are the best in terms of accuracy in most of the images. We elaborate these in Figure 2. It is very important that the error values are more concentrated towards zero. Figure 2(a) illustrates that for λ=0.75 and λ=1.00, around 69% of the images have GCE value of 0.5 or less. Similarly, Figure 2(b) illustrates that for λ=0.75 and λ=1.00, around 70% of the images have BCE∗ value of 0.5 or less.Figure 2
Error metrics for different values ofλ for 100 images.
(a)
(b)Figure1 shows the original image, the three human segmentations comprising the ground truth, the initial segmentation obtained by the k-means clustering technique, and the result of the Douglas-Rachford algorithm for two selected images from the dataset. It is clear that when λ is large, the minimization of E leads to fewer boundaries and coarser segmentation and when λ is small, we obtain a finer segmentation. Minimization of E is achieved when the difference of E between two iterations is less than 0.1 tolerance. We vary λ from 0.1 to 1.00 to have segmentation from finer to coarser, as established in [15]. Figure 3 compares the performance and efficiency of the iterations with different λ values. It shows that the value of the minimization of energy functional keeps decreasing when iterated long enough. We observe that for λ values of 0.75 and 1.00, minimization of energy functional reaches the lowest values. At the same time for λ value of 0.75, it converges faster.Figure 3
Behaviour ofE for different values of λ with tolerance=0.1.
(a)
(b)Further, for comparison with other methods, we select three images from the dataset such that their segmentation is not an easy task. Figure4 shows these three selected images and their different ground truths. For the church image, it is difficult to determine which is the object, the church, or the signboard. This confusion is also seen in the ground truths. The same is observed in the case of a lion, where either the lion or the ball, or both, is considered as the object. The segmentation of caterpillar is the most difficult because it has stripes which can act as an artificial boundary.Figure 4
Images and their ground truths.Figure5 shows segmentations obtained by the reg-KM algorithm [14], Chan-Vese algorithm [8], and the proposed DR algorithm. The results are optimized as follows. Results of reg-KM are with the regularization parameter 0.5. The Chan-Vese algorithm has been used with various initial contours like a small circle, large circle, and user-defined box. We tested all these and chose to use small circles (see Figure 6) as it provided the best segmentation for these images. In our algorithm, we have used the k-means clustering method for obtaining initial segmentation and then Mumford-Shah functional minimized by Douglas-Rachford using the parameter λ=0.5.Figure 5
Segmentation results using reg-KM, Chan-Vese, and Douglas-Rachford (DR) algorithms.Figure 6
Result of the Chan-Vese algorithm.Table1 shows the numerical comparison of reg-KM, Chan-Vese, and the proposed DR algorithm. The results show that the reg-KM algorithm has worse accuracy in all three cases which we obtained through the error evaluation criteria mentioned in Section 4.2. The results of the Chan-Vese method and the proposed DR method have similar accuracy; Chan-Vese is better with caterpillar while DR is better with lion.Table 1
Comparison of reg-KM, Chan-Vese algorithm, and proposed DR algorithm in terms of BCE∗ error.
Image
reg-KM
Chan-Vese
Proposed DR
Church
0.2732
0.1563
0.1552
Lion
0.2682
0.1893
0.1435
Caterpillar
0.3919
0.3327
0.3723Table2 summarizes the SSIM values for the three methods: reg-KM, Chan-Vese, and the proposed DR. In case of all images, the proposed DR algorithm achieves significantly higher SSIM values which indicates its superiority over the compared methods.Table 2
Comparison of reg-KM, Chan-Vese algorithm, and proposed DR algorithm in terms of SSIM.
Image
reg-KM
Chan-Vese
Proposed DR
Church
0.6864
0.6878
0.7334
Lion
0.7026
0.6880
0.7948
Caterpillar
0.1238
0.1187
0.2059Table3 summarizes the processing times of these algorithms. We can see that the proposed DR algorithm is an order of magnitude faster than Chan-Vase. The three images are segmented in about 1.5 seconds, on average, while Chan-Vese takes 40 seconds per image, on average. The reg-KM technique does not give very good segmentation but is definitely fast in terms of execution time. However, the DR algorithm uses a minimization technique and therefore it takes more time but gives much better results than reg-KM and nearly the same results when compared to the popular Chan-Vese method.Table 3
Comparison of reg-KM, Chan-Vese algorithm, and proposed DR algorithm in terms of execution time (seconds).
Image
reg-KM
Chan-Vese
Proposed DR
Church
<1
60.00
2.02
Caterpillar
48.42
1.42
Lion
10.63
1.38


## 5. Conclusions

In this paper, we proposed a reduced two-phase Mumford-Shah model to segment images with one object. We obtained initial segmentation by the k-means clustering technique and further minimized it by our implementation of the Douglas-Rachford algorithm. In k-means,k is a user-selected parameter about how many segments. If k needs to be automatically selected, then it should also be included in the Mumford-Shah functional in (2). Possible solutions how to do it, we refer to solutions in the clustering context [40]. Experiments with various error metrics LCE, GCE, and BCE∗ show that 70 percent of the segmentations keep the error values below 0.5. We also compare our algorithm to Chan-Vese implementation which is one of the classical methods to optimize the Mumford-Shah functional. Our segmentation is slightly better in two of the three selected images, while the main benefit is that our method is significantly faster than the Chan-Vese algorithm. Secondly, we compare our algorithm with a recent k-means variant, reg-KM. According to the BCE∗ error metric, our segmentations are significantly better. However, reg-KM is much faster as it has a statistical approach of clustering but it does not give smooth segmentation with convex boundaries. Our algorithm balances time as well as quality of segmentation. As a future research, we consider adopting the Mumford-Shah model in k-means context.



---

*Source: 6618505-2021-04-13.xml*