\documentclass[fleqn,10pt]{wlscirep}
\usepackage{cite}
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{algorithmic}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{textcomp}
\usepackage{comment}
\def\BibTeX{{\rm B\kern-.05em{\sc i\kern-.025em b}\kern-.08em
    T\kern-.1667em\lower.7ex\hbox{E}\kern-.125emX}}
\usepackage{geometry}
\usepackage{xcolor,colortbl}
\geometry{a4paper, top=3.0cm, left=3.0cm, right=3.0cm, bottom=3.0cm}
\usepackage{url}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\urldef{\mailsa}\path|{afzal85, dr.haider}@uop.edu.pk, noor2knoor@gmail.com, {mohd.zakarya, izaz, aftab.ahmed.khan}@awkum.edu.pk, dr.hameed@ubuner.edu.pk, m.haleem@kardan.edu.af |

\author[1]{Afzal Rahman}
\author[1]{Haider Ali}
\author[2,+]{Noor Badshah}
\author[3,+]{Muhammad Zakarya}
\author[4]{Hameed Hussain}
\author[3,+]{Izaz Ur Rahman}
\author[3]{Aftab Ahmed}
\author[5,+,*]{Muhammad Haleem}
\affil[1]{Department of Mathematics, University of Peshawar, Pakistan}
\affil[2]{Department of Basic Sciences, University of Engineering and Technology Peshawar, Pakistan}
\affil[3]{Department of Computer Science, Abdul Wali Khan University, Mardan, Pakistan}
\affil[4]{Department of Computer Science, University of Buner, Pakistan}
\affil[5]{Department of Computer Science, Kardan University, Kabul, Afghanistan}

\affil[*]{m.haleem@kardan.edu.af}

\affil[+]{these authors contributed equally to this work}

\title{Power mean based image segmentation in the presence of noise}
\keywords{Convex functional, Fuzzy set, Image segmentation, Outliers, Power mean}
\begin{abstract}
	In image segmentation and in general in image processing, noise and outliers distort contained information posing in this way a great challenge for accurate image segmentation results. To ensure a correct image segmentation in presence of noise and outliers, it is necessary to identify the outliers and isolate them during a denoising pre-processing or impose suitable constraints into a segmentation framework. In this paper, we impose suitable removing outliers constraints supported by a well-designed theory in a variational framework for accurate image segmentation. We investigate a novel approach based on the power mean function equipped with a well established theoretical base. \color{red}The power mean function has the capability to distinguishes between true image pixels and outliers and, therefore, is robust against outliers\color{black}. To deploy the novel image data term and to guaranteed unique segmentation results, a fuzzy-membership function is employed in the proposed energy functional. Based on qualitative and quantitative extensive analysis on various standard data sets, it has been observed that the proposed model works well in images having multi-objects with high noise and in images with intensity inhomogeneity in contrast with the latest and \color{red}state-of-the-art \color{black}models.
\keywords{Convex functional, Fuzzy set, Image segmentation, Outliers, Power mean.}	
\end{abstract}

\begin{document}
\flushbottom
\begin{comment}
	
\history{Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.}
\doi{10.1109/ACCESS.2017.DOI}

\title{Power Mean Based Image Segmentation in the Presence of Noise}
\author{\uppercase{Afzal Rahman}\authorrefmark{1}, %\IEEEmembership{Fellow, IEEE},
	\uppercase{Haider Ali\authorrefmark{1}, Noor Badshah}\authorrefmark{2},
	%\IEEEmembership{Member, IEEE}
	\uppercase{Muhammad Zakarya}\authorrefmark{3}, \uppercase{Aftab Ahmed}\authorrefmark{3}, \uppercase{Hameed Hussain}\authorrefmark{4}, \uppercase{Izaz Ur Rahman}\authorrefmark{3}, \uppercase{Muhammad Haleem}\authorrefmark{5}
}
\address[1]{Department of Mathematics, University of Peshawar, Pakistan (e-mail: \{afzal85,  dr.haider\}@uop.edu.pk)}
\address[2]{Department of Basic Sciences, University of Engineering and Technology Peshawar, Pakistan (e-mail: noor2knoor@gmail.com)}
\address[3]{Department of Computer Science, Abdul Wali Khan University, Mardan, Pakistan
	(e-mail: \{mohd.zakarya, aftab.ahmed.khan, izaz\}@awkum.edu.pk)}
\address[4]{Department of Computer Science, University of Buner, Pakistan (e-mail: ham.hamdard@gmail.com)}
\address[5]{Department of Computer Science, Kardan University, Kabul, Afghanistan (e-mail: m.haleem@kardan.edu.af)}

\end{comment}

\maketitle

%\author{\uppercase{First A. Author}\authorrefmark{1}, \IEEEmembership{Fellow, IEEE},
%\uppercase{Second B. Author\authorrefmark{2}, and Third C. Author,
%Jr}.\authorrefmark{3},
%\IEEEmembership{Member, IEEE}}
%\address[1]{National Institute of Standards and 
%Technology, Boulder, CO 80305 USA (e-mail: author@boulder.nist.gov)}
%\address[2]{Department of Physics, Colorado State University, Fort Collins, 
%CO 80523 USA (e-mail: author@lamar.colostate.edu)}
%\address[3]{Electrical Engineering Department, University of Colorado, Boulder, CO 
%80309 USA}
%\tfootnote{This paragraph of the first footnote will contain support 
%information, including sponsor and financial support acknowledgment. For 
%example, ``This work was supported in part by the U.S. Department of 
%Commerce under Grant BS123456.''}

%\markboth
%{A. Rahman \headeretal: Power Mean Based Image Segmentation in the Presence of Noise}
%{A. Rahman \headeretal: Power Mean Based Image Segmentation in the Presence of Noise}

%\corresp{Corresponding author: Muhammad Haleem (e-mail: m.haleem@kardan.edu.af).}


%\titlepgskip=-15pt

%\maketitle

\section{Introduction}\label{sec:introduction}
Image segmentation is a fundamental step in computer vision and in digital image processing. The main idea of image segmentation is to visualize meaningful objects in a given scene or image \cite{8,7} linked to many important fields such as medical imaging, object detection, video, traffic control systems, surveillance, automated surgeries, and so on \cite{1,10,11}.
%Few applications include content-based image retrieval, machine vision, medical imaging, automated surgeries, object detection, recognition tasks, traffic control systems, video surveillance and so on \cite{1,10,11}. 
\color{red}Several state-of-the-art approaches for image segmentation exist; and some well-known methods include clustering \cite{1}, thresholding \cite{2}, edge detection and region-based models \cite{18,27,28,29}, Markov random fields \cite{18,27}, and stochastic methods \cite{28,29}, etc. \color{black}However, images are diverse in nature, and frequently happens that one model working for a particular class or type of images \color{red}but \color{black} may not properly work for other types. Some well-known factors which normally affect the performance of the segmentation models are noise and outliers, intensity in-homogeneity and low contrast between background and foreground in a given image \cite{3,10}. \color{red}To cope with these issues, and to achieve accurate image segmentation, the active contour framework is a very popular technique due to its flexibility of allowing and imposing the desired constraints. \color{black}Moreover, the availability of robust implementations, such as efficient optimization, and fast numerical methods is crucial. The main idea in active contour methods \cite{8,20,18,27,28,29,31,33,34,36,Wang} is to allow dynamical curves to move autonomously on a given image which locates boundaries of the objects/regions therein. 
\par
The active contour models mainly use the concept of variational calculus \color{red}\cite{5,8,9,13,15,24}\color{black}, that is functional optimization. It can be easily observed that the variational models for image processing in general and for segmentation, in particular, derive an energy functional which is minimized to get the desired results. The active contour models can be mainly divided into three categories, region-based \cite{8,20,31,33,34,36,Joshi}, edge-based models \cite{18,27,28,29} and region and edge based models \cite{selec}. A milestone variational model for segmentation purposes introduced by Mumford-Shah (MS ) \cite{8} aims to obtain a smooth cartoon image that leads to edge detection. It is important to note that the design of the MS model is for ideal images, that is, images without noise, outliers and intensity in-homogeneity.  Moreover, the direct implementation of this model is not feasible \cite{9}. To easily implement the MS model, the Chan {\it et al.} (CV) \cite{28} restricted the MS model to a piece-wise function reconstruction (two phases). By phase here we mean the set of homogeneous intensity pixels which can be easily distinguished from other sets of homogeneous intensity pixels in a given image. We should emphasize that the CV model ignores the presence of noise and other factors such as intensity in-homogeneity \cite{13}. To improve CV model, many techniques has been introduced in the last decades. In contrast with CV model, Li {\it et al.} \cite{13a} proposed the Local Binary Fitting (LBF) model which performs much better and resembles \color{red} than \color{black} the MS model by carrying out the segmentation via approximating given image with two locally smooth functions.
% The (LBF) model Li {\it et al.} tackles intensity in-homogeneity but not noise. This phenomenon can be observed by analyzing the data fitting term 
%\begin{equation}\label{3}
%	\displaystyle{\int_{\Omega} K_{\sigma}\ast|\mathbf{u}(x,y)-c_1|^{2}dx dy+\int_{\Omega}K_{\sigma}\ast|\mathbf{u}(x,y)-c_2|^{2}dx dy}, 
%\end{equation}
%where $K_{\sigma}$ is Gaussian kernel. In a discrete and local sense this data term is also based on the least square based objective function
%\begin{equation}\label{4}
%	\displaystyle\sum_{N_x}(\mathbf{u}(i,j)-{\textbf{c}})^2, 
%\end{equation}
%in local neighborhood $N_x$ whose minimum is also the sample mean $\textbf{c}=\bar{\mathbf{x}}$. In other words, the LBF model uses the concept of the CV model but in local neighborhoods throughout image domain $\Omega$. This leads to wider image intensity variation in small patches but on the other hand, it is more prone to noise and outliers as compared to the CV model. In this way, the fitting term takes into account the image intensity variance in small patches but on the other hand, it is more prone to noise and outliers as compared to the CV model. Moreover, this model is not convex so consequently one need tune several times the position of initial guess for the desired results. 

%Moreover, this model is not convex, consequently, the initial guess position requires several times tuning to obtain desired results.
\par 
Mondal {\it et al.}  \cite{36} introduced a new model based on a trade-off local and global information. The model shows good performance for images with intensity in-homogeneity and presence of noise and outliers. The experimental results shown in different data sets indicate that the model can deal with intensity in-homogeneity images, blurred boundary or discontinuous edges, and presence of moderated noise. Chuang {\it et al.} \cite{Chuang} and Tripathy {\it et al.} \cite{Tripathy} introduced the models which work better for the noisy MRI images, but may not work well for the images with intensity in-homogeneity. This model is not convex so consequently one need tune several times the position of initial guess for the desired results. Recently, \color{red}Wu {\it et al.} \cite{15}\color{black} proposed a convex variational segmentation model based on the concept of coefficient of variation (CoV) by ignoring the factor and presence of noise and outliers in given images. This argument is validated by the fact that their CoV based image data fitting term is sum of squares divided by sum of image intensity. 
%$\displaystyle{\int_{\textrm{inside}(\Gamma)}\frac{|\mathbf{u}(x,y)-c_1|^{2}}{{a_1}^2}dx dy+\int_{\textrm{outside}(\Gamma)}\frac{|\mathbf{u}(x,y)-c_2|^{2}}{{c_2}^{2}}dx dy}$ In discrete sense this image data term is based on the concept of squared CoV, ${CoV}^2=\displaystyle{\sum_{i,j}\frac{(\mathbf{u}(i,j)-{a})^2}{a^2}}$ whose minimum turns out to be $a= \frac{\displaystyle\sum_{i,j}{\mathbf{u}(i,j)}^2}{\displaystyle\sum_{i,j}{\mathbf{u}(i,j)}}$. 
Similar to the average fitting term in the CV model the value of the average fitting term in the \color{red}Wu {\it et al.} \cite{15}\color{black} model shows sensitivity to noise and outliers \cite{41}. To further improve this model, Wu {\it et al.} \cite{20} proposed an active contour model incorporating a kernel metric, which is robust, stable, and works well for images with low noise and outliers. Ali et al. \cite{35} introduced Lehmer's type generalized mean
\begin{equation}\label{5}
	\frac{\sum_{i,j}{\mathbf{u}(i,j)}^p}{\displaystyle\sum_{i,j}{\mathbf{u}(i,j)^{p-1}}}, 
\end{equation}
\color{red}in a segmentation \color{black}framework, where $p$ is any real number. Although this average is very effective in multi-region segmentation and suitable to different image intensity backgrounds, \color{red}however \color{black}it requires further analysis to tackle noise and outliers. 
\par
\color{red}From the above discussion, we can observe \color{black}that most of the variational region-based image segmentation models, \color{red}in the existing literature, \color{black}are based on the least square function. \color{red}In fact, this \color{black}forces the fit of the data to a piece-wise function of the mean intensity values of the foreground and background. \color{red}Moreover, they are unable to fully discriminate the noise and intrinsic intensities in the images. \color{black}This is one the main reasons that the aforementioned models and \color{red}similar frameworks \color{black}are unable to correctly and appropriately segment noisy and outliers affected images \cite{35}. 
\par
In this article, we mainly focus to design an efficient image data fitting term based on a novel objective \color{red}function, as given by (Eq. \ref{6}).\color{black}
\begin{equation}\label{6}
	\displaystyle\left(\frac{1}{\mid\Omega\mid}\sum_{i,j}(||\mathbf{u}(i,j)-a||_{2}^{2})^{p}\right)^{\frac{1}{p}}. 
\end{equation}
As further will be explained in Section \ref{proposed}, \color{red}this term is robust against the outliers by giving very fewer weights to outliers and noise in contrast compare with the traditional and old objective function which gives equal or almost equal weights to outliers and true image pixels \cite{37}. Moreover, besides the new data fitting term of the proposed model, a fuzzy level set function is employed which has two main benefits over the traditional level set function. Firstly, a single fuzzy function can capture more than one phase or objects of different intensities at the same time \cite{4,313}. Secondly, it plays an important role in efficiently imposing constraints for implementing convexity. This lead to non-dependence of the initial guess. Furthermore,  for a deeper understanding of the proposed model, the mathematical analysis is presented. For the regularization of the fuzzy membership function, the Gaussian smoothing filtering is employed \color{black}. Following are the major contributions of this work:
\begin{itemize}
	\item we impose suitable removing outliers constraints supported by a well-designed theory in a variational framework for accurate image segmentation;
	\item we investigate a novel approach based on the power mean function equipped with a well established theoretical base; 
	\item to guarantee unique segmentation results, a fuzzy-membership function is employed in the proposed energy functional; and
	\item extensive analysis on various standard data sets, it has been observed that the proposed model works well in images having multi-objects with high noise.
\end{itemize}
\par
The rest of the paper is organized as follows. In Section \ref{relwork}, we give a brief review of related segmentation models. The design and analysis of the proposed novel model are presented in Section  \ref{proposed}. In Section  \ref{exp}, a comprehensive experimental analysis is carried out both qualitative and quantitatively for types of outdoor natural, synthetic and medical images compared to existing and latest state-of-the-art segmentation techniques. %For validated comparisons the test images are taken from Berkeley’s data of various nature. It has been observed that the proposed novel model performs far better than existing and latest state-of-the-art segmentation techniques.
Final remarks and conclusions are made in Section  \ref{con}. 

\section{Related works}\label{relwork}
%say something

\subsection{Active contours without edges (CV)} 
To easily implement the MS model \cite{8}, Chan {\it et al.} (CV) \cite{28} restricted the MS model to a piecewise function reconstruction (two phases). Chan {\it et al.} \cite{28} \color{red}considered \color{black}a piecewise constant function which divides the image into different homogeneous regions representing the foreground and background \cite{WuLGC21}. %The Chan and Vese (CV) model is initially derived from the functional (MS)%it is region based model and famous for the image of two regions having different intensities,this model is based on the technique of curve evolution, mumford-shah functional for segmentation and level set.
For the image $\mathbf{u}$, %an image defined on $\Omega\subset\Re^{2}$ %image domain. let C be the evolving curve .In \cite{ACV}  a region based model is given whose
the minimization energy functional is given by:
\begin{align} \label{7}
	F^{CV}(a_1,a_2,\Gamma)&={\mu}\ length(\Gamma)\nonumber\\&+\lambda_1\int_{inside(\Gamma)}|\mathbf{u}(x,y)-c_1|^2dxdy \nonumber\\&+
	\lambda_2\int_{outside(\Gamma)}|\mathbf{u}(x,y)-c_2|^2dxdy,
\end{align}
where $\lambda_1,\ \lambda_2,\ \mu \geq 0$ are constants which tune the weight between the smoothing and the fitting terms. $\Gamma$ is the contour, and $c_1$, $c_2$ are average intensities of given image $I_{0}(x,y)$ for foreground and background, respectively. This is a non-convex model, so consequently one need to tune several times the position of initial guess for the desired results \cite{WuH16}.

%%%%%%%exclude
\subsection{Fuzzy energy-based minimization (FEBM)}
Given an image, $\mathbf{u}(x,y)$ in a spacial domain $\Omega$ Krinidis et al. \cite{313} proposed a segmentation model based on fuzzy function embedded in active contour variational framework:
\begin{align}\label{8}
	&F(\Gamma,c_1,c_2,v)=\mu  length(\Gamma)\nonumber\\&+\eta_{1}\int_{\Omega}[\mathbf{z}(x,y)]^{m}|\mathbf{u}(x,y)-c_1|^{2}dx dy \nonumber\\ &+\eta_{2} \int_{\Omega}[1-\mathbf{z}(x,y)]^{m}|\mathbf{u}(x,y)-c_2|^{2}dx dy
\end{align}
where the constants $c_1$, $c_2$  stand for average values inside and outside  the contour $\Gamma$, respectively, $m$ is the weight exponent (normally taking the value $2$), $\eta_1,\eta_2>0$ and $\mu\geq0$ are constants. The function $\mathbf{z}(x,y)\in[0,1]$ is the fuzzy membership function representing the membership degree of $\mathbf{u}(x,y)$  inside the $\Gamma$ and $1-\mathbf{z}(x,y)$ is the membership degree of $\mathbf{u}(x,y)$ outside the $\Gamma$.
%\par
%Pseudo level set $\mathbf{z}(x,y)$ formulation which based on the membership values and active contour related as follow:
%\begin{align} \nonumber
%\left\{ \begin{array}{lll}
%\mathcal{C}=\{(x,y)\in\Omega: \mathbf{z}(x,y)=0.5\},\\
%\textit{inside}(\mathcal{C})=\{(x,y)\in\Omega: \mathbf{z}(x,y)>0.5\},\\
%\textit{outside}(\mathcal{C})=\{(x,y)\in\Omega: \mathbf{z}(x,y)<0.5\}.
%\end{array}\right.
%\end{align}
%From Eq. (\ref{FEBM}), through minimization with respect to $c_1$ and $c_2$, we obtained:
%\begin{align}\nonumber
%a_1={\int_\Omega (\mathbf{z}(x,y))^{m} \mathbf{u}(x,y)dxdy\over \int_\Omega (\mathbf{z}(x,y))^{m}dxdy}
%\end{align}
%\begin{align}\nonumber
%a_2={\int_\Omega (1-\mathbf{z}(x,y))^{m} \mathbf{u}(x,y) dxdy\over \int_\Omega (1-\mathbf{z}(x,y))^{m}dxdy}.
%\end{align}
%The updated fuzzy membership function $\mathbf{z}(x,y)$ is:
%\begin{align}\label{fuzzy1}
%\mathbf{z}(x,y)=\frac{1}{1+(\frac{\eta_{1(\mathbf{u}(x,y)-c_1)^{2}}}{\eta_2(\mathbf{u}(x,y)-c_2)^{2}})^{\frac{1}{m-1}}}.
%\end{align}
For a fast convergence of the minimization problem in Eq. (\ref{8}) the authors use a fast algorithm as proposed by Song and Chan \cite{16}. This model can segment images with multi-objects, different intensity variations objects, and hazy boundaries, however, it may not properly segment noisy images. The reason is that this model uses the same conventional least square objective function which fits the data to the mean value of the foreground and background. % ($\displaystyle\frac{1}{\mid\Omega\mid}\displaystyle\sum_{i,j}(\mathbf{u}(i,j)-{a})^2$ ) based image data terms as used by earlier most of the similar models.
%%Initial guess 1

\subsection{A convex variational level set model for image segmentation (CVMS)}
\color{red}Wu {\it et al.} \cite{15}\color{black} proposed a strictly convex functional for two-phase image segmentation:
\begin{align}\label{9}
	F_{WH}(\psi)&=\eta\int_\Omega\frac{(\mathbf{u}(x,y)-c_1)^{2}}{{c_1}^{2}}(\psi(x,y)+1)^{2}dx dy \nonumber\\
	&+\int_\Omega\frac{(\mathbf{u}(x,y)-c_2)^{2}}{{c_2}^{2}}(\psi(x,y)-1)^{2}dx dy
\end{align}
where $\psi$ denote the level set function \cite{11}, $\eta>0$ is a parameter. 
%If the level set function $\psi$ is non-empty inside and outside, then the constant samples $c_1$ and $c_2$ are defined by:
%%%%%%%%%figure 1%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%
%\begin{align}\nonumber
%c_1=\frac{\int_\Omega I^{2}(x,y) H(\ps\mathbf{u}(x,y))dxdy}{\int_\Omega \mathbf{u}(x,y) H(\ps\mathbf{u}(x,y))dxdy}
%\end{align}
%\begin{align}\nonumber
%c_2=\frac{\int_\Omega I^{2}(x,y) (1-H(\ps\mathbf{u}(x,y)))dxdy}{\int_\Omega \mathbf{u}(x,y) (1-H(\ps\mathbf{u}(x,y)))dxdy}
%\end{align}
%Otherwise
%\begin{align}\nonumber
%c_1=c_2=\frac{\int_\Omega I^{2}(x,y) dxdy}{\int_\Omega \mathbf{u}(x,y) dxdy}.
%\end{align}
%The following descent flow gives the minimization of energy functional:
%\begin{align}\nonumber
%\frac{\partial \psi}{\partial t}&=-\eta\frac{(\mathbf{u}(x,y)-c_1)^{2}}{{c_1}^{2}}(\ps\mathbf{u}(x,y)+1) \nonumber\\
%&-\frac{(\mathbf{u}(x,y)-c_2)^{2}}{{c_2}^{2}}(\ps\mathbf{u}(x,y)-1).
%\end{align}
Eq. (\ref{9}) is strictly convex, and it is flexible to its initial contour place, but it may not work for the noisy images as we can see in Fig. \ref{jac14} and Fig. \ref{jac15}. \color{red}Note that Fig. \ref{jac14} has been taken from the from Berkeley’s data set and is publicly available online [\textit{https://www2.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/resources.html}]. \color{black}In the theoretical aspect, the image data fitting term in a discrete sense is based on the concept of squared CoV, ${CoV}^2=\displaystyle{\sum_{i,j}\frac{(\mathbf{u}(i,j)-\mathbf{a})^2}{\mathbf{a}^2}}$ whose minimum turns out to be $\mathbf{a}= \frac{\displaystyle\sum_{i,j}{\mathbf{u}(i,j)}^2}{\displaystyle\sum_{i,j}{\mathbf{u}(i,j)}}$. The value of this average and the objective function both are sensitive to noise and outliers \cite{41} similar to \color{red}the \color{black}CV model. That is the main reason why the Wu \textit{et al.} \cite{20} model is unable to work in noisy images and performs \color{red}even the \color{black}worst than the CV model.
\subsection{Fuzzy active contour (FAC) model }
In contrast with the traditional $L_2$ norm fidelity term based models, a fuzzy active contour model with kernel metric is proposed by Wu {\it et al.} \cite{20}, which is based on the following fuzzy function:
\begin{align}\label{10}
	&F({\Gamma},c_1,c_2,z)=\mu length({\Gamma}) \nonumber\\
	&+\eta_1\int_{\Omega}[\mathbf{z}(x,y)]^{m}(1-\hat{k}(\mathbf{u}(x,y),c_1))dx dy \nonumber\\ &+\eta_2\int_{\Omega}[1-\mathbf{z}(x,y)]^{m}(1-\hat{k}(\mathbf{u}(x,y),c_2))dxdy,
\end{align}
%The implicit representation of the active contour is $0.5$ or pseudo level set of $\mathbf{z}(x,y)$, that's why the length-term is defined by:
%\begin{align}\nonumber
%\mu.lenth(\mathcal{C})&=\mu\int_{\Omega}|\nabla H(\mathbf{z}(x,y)-0.5)|dx dy \nonumber\\
%&=\mu \int_{\Omega} {\delta}(z-0.5)|\nabla (z-0.5)|dx dy,
%\end{align}
%where $\delta$, $H$ are the Dirac delta and the Heaviside function respectively. So that the energy functional given by:
%\begin{align}\label{KFAC}
%F(c_1,c_2,z)&=\mu.\int_{\Omega} {\delta}(z-0.5)|\nabla (z-0.5)|dx dy \nonumber\\ &+\eta_1\int_{\Omega}[z]^{m}(1-\hat{k}(\mathbf{u}(x,y),c_1))dx dy \nonumber\\ &+\eta_2\int_{\Omega}[1-z]^{m}(1-\hat{k}(\mathbf{u}(x,y),c_2))dx dy \nonumber
%\end{align}
%where $\eta_1,\eta_2>0$ and $\mu\geq0$ are the parameters. During the evolution $z$ must satisfy the condition $0\leq{z}\leq1$. The constants $c_1,c_2$ can be obtained as:
%\begin{eqnarray}\nonumber
%c_1=\frac{\int_{\Omega}[z]^{m}\hat{k}(\mathbf{u}(x,y),c_1)\mathbf{u}(x,y)dx dy}{\int_{\Omega}[z]^{m}\hat{k}(\mathbf{u}(x,y),c_1)dxdy}
%\end{eqnarray}
%\begin{eqnarray}\nonumber
%c_2=\frac{\int_{\Omega}[1-z]^{m}\hat{k}(\mathbf{u}(x,y),c_1)\mathbf{u}(x,y)dx dy}{\int_{\Omega}[1-z]^{m}\hat{k}(\mathbf{u}(x,y),c_1)dxdy}
%\end{eqnarray}
where the kernel metric $\hat{k}(\xi_{1},\xi_{2})=\langle{\chi(\xi_{1}),\chi(\xi_{2})}\rangle$ and
%\begin{eqnarray}\label{kernel}
%\hat{k}(\xi_{1},\xi_{2})=\langle{\chi(\xi_{1}),\chi(\xi_{2})}\rangle
%\end{eqnarray}
given  $\xi_{1}$, $\xi_{2}$  vectors and $\chi(.)$  a nonlinear map. Here $\langle{\chi(\xi_{1}),\chi(\xi_{2})}\rangle$ is the inner product of $\chi(\xi_{1})$ and $\chi(\xi_{2})$. Gaussian radial basis function $\hat{k}(\xi_{1},\xi_{2})$ is given by:
\begin{eqnarray}\label{11}
	\hat{k}(\xi_{1},\xi_{2})=\exp\left(-\frac{(\xi_{1}-\xi_{2})^{2}}{\rho}\right)
\end{eqnarray}
where  $\rho$ is the parameter.
%
%(called bandwidth parameter) of the kernel function (special case of distance covariance) and is the distance variance of all the data points of the image, and can be calculated as follow:
%\begin{eqnarray}\nonumber
%\rho=\left (\frac{1}{m-1}\sum_{i=1}^{m}(d_{i}-d)\right )^{\frac{1}{2}}
%\end{eqnarray}
%where $m$ denotes the total number of data points, $d_{i}=|I(x_{i})-\hat{I}(x_{i})|$ is the distance from the image intensity (data points) to average intensity (data center). More information to choose $\rho$ can be found in \cite{21,20}.\\
%For the minimization of Eq. (\ref{FAC}), instead of fast algorithm developed in \cite{16} Euler Lagrange equations are used. The gradient descent flow of the model is given by:
%\begin{align}\nonumber
%\frac{\partial z}{\partial t} &=-\eta_1.m.[z]^{m-1}(1-\hat{k}(\mathbf{u}(x,y),c_1)) \nonumber\\
%&+\eta_2.m.[1-z]^{m-1}(1-\hat{k}(\mathbf{u}(x,y),c_2)) \nonumber\\
%&+\mu.{\delta}(z-0.5)div \left(\frac{\nabla (z-0.5)}{|\nabla (z-0.5)|}\right). \nonumber
%\end{align}
From Figs. \ref{jac14}-\ref{jac24} in the experimental section it is clear that Eq. (\ref{10})  may not work for noisy images although $\hat{k}$ serves as a weight function which is supposed to assign suitable weights to image true pixel and outliers \cite{LiuHGWR19}. In the model implementation, the outliers get enough weights to affect the segmentation performance of this model in noisy images.

\subsection{Unconditional stable method for bimodal (USMB) image segmentation}
Li {\it et al.} \cite{34}, \color{red}proposed the following \color{black}energy functional which is based on Lee {\it et al.} \cite{39} idea \color{red}of a \color{black}stationary global minimum:
\begin{eqnarray} \label{12}
	&F(c_1,c_2,\psi)\nonumber\\
	&=\eta_1\int_{\Omega}(\mathbf{u}(x,y)-c_1)^{2}\psi(x,y)H(1+\psi(x,y))dxdy \nonumber\\
	&-\eta_2\int_{\Omega}(\mathbf{u}(x,y)-c_2)^{2}\psi(x,y)H(1-\psi(x,y))dxdy,
\end{eqnarray}
\color{red}where $H$ denotes the Heaviside function and $c_1$, $c_2$ are constants.\color{black}
% given by
%\begin{eqnarray}
%c_1=\frac{\int_{\Omega}I(x)H(\psi(x))dx}{\int_{\Omega}H(\psi(x))dx}, \nonumber
%%\end{eqnarray}
%%and
%%\begin{eqnarray}
%c_2=\frac{\int_{\Omega}I(x)(1-H(\psi(x)))dx}{\int_{\Omega}(1-H(\psi(x)))dx}. \nonumber
%\end{eqnarray}
%The parabolic equation can be obtained from Eq. (\ref{Us}) as
%\begin{align}
%&\psi_{t} =-\eta_1(I-c_1)^{2}(1+\psi) + \eta_2(I-c_2)^{2}(1-\psi)
%= -[\eta_1 \nonumber \\&(I-c_1)^{2} + \eta_2(I-c_2)^{2}]\psi-[\eta_1(I-c_1)^{2} - \eta_2(I-c_2)^{2}]. \nonumber
%\end{align}
\par
\color{red}The fact that the CV model ignores the presence of noise and other factors such as intensity in-homogeneity \cite{13}, can be easily observed from the fitting data term used in the CV model.
\begin{equation}\label{1}
	\displaystyle{\int_{\textrm{inside}(\Gamma)}|\mathbf{u}(x,y)-c_1|^{2}dxdy+\int_{\textrm{outside}(\Gamma)}|\mathbf{u}(x,y)-c_2|^{2}dxdy},
\end{equation}
where $\mathbf{u}(x,y)$ is the given image with $(x,y)\in \Omega$ a rectangular domain, $c_1$, $c_2$ are constants, and $\Gamma$ denotes the boundary of the objects. In discrete sense, this data term is based on the least square method and the objective function given by:
\begin{equation}\label{2}
	\displaystyle\frac{1}{\mid\Omega\mid}\displaystyle\sum_{i,j}(\mathbf{u}(i,j)-{\mathbf{c}})^2, 
\end{equation}
whose minimum is the sample mean $\mathbf{c}=(c_1,c_2)=\bar{\mathbf{x}}$ inside and outside the $\Gamma$. From the formula, this can be easily observed that the sample mean is largely affected by the outliers. \color{black}
\color{red}
\par
The (LBF) model Li {\it et al.} tackles intensity in-homogeneity but not noise. This phenomenon can be observed by analyzing the data fitting term:
\begin{equation}\label{3}
	\displaystyle{\int_{\Omega} K_{\sigma}\ast|\mathbf{u}(x,y)-c_1|^{2}dx dy+\int_{\Omega}K_{\sigma}\ast|\mathbf{u}(x,y)-c_2|^{2}dx dy}, 
\end{equation}
where $K_{\sigma}$ is Gaussian kernel. In a discrete and local sense this data term is also based on the least square method and the corresponding objective function:
\begin{equation}\label{4}
	\displaystyle\sum_{N_x}(\mathbf{u}(i,j)-{\textbf{c}})^2, 
\end{equation}
The above is true, particularly, in local neighborhood $N_x$ whose minimum is also the sample mean $\textbf{c}=\bar{\mathbf{x}}$. In other words, the LBF model uses the concept of the CV model but in local neighborhoods throughout image domain $\Omega$. This leads to wider image intensity variation in small patches but on the other hand, it is more prone to noise and outliers as compared to the CV model \cite{LiuHGWR19, LiuHW18, LiuHWR18}. In this way, the fitting term takes into account the image intensity variance in small patches but on the other hand, it is more prone to noise and outliers as compared to the CV model. Moreover, this model is not convex so consequently one need tune several times the position of initial guess for the desired results.\color{black}
\par
Li {\it et al.} \cite{34} showed that for any time step the proposed scheme is unconditionally stable. Moreover, with the assumption that $|\psi^{n}| \leq  1$  it is easy to show that $|\psi^{n+1}| \leq  1$, which leads to  a straightforward update of $\psi^{n+1}$ from given $\psi^{n}$. Although the method shows stability for image segmentation of synthetic and real images with moderated noise the method, similar to the above ideas were the least square fit directs to the mean of the foreground and background, shows sensitivity to high noise and outliers \cite{37}. %falls into a which .  % as in $Fig.2$.
% Bregmann Camramaan
\par
Ali \textit{et al.} \cite{43} introduced Lehmer's type generalized mean in an segmentation framework. Although this average is very effective in multi-region segmentation and suitable to different image intensity backgrounds it requires further analysis to tackle noise and outliers. Goldstein {\it et al.} \cite{40} used Bregman-split method which is well known for its speed, but may not work very well for images with intensity inhomogeneous. Furat {\it et al.} \cite{Furat} proposed Techniques for the Segmentation of Tomographic Image Data of Functional Materials by combining machine learning methods and conventional image processing steps. This approach produced good segmentation results specially for tomographic images.

\section{Proposed model}\label{proposed}
As mentioned above, most of the active contour region-based variational segmentation models consider ideal image while constructing the energy functional(s). This can be very easily observed by investigating the utilized image statistical information incorporated in objective functions, such as averages, the measures of dispersions, statistical variance, and standard deviation. \color{red}In the literature, \color{black}most of the variational region-based image segmentation models are based on the CV model fitting term idea which is sensitive to noise and outliers \cite{37}, or similar to the works in \cite{15,selec}. \color{red}Albeit, these methods are demonstrated to be very robust and effective when detecting edges and boundaries in images of low contrast; however, these methods can be very sensitive when there exists noise and outliers \cite{41}. Therefore, other methods or, at least, improvements to the classical CV model should be made in order to ensure detection of noise and outliers in low contrast images. \color{black}
%objective function  $\displaystyle\frac{1}{\mid\Omega\mid}\displaystyle\sum_{i,j}(\mathbf{u}(i,j)-{a})^2$ which is sensitive to noise and outliers \cite{37}. Few models have used  $\displaystyle\frac{1}{\mid\Omega\mid}{\displaystyle\sum_{i,j}\frac{(\mathbf{u}(i,j)-{a})^2}{a^2}}$ \cite{15,selec}, which is very robust and effective when detecting edges and boundaries of low contrast \cite{15,selec} but very sensitive and vigilant when there is noise and outliers \cite{41}.
\par
To improve the \color{red}state-of-the-art \color{black}models mentioned in related works, \color{red}we propose a new method which incorporates the power mean into the robust discrete objective function by replacing the traditional models where the arithmetic mean has been used. The sate-of-the-art with in the domain of the power mean indicates that it has the capability to discriminate the noise and intrinsic intensity \cite{37}. To handle a noisy image one can design a formulation in the continuous framework based on averages and measure of dispersion's. Furthermore, the employment of a fuzzy membership function has its advantages over the traditional level set function, \color{red}as this allows, the involvement \color{black}of less number of functions to capture many objects of different intensities \cite{4,313}.\color{black}
 %, as natural images may involve many objects of different intensities \cite{4,313}. 
\par
Initially, we discuss the power mean function and its property of canceling the negative effect of outliers. We continue in the second subsection with the presentation of the proposed model guided by a fuzzy function based formulation. The rest of the section analyses the convexity of the energy functional, its semi-continuity and coercivity. 
\subsection{Power mean}%{Principal component}
\par{\textbf{Definition.}}
For a given gray scale image $\mathbf{u}(x,y)\in \Omega$ of size $N\times M$, power mean can be defined in discrete form  as follows \cite{37,43}:
\begin{align} \label{13}
	M_{p}(\textbf{I})=\left(\frac{1}{NM}\sum_{i=1,j=1}^{N,M}{\mathbf{u}(i,j)}^{p}\right)^{\frac{1}{p}},
\end{align}
where $p\neq0$, and $\mathbf{u}(i,j)>0$ %$\textbf{x}=\mathbf{u}(i,j)>0$
is the intensity value at a certain pixel $(i,j)$. For different value of $p$, such as $p=1,0,-1$, the general mean represents specific mean variations such as arithmetic, geometric or  harmonic mean. \\
\color{red}The parameter $p$ controls the contribution of each sample's element by handling each of them differently according to their significance. \color{black}Oh \textit{et al.} \cite{37} and Ali \textit{et al.} \cite{43} has been introduced an implementation of such feature. \color{red}The authors have \color{black}expressed the general power mean as a linear combination of the elements in the set and its simplification form as:
%In those works, 
%the general mean of a set of positive numbers is expressed by a non-negative linear combination of the elements in the set and its simplification as in \cite{37,43}:
\begin{eqnarray}\label{14}
	\sum_{i=1,j=1}^{N,M}{\mathbf{u}(i,j)}^{p}=\sum_{i=1,j=1}^{N,M}y(i,j)\mathbf{u}(i,j),  \\
	y(i,j)=\mathbf{u}(i,j)^{p-1}, \textrm{ for } i=1,2,...,N\ \textrm{and}\ j=1,2,...,M.  \nonumber
\end{eqnarray}
The employment of the generalized mean controls \color{red}the existing trade-off \color{black}between the negativeness of outliers in the observed set. It is easy to observe that, the generalized mean in Eq. (\ref{14}) is an arithmetic mean if $p = 1$. The weight $y(i,j)$ decreases (increases) as $\mathbf{u}(i,j)$ increases (decreases) if $p\leq 1$. This indicates that Eq. (\ref{14}) is more affected by the small intensity values in the given image $\{\mathbf{u}(i,j)\}_{i=1,j=1}^{N, M}$ and if $p$ decreases, the extent of the \color{red}effeteness \color{black}increases. In %a generalized principal component analysis framework
\cite{21}, this information played a key role in applying the generalized mean approach. To develop the ancient models, %of principle component analysis, 
Oh \textit{et al.} \cite{37} exchange the conventional least square sample mean fitting term with the generalized mean fitting term as below: %Based on these observations, we can write and proceed as 
\begin{eqnarray}\label{15}
	m_G(\textbf{I})={\arg\min}_\mathbf{a}\left(\frac{1}{NM}\sum_{i=1,j=1}^{N,M}\left(||\mathbf{u}(i,j)-\mathbf{a}||_{2}^{2}\right)^{p}\right)^{\frac{1}{p}},
\end{eqnarray}
where $\mathbf{a}$ is any arbitrary value in a given image intensity values. \color{red}We observe that Eq. (\ref{15}) converts to the traditional CV objective function for $p=1$ which is based on the conventional arithmetic mean of the squared distance \cite{37,43}. One can choose $p < 1$ \cite{37,43} to reduce the negative effects of outliers. In such a way, as $p$ decreases the contribution of a large number to the objective function decreases. This means that the power mean can discriminate the noise and intrinsic intensity. \color{black}Furthermore, Eq. (\ref{15}) can be written as \cite{37,43}:
\begin{eqnarray}\label{16}
	m_G(\textbf{I})={\arg \min}_\mathbf{a}\sum_{i=1,j=1}^{N,M}\left(||\mathbf{u}(i,j)-\mathbf{a}||_{2}^{2}\right)^{p}.
\end{eqnarray}
The basic condition for the generalized sample mean $m_G$ to be a local minimum of the objective function (\ref{15}) is that the gradient of this function with respect to $\mathbf{a}$ is equal to zero \cite{37,43}, \color{red}that is mathematically described as given by Eq. (\ref{17}).\color{black}
\begin{eqnarray}\label{17}
	\frac{\partial}{\partial \mathbf{a}}\Big(\sum_{i=1,j=1}^{N,M}\left(||\mathbf{u}(i,j)-\mathbf{a}||_{2}^{2}\right)^{p}\Big)=0.
\end{eqnarray}
Similar to the expectation-maximization algorithm scheme, Oh \textit{et al.} \cite{37} develop an iterative form for easily \color{red}solving \color{black}Eq. (\ref{16}). First, rewriting Eq. (\ref{16}) in the form of Eq. (\ref{14}) and then approximated by a quadratic function $$||\mathbf{u}(i,j)-\mathbf{a}||_{2}^{2}$$ which
can be optimized as:
\begin{eqnarray}\label{18}
	\sum_{i=1,j=1}^{N,M}\left(||\mathbf{u}(i,j)-\mathbf{a}||_{2}^{2}\right)^{p}\approx\sum_{i=1,j=1}^{N,M}{\beta(i,j)}^{(k)}||\mathbf{u}(i,j)-\mathbf{a}||_{2}^{2} 
\end{eqnarray}
where
\begin{eqnarray}\label{19}
	{\beta^{(k)}(i,j)}=\left(||\mathbf{u}(i,j)-\mathbf{a}^{(k)}||_{2}^{2}\right)^{p-1},
\end{eqnarray}
for $k$ number of the iterations. The approximation is exact when $\mathbf{a} = \mathbf{a}^{(k)}$. 
%Without computing the Hessian matrix of the objective function nearby
Here, $\mathbf{a}^{(k)}$ can be updated based on the computed ${\beta(i,j)}$  in Eq. (\ref{18}). The approximated function based on computed  ${\beta(i,j)}$ is:
\begin{eqnarray}\label{20}
	\frac{\partial}{\partial \mathbf{a}}{\beta^{(k)}(i,j)}||\mathbf{u}(i,j)-\mathbf{a}||_{2}^{2}=0.
\end{eqnarray}
Then, as a weighted average of the samples $\mathbf{a}^{(k+1)}$ can be computed by (\ref{19}) and it gives:
\begin{eqnarray}\label{21}
	\mathbf{a}^{(k+1)}=\frac{1}{\sum_{i=1,j=1}^{N,M}{\beta^{(k)}(i,j)}}\Big(\sum_{i=1,j=1}^{N,M}{\beta^{(k)}(i,j)}\mathbf{u}(i,j)\Big).
\end{eqnarray}
It is important to point out that the function $\beta$ serves as a weight function that assigns suitable weights to the true image pixels and outliers \cite{37,43}. The parameter $p$ controls the function $\beta$ and its optimal tuning value has been shown to be in the range $0.6\leq p\leq 0.8$ \cite{37,43}. In the following, we show a new implementation of generalized mean in fuzzy membership variational segmentation framework, which has been fully studied in the work of Oh \textit{et al.} \cite{37}.  %and H. Ali et al. \cite{43}

\begin{figure*}
	\centering
	\hspace*{-1cm}
	\begin{tabular}{cc}
		%\underline{Contour evaluation for different size of initialization:}
		\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{Initial1.eps}}
		\subfigure{\label{wu1}\includegraphics[width=2cm,height=2cm]{Initial2.eps}}
		\subfigure{\label{wu1}\includegraphics[width=2cm,height=2cm]{Initial4.eps}}
		\subfigure{\label{wu2}\includegraphics[width=2cm,height=2cm]{Initial3.eps}}
		\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{holeinitial1.eps}}
		\subfigure{\label{wu1}\includegraphics[width=2cm,height=2cm]{holeinitial2.eps}}
		\subfigure{\label{wu1}\includegraphics[width=2cm,height=2cm]{holeinitial3.eps}}
		\subfigure{\label{wu2}\includegraphics[width=2cm,height=2cm]{holeinitial4.eps}}\\
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{InitialResult1.eps}}
		\subfigure{\label{Li2}\includegraphics[width=2cm,height=2cm]{InitialResult2.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{Initial4result.eps}}
		\subfigure{\label{gspf1}\includegraphics[width=2cm,height=2cm]{InitialResult3.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{holefinalc1.eps}}
		\subfigure{\label{Li2}\includegraphics[width=2cm,height=2cm]{holefinalc2.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{holefinalc3.eps}}
		\subfigure{\label{gspf1}\includegraphics[width=2cm,height=2cm]{holefinalc4.eps}}\\
		\subfigure{\label{gspf2}\includegraphics[width=2cm,height=2cm]{InitialSeg1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=2cm,height=2cm]{InitialSeg2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=2cm,height=2cm]{Initial4seg.eps}}
		\subfigure{\label{L3S}\includegraphics[width=2cm,height=2cm]{InitialSeg3.eps}}
		\subfigure{\label{gspf2}\includegraphics[width=2cm,height=2cm]{hole1seg.eps}}
		\subfigure{\label{L3S}\includegraphics[width=2cm,height=2cm]{hole2seg.eps}}
		\subfigure{\label{L3S}\includegraphics[width=2cm,height=2cm]{hole3seg.eps}}
		\subfigure{\label{L3S}\includegraphics[width=2cm,height=2cm]{hole4seg.eps}}
	\end{tabular}
	\caption{Illustration of the proposed segmentation method for the first ever black hole image and noisy images taken from \cite{38,313} papers. First row: different initial contours. Second row: final contours. Third row: segmented results with $p=0.5$, $\mu$=0.7, and $\sigma$=3 .}\label{exp1}\label{hole}
\end{figure*}
%for p

\begin{figure*}
	\centering
	\hspace*{-1cm}
	\color{red}
	\hspace*{0.3cm}
	\begin{tabular}{cc}
		%\underline{Contour evaluation for different size of initialization:}
		%\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{p=0.eps}}
		\subfigure[p=0.1]{\label{wu1}\includegraphics[width=3cm,height=3cm]{pointone.eps}}
		\subfigure[p=0.2]{\label{wu1}\includegraphics[width=3cm,height=3cm]{pointtwo.eps}}
		\subfigure[p=0.3]{\label{wu2}\includegraphics[width=3cm,height=3cm]{pointthree.eps}}
		\subfigure[p=0.4]{\label{circle}\includegraphics[width=3cm,height=3cm]{p4.eps}}
		\subfigure[p=0.5]{\label{wu1}\includegraphics[width=3cm,height=3cm]{InitialResult1.eps}}\\
		\subfigure[p=0.6]{\label{wu1}\includegraphics[width=3cm,height=3cm]{p6.eps}}
		\subfigure[p=0.7]{\label{wu2}\includegraphics[width=3cm,height=3cm]{p7.eps}}
		\subfigure[p=0.8]{\label{Li1}\includegraphics[width=3cm,height=3cm]{p8.eps}}
		\subfigure[p=0.9]{\label{Li2}\includegraphics[width=3cm,height=3cm]{p9.eps}}
		\subfigure[p=0.10]{\label{Li1}\includegraphics[width=3cm,height=3cm]{p10.eps}}
	\color{black}
	\end{tabular}
	\caption{\color{red} The performance of the proposed method for noisy image taken from \cite{38,313} papers for different values of $p$, and fixed $\mu$=0.7, $\sigma$=3, iteration$=50$.\color{black}} \label{pvalue}
\end{figure*}

%%%%% 1st comparison ==Circle fig..3
\begin{figure*}
	\centering
	\hspace*{-0.4cm}
	\begin{tabular}{cc}
		%\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{Img1.eps}}
		\subfigure{\label{circle}\includegraphics[width=3cm,height=2cm]{Image1initial.eps}}
		\subfigure{\label{circle}\includegraphics[width=3cm,height=2cm]{Image1initial.eps}}
		\subfigure{\label{circle}\includegraphics[width=3cm,height=2cm]{Image1initial.eps}}
		\subfigure{\label{circle}\includegraphics[width=3cm,height=2cm]{Image1initial.eps}}
		\subfigure{\label{circle}\includegraphics[width=3cm,height=2cm]{Image1initial.eps}}\\
		\subfigure{\label{wu1}\includegraphics[width=3cm,height=2cm]{Img1wu1.eps}}
		\subfigure{\label{Li1}\includegraphics[width=3cm,height=2cm]{Img1li1.eps}}
		\subfigure{\label{gspf1}\includegraphics[width=3cm,height=2cm]{Img1gspf1.eps}}
		\subfigure{\label{gspf1}\includegraphics[width=3cm,height=2cm]{Img1krin1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img1our1.eps}}\\	
		\subfigure{\label{wu2}\includegraphics[width=3cm,height=2cm]{Img1wu2.eps}}
		\subfigure{\label{Li2}\includegraphics[width=3cm,height=2cm]{Img1li2.eps}}
		\subfigure{\label{gspf2}\includegraphics[width=3cm,height=2cm]{Img1gspf2.eps}}
		\subfigure{\label{gspf2}\includegraphics[width=3cm,height=2cm]{Img1krin2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img1our2.eps}}
	\end{tabular}
	\caption{This image is taken from the Berkeley's data set$^1$. First row shows the given image and the zero level set initialization, the second row shows the segmentation contour and the last row \color{red}denotes \color{black}the binary image resulting from the reconstruction of each method involved in the comparison. First, second, third and fourth \color{red}columns \color{black}are the segmentation results of Wu {\it et al.} \cite{15}, Li  {\it et al.} \cite{34}, Wu {\it et al.} \cite{20} and Krinidis {\it et al.} \cite{313}, respectively. \color{red}The fifth column \color{black}illustrates the result of proposed model with $p=0.5$, $\mu$=0.7, $\sigma$=3, with speckle noise = 0.2.}\label{jac14}
\end{figure*}
%\FloatBarrier
%
%%second comparison airplane fig..4
\par
\begin{figure*}
	\centering
	\hspace*{-0.4cm}
	\begin{tabular}{cc}
		%\subfigure{\label{LBF1}\includegraphics[width=2cm,height=2cm]{Img2.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Image2initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Image2initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Image2initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Image2initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Image2initial.eps}}
		\\
		\subfigure{\label{LBF2}\includegraphics[width=3cm,height=2cm]{Img2wu1.eps}}
		\subfigure{\label{LBF1S}\includegraphics[width=3cm,height=2cm]{Img2li1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img2gspf1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img2krin1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img2our1.eps}}
		\\
		\subfigure{\label{L3}\includegraphics[width=3cm,height=2cm]{Img2wu2.eps}}
		\subfigure{\label{LS}\includegraphics[width=3cm,height=2cm]{Img2li2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img2gspf2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img2krin2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img2our2.eps}}
	\end{tabular}
	\caption{First row shows the given image and the zero level set initialization, the second row shows the segmentation contour and the last row the binary image resulting from the reconstruction of each method involved in the comparison. First, second, third and fourth \color{red}columns \color{black}are the segmentation results of Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34}, Wu {\it et al.} \cite{20} and Krinidis {\it et al.} \cite{313}, respectively. \color{red}The Fifth column \color{black}illustrates the result of proposed model with $p=0.5$, $\mu$=0.7, $\sigma$=3, with speckle noise = 0.2.} \label{jac15}
\end{figure*}
%\FloatBarrier
%%%%%%%%%%%% 3rd comparison Foot fig.5
\begin{figure*}
	\centering
	\hspace*{-0.4cm}
	\begin{tabular}{cc}
		%\subfigure{\label{LBF1}\includegraphics[width=2cm,height=2cm]{Foot.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Foot3initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Foot3initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Foot3initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Foot3initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{Foot3initial.eps}}
		\\
		\subfigure{\label{LBF2}\includegraphics[width=3cm,height=2cm]{WuheresultF.eps}}
		\subfigure{\label{LBF1S}\includegraphics[width=3cm,height=2cm]{LikimresultF.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{GspfresultF.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{KrinresultF.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{oursegF.eps}}
		\\
		\subfigure{\label{L3}\includegraphics[width=3cm,height=2cm]{WuhesegF.eps}}
		\subfigure{\label{LS}\includegraphics[width=3cm,height=2cm]{LikimsegF.eps}}	\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{GspfsegF.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{KrinsegF.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{OurresultF.eps}}
	\end{tabular}
	\caption{First, second, third and fourth \color{red}columns \color{black}are the segmentation results of Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34}, Wu {\it et al.} \cite{20} and Krinidis {\it et al.} \cite{313}, respectively. \color{red}The \color{black}fifth column is the segmentation result of proposed model with p=0.6, $\mu$=0.7, $\sigma$=3, and noise = 0.1.} \label{foot}
\end{figure*}

\subsection{A new fuzzy function segmentation model led by data-guided outliers avoidance}
Defining the image $\mathbf{u}$ on $\Omega \subset \Re^{2}$, and $\Omega_i \subseteq \Omega$ are disjoint connected open subsets with a piecewise smooth boundary $C$ ($\cup_{i}\Omega_i$). $\{C_i\in \Re^2|^n_{i=1}\}$ are the curves of the samples to be segmented and $\{c_i\in \Re^2|^n_{i=1}\}$ are their homogeneous associated means. The task of image segmentation is to divide an image into $n$ group of data samples $\{C_i|^n_{i=1}\}$. To improve the segmentation accuracy in the presence of outliers we investigate a novel approach where the quality of generalized sample mean is taken into account and incorporated as a fitting term in a minimization functional. Concretely, we use the Euclidean distance of an input sample $C_i$ to representative samples $c_i$  by allowing in this way some pixels of $C_i$ to be recognized as outliers. In this case, the  $c_i$ value not necessary must be near to these samples which consequently brings robustness to outliers. Based on this simple idea we \color{red}can \color{black}modify the Chan-Vese active contour model, \color{red}and described it \color{black}in the following form:
\begin{eqnarray}\label{22}
	&&F(c_{1},c_{2},C)=\mu length(C)\\
	&+&\Big(\frac{1}{N_1(C)}\int_{inside(C)}(\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}dxdy)^{p}\Big)^\frac{1}{p}\nonumber\\
	&+&\Big(\frac{1}{N_2(C)}\int_{outside(C)}(\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}dxdy)^{p}\Big)^\frac{1}{p},\nonumber
\end{eqnarray}
here we have three terms, first on is the curve length term, the second and third terms, which we will further refer as $F_1(C)$ and $ F_2(C)$, are the new %generalized segmentation
fitting-terms with $N_1(C)$ and $N_2(C)$ the number of points inside and outside the curve $C$, respectively.
By the same argument as we get Eq. (\ref{16}):
\begin{eqnarray} \label{23,24}
	F_{1}(C)=\int_{inside(C)}(\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}dxdy)^{p},\\
	F_{2}(C)=\int_{outside(C)}(\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}dxdy)^{p},
\end{eqnarray}
and by the same way to Eq. (\ref{18}), we get:
\begin{eqnarray}\label{25,26}
	F_{1}(C) \approx \int_{inside(C)}\alpha(x,y)\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}dxdy, \\
	F_{2}(C) \approx \int_{outside(C)}\beta(x,y)\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}dxdy, 
\end{eqnarray}
where
\begin{eqnarray} \label{27}
	\alpha(x,y)=(\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2})^{p-1},
\end{eqnarray}
and
\begin{eqnarray} \label{28}
	\beta(x,y)=(\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2})^{p-1}. 
\end{eqnarray}
Incorporating the fuzzy membership function $\mathbf{z}(x,y)$, Eq. (\ref{16}) can be rewritten as: %the above terms as: %each term of Eq. (\ref{eq5}) in the following form: %and the level set formulation of our model we replace unknown variable C by the unknown variable $\phi$ thus we have :
\begin{eqnarray} \label{29}
	\int_{z>0.5}&\alpha(x,y)&\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}dxdy\\
	&=& \int_{\Omega}\alpha(x,y){\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}}[\mathbf{z}(x,y)]^{m}dxdy,\nonumber
\end{eqnarray}
and
\begin{eqnarray} \label{30}
	\int_{z<0.5}&&\beta(x,y)\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}dxdy\\\
	&=&\int_{\Omega}\beta(x,y){\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}}[1-\mathbf{z}(x,y)]^{m}dxdy,\nonumber
\end{eqnarray}
%Therefor, the fitting energy functional is written as:
%\begin{eqnarray}
%F_{ft}(c_{1},c_{2},\alpha,\beta,\phi)&=&\int_{\Omega}\alpha(x,y){\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}}(H(\Ph\mathbf{u}(x,y)))dxdy\nonumber\\
%&+&\int_{\Omega}\beta(x,y){\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}}(1-H(\Ph\mathbf{u}(x,y)))dxdy\label{fit1}
%\end{eqnarray}
where $\alpha$, $\beta$ are updated through $c_{1}$ and $c_{2}$ in each step using Eq. (\ref{27}), (\ref{28}) and $\mathbf{z}$ is the fuzzy membership function.
%As we know the influence of noise depend on the value of p , where p is the parameter so it must satisfy $0<p<1$ ,now
%In order to controls size of jumps and the geometry of boundaries in an image segmentation process it has been shown that the alternation between a smoothing and a fitting term is beneficial \cite{cv}. To measure the variation of the curve with penalizing discontinuities Chan-Vese \cite{cv} hired the length of the curve into the minimization functional leading to a Total Variation (TV) term. In similar way we use this term as a smoothing term in our model.
We propose the following minimization functional:
\begin{align}\label{31}
	&F(\mathbf{z},c_{1},c_{2})= \mu\int_{\Omega}|\nabla \mathbf{z}(x,y)|dxdy \nonumber\\ &+\int_{\Omega}\alpha(x,y){\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}}[\mathbf{z}(x,y)]^{m}dxdy\nonumber\\
	&+\int_{\Omega}\beta(x,y){\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}}[1-\mathbf{z}(x,y)]^{m}dxdy.
\end{align}
%where $p\in\Re$.
%The Euler-Lagrange equation acts locally on all level curves around ${\ph\mathbf{u}(x,y)=0}$ when we use $H_{\mu}$ and $\delta_{\mu}$, so we can obtain a global minimizer, independently of the position of the initial curve; moreover, this allows to automatically detect interior contours. To remain close to the initial minimization problem we prefer to work with $\delta(\phi)$ and the problem of extending the evolution to all level sets of $\phi$ was solved here using the approximation $\delta_{\mu}$ of $\delta$.

%In order to compute the associated Euler-Lagrange equation for $z$ in equation (\ref{eq13}) we consider the regularized version of the Heaviside and Delta functions as $H_{\mu}$ and $\delta_{\mu}$, respectively. In this paper, we use the following regularization for the heaviside function  \cite{VLS}:
%\begin{equation}
% H_{\mu}(z)=\frac{1}{2}(1+\frac{2}{\Pi}\arctan(\frac{z}{\mu}))
% \end{equation} % $\mu\rightarrow 0$ where $H_{\mu}$ is the regularization of H and
% and its corresponding Delta function $\delta_{\mu}=H_{\mu}^{'}$.
%
%Then, equation (\ref{eq13}) is rewritten as:
%\begin{eqnarray}
%F_{\mu}(\phi,c_{1},c_{2})&=&\mu\int_{\Omega}\delta_\mu(\ph\mathbf{u}(x,y))|\nabla\ph\mathbf{u}(x,y)\|dxdy+\int_{\Omega}\alpha{\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2}}(H_{\mu}(\ph\mathbf{u}(x,y)))dxdy\nonumber\\
%&+&\int_{\Omega}\beta{\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}}(1-H_{\mu}(\ph\mathbf{u}(x,y)))dxdy.\label{eq14}
%\end{eqnarray}
Keeping $c_{1}$ and $c_{2}$ fixed in Eq. (\ref{31}), then minimizing $F(\mathbf{z},c_{1},c_{2})$ with respect to $\mathbf{z}$, we get the associated Euler-Lagrange equation for $\mathbf{z}$, $t$ is an artificial time parameterizing the descent direction as:
\begin{eqnarray} \label{32}
	\frac{\partial \mathbf{z}}{\partial t}&=&\mu\nabla\Big({\nabla \mathbf{z} \over |\nabla \mathbf{z}|}\Big) \nonumber \\
	&-& m\alpha(x,y) [\mathbf{z}(x,y)]^{m-1}\|\mathbf{u}(x,y)-c_{1}\|_{2}^{2} \nonumber \\
	&+&m\beta(x,y) [1-\mathbf{z}(x,y)]^{m-1}\|\mathbf{u}(x,y)-c_{2}\|_{2}^{2}=0 \nonumber \\ &&   \text{in} \hspace{0.5cm} (0,\infty)\times\Omega, 
\end{eqnarray}
with
\begin{eqnarray} \label{33}
	\mathbf{z}(0,x,y)=\mathbf{z}(x,y) \nonumber \hspace{0.5cm}\text{in} \hspace{0.5cm} \Omega \\
	\frac{\mathbf{z}}{|\nabla \mathbf{z}|}\frac{\partial \mathbf{z}}{\partial\overrightarrow{n}} \hspace{0.5cm} \text{on} \hspace{0.5cm} \partial\Omega,
\end{eqnarray}
where $\overrightarrow{n}$ is the normal to the boundary $\partial\Omega$ in exterior, $\alpha(x,y)$ is defined in Eq. (\ref{27}), $\beta(x,y)$ is defined in Eq. (\ref{28}) and $\frac{\partial \mathbf{z}}{\partial\overrightarrow{n}}$ is the normal derivative of $\mathbf{z}$ at $\partial\Omega$.

%In order to prevent the level set function to become too steep or too flat, different level set regularization techniques are applied forcing the level set to be a signed distance function (SDF). Re-scaling, re-initialize or fitting constraints are applied in order to keep $\phi$ close to signed distance function \cite{CSI, VLS}. %it act as a regularization or rescaling.
%As re-initialization prevents interior contour from growing so it should not be too strong. Only for a few numerical results we have applied the reinitialization, solving the following evolution Eq. \cite{CSI} .
%\begin{equation}
%\{\Psi_{r}=sign(\Phi(t))(1-|\nabla|)\\
%\\ \Psi(0,.)=\Phi(t,.)
%\end{equation}
%The solution of $\phi$ at time t is $\Phi(t,.)$  ,then new $\phi(t,.)$ will be $\Psi$ ,such that $\Psi$ is obtained at the steady state of (29).The solution $\Psi(t,.)$ of (29) will have the same zero level set as $\phi(t,.)$ and away from this
%set.$|\nabla\Psi|$ will converges to 1.
%Such techniques are the curvature-based term $div(\nabla\phi/|\nabla\phi|)$.
%Given $z$ is an SDF satisfying $|\nabla z| = 1$, the regularized term can be rewritten
%as $\triangle z$, which is the Laplacian of the fuzzy membership function $z$. As pointed
%out in \cite{Shi:05} and based on the theory of scale-space \cite{Perona:90}, the evolution
%of a function with its Laplacian is equivalent to a Gaussian kernel
%filtering the initial condition of the function \cite{ZHANG2010668}. Thus one can use a
%Gaussian filtering process to further regularize the level set function.
%The standard deviation of the Gaussian filter can control the regularization strength, just as the parameter $\mu$.
%In our work we utilize a Gaussian filter to smooth the level set function
%in order to guaranty a regular  interface by discarding the regularization  term $div \frac{\nabla z}{|\nabla z|}$ and update $z$ using the following formula: 
%\begin{eqnarray}
%\frac{\partial z}{\partial t}&=\mathbf{z}(x,y)\Big[-\alpha\|\mathbf{u}(x,y)-c_{1}(\mathbf{z}(x,y))\|_{2}^{2}\nonumber\\ &+\beta\|\mathbf{u}(x,y)-c_{2}(\mathbf{z}(x,y))\|_{2}^{2}\Big].\label{eq15}
%\end{eqnarray}
It is important to note that $c_1$ and $c_2$ are updated through $\alpha$ and $\beta$ in each step using Eq. (\ref{21}).
\begin{equation} \label{34}
	c^{}_{1}=\frac{1}{\int_{\Omega}\alpha^{}(x,y)[\mathbf{z}(x,y)]^{m}}\int_{\Omega}\alpha^{}(x,y)\mathbf{u}(x,y)[\mathbf{z}(x,y)]^{m},
\end{equation}
\begin{equation} \label{35}
	c^{}_{2}=\frac{1}{\int_{\Omega}\beta^{}(x,y)[1-\mathbf{z}(x,y)]^{m}}\int_{\Omega}\beta^{}(x,y)\mathbf{u}(x,y)[1-\mathbf{z}(x,y)]^{m}.
\end{equation}
Keeping $c_{1}$, $c_{2}$ fixed and $\mu=0$, then minimizing the energy functional (\ref{30}) with respect to the fuzzy membership function $\mathbf{z}$, as in \cite{313} we get:
\begin{equation}\label{36}
	z=\frac{1}{1+\left(\frac{\alpha(x,y)||\mathbf{u}(x,y)-c_1||_{2}^{2}}{\beta(x,y)||\mathbf{u}(x,y)-c_2||_{2}^{2}}\right)^{\frac{1}{m-1}}}.
\end{equation}
Moreover, this updated  value is used in the numerical explicit solution of the following Euler Lagrange's equation:
\begin{eqnarray}\label{37}
	\frac{\partial \mathbf{z}(x,y)}{\partial t}&=&\mu \nabla \frac{\nabla \mathbf{z}(x,y)}{|\nabla \mathbf{z}(x,y)|}\nonumber\\
	&+& m[\mathbf{z}(x,y)]^{m-1}||\mathbf{u}(x,y)-c_1||_{2}^{2}\nonumber\\
	&+& m[1-\mathbf{z}(x,y)]^{m-1}||\mathbf{u}(x,y)-c_2||_{2}^{2}.
\end{eqnarray}
Introduction a time step $\Delta t $ the above equation can be solved with the time marching method as following:
\begin{eqnarray}\label{38}
	z^{k+1}(x,y)&=&z^k(x,y)+\Delta t\Big[\mu \nabla \frac{\nabla \mathbf{z}(x,y)}{|\nabla \mathbf{z}(x,y)|}\nonumber\\
	&+& m[\mathbf{z}(x,y)]^{m-1}||\mathbf{u}(x,y)-c_1||_{2}^{2}\nonumber\\
	&+& m[1-\mathbf{z}(x,y)]^{m-1}||\mathbf{u}(x,y)-c_2||_{2}^{2}\Big].
\end{eqnarray}

%The principal steps of the proposed model are shown in an algorithm form, as given by Algorithm \ref{}:\\
%\begin{algorithm}
%	\caption{ 
%		$ z\leftarrow PMBIS(\mathbf{u},\mu,\Delta t, \sigma)$}
%	\begin{algorithmic}
%		\STATE  Step 0. Choose $\mu$,  $\Delta t$, $\mathbf{a}$ and $\sigma$; 
%		\STATE Step 1. Initialize the function $\mathbf{z}(x,y)$ as  binary function, and $\alpha=\beta=I$.\\
%		\STATE    Step 2.
%		\FOR {$k=0: maxit$} 
%		\hspace{2.cm} $\bullet$ Compute $\alpha$ using Eq. (\ref{alph1}),\newline
%		\hspace{2.cm}$\bullet$  Compute $\beta$ using Eq. (\ref{alph2}), \newline
%		\hspace{2.cm}$\bullet$  Compute  $c_1$ and $c_2$ using Eqs.  (\ref{eq16}) and (\ref{eq17}),\newline
%		\hspace{2.cm}$\bullet$  Update  $\mathbf{z}(x,y)$ from Eq. (\ref{fuzzy}),\newline
%		\hspace{2.cm}$\bullet$ Evolve $\mathbf{z}(x,y)$ according to Eq. (\ref{EL}).\newline
%		\hspace{2.cm}$\bullet$   Smooth the function $\mathbf{z}(x,y)$ by Gaussian kernel.\newline
%		\hspace{2.cm}$\bullet$   If the stationary evolution criterion is satisfied, the algorithm stops. Otherwise go to $k=k+1$.%\newline
%				\STATE \ENDFOR
%	\end{algorithmic}
%	\label{Afzalalg}
%\end{algorithm}
\par
In the following section, we explore some mathematical properties \color{red}that are related \color{black}to the convexity of the proposed \color{red}functional measurements, as determined in Eq. (\ref{30}), \color{black}which are important to obtain the global minimum. 
\color{red}
\subsection{Convexity, Lower semi-continuity and coercivity of the energy functional} 
\par
{\textbf{Theorem 1.}} The energy functional (\ref{31}) is convex. The mathematical discussion over the proof of this theorem can be found in the Appendix.
\par
{Proof:} (\ref{Appendix}).\\
{\textbf{Theorem 2.}}
For the energy functional Eq. (\ref{31}) and for fixed $\alpha, \beta, c_{1}, c_{2}$, there exists at least one solution $\mathbf{z}^{\ast}$ in the admissible set $\Lambda=\{\mathbf{z}: \mathbf{z}\in BV(\Omega), 0\leq \mathbf{z} \leq 1\}$. The mathematical discussion over the proof of this theorem can be found in the Appendix.
\par{Proof:} (\ref{Theorem 2}).
\color{black}
%%%%%%%%%%%%%%%  4th comparison child fig.6
\begin{figure*}[]
	\centering
	\hspace*{-0.4cm}
	\begin{tabular}{cc}
		%\subfigure[]{\label{LBF1}\includegraphics[width=2cm,height=2cm]{48.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{child4initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{child4initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{child4initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{child4initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{child4initial.eps}}
		\\
		\subfigure{\label{LBF2}\includegraphics[width=3cm,height=2cm]{wuhe48.eps}}
		\subfigure{\label{LBF1S}\includegraphics[width=3cm,height=2cm]{Likim48.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Gspf48.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{krin48.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{our48.eps}}
		\\
		\subfigure{\label{L3}\includegraphics[width=3cm,height=2cm]{wuhe48seg.eps}}
		\subfigure{\label{LS}\includegraphics[width=3cm,height=2cm]{Likim48seg.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Gspf48seg.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{krinidis48seg.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{our48seg.eps}}
	\end{tabular}
	\caption{First, second, third and fourth \color{red}columns \color{black}are the segmentation results of Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34}, Wu {\it et al.} \cite{20} and Krinidis {\it et al.} \cite{313}, respectively. \color{red}The fifth column \color{black}is the segmentation result of proposed model with $p=0.6$, $\mu$=0.7, and $\sigma$=3.} \label{child}
\end{figure*}
\begin{figure*}[]
	\centering
	\hspace*{-0.4cm}
	\begin{tabular}{cc}
		%\subfigure{\label{LBF1}\includegraphics[width=2cm,height=2cm]{Img5.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{sm5initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{sm5initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{sm5initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{sm5initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{sm5initial.eps}}
		\\
		\subfigure{\label{LBF2}\includegraphics[width=3cm,height=2cm]{Img5wu1.eps}}
		\subfigure{\label{LBF1S}\includegraphics[width=3cm,height=2cm]{Img5li1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img5gspf1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img5krin1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img5our1.eps}}
		\\
		\subfigure{\label{L3}\includegraphics[width=3cm,height=2cm]{Img5wu2.eps}}
		\subfigure{\label{LS}\includegraphics[width=3cm,height=2cm]{Img5li2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img5gspf2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img5krin2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img5our2.eps}}
	\end{tabular}
	\caption{First, second, third and fourth \color{red}columns \color{black}are the segmentation results of Wu {\it et al.} \cite{15}, Li-kim \cite{34}, Wu {\it et al.} \cite{20} and Krinidis {\it et al.} \cite{313}, respectively. \color{red}The fifth column \color{black}is the segmentation result of proposed model with $p=0.5$, $\mu$=0.7, and $\sigma$=3.}\label{moon and star}
\end{figure*}
\begin{figure*}[]
	\centering
	\hspace*{-0.4cm}
	\begin{tabular}{cc}
		%\subfigure[]{\label{LBF1}\includegraphics[width=2cm,height=2cm]{Img12.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{t12initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{t12initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{t12initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{t12initial.eps}}
		\subfigure{\label{LBF1}\includegraphics[width=3cm,height=2cm]{t12initial.eps}}
		\\
		\subfigure{\label{LBF2}\includegraphics[width=3cm,height=2cm]{Img12wu1.eps}}
		\subfigure{\label{LBF1S}\includegraphics[width=3cm,height=2cm]{Img12li1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img12gspf1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img12krin1.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img12our1.eps}}
		\\
		\subfigure{\label{L3}\includegraphics[width=3cm,height=2cm]{Img12wu2.eps}}
		\subfigure{\label{LS}\includegraphics[width=3cm,height=2cm]{Img12li2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img12gspf2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img12krin2.eps}}
		\subfigure{\label{L3S}\includegraphics[width=3cm,height=2cm]{Img12our2.eps}}
	\end{tabular}
	\caption{First, second, third and fourth \color{red}columns \color{black}are the segmentation results of Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34}, Wu {\it et al.} \cite{20} and Krinidis {\it et al.} \cite{313}, respectively. \color{red}The fifth column \color{black}is the segmentation result of proposed model with $p=0.5$,  $\mu$=0.7, $\sigma$=3, and speckle noise = 0.2.}\label{jac24}
\end{figure*}
%\FloatBarrier
\begin{figure*}
	%\underline{Inhomogeneous case:}
	\centering	
	\begin{tabular}{cc}
		%\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{camragiven.eps}}
		%\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{breginitial.eps}}\\
		  %ZAKKY%\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{camragiven.eps}}
		  %ZAKKY%\subfigure{\label{plane}\includegraphics[width=2cm,height=2cm]{bregfcont1.eps}}
		  %ZAKKY%\subfigure{\label{wu1}\includegraphics[width=2cm,height=2cm]{bregseg1.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{givenplane.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{Bregplanef.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{Bregplaneseg.eps}}\\
		  %ZAKKY%\subfigure{\label{circle}\includegraphics[width=2cm,height=2cm]{breginitial.eps}}
		  %ZAKKY%\subfigure{\label{wu2}\includegraphics[width=2cm,height=2cm]{camrafcontmy.eps}}
		%%\subfigure[]{\label{Li1}\includegraphics[width=2cm,height=2cm]{J1initial.eps}}
		  %ZAKKY%\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{camrasegmy.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{Bregplaneinitial.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{Bregplaneourf.eps}}
		\subfigure{\label{Li1}\includegraphics[width=2cm,height=2cm]{Bregplaneourseg.eps}}
	\end{tabular}
	\caption{\color{red}These inhomogeneous images are taken from Goldstein {\it et al.} \cite{40} and Rodriguez et al. \cite{rodriguez2011nonconvex} and are publicly available: \color{black}As Bregman-split method is well known for its speed. But it can be seen that the \color{red}performance \color{black}of our proposed model is better than Goldstein {\it et al.} \cite{40}. The first row is the segmentation result of Goldstein {\it et al.} \cite{40} and second row is the segmentation result of our proposed model with $p=0.5$, $\mu$=0.7, and $\sigma$=0.6.}\label{camraman}
\end{figure*}
% comparison with machine learning Tomographic image
\begin{figure*}
	%\underline{Inhomogeneous case:}
	\centering	
	\begin{tabular}{cc}
		\subfigure[]{\label{circle}\includegraphics[width=3cm,height=2cm]{Tomographic2given.eps}}
		\subfigure[]{\label{plane}\includegraphics[width=3cm,height=2cm]{Tomographic2initial.eps}}
		\subfigure[]{\label{wu1}\includegraphics[width=3cm,height=2cm]{Tomographic2-1seg.eps}}
		\subfigure[]{\label{Li1}\includegraphics[width=3cm,height=2cm]{Tomographic2-1segFeat.eps}}
	\end{tabular}
	\caption{The segmentation results of the proposed model for the image taken from Furat {\it et al.} \cite{Furat}: (a) 2D cut-out of tomographic image data of ore particles, (b) Initial contour, (c) Final contour and (d) Segmentated result of our proposed model with $p=0.6$, $\mu$=3,  and $\sigma$=0.5.}\label{Furat}
\end{figure*}

%%%%%%%fig.15    Jaccard 1  without noise
\begin{figure*}
	\centering
	\begin{tabular}{cc}
		\subfigure%[Clear image]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{GivenGC1.eps}}
		\subfigure%[Ground truth]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{gtruth1.eps}}
		\subfigure%[Wu {\it et al.}]
		{\label{wu1}\includegraphics[width=2cm,height=2cm]{jwuhec1.eps}}
		\subfigure
		%[Li {\it et al.}]
		{\label{wu2}\includegraphics[width=2cm,height=2cm]{jlikimc1.eps}}
		%%\subfigure[]{\label{Li1}\includegraphics[width=2cm,height=2cm]{J1initial.eps}}
		\subfigure%[Krinidis {\it et al.}]
		{\label{Li1}\includegraphics[width=2cm,height=2cm]{jkrinc2.eps}}
		\subfigure%[Our]
		{\label{Li2}\includegraphics[width=2cm,height=2cm]{OURGC8.eps}}
		\\
		\subfigure%[Noisy image]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{GivenGC2.eps}}
		\subfigure%[Ground truth]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{gtruth2.eps}}
		\subfigure%[Wu {\it et al.}]
		{\label{wu1}\includegraphics[width=2cm,height=2cm]{jwuhec2.eps}}
		\subfigure%[Li {\it et al.}]
		{\label{wu2}\includegraphics[width=2cm,height=2cm]{jlikimc2.eps}}
		%\subfigure[]{\label{Li1}\includegraphics[width=2cm,height=2cm]{J2initial.eps}}
		\subfigure%[Krinidis {\it et al.}]
		{\label{Li1}\includegraphics[width=2cm,height=2cm]{jkrinc1.eps}}
		\subfigure%[Our]
		{\label{Li2}\includegraphics[width=2cm,height=2cm]{OURGC2.eps}}
	\end{tabular}
	\caption{First row image clear image, second row nosy image. Second column shows the clean image. Jaccard similarity coefficient for Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34}, Krinidis {\it et al.} \cite{313} and the proposed model with $p=0.5$,$\mu$=0.7,  and $\sigma$=3 is shown in the third, forth, fifth and sixth column. x-axis denote the iterations and y-axis is the Jaccard accuracy in each time step iteration.}\label{jac1}\label{jac2}
\end{figure*}
%\FloatBarrier

%%%%%%% 19 plane %%%%   Jaccard without noise
\begin{figure*}
	\centering
	\begin{tabular}{cc}
		\subfigure%[Clear image]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{GC8.eps}}
		\subfigure%[Ground truth]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{gtruth4.eps}}
		\subfigure%[Wu {\it et al.}]
		{\label{wu1}\includegraphics[width=2cm,height=2cm]{jwuhep2.eps}}
		\subfigure%[Li {\it et al.}]
		{\label{wu2}\includegraphics[width=2cm,height=2cm]{jlikimp2.eps}}
		%\subfigure[]{\label{Li1}\includegraphics[width=2cm,height=2cm]{J6initial.eps}}
		\subfigure%[Krinidis {\it et al.}]
		{\label{Li1}\includegraphics[width=2cm,height=2cm]{jkrinp2.eps}}
		\subfigure%[Our model]
		{\label{Li2}\includegraphics[width=2cm,height=2cm]{OURGC8.eps}}
		\\
		%\subfigure[]{\label{circle}\includegraphics[width=2cm,height=2cm]{GC7.eps}}
		\subfigure%[Noisy image]
		{\label{Li1}\includegraphics[width=2cm,height=2cm]{J5initial.eps}}
		\subfigure%[Ground truth]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{gtruth3.eps}}
		\subfigure%[Wu {\it et al.}]
		{\label{wu1}\includegraphics[width=2cm,height=2cm]{jwu-hep.eps}}
		\subfigure%[Li {\it et al.}]
		{\label{wu2}\includegraphics[width=2cm,height=2cm]{jli-kim.eps}}
		%%\subfigure[]{\label{Li1}\includegraphics[width=2cm,height=2cm]{J5initial.eps}}
		\subfigure%[Krinidis {\it et al.}]
		{\label{Li1}\includegraphics[width=2cm,height=2cm]{jkrinp1.eps}}
		\subfigure%[Our model]
		{\label{Li2}\includegraphics[width=2cm,height=2cm]{OURGC7.eps}}
		\\
		\subfigure%[Given Image]
		{\label{Eagle}\includegraphics[width=2cm,height=2cm]{eaglegivennoise.eps}}
		\subfigure%[Ground truth]
		{\label{circle}\includegraphics[width=2cm,height=2cm]{eaglegtruthnoise.eps}}
		\subfigure%[Wu {\it et al.}]
		{\label{wu1}\includegraphics[width=2cm,height=2cm]{jwuheE1.eps}}
		\subfigure%[Li {\it et al.}]
		{\label{wu2}\includegraphics[width=2cm,height=2cm]{jlikimE2.eps}}
		\subfigure%[Krinidis {\it et al.}]
		{\label{Li1}\includegraphics[width=2cm,height=2cm]{jkrinE1.eps}}
		\subfigure%[Our model]
		{\label{Li2}\includegraphics[width=2cm,height=2cm]{KrinGC6.eps}}
	\end{tabular}
	\caption{First row a clear image, second row nosy image, and third row is an is taken from Berkeley's data set. Second column shows the clean image. Jaccard similarity coefficient for Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34}, Krinidis {\it et al.} \cite{313} and the proposed model with $p=0.5$, $\mu=0.7$, and $\sigma$=3 is shown in the third, forth, fifth and sixth column. x-axis denote the iterations and y-axis is the Jaccard accuracy in each time step iteration.}  \label{jplane1}\label{jplane2} \label{jaceagle1}
\end{figure*}

%\subsection{Principal steps of the algorithm:}

%In the next section, we discuss some experimental results of our proposed model compared with other existing models.

\section{Experimental results}\label{exp}
In this section, we present experiments for real and synthetic image compare the performance of our method to other existing models such as Wu {\it et al.} \cite{15}, Krinidis {\it et al.} \cite{313}, Li {\it et al.} \cite{34}, and Wu {\it et al.} \cite{20}. The images used in our experiments are of a wide range including medical and real-world images having different sizes and different noise level. The proposed model is also tested for images with intensity in-homogeneity and compared with Goldstein {\it et al.} \cite{40}. Moreover, different initial guesses have been applied to show the proposed model does not depend on the initialization and stuck in local minima. \color{red} In our experiments the parameters $\mu=0.7$, $p=0.5$ ($p=0.6$) has been fixed through the experiments. Through the experiments we observed that for the parameter $p$ in the range  $0.5\leq p \leq 0.9$ the new model work, but from Fig. \ref{pvalue} it is clear that the best value for $p$ is $0.5$\color{black}. All the experiments were performed on a $1.61$ GHz Core $m3-7y30$ CPU $@ 1.00$ GHz with $8$ GB memory. The algorithm was implemented and carried out using Matlab 9.4, in Windows 10 environment. The image size varies from $100 \times 100$ to $256 \times 256$. \color{red}The datasets and images used during the experimental study are publicly available in the kaggle repository, and can be accessed at  [\color{red}\textit{https://www.kaggle.com/datasets/mnavaidd/image-segmentation-dataset}\color{black}]. \color{black}
%\subsection{Performance of the proposed model}
%Here we show the successful performance of the proposed model on clear and noisy images, also on images with intensity inhomogeneity.
\subsection*{Test Set 1: Global minima achievement of the new model.} 
To show the global minima achievement of the proposed model due to its convexity property we run experiments with diferent initialization.  Fig. \ref{exp1} consist of two images (noisy image with three objects and first-ever black hole image) with different initial guess as shown in the first row. As clearly the proposed model does not depend on the initial guesses to archive the same segmentation results. This indicates that the method is independent on the initialization and that there \color{red}is no need to \color{black}check several times for different initial points. \color{red} Fig. \ref{pvalue} is the experimental results of our proposed model for different value of $p$, it shows that the best value for $p$ is $0.5$. \color{black} 
\subsection*{Test Set 2: Robustness and accuracy of the new model}
This test set consist of showing the successful performance of the proposed model on noisy images with a single and multiple objects in comparison with well-known models, such as Wu {\it et al.} \cite{15}, Krinidis {\it et al.} \cite{313}, Li {\it et al.} \cite{34}, and Wu {\it et al.} \cite{20}. Fig. \ref{jac14} (Berkeley's data set) and \ref{jac15} are images in presence of high noise and outlier, Fig. \ref{foot} and Fig. \ref{child} are medical images, and Fig. \ref{moon and star} and Fig. \ref{jac24} are noisy images with multi-objects. From all this experiments it can be observed that Wu {\it et al.} \cite{15}, Krinidis {\it et al.} \cite{313}, Li {\it et al.} \cite{34}, and Wu {\it et al.} \cite{20} fail or partially fail to properly segment the objects in the given images \color{red}whereas the proposed \color{black}method gives satisfactory results. 
\subsection*{Test Set 3: Comparison of the proposed model on images with intensity inhomogeneity} Fig. \ref{camraman} shows the comparison of the proposed model and Goldstein {\it et al.} \cite{40}. \color{red}The images with intensity inhomogeneity are also taken from the Goldstein {\it et al.} \cite{40}. The images are publicly available online [\color{red}\textit{https://sites.google.com/a/istec.net/prodrig/Home}\color{black}]. This \color{black}can be observed that the proposed method gives satisfactory results as compare to  Goldstein {\it et al.} \cite{40}. Similarly, Fig. \ref{Furat} offers a comparison of obtained results using the proposed approach and the model demonstrated in \cite{Furat}. Note that, the image were taken from Furat {\it et al.} \cite{Furat}. 
\subsection*{Test Set 4: Accuracy analysis through Jaccard similarity coefficient \color{red}and Sørensen-Dice similarity\color{black}}
We evaluate the accuracy of the proposed model using the Jaccard similarity coefficient \cite{Afzal22}. One can quantifying the similarities between the obtained image $X$ and the ground truth $Y$ using the Jaccard index that is mathematically defined by:
\begin{align}\label{51}	
	J(X,Y)=\frac{|X \cap Y|}{|X\cup Y|}
\end{align}
In Figs. \ref{jac1} and \ref{jaceagle1} we show the quantitative comparison of our proposed model compared to the other existing models such as Wu {\it et al.} \cite{15}, Krinidis {\it et al.} \cite{313}, and  Li {\it et al.} \cite{34} for $5$ different images with or without noise.  It can be observed that Krinidis {\it et al.}   \cite{313} produced relatively better results compared to  Wu {\it et al.} \cite{15} and Li {\it et al.} \cite{34}, but the results of the proposed model are better than Krinidis {\it et al.} \cite{313}  as clearly seen in the last column of those figures. From the quantitative comparisons, it can be seen that the proposed model performs better than other existing models \cite{15,34,313}. Table \ref{tbl:afzal1} shows the JS coefficients comparison of our model with other competing models. The results of this table show 10 images from Berkeley's data set. It can be observed that  in terms of accuracy the proposed model is performing better than the competing three other models in almost each image. %Table \ref{tbl:afzal1} shows rather good speed performance of our proposed segmentation algorithm in comparison with the competing models. 
\begin{table*}[!htbp]
	%\subsection{$DMD_{AOS}$ vs. $DMD_{FD}$ in Terms of CPU Performance}
	\begin{center}
		\color{red}
		\begin{tabular}{|p{.8cm}|p{.5cm}|p{.8cm}|p{.8cm}|p{.5cm}|p{.8cm}|p{.8cm}|p{.5cm}|p{.8cm}|p{.8cm}|p{.5cm}|p{.8cm}|p{.8cm}|} \hline
			Image  
			&\multicolumn{3}{|c|}{proposed model}  &\multicolumn{3}{|c|}{Krinidis et al. \cite{313}} &\multicolumn{3}{|c|}{Wu et al. \cite{15}}&\multicolumn{3}{|c|}{Li et al. \cite{34}}\\
			\cline{2-13}
			&Iter.  & JS & Time & Iter.  & JS & Time   &Iter. & JS & Time  &Iter.  & JS & Time \\
			%\hline
			%   $100\times 100$ &  5&0.8 & 2&0.3 & 2&0.2\\
			%\hline
			%  $200\times200$ &  5&0.9 & 2&0.4 & 2&0.2\\
			\hline
			$1$ & 50  &0.9999 &0.0674 &50  &0.9912 &0.0935 &300  &0.7525 &0.1595 &1000  &0.2806 &0.1085\\ \hline
			$2$ & 100  &0.8872 &0.0850 &100  &0.8720 &0.1395 &300 & 0.7747 &0.1705 &1000  &0.4780 &0.0975\\ \hline
			$3$ &100  &0.9210 &0.0873 &100  &0.8890 &0.1203 &300  &0.4976 &0.1773 &800  &0.3003 &0.0935 \\ \hline
			$4$ &100  &0.7989 &0.0845 &100 & 0.7922 &0.1253 &250  &0.7227 &0.1543 &1000  &0.4043 &0.0971\\ \hline
			$5$ & 100 & 0.7210 &0.0715 &100  &0.6709 &0.1257 &300  &0.5606 &0.1784 &900  &0.5437 &0.0926\\ \hline
			$6$ & 80  &0.9194 &0.0656 &100  &0.5488 &0.0936 &300  &0.5427 &0.1475 &1000  &0.3563 &0.0981  \\ \hline
			$7$ & 100 &0.7273 &0.0881 &100  &0.5342 &0.0982  &300  &0.3440 &0.1421  &1000 & 0.2863 &0.1103 \\ \hline
			$8$ & 100  &0.7840 &0..0752 &100& 0.7839 &0.1052  &300 & 0.7372 &0.1522  &1000  &0.5321 &0.0973\\ \hline
			$9$ &100  &0.8458 &0.0661 &100  &0.8325 &0.1203 &300 &0.6735 &0.1253 &1000  &0.3425 &0.0953 \\ \hline	
			$10$ &100  &0.8053 &0.0771 & 100 &0.7832 &0.1283 &300 &0.7452 &0.1523 &1000  &0.4731 &0.1150 \\ \hline
		\end{tabular}
	\color{black}
	\end{center}
	\caption{Jaccard similarity measure, number of iterations \color{red} and CPU time (second) \color{black}of Krinidis {\it et al.} \cite{313}, Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34} and of our proposed model on 10 images from Berkeley's data set, image size $110\times110$.}\label{tbl:afzal1}
\end{table*}
\color{red}
\paragraph{Sørensen-Dice similarity}
The Sørensen-Dice similarity is computed using Eq. \ref{522}:  
\begin{align}\label{522}	
	D(X,Y)=\frac{2|X \cap Y|}{|X|+|Y|}.
\end{align}
The Sørensen-Dice similarity values are normalized and given with in the range of $[0, 1]$. The higher Dice value shows better segmentation results and vice versa.
\begin{table*}[!htbp]
	%\setlength{\tabcolsep}{3pt}
	\color{red}
	\begin{center}
		\begin{tabular}{|p{2.5cm}|p{2.5cm}|p{2.5cm}|p{2.5cm}|} \hline 
			{Krinidis et al. \cite{313}}& {Wu et al. \cite{15}} &{Li et al. \cite{34}} &{proposed model}
			\\ 
			\hline
			 0.95 ± 0.082 &  0.91 ± 0.058 &  0.80 ± 0.069 & 0.98 ± 0.085 \\ 
			\hline
			
		\end{tabular}
	\end{center}
\caption{\color{red} Sørensen-Dice similarity for Krinidis {\it et al.} \cite{313}, Wu {\it et al.} \cite{15}, Li {\it et al.} \cite{34} and of our proposed model on $10$ different images.}\label{Dice}
\color{black}	
\end{table*}
\\Table \ref{Dice} shows Sørensen-Dice coefficients comparison of our model with other competing models Krinidis {\it et al.} \cite{313}, Wu {\it et al.} \cite{15}, and  Li {\it et al.} \cite{34}. The results are obtained from experiments on $10$ different images suitable for interactive segmentation with a pre-labeled ground truth consisting of means of the labeled ground truth. It can be observed that Krinidis {\it et al.} produced relatively better results compared to Wu {\it et al.} and Li {\it et al.}, but for a high noisy or low intensity image it loses the details. From the results it is clear that the proposed model performs better than the other competing models.
\color{black}

\section{Conclusions and future work}\label{con}
This article mainly focuses to design an efficient image data term based on an unconventional and novel objective function - as given by Equation \ref{6}. The reason is that this metric is robust against the outliers by giving \color{red}fewer weights \color{black}to outliers and noise in contrast with the conventional and old objective function, given by Equation \ref{2}, which give importance to outliers. Besides this a fuzzy level set function is employed with two main benefits over the conventional level set function:  capturing  more than one phase or objects of different intensities  plays an important role while designing a convex functional. In this way, one can impose constraints for convexity, which can be efficiently implemented, avoiding the initial guess tuning. For a deeper understanding of the properties of the proposed model, a mathematical analysis is presented and demonstrated. Moreover, the Gaussian smoothing filtering is employed for the regularization of the fuzzy membership function. Furthermore, for comprehensive analysis of the performance of the proposed model qualitative and quantitative measures are performed on various images. It has been observed that the proposed novel model performs far, and much, better than the existing and latest state-of-the-art segmentation techniques. \\
Selective image segmentation is one of the most important topics in medical imaging and real applications. In the future, we will work and propose a robust selective segmentation model using a dual-level set variational formulation model that should be based on the local spatial distance. A similar model should aim to segment all objects with one level set function (global) and the selected object with another level set function (local). Furthermore, the combination of marker distance function, edge detection, local spatial distance, and active contour without edges should be considered in the future. \color{red}Outliers must be discovered and segregated during the denoising pre-processing or suitable limits must be put on the segmentation framework to ensure correct and the most appropriate image segmentation in the presence of noise and outliers. In the future, we will use suitable removing outliers criteria backed by a well-designed theory in a variational framework for accurate and appropriate image segmentation. \color{black}
%\appendices

%Appendixes, if needed, appear before the acknowledgment.
%
\section*{Acknowledgment}
The code of the proposed method will be provided for a research purpose if requested from the principal author through email. The work is supported, in parts, by the University of Peshawar and, in parts, by the Abdul wali Khan University, Pakistan.
%The authors would like to express their gratitude to the anonymous referees for their valuable comments.

\section*{Data availability}
The datasets generated and/or analysed during the current study are publicly available in the kaggle repository, and can be accessed at [\textit{https://www.kaggle.com/datasets/mnavaidd/image-segmentation-dataset}]. Moreover, various images used within the experimental work are publicly available online.

\begin{thebibliography}{00}

\bibitem{8}
D. Mumford, J. Shah, ``Optimal approximations by piecewise smooth functions and associated variational problems,'' \emph{Commun. Pure Appl. Math}., vol. 42, no. 5, pp. 577-685, Jul. 1989.

\bibitem{7}
V. Caselles, R. Kimmel, G. Sapiro, ``Geodesic active contours,''
\emph{Int. J. Comput. Vis}., vol. 22, no. 1, pp. 61-79, Feb. 1997.

\bibitem{1}
C. Li, L. Liu, X. Sun, J. Zhao, J. Yin, ``Image segmentation based on fuzzy clustering with cellular automata and features weighting,'' \emph{EURASIP J.Image Vide}., vol. 1, pp. 1-11, Dec. 2019.

\bibitem{10}
M. Kass, A. Witkin, D. Terzopoulos, ``Snakes: Active contour models,''
\emph{Int. J. Comput. Vis}., vol. 1, no. 4, pp. 321-331, Jan. 1988.

\bibitem{11}
S. Osher, J. A. Sethian, ``Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,''
\emph{J. Comput. Phys}., vol. 79, no. 1, pp. 12-49, Nov. 1988.

\bibitem{2}
S. Pare, A. Kumar, G. K. Singh, V. Bajaj, ``Image Segmentation Using Multilevel Thresholding: A Research Review,'' \emph{Iran. J. Sci. Technol. Trans. Electr. Eng}., vol. 44, pp. 1-29, Mar. 2020.

\bibitem{18}
B. Appleton, H. Talbot, ``Globally optimal geodesic active contours,''
\emph{J. Math. Imaging. Vis}., vol. 23, no. 1, pp. 67-86, Jul. 2005.

\bibitem{27}
A. Yezzi, S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, ``A geometric snake model for segmentation of medical imagery,''
\emph{IEEE Trans. med. imaging}., vol. 16, no. 2, pp. 199-209, Apr. 1997.

\bibitem{28}
T. F. Chan, L. A. Vese, ``Active contours without edges,''
\emph{IEEE Trans. Image Process}., vol. 10, no. 2, pp. 266-277, Feb. 2001.

\bibitem{29}
T. F. Chan, S. Esedoglu, M. Nikolova, ``Algorithms for finding global minimizers of image segmentation and denoising models,''
\emph{SIAM j. Appl. math}., vol. 66, no. 5, pp. 1632-1648, 2006.

\bibitem{3}
L. D. Cohen, ``On active contour models and balloons,'' \emph{CVGIP: Image Und}., vol. 53, no. 2, pp. 211-218, Mar. 1991.

\bibitem{20}
Y. Wu, W. Ma, M. Gong, H. Li, L. Jiao, ``Novel fuzzy active contour model with kernel metric for image segmentation,''
\emph{Appl. Soft Comput}., vol. 34, pp. 301-311, Sep. 2015.

\bibitem{31}
J. Lie, M. Lysaker, X. C. Tai, ``A binary level set model and some applications for Mumford {\it et al.} image segmentation,''
\emph{IEEE Trans. Image Process}., vol. 15, no. 4, pp. 1171-1181, Apr. 2006.
\color{red}
\bibitem{33}
K. Zhang, L. Zhang, H. Song, W. Zhou, ``Active contours with selective local or global segmentation: a new formulation and level set method,''
\emph{Image. Vision. comput}., vol. 28, no. 4, pp. 668-676, Apr. 2010.
\color{black}

\bibitem{34}
Y. Li, J. Kim, ``An unconditionally stable numerical method for bimodal image segmentation,''
\emph{Appl. Math. Comput}., vol. 219, no. 6, pp. 3083-3090, Nov. 2012.

\bibitem{36}
A. Mondal, S. Ghosh, A. Ghosh, ``Robust global and local fuzzy energy based active contour for image segmentation,''
\emph{Appl. Soft Comput}., vol. 47, pp. 191-215, Oct. 2016.

\bibitem{Wang}
C. Wang and B. Zhu, ``Image Segmentation and Adaptive Contrast Enhancement for Haze Removal,''
\emph{IEEE 63rd International Midwest Symposium on Circuits and Systems (MWSCAS)}., pp. 1036-1039, 2020, doi: 10.1109/MWSCAS48704.2020.9184525.

\bibitem{Joshi}
A. Joshi, M. S. Khan, S. Soomro, A. Niaz, B. S. Han and K. N. Choi, ``SRIS: Saliency-Based Region Detection and Image Segmentation of COVID-19 Infected Cases,''
\emph{IEEE Access}., vol. 8, pp. 190487-190503, 2020, doi: 10.1109/ACCESS.2020.3032288.

\bibitem{Chuang}
K. S. Chuang, H. L. Tzeng, S. Chen, J.Wu, T. J. Chen, ``Fuzzy c-means clustering with spatial information for image segmentation,''
\emph{computerized medical imaging and graphics}., vol. 30, no. (1) pp. 9-15, Jan. 2006.

\bibitem{Tripathy}
B. K. Tripathy, A. Basu, S. Govel, ``Image segmentation using spatial intuitionistic fuzzy C means clustering,''
\emph{IEEE International Conference on Computational Intelligence and Computing Research}., pp. 1-5, Dec. 2014.

\bibitem{5}
Y. Chen, H. D. Tagare, S. Thiruvenkadam, F. Huang, D. Wilson, K. S. Gopinath, E. A. Geiser, ``Using prior shapes in geometric active contours in a variational framework,''
\emph{Int. J. Comput. Vis}., vol. 50, no. 3, pp. 315-328, Dec. 2002.

\bibitem{9}
H. Ali, N. Shah, K. Chen, G. A. Khan, ``A variational model with hybrid images data fitting energies for segmentation of images with intensity inhomogeneity,'' \emph{Pattern Recognit}., vol. 51, pp. 27-42, Mar. 2016.

\bibitem{13}
C. Li, C. Xu, C. Gui, M. D. Fox, ``Level set evolution without re-initialization: a new variational formulation,''
\emph{Proc. IEEE Conf. Comput. Vis. Pattern Recognit}., vol. 1, pp. 430-436, Jun. 2005.
\color{red}
\bibitem{15}
Y. Wu, C. He, ``A convex variational level set model for image segmentation,''
\emph{Signal Process}., vol. 106, pp. 123-133, Jan. 2015.
\color{black}
\bibitem{24}
L. I. Rudin, S. Osher, E. Fatemi, ``Nonlinear total variation based noise removal algorithms,''
\emph{Physica D: nonlinear phenomena}., vol. 60, no. 1-4, pp. 259-268, Nov. 1992.

\bibitem{selec}
N. Badshah, K. Chen2, H. Ali and G. Murtaza. Coefficient of Variation Based Image Selective Segmentation Model Using Active Contours. \emph{East Asian Journal on Applied Mathematics}., vol. 2, no. 2, pp. 150-169, May 2012. doi: 10.4208/eajam.090312.190412a.

\bibitem{41}
Y. Yu, S. T. Acton, ``Edge Detection in Ultrasound Imagery Using the Instantaneous Coefficient of Variation,''
\emph{IEEE Trans. Image Process}., vol. 13, no. 12, pp. 1640 - 1655J, Jan 2005.
DOI: 10.1109/TIP.2004.836166.

\bibitem{35}
H. Ali, N. Shah, K. Chen, G. A. Khan, N. Zikria, ``Multiphase segmentation based on new signed pressure force functions and one level set function,''
\emph{Turk. J. Electr. Eng. Co}., vol. 25, no. 4, pp. 2943-2955, Jul. 2017.

\bibitem{37}
J. Oh, N. Kwak, ``Generalized mean for robust principal component analysis,''
\emph{Pattern Recognit}., vol. 54, pp. 116-127, Jun. 2016.

\bibitem{4}
S. Balla-Arabé, X. Gao, B. Wang, ``A fast and robust level set method for image segmentation using fuzzy clustering and lattice Boltzmann method,''
\emph{IEEE T. Cybernetics}., vol. 43, no. 3, pp. 910-920, Mar. 2013.

\bibitem{313}
S. Krinidis, V. Chatzis, ``Fuzzy energy-based active contours,''
\emph{IEEE Trans. Image Process}., vol. 18, no. 12, pp. 2747-2755, Aug. 2009.

\bibitem{16}
B. Song, T. Chan, ``A fast algorithm for level set based optimization,''
\emph{UCLA CamReport}., vol. 2, pp. 68, Dec. 2002.

\bibitem{39}
S. H. Lee, J. K. Seo, ``Level set-based bimodal segmentation with stationary global minimum,''
\emph{IEEE Trans. Image Process}., vol. 15, no. 9, pp. 2843-2852, Aug. 2006.

\bibitem{43}
H. Ali, A. Shujjahuddin, L. Rada, ``A New Active Contours Image Segmentation Model Driven by Generalized Mean with Outlier Restoration Achievements''
\emph{International J. of Pattern Recognit and Artificial Intelligence}., 2019.  DOI: 10.1142/S0218001420540269.

\bibitem{21}
M. Gong, Y. Liang, J. Shi, W. Ma, J. Ma, ``Fuzzy c-means clustering with local information and kernel metric for image segmentation,''
\emph{IEEE Trans. Image Process}., vol. 22, no. 2, pp. 573-584, Sep. 2013.

\bibitem{38}
X. Cai, ``Variational image segmentation model coupled with image restoration achievements,''
\emph{Pattern Recognit}., vol. 48, no. 6, pp. 2029-2042, Jun. 2015.

\bibitem{40}
T. Goldstein, X. Bresson, S. Osher, ``Geometric applications of the split Bregman method: segmentation and surface reconstruction,''
\emph{J. Sci. Comput}., vol. 45, no. 1-3, pp. 272-93, Oct. 2010.

\bibitem{Afzal22}
A. Rahman, H. Ali, N. Badshah, L. Rada, A. A. Khan, H. Hussain, M. Zakarya, A. Ahmed, I. U. Rahman, M. Raza, M. Haleem,  
``A Selective Segmentation Model Using Dual-Level Set Functions and Local Spatial Distance,'' 
\emph{{IEEE} Access,} vol. 10, pp. 22344-22358, Mar 2022.

\bibitem{13a}
C. Li,  C. Y. Kao,  J. C. Gore, Z. Ding,
``Implicit active contours driven by local binary fitting energy,''
\emph{in proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR).,} vol. 42, pp. 1-7, Jun. 2007.

\bibitem{Furat}
O. Furat, M. Wang, M. Neumann, L. Petrich, M. Weber, C. E. Krill III, V. Schmidt, 
``Machine learning techniques for the segmentation of tomographic image data of functional materials,'' 
\emph{Frontiers in Materials.,} vol. 6, pp. 145, Jun. 2019.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\color{red}
\bibitem{WuLGC21}
Yongfei Wu, Xilin Liu, Peiting Gao, Zehua Chen, 
``A variational level set model with closed-form solution for bimodal image segmentation,'' 
\emph{Multim. Tools Appl.,} vol. 80, pp. 25943-25963, Jul 2021.

\bibitem{LiuHGWR19}
Yang Liu, Chuanjiang He, Peiting Gao, Yongfei Wu, Zemin Ren,  
``A binary level set variational model with \emph{L}\({}^{\mbox{1}}\) data term for image segmentation,'' 
\emph{Signal Process.,} vol. 155, pp. 193-201, Feb 2020.

\bibitem{LiuHWR18}
Yang Liu, Chuanjiang He, Yongfei Wu, Zemin Ren, 
``The \emph{L}\({}_{\mbox{0}}\)-regularized discrete variational level set method for image segmentation,'' 
\emph{Image Vis. Comput.,} vol. 75, pp. 32-43, Feb 2020.

\bibitem{LiuHW18}
Yang Liu, Chuanjiang He, Yongfei Wu, 
``Variational model with kernel metric-based data term for noisy image segmentation,'' 
\emph{Digit. Signal Process.,} vol. 78, pp. 42-55, Feb 2020.

\bibitem{WuH16}
Yongfei Wu, Chuanjiang He, 
``Indirectly regularized variational level set model for image segmentation,'' 
\emph{Neurocomputing,} vol. 171, pp. 194-208, May 2017.

\bibitem{rodriguez2011nonconvex}
Rodriguez, Paul, 
``Nonconvex total variation speckled image restoration via nonnegative quadratic programming algorithm,'' 
\emph{19th European Signal Processing Conference, IEEE,} pp. 288-292, 2011.

\color{black}
\end{thebibliography}


\newpage
\section*{Appendix}\label{Appendix}
\color{red}
\par{\bf{Proof of Theorem 1:}}
For simplicity, consider the energy functional in Eq. (\ref{31}) as follow:
\begin{align}\label{39}
	f(\zeta)=\mu f_{1}(\zeta)+f_{2}(\zeta)+f_{3}(\zeta)
\end{align}
where $\zeta=(x,y)$,
\begin{align}
	&f_{1}(\zeta)=\int_{\Omega}|\bigtriangledown \mathbf{z}(\zeta)|d\zeta, \label{40}
	%\end{align}
	%and
	%\begin{align}\label{f2}
	\\ &f_{2}(\zeta)=\int_{\Omega}\alpha(\zeta)||\mathbf{u}(\zeta)-c_1||_{2}^{2}[\mathbf{z}(\zeta)]^{m}d\zeta, \label{41}
	%\end{align}
	%\text{and}
	%\begin{align}\label{simple}
	\\ &f_{3}(\zeta)=\int_{\Omega}\beta(\zeta)||\mathbf{u}(\zeta)-c_2||_{2}^{2}[1-\mathbf{z}(\zeta)]^{m}d\zeta. \label{42}
\end{align}
First of all, the domain $\Omega$ is convex, because it is a rectangle. The function $f_{1}(\zeta)$ is convex as in \cite{24}.
Consider
\begin{align}\label{43}
	f_{2}(\zeta)=\int_{\Omega}\alpha(\zeta)||\mathbf{u}(\zeta)-c_1||_{2}^{2}[\mathbf{z}(\zeta)]^{m}d\zeta.
\end{align}
%%%%%%%%%%kam ngec
Taking $F_{2}(\zeta)=\alpha(\zeta)||\mathbf{u}(\zeta)-c_1||_{2}^{2}[\mathbf{z}(\zeta)]^{m}$,
%\begin{eqnarray}\label{simple}
%F_{2}(\zeta)=\alpha(x,y)||I-c_1||_{2}^{2}[z]^{m}
%\end{eqnarray}
where $F_{2}:\Omega\rightarrow R$ such that
\begin{align}\label{44}
	f_{2}(\zeta)=\int_{\Omega}F_{2}(\zeta)d\zeta.
\end{align}
Let ${\zeta_1=(x_{1},y_{1}), \zeta_2=(x_{2},y_{2})}\in \Omega$ and $\kappa \in [0,1]$, since $\Omega$ is convex,  we can write:
\begin{align}\label{45}
	\kappa\zeta_1+&(1-\kappa)\zeta_2=(\kappa(x_{1},y_{1})+(1-\kappa)(x_{2},y_{2})) \nonumber \\
	&=(\kappa(x_{1}-x_{2})+x_{2},\kappa(y_{1}-y_{2})+y_{2})\in \Omega.
\end{align}
%$\Omega$ is convex because $\kappa(x_{1}-x_{2})+x_{2} \in R$ as $\kappa \in [0,1]$ and $x_{1}-x_{2} \in R$, also $\kappa(y_{1}-y_{2})+y_{2} \in R$ as $y_{1}-y_{2} \in R$, this implies that $\kappa\zeta_1+(1-\kappa)\zeta_2 \in \Omega$.\\
%\par{}
Taking the derivative of $F_{2}(\zeta)$ with respect to the function $\mathbf{z}(\zeta)$, we get
\begin{eqnarray}\label{46}
	\frac{\partial F_{2}}{\partial \mathbf{z}}=m[\mathbf{z}(\zeta)]^{m-1}\alpha(\zeta)||\mathbf{u(\zeta)}-c_1||_{2}^{2}.
\end{eqnarray}
Differentiating again with respect to $\mathbf{z}(\zeta)$, we get
\begin{eqnarray}\label{47}
	\frac{\partial^{2} F_{2}}{\partial \mathbf{z}^{2}}=m(m-1)[\mathbf{z}(\zeta)]^{m-2}\alpha(x,y)||\mathbf{u(\zeta)}-c_1||_{2}^{2}.
\end{eqnarray}
\emph{$\frac{\partial^{2} F_{2}}{\partial \mathbf{z}^{2}} \geq 0$}, as $\mathbf{z}(\zeta) \in [0,1]$, $m=2$, $\alpha(\zeta) \geq 0$ and $||\mathbf{u}-c_1||_{2}^{2} \geq 0$, also $\Omega$ is convex. Thus $F_{2}(\zeta)$ is convex and for all $\zeta_1, \zeta_2 \in \Omega$ and $\kappa \in [0,1]$ the inequality
\begin{eqnarray}\label{48}
	F_{2}(\kappa\zeta_1+(1-\kappa)\zeta_2) \leq \kappa F_{2}(\zeta_1)+(1-\kappa)F_{2}(\zeta_2)
\end{eqnarray}
holds. From Eq. (\ref{48}), we have
\begin{align}\label{49}
	\int_{\Omega}F_{2}(\kappa\zeta_1&+(1-\kappa)\zeta_2)d\zeta \nonumber\\ & \leq \kappa \int_{\Omega}F_{2}(\zeta_1)d\zeta+(1-\kappa)\int_{\Omega}F_{2}(\zeta_2)d\zeta.
\end{align}
Using Eq. (\ref{39}), we get
\begin{eqnarray}\label{50}
	f_{2}(\kappa\zeta_1+(1-\kappa)\zeta_2) \leq \kappa f_{2}(\zeta_1)+(1-\kappa)f_{2}(\zeta_2),
\end{eqnarray}
which means that $f_{2}$ is convex. In the same way, one can prove the convexity of $f_{3}$. \color{red}Thus, \color{black}$f(\zeta)$ is convex with respect to $\mathbf{z}(\zeta)$ being the sum of convex functions.\color{black}\\
%proved\\
\par{\bf{Proof of Theorem 2:}} \label{Theorem 2}
Let $\{\mathbf{z}^{n}\}$ be a minimizing sequence of the energy functional Eq. (\ref{31}), then there exists a constant $M$, such that $F(\mathbf{z}^{n},c_1,c_2,\alpha,\beta) \leq M$. This implies that
\begin{align}\label{52}
	&\mu \int_{\Omega}|\bigtriangledown \mathbf{z}^{n}(x,y)|dxdy \nonumber\\ &+\int_{\Omega}\alpha(x,y)||I-c_1||_{2}^{2}[\mathbf{z}^{n}(x,y)]^{m}dxdy \nonumber \\
	&+ \int_{\Omega}\beta(x,y)||I-c_2||_{2}^{2}[1-\mathbf{z}^{n}(x,y)]^{m}dxdy \leq M.
\end{align}
The constraint $0\leq \mathbf{z} \leq 1$, ensure that $\{\mathbf{z}^{n}\}$ is uniformly bounded in $BV(\Omega)$. \color{red}Moreover, \color{black}$BV(\Omega)$ is compact w.r.t $BV_{w}^{\ast}(\Omega)$ topology, then for the subsequence which we also denote by $\{\mathbf{z}^{n}\}$, $\exists$ $\{\mathbf{z}^{\ast}\} \in BV(\Omega)$ such that
%\begin{center}
$\mathbf{z}^{n} \xrightarrow{L^{1}(\Omega)} \mathbf{z}^{\ast}$ \text{and} $\mathbf{z}^{n} \rightarrow \mathbf{z}^{\ast}$ \text{a.e} $x \in \Omega$
%\end{center}
and by convergence result the constraint $0\leq \mathbf{z}^{\ast} \leq 1$ also holds. Since $\mathbf{z}^{n} \rightarrow \mathbf{z}^{\ast}$ implies that $[\mathbf{z}^{n}]^{m} \rightarrow [\mathbf{z}^{\ast}]^{m}$ and this implies that $\alpha ||\mathbf{u}-c_1||_{2}^{2}[\mathbf{z}^{n}]^{m} \rightarrow \alpha ||\mathbf{u}-c_1||_{2}^{2}[\mathbf{z}^{\ast}]^{m}$. By Fatou's lemma,
\begin{align}\label{53}
	\int_{\Omega} \alpha(\zeta) &||\mathbf{u}-c_1||_{2}^{2}[\mathbf{z}^{\ast}(\zeta)]^{m} d\zeta \nonumber \\&\leq \lim_{n\to\infty} \inf \int_{\Omega}\alpha(\zeta) ||\mathbf{u}-c_1||_{2}^{2}[\mathbf{z}^{n}(\zeta)]^{m} d\zeta
\end{align}
similarly we can write
\begin{align}\label{54}
	\int_{\Omega} \beta(\zeta) &||\mathbf{u}-c_2||_{2}^{2}[1-\mathbf{z}^{\ast}(\zeta)]^{m} d\zeta \nonumber \\&\leq \lim_{n\to\infty} \inf \int_{\Omega}\beta(\zeta) ||\mathbf{u}-c_2||_{2}^{2}[1-\mathbf{z}^{n}(\zeta)]^{m} d\zeta
\end{align}
also by lower semi-continuity of total variation, we have
\begin{eqnarray}\label{55}
	\int_{\Omega} |\nabla \mathbf{z}^{\ast}(\zeta)| d\zeta \leq \lim_{n\to\infty} \inf \int_{\Omega} |\mathbf{z}^{n}(\zeta)| d\zeta.
\end{eqnarray}
From Eqs. (\ref{41}), (\ref{42}) and (\ref{47}), we have
\begin{eqnarray}\label{56}
	F(\mathbf{z}^{\ast},c_1,c_2,\alpha,\beta) \leq \lim_{n\to\infty} \inf F(\mathbf{z}^{n},c_1,c_2,\alpha,\beta)
\end{eqnarray}
thus $\mathbf{z}^{\ast} \in \Lambda$ (this complete the proof).\\
\color{red}Therefore, the minimizer of our proposed model has a global minimum.\color{black}


\begin{comment}
\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{afzalPhoto.JPG}}]{Afzal Rahman} %(M'76--SM'81--F'87)
	received the M.sc degree from
	the Department of Mathematics University of Peshawar, in 2010, the M.S. degree in Applied Mathematics from Trento University, in 2015, and currently he is a Ph.D. scholar at university of Peshawar, Pakistan. His research interests include variational models and computational algorithms, image processing, machine learning, image segmentation, image inpainting, medical image analysis.
\end{IEEEbiography}


\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{haiderphpto.JPG}}]{Haider Ali} 
	received his M.Sc. degree in mathematics from the
	University of Peshawar in 2007,
	M.S. degree in applied mathematics (variational models in image processing) from University of
	Engineering and Technology Peshawar, in 2012 and Ph.D. degree in applied mathematics (variational
	models in image processing), from
	the University of Peshawar, Pakistan, in 2016. He has been a lecturer and researcher at the University of Peshawar, Pakistan since 2010. His research interests
	include image processing, medical image analysis, robotics, machine learning, segmentation, registration using PDEs and variational methods.
	
\end{IEEEbiography}

\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{noor.JPG}}]{Noor Badshah} 
	received the B.Sc. and M.Sc. degrees
	in mathematics from the University of Peshawar
	and the Ph.D. degree in applied mathematics (fast
	algorithms for image segmentation and variational
	models) from the Department of Mathematical Sciences, University of Liverpool, Liverpool, U.K. His
	research interests include fast algorithms for variational image processing models and medical image
	analysis.
\end{IEEEbiography}

\end{comment}

%\EOD

\end{document}
