# Difference of Anisotropic and Isotropic TV for Segmentation under Blur and
Poisson Noise

Kevin Bui Department of Mathematics; University of California, Irvine; Irvine,
CA 92697, United States<EMAIL_ADDRESS>Yifei Lou Department of Mathematical
Sciences; University of Texas, Dallas; Richardson, TX 75080, United States;
<EMAIL_ADDRESS>Fredrick Park Department of Mathematics & Computer
Science; Whittier College; Whittier, CA 90602, United States;
<EMAIL_ADDRESS>Jack Xin Department of Mathematics; University of
California, Irvine; Irvine, CA 92697, United States<EMAIL_ADDRESS>
###### Abstract

In this paper, we aim to segment an image degraded by blur and Poisson noise.
We adopt a smoothing-and-thresholding (SaT) segmentation framework that finds
a piecewise-smooth solution, followed by $k$-means clustering to segment the
image. Specifically for the image smoothing step, we replace the least-squares
fidelity for Gaussian noise in the Mumford-Shah model with a maximum posterior
(MAP) term to deal with Poisson noise and we incorporate the weighted
difference of anisotropic and isotropic total variation (AITV) as a
regularization to promote the sparsity of image gradients. For such a
nonconvex model, we develop a specific splitting scheme and utilize a proximal
operator to apply the alternating direction method of multipliers (ADMM).
Convergence analysis is provided to validate the efficacy of the ADMM scheme.
Numerical experiments on various segmentation scenarios (grayscale/color and
multiphase) showcase that our proposed method outperforms a number of
segmentation methods, including the original SaT.

##  1 Introduction

Image segmentation partitions an image into multiple, coherent regions, where
pixels of one region share similar characteristics such as colors, textures,
and edges. It remains an important yet challenging problem in computer vision
that has various applications, including magnetic resonance imaging [26, 40,
62] and microscopy [7, 81]. One of the most fundamental models for image
segmentation is the Mumford-Shah model [50] because of its robustness to
noise. Given an input image $f:\Omega\rightarrow\mathbb{R}$ defined on an
open, bounded, and connected domain $\Omega\subset\mathbb{R}^{2}$, the
Mumford-Shah model is formulated as

$\displaystyle\min_{u,\Gamma}E_{MS}(u,\Gamma)\coloneqq$
$\displaystyle\frac{\lambda}{2}\int_{\Omega}(f-u)^{2}\;dx+\frac{\mu}{2}\int_{\Omega\setminus\Gamma}|\nabla
u|^{2}\;dx+\text{Length}(\Gamma),$ (1)

where $u:\Omega\rightarrow\mathbb{R}$ is a piecewise-smooth approximation of
the image $f$, $\Gamma\subset\Omega$ is a compact curve representing the
region boundaries, and $\lambda,\mu>0$ are the weight parameters. The first
term in (1) is the fidelity term that ensures that the solution $u$
approximates the image $f$. The second term enforces $u$ to be piecewise
smooth on $\Omega\setminus\Gamma$. The last term measures the perimeter, or
more mathematically the one-dimensional Haussdorf measure in $\mathbb{R}^{2}$
[4], of the curve $\Gamma$. However, (1) is difficult to solve because the
unknown set of boundaries needs to be discretized. One common approach
involves approximating the objective function in (1) by a sequence of elliptic
functionals [1].

Alternatively, Chan and Vese (CV) [18] simplified (1) by assuming the solution
$u$ to be piecewise constant with two phases or regions, thereby making the
model easier to solve via the level-set method [51]. Let the level-set
function $\phi$ be Lipschitz continuous and be defined as follows:

$\displaystyle\begin{cases}\phi(x)>0&\text{ if }x\text{ is inside }\Gamma,\\\
\phi(x)=0&\text{ if }x\text{ is at }\Gamma,\\\ \phi(x)<0&\text{ if }x\text{ is
outside }\Gamma.\end{cases}$

By the definition of $\phi$, the curve $\Gamma$ is represented by $\phi(x)=0.$
The image region can be defined as either inside or outside the curve
$\Gamma$. In short, the CV model is formulated as

$\displaystyle\min_{c_{1},c_{2},\phi}E_{CV}(c_{1},c_{2},\phi)$
$\displaystyle\coloneqq\lambda\left(\int_{\Omega}|f-c_{1}|^{2}H(\phi)\;dx+\int_{\Omega}|f-c_{2}|^{2}(1-H(\phi))\;dx\right)+\nu\int_{\Omega}|\nabla
H(\phi)|\;dx,$ (2)

where $\lambda,\nu$ are weight parameters, the constants $c_{1},c_{2}$ are the
mean intensity values of the two regions, and $H(\phi)$ is the Heaviside
function defined by $H(\phi)=1$ if $\phi\geq 0$ and $H(\phi)=0$ otherwise. A
convex relaxation [17] of (2) was formulated as

$\displaystyle\min_{c_{1},c_{2},u\in[0,1]}\lambda\left(\int_{\Omega}|f-c_{1}|^{2}u\;dx+\int_{\Omega}|f-c_{2}|^{2}(1-u)\;dx\right)+\nu\int_{\Omega}|\nabla
u|\;dx,$

where an image segmentation $\tilde{u}$ is obtained by thresholding $u$, that
is

$\displaystyle\tilde{u}(x)=\begin{cases}1&\text{ if }u(x)>\tau,\\\ 0&\text{ if
}u(x)\leq\tau,\end{cases}$

for some value $\tau\in(0,1)$. It can be solved efficiently by convex
optimization algorithms, such as the alternating direction method of
multipliers (ADMM) [6] and primal-dual hybrid gradient [15]. A multiphase
extension of (2) was proposed in [64], but it requires that the number of
regions to be segmented is a power of 2. For segmenting into an arbitrary
number of regions, fuzzy membership functions were incorporated [38].

Cai et al. [12] proposed the smoothing-and-thresholding (SaT) framework that
is related to the model (1). In the smoothing step of SaT, a convex variant of
(1) is formulated as

$\displaystyle
u^{*}=\arg\min_{u}\frac{\lambda}{2}\int_{\Omega}(f-Au)^{2}\;dx+\frac{\mu}{2}\int_{\Omega}|\nabla
u|^{2}\;dx+\int_{\Omega}|\nabla u|\;dx,$ (3)

yielding a piecewise-smooth solution $u^{*}$. The blurring operator $A$ is
included in the case when the image $f$ is blurred. The total variation (TV)
term $\int_{\Omega}|\nabla u|\;dx$ is a convex approximation of the length
term in (2) by the coarea formula [17]. After the smoothing step, a
thresholding step is applied to the smooth image $u^{*}$ to segment it into
multiple regions. The two-stage framework has many advantages. First, the
smoothing model (3) is strongly convex, so it can be solved by any convex
optimization algorithm to obtain a unique solution $u^{*}$. Second, the user
can adjust the number of thresholds to segment $u^{*}$ and the threshold
values to obtain a satisfactory segmentation result, thanks to the flexibility
of the thresholding step. Furthermore, the SaT framework can be adapted to
color images by incorporating an intermediate lifting step [10]. Before
performing the thresholding step, the lifting step converts the RGB space to
Lab (perceived lightness, red- green and yellow-blue) color space and
concatenates both RGB and Lab intensity values into a six-channel image. The
multi-stage framework for color image segmentation is called smoothing,
lifting, and thresholding (SLaT).

One limitation of (3) lies in the $\ell_{2}$ fidelity term that is
statistically designed for images corrupted by additive Gaussian noise, and as
a result, the smoothing step is not applicable to other types of noise
distribution. In this paper, we aim at Poisson noise, which is commonly
encountered when an image is taken by photon-capturing devices such as in
positron emission tomography [63] and astronomical imaging [36]. By using the
data fidelity term of $Au-f\log Au$ [37], we obtain a smoothing model that is
appropriate for Poisson noise [16]:

$\displaystyle\min_{u}\lambda\int_{\Omega}(Au-f\log
Au)\;dx+\frac{\mu}{2}\int_{\Omega}|\nabla u|^{2}\;dx+\int_{\Omega}|\nabla
u|\;dx.$ (4)

As a convex approximation of the length term in (1), the TV term in (4) can be
further improved by nonconvex regularizations. The TV regularization is
defined by the $\ell_{1}$ norm of the image gradient. Literature has shown
that nonconvex regularizations often yield better performance than the convex
$\ell_{1}$ norm in identifying sparse solutions. Examples of nonconvex
regularization include $\ell_{p},0<p<1,$ [13, 20, 77],
$\ell_{1}-\alpha\ell_{2},\alpha\in(0,1]$ [25, 29, 42, 45, 47],
$\ell_{1}/\ell_{2}$ [56, 67, 75], and an error function [31]. Lou et al. [48]
designed a TV version of $\ell_{1}-\alpha\ell_{2}$ called the weighted
difference of anisotropic–isotropic total variation (AITV), which outperforms
TV in various imaging applications, such as image denoising [48], image
reconstruction [48, 42], and image segmentation [8, 9, 74].

In this paper, we propose an AITV variant of (4) to improve the smoothing step
of the SaT/SLaT framework for images degraded by Poisson noise and/or blur.
Incorporating AITV regularization is motivated by our previous works [8, 9,
53], where we demonstrated that AITV regularization is effective in preserving
edges and details, especially under Gaussian and impulsive noise. To maintain
similar computational efficiency as the original SaT/SLaT framework, we
propose an ADMM algorithm that utilizes the $\ell_{1}-\alpha\ell_{2}$ proximal
operator [46]. The main contributions of this paper are as follows:

  * •

We propose an AITV-regularized variant of (4) and prove the existence of a
minimizer for the model.

  * •

We develop a computationally efficient ADMM algorithm and provide its
convergence analysis under certain conditions.

  * •

We conduct numerical experiments on various grayscale/color images to
demonstrate the effectiveness of the proposed approach.

The rest of the paper is organized as follows. Section 2 describes the
background information such as notations, Poisson noise, and the SaT/SLaT
framework. In Section 3, we propose a simplified Mumford-Shah model with AITV
and a MAP data fidelity term for Poisson noise. In the same section, we show
that the model has a global minimizer and develop an ADMM algorithm with
convergence analysis. In Section 4, we evaluate the performance of the AITV
Poisson SaT/SLaT framework on various grayscale and color images. Lastly, we
conclude the paper in Section 5.

##  2 Preliminaries

###  2.1 Notation

Throughout the rest of the paper, we represent images and mathematical models
in discrete notations (i.e., vectors and matrices). An image is represented as
an $M\times N$ matrix, and hence the image domain is denoted by
$\Omega=\\{1,2,\ldots,M\\}\times\\{1,2,\ldots,N\\}$. We define two inner
product spaces: $X\coloneqq\mathbb{R}^{M\times N}$ and $Y\coloneqq X\times X$.
Let $u\in X$. For shorthand notation, we define $u\geq 0$ if $u_{i,j}\geq 0$
for all $(i,j)\in\Omega$. The discrete gradient operator $\nabla:X\rightarrow
Y$ is defined by $(\nabla
u)_{i,j}=\left((\nabla_{x}u)_{i,j},(\nabla_{y}u)_{i,j}\right)$, where

$\displaystyle(\nabla_{x}u)_{i,j}=\begin{cases}u_{i,j}-u_{i,j-1}&\text{ if
}2\leq j\leq N,\\\ u_{i,1}-u_{i,N}&\text{ if }j=1,\end{cases}\mbox{\ and\
}(\nabla_{y}u)_{i,j}=\begin{cases}u_{i,j}-u_{i-1,j}&\text{ if }2\leq i\leq
M,\\\ u_{1,j}-u_{M,j}&\text{ if }i=1.\end{cases}$

The space $X$ is equipped with the standard inner product
$\langle\cdot,\cdot\rangle_{X}$ and Euclidean norm $\|\cdot\|_{2}$. The space
$Y$ has the following inner product and norms: for $p=(p_{1},p_{2})\in Y$ and
$q=(q_{1},q_{2})\in Y$,

$\displaystyle\langle p,q\rangle_{Y}$ $\displaystyle=\langle
p_{1},q_{1}\rangle_{X}+\langle p_{2},q_{2}\rangle_{X},$
$\displaystyle\|p\|_{1}$
$\displaystyle=\sum_{i=1}^{M}\sum_{j=1}^{N}|(p_{1})_{i,j}|+|(p_{2})_{i,j}|,$
$\displaystyle\|p\|_{2}$
$\displaystyle=\sqrt{\sum_{i=1}^{M}\sum_{j=1}^{N}|(p_{1})_{i,j}|^{2}+|(p_{2})_{i,j}|^{2}},$
$\displaystyle\|p\|_{2,1}$
$\displaystyle=\sum_{i=1}^{M}\sum_{j=1}^{N}\sqrt{(p_{1})_{i,j}^{2}+(p_{2})_{i,j}^{2}}.$

For brevity, we omit the subscript $X$ or $Y$ in the inner product when its
context is clear.

###  2.2 AITV Regularization

There are two popular discretizations of total variation: the isotropic TV
[58] and the anisotropic TV [21], which are defined by

$\displaystyle\|\nabla
u\|_{2,1}=\sum_{i=1}^{M}\sum_{j=1}^{N}\sqrt{|(\nabla_{x}u)_{i,j}|^{2}+|(\nabla_{y}u)_{i,j}|^{2}},$
$\displaystyle\|\nabla
u\|_{1}=\sum_{i=1}^{M}\sum_{j=1}^{N}|(\nabla_{x}u)_{i,j}|+|(\nabla_{y}u)_{i,j}|,$

respectively. This work is based on the weighted difference between
anisotropic and isotropic TV (AITV) regularization [48], defined by

$\displaystyle\|\nabla u\|_{1}-\alpha\|\nabla
u\|_{2,1}=\sum_{i=1}^{M}\sum_{j=1}^{N}\left(|(\nabla_{x}u)_{i,j}|+|(\nabla_{y}u)_{i,j}|-\alpha\sqrt{|(\nabla_{x}u)_{i,j}|^{2}+|(\nabla_{y}u)_{i,j}|^{2}}\right),$
(5)

for a weighting parameter $\alpha\in[0,1].$ The range of $\alpha$ ensures the
non-negativity of the AITV regularization. Note that anisotropic TV is defined
as the $\ell_{1}$ norm of the image gradient
$((\nabla_{x}u)_{i,j},(\nabla_{y}u)_{i,j})$ at the pixel location
$(i,j)\in\Omega$, while isotropic TV is the $\ell_{2}$ norm on the gradient
vector. As a result, AITV can be viewed as the $\ell_{1}-\alpha\ell_{2}$
regularization on the gradient vector at every pixel, thereby enforcing
sparsity individually at each gradient vector.

###  2.3 Poisson Noise

Poisson noise follows the Poisson distribution with mean and variance $\eta$,
whose probability mass function is given by

$\displaystyle\mathbb{P}_{\eta}(n)=\frac{e^{-\eta}\eta^{n}}{n!},\;n\geq 0.$
(6)

For a clean image $g\in X$, its intensity value at each pixel $g_{i,j}$ serves
as the mean and variance for the corresponding noisy observation $f\in X$
defined by

$\displaystyle f_{i,j}\sim\text{Poisson}(g_{i,j})\;\forall(i,j)\in\Omega.$

To recover the image $g$ from the noisy image $f$, we find its maximum a
posteriori (MAP) estimation $u$, which maximizes the probability
$\mathbb{P}(u|f)$. By Bayes’ theorem, we have

$\displaystyle\mathbb{P}(u|f)=\frac{\mathbb{P}(f|u)\mathbb{P}(u)}{\mathbb{P}(f)}.$

It further follows from the definition (6) that

$\displaystyle\mathbb{P}(f_{i,j}|u_{i,j})\mathbb{P}(u_{i,j})=\mathbb{P}_{u_{i,j}}(f_{i,j})\mathbb{P}(u_{i,j})=\frac{e^{-u_{i,j}}u_{i,j}^{f_{i,j}}}{(f_{i,j})!}\mathbb{P}(u_{i,j}).$

Since Poisson noise is i.i.d. pixelwise, we have

$\displaystyle\mathbb{P}(u|f)=\prod_{(i,j)\in\Omega}\mathbb{P}(f_{i,j}|u_{i,j})\frac{\mathbb{P}(u_{i,j})}{\mathbb{P}(f_{i,j})}=\prod_{(i,j)\in\Omega}\frac{e^{-u_{i,j}}u_{i,j}^{f_{i,j}}}{(f_{i,j})!}\frac{\mathbb{P}(u_{i,j})}{\mathbb{P}(f_{i,j})}.$

The MAP estimate of $\mathbb{P}(u|f)$ is equivalent to its negative logarithm,
thus leading to the following optimization problem:

$\displaystyle\min_{u\geq 0}\sum_{(i,j)\in\Omega}u_{i,j}-f_{i,j}\log
u_{i,j}-\log\mathbb{P}(u_{i,j}).$ (7)

The last term $-\log\mathbb{P}(u_{i,j})$ can be regarded as an image prior or
a regularization. For example, Le et al. [37] considered the isotropic total
variation as the image prior and proposed a Poisson denoising model

$\displaystyle\min_{u\geq 0}\langle u-f\log u,\mathbbm{1}\rangle+\|\nabla
u\|_{2,1},$ (8)

where $\log$ is applied pixelwise and $\mathbbm{1}$ is the matrix whose
entries are all 1’s. The first term in (8) is a concise notation that is
commonly used as a fidelity term for Poisson denoising in various imaging
applications [16, 19, 22, 23, 37, 70].

###  2.4 Review of Poisson SaT/SLaT

A Poisson SaT framework [16] consists of two steps. Given a noisy grayscale
image $f\in X$ corrupted by Poisson noise, the first step is the smoothing
step that finds a piecewise-smooth solution $u^{*}$ from the optimization
model:

$\displaystyle u^{*}=\operatorname*{arg\,min}_{u\geq 0}\lambda\langle Au-f\log
Au,\mathbbm{1}\rangle+\frac{\mu}{2}\|\nabla u\|_{2}^{2}+\|\nabla u\|_{2,1}.$
(9)

Then in the thresholding step, $K-1$ threshold values
$\tau_{1}\leq\tau_{2}\leq\ldots\leq\tau_{K-1}$ are appropriately chosen to
segment $u^{*}$ into $K$ regions, where the $k$th region is given by

$\displaystyle\Omega_{k}=\\{(i,j)\in\Omega:\tau_{k-1}\leq
u^{*}_{i,j}<\tau_{k}\\},$

with $\tau_{0}\coloneqq\inf_{x\in\Omega}u^{*}(x)$. The thresholding step is
typically performed by $k$-means clustering.

The Poisson smoothing, lifting, and thresholding (SLaT) framework [10] extends
the Poisson SaT framework to color images. For a color image
$f=(f_{1},f_{2},f_{3})\in X\times X\times X$, the model (9) is applied to each
color channel $f_{i}$ for $i=1,2,3$, thus leading to a smoothed color image
$u^{*}=(u_{1}^{*},u_{2}^{*},u_{3}^{*})$. An additional lifting step [49] is
performed to transform $u^{*}$ to
$(u_{1}^{\prime},u_{2}^{\prime},u_{3}^{\prime})$ in the Lab space (perceived
lightness, red-green, and yellow-blue). The channels in Lab space are less
correlated than in RGB space, so they may have useful information for
segmentation. The RGB image and the Lab image are concatenated to form the
multichannel image
$\hat{u}\coloneqq(u_{1}^{*},u_{2}^{*},u_{3}^{*},u_{1}^{\prime},u_{2}^{\prime},u_{3}^{\prime})$,
followed by the thresholding stage. Generally, $k$-means clustering yields $K$
centroids $c_{1},\ldots,c_{K}$ as constant vectors, which are used to form the
region

$\displaystyle\Omega_{k}=\left\\{(i,j)\in\Omega:\|\hat{u}_{i,j}-c_{k}\|_{2}=\min_{1\leq\kappa\leq
K}\|\hat{u}_{i,j}-c_{\kappa}\|_{2}\right\\}$

for $k=1,\ldots,K$ such that $\Omega_{k}$’s are disjoint and
$\bigcup_{k=1}^{K}\Omega_{k}=\Omega$.

After the thresholding step for both SaT/SLaT, we define a piecewise-constant
approximation of the image $f$ by

$\displaystyle\tilde{f}=(\tilde{f}_{1},\ldots,\tilde{f}_{d})\text{ such that
}\tilde{f}_{\ell}=\sum_{k=1}^{K}c_{k,\ell}\mathbbm{1}_{\Omega_{k}}\;\forall\ell=1,\ldots,d,$
(10)

where $c_{k,\ell}$ is the $\ell$th entry of the constant vector $c_{k}$ and

$\displaystyle\mathbbm{1}_{\Omega_{k}}=\begin{cases}1&\text{ if
}(i,j)\in\Omega_{k},\\\ 0&\text{ if }(i,j)\not\in\Omega_{k}.\end{cases}$

Recall that $d=1$ when $f$ is grayscale, and $d=3$ when $f$ is color.

##  3 Proposed Approach

To improve the Poisson SaT/SLaT framework, we propose to replace the isotropic
TV in (9) with AITV regularization. In other words, in the smoothing step, we
obtain the smoothed image $u^{*}$ from the optimization problem

$\displaystyle u^{*}=\operatorname*{arg\,min}_{u}F(u)\coloneqq\lambda\langle
Au-f\log Au,\mathbbm{1}\rangle+\frac{\mu}{2}\|\nabla u\|_{2}^{2}+\|\nabla
u\|_{1}-\alpha\|\nabla u\|_{2,1},$ (11)

for $\alpha\in[0,1]$. We establish that this model admits a global solution.
We then develop an ADMM algorithm to find a solution and provide the
convergence analysis. The overall segmentation approach is described in
Algorithm 1.

1 Input:

  * •

image $f=(f_{1},\ldots,f_{d})$

  * •

blurring operator $A$

  * •

fidelity parameter $\lambda>0$

  * •

smoothing parameter $\mu\geq 0$

  * •

AITV parameter $\alpha\in[0,1]$

  * •

the number of regions in the image $K$

Output: Segmentation $\tilde{f}$

2 Stage one: Compute $u_{\ell}$ by solving (11) separately for
$\ell=1,\ldots,d$.

3 Stage two: if  _$f$ is a grayscale image, i.e., $d=1$_ then

4 Go to stage three.

5else if  _$f$ is a color image, i.e., $d=3$_ then

6 Transfer the solution $u^{*}=(u_{1}^{*},u_{2}^{*},u_{3}^{*})$ into Lab space
to obtain $(u_{1}^{\prime},u_{2}^{\prime},u_{3}^{\prime})$ and concatenate to
form
$\hat{u}=(u_{1}^{*},u_{2}^{*},u_{3}^{*},u_{1}^{\prime},u_{2}^{\prime},u_{3}^{\prime})$.

Stage three: Apply $k$-means to obtain $\\{(c_{k},\Omega_{k})\\}_{k=1}^{K}$
and compute $\tilde{f}$ by (10).

Algorithm 1 AITV Poisson SaT/SLaT

###  3.1 Model Analysis

To establish the solution’s existence of the proposed model (11), we start
with Lemma 1, a discrete version of Poincaré’s inequality [27]. In addition,
we prove Lemma 2 and Proposition 3, leading to the global existence theorem
(Theorem 4).

######  Lemma 1.

There exists a constant $C>0$ such that

$\displaystyle\left\|u-\bar{u}\mathbbm{1}\right\|_{2}\leq C\|\nabla
u\|_{2,1},$ (12)

for every $u\in X$ and
$\bar{u}\coloneqq\displaystyle\frac{1}{MN}\sum_{i=1}^{M}\sum_{j=1}^{N}u_{i,j}$.

###### Proof.

We prove it by contradiction. Suppose there exists a sequence
$\\{u_{k}\\}_{k=1}^{\infty}$ such that

$\displaystyle\left\|u_{k}-\bar{u}_{k}\mathbbm{1}\right\|_{2}>k\|\nabla
u_{k}\|_{2,1},$ (13)

where
$\bar{u}_{k}=\displaystyle\frac{1}{MN}\sum_{i=1}^{M}\sum_{j=1}^{N}(u_{k})_{i,j}$.
For every $k$, we normalize each element in the sequence by
$v_{k}=\frac{u_{k}-\bar{u}_{k}\mathbbm{1}}{\left\|u_{k}-\bar{u}_{k}\mathbbm{1}\right\|_{2}}.$
It is straightforward that

$\displaystyle\bar{v}_{k}=\frac{1}{MN}\sum_{i=1}^{M}\sum_{j=1}^{N}(v_{k})_{i,j}=0,\quad\|v_{k}\|_{2}=1\quad\forall
k\in\mathbb{N}.$ (14)

By (13), we have

$\displaystyle\|\nabla v_{k}\|_{2,1}<\frac{1}{k}.$ (15)

As $\\{v_{k}\\}_{k=1}^{\infty}$ is bounded, there exists a convergent
subsequence $\\{v_{k_{j}}\\}_{j=1}^{\infty}$ such that $v_{k_{j}}\rightarrow
v^{*}$ for $v^{*}\in X$. It follows from (15) that $\|\nabla v^{*}\|_{2,1}=0.$
Since $\text{ker}(\nabla)=\\{c\mathbbm{1}:c\in\mathbb{R}\\}$, then $v^{*}$ is
a constant vector. However, (14) implies that $\bar{v}^{*}=0$ and
$\|v^{*}\|_{2}=1$. This contradiction proves the lemma. ∎

######  Lemma 2.

Suppose $\|f\|_{\infty}<\infty$ and $\min_{i,j}f_{i,j}>0$. There exists a
scalar $u_{0}>0$ such that we have $2(x-f_{i,j}\log x)\geq x$ for any $\ x\geq
u_{0}$ and $(i,j)\in\Omega$.

###### Proof.

For each $(i,j)\in\Omega$, we want to show that there exists $u_{i,j}>0$ such
that $H(x)\coloneqq x-2f_{i,j}\log x\geq 0$ for $x\geq u_{i,j}$. Since $H(x)$
is strictly convex and it attains a global minimum at $x=2f_{i,j}$, it is
increasing on the domain $x>2f_{i,j}$. Additionally as $x$ dominates $\log(x)$
as $x\rightarrow+\infty$, there exists $u_{i,j}>2f_{i,j}>0$ such that
$\frac{u_{i,j}}{\log u_{i,j}}\geq 2f_{i,j},$ which implies that
$H(u_{i,j})=u_{i,j}-2f_{i,j}\log u_{i,j}\geq 0.$ As a result, for $x\geq
u_{i,j}>2f_{i,j}$, we obtain $x-2f_{i,j}\log x=H(x)\geq H(u_{i,j})\geq 0.$
Define $u_{0}\coloneqq\max_{i,j}u_{i,j}$, and hence we have $2(x-f_{i,j}\log
x)\geq x$ for $x\geq u_{0}\geq u_{i,j},\forall(i,j)\in\Omega$. ∎

######  Proposition 3.

Suppose $\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$ and
$\\{u_{k}\\}_{k=1}^{\infty}\subset X$. If $\\{(Au_{k},\nabla
u_{k})\\}_{k=1}^{\infty}$ is bounded, then $\\{u_{k}\\}_{k=1}^{\infty}$ is
bounded.

###### Proof.

Since $\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$, we have $A\mathbbm{1}\neq
0$. Simple calculations lead to

$\displaystyle\begin{aligned}
|\bar{u}_{k}|\|A\mathbbm{1}\|_{2}&=\|A(\bar{u}_{k}\mathbbm{1})\|_{2}\leq\|A(\bar{u}_{k}\mathbbm{1}-u_{k})\|_{2}+\|Au_{k}\|_{2}\\\
&\leq\|A\|\|u_{k}-\bar{u}_{k}\mathbbm{1}\|_{2}+\|Au_{k}\|_{2}\\\ &\leq
C\|A\|\|\nabla u_{k}\|_{2,1}+\|Au_{k}\|_{2},\end{aligned}$ (16)

where the last inequality is due to Lemma 1. The boundedness of
$\\{Au_{k}\\}_{k=1}^{\infty}$ and $\\{\nabla u_{k}\\}_{k=1}^{\infty}$ implies
that $\\{\bar{u}_{k}\\}_{k=1}^{\infty}$ is also bounded by (16). We apply
Lemma 1 to obtain

$\displaystyle\|u_{k}\|_{2}\leq\|u_{k}-\bar{u}_{k}\mathbbm{1}\|_{2}+\|\bar{u}_{k}\mathbbm{1}\|_{2}<C\|\nabla
u_{k}\|_{2,1}+\|\bar{u}_{k}\mathbbm{1}\|_{2}<\infty,$

which thereby proves that $\\{u_{k}\\}_{k=1}^{\infty}$ is bounded. ∎

Finally, we adapt the proof in [16] to establish that $F$ has a global
minimizer.

######  Theorem 4.

Suppose $\|f\|_{\infty}<\infty$ and $\min_{i,j}f_{i,j}>0$. If
$\lambda>0,\mu\geq 0,\alpha\in[0,1)$, and
$\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$, then $F$ has a global
minimizer.

###### Proof.

It is straightforward that $\|\nabla u\|_{2,1}\leq\|\nabla u\|_{1}$, thus
$\|\nabla u\|_{1}-\alpha\|\nabla u\|_{2,1}\geq 0$ for $\alpha\in[0,1)$. As a
result, we have

$\displaystyle F(u)\geq\lambda\langle Au-f\log
Au,\mathbbm{1}\rangle=\lambda\sum_{i=1}^{M}\sum_{j=1}^{N}(Au)_{i,j}-f_{i,j}\log(Au)_{i,j}.$

Given a scalar $f>0,$ the function $G(x)=x-f\log(x)$ attains its global
minimum at $x=f.$ Therefore, we have $x-f_{i,j}\log x\geq f_{i,j}-f_{i,j}\log
f_{i,j}$ for all $x>0$ and $(i,j)\in\Omega$, which leads to a lower bound of
$F(u)$, i.e.,

$\displaystyle
F(u)\geq\lambda\sum_{i=1}^{M}\sum_{j=1}^{N}(Au)_{i,j}-f_{i,j}\log(Au)_{i,j}\geq\lambda\sum_{i=1}^{M}\sum_{j=1}^{N}f_{i,j}-f_{i,j}\log
f_{i,j}\eqqcolon F_{0}.$ (17)

As $F(u)$ is lower bounded by $F_{0}$, we can choose a minimizing sequence
$\\{u_{k}\\}_{k=1}^{\infty}$ and hence $F(u_{k})$ has a uniform upper bound,
denoted by $B_{1}$, i.e., $F(u_{k})<B_{1}$ for all $k\in\mathbb{N}$. It
further follows from (17) that

$\displaystyle B_{1}\geq F(u_{k})\geq\lambda\langle Au_{k}-f\log
Au_{k},\mathbbm{1}\rangle\geq F_{0},$

which implies that $\\{|\langle Au_{k}-f\log
Au_{k},\mathbbm{1}\rangle|\\}_{k=1}^{\infty}$ is uniformly bounded, i.e.,
there exists a constant $B_{2}>0$ such that $\left|\langle Au_{k}-f\log
Au_{k},\mathbbm{1}\rangle\right|<B_{2},\ \forall k$. Using these uniform
bounds, we derive that

$\displaystyle(1-\alpha)\|\nabla u_{k}\|_{1}\leq\frac{\mu}{2}\|\nabla
u_{k}\|_{2}^{2}+\|\nabla u_{k}\|_{1}-\alpha\|\nabla
u_{k}\|_{2,1}=F(u_{k})-\lambda\langle Au_{k}-f\log
Au_{k},\mathbbm{1}\rangle\leq B_{1}+\lambda B_{2}.$

As $\alpha<1$, the sequence $\\{\nabla u_{k}\\}_{k=1}^{\infty}$ is bounded.

To prove the boundedness of $\\{Au_{k}\\}_{k=1}^{\infty}$, we introduce the
notations of $x^{+}=\max(x,0)$ and $x^{-}=-\min(x,0)$ for any
$x\in\mathbb{R}$. Then $x=x^{+}-x^{-}$. By Lemma 2, there exists $u_{0}>0$
such that $2(x-f_{i,j}\log x)\geq x,\ \forall x\geq u_{0}$ and
$(i,j)\in\Omega$. We observe that

$\displaystyle\begin{aligned}
\|Au_{k}\|_{1}&=\sum_{i=1}^{M}\sum_{j=1}^{N}|(Au_{k})_{i,j}|\leq\sum_{i=1}^{M}\sum_{j=1}^{N}\max\\{2((Au_{k})_{i,j}-f_{i,j}\log(Au_{k})_{i,j}),u_{0}\\}\\\
&\leq
2\sum_{i=1}^{M}\sum_{j=1}^{N}\left((Au_{k})_{i,j}-f_{i,j}\log(Au_{k})_{i,j}\right)^{+}+MNu_{0}\\\
&=2\sum_{i=1}^{M}\sum_{j=1}^{N}\left[\left((Au_{k})_{i,j}-f_{i,j}\log(Au_{k})_{i,j}\right)+\left((Au_{k})_{i,j}-f_{i,j}\log(Au_{k})_{i,j}\right)^{-}\right]+MNu_{0}\\\
&=2\langle Au_{k}-f\log
Au_{k},\mathbbm{1}\rangle+2\sum_{i=1}^{M}\sum_{j=1}^{N}\left((Au_{k})_{i,j}-f_{i,j}\log(Au_{k})_{i,j}\right)^{-}+MNu_{0}\\\
&\leq 2B_{2}+2\sum_{i=1}^{M}\sum_{j=1}^{N}\left|f_{i,j}-f_{i,j}\log
f_{i,j}\right|+MNu_{0}<\infty.\end{aligned}$ (18)

This shows that $\\{Au_{k}\\}_{k=1}^{\infty}$ is bounded.

Since both $\\{\nabla u_{k}\\}_{k=1}^{\infty}$ and
$\\{Au_{k}\\}_{k=1}^{\infty}$ are bounded, then $\\{u_{k}\\}_{k=1}^{\infty}$
is bounded due to Proposition 3. Therefore, there exists a subsequence
$\\{u_{k_{n}}\\}_{n=1}^{\infty}$ that converges to some $u^{*}\in X$. As $F$
is continuous and thus lower semicontinuous, we have

$\displaystyle F(u^{*})\leq\liminf_{n\rightarrow\infty}F(u_{k_{n}}),$

which means that $u^{*}$ minimizes $F$. ∎

###  3.2 Numerical Algorithm

To minimize (11), we introduce two auxiliary variables $v\in X$ and
$w=(w_{x},w_{y})\in Y$, leading to an equivalent constrained optimization
problem:

$\displaystyle\min_{u,v,w}$ $\displaystyle\quad\lambda\langle v-f\log
v,\mathbbm{1}\rangle+\frac{\mu}{2}\|\nabla
u\|_{2}^{2}+\|w\|_{1}-\alpha\|w\|_{2,1}$ (19) s.t. $\displaystyle\quad
Au=v,\quad\nabla u=w.$

The corresponding augmented Lagrangian is expressed as

$\displaystyle\begin{aligned}
\mathcal{L}_{\beta_{1},\beta_{2}}(u,v,w,y,z)=&\lambda\langle v-f\log
v,\mathbbm{1}\rangle+\frac{\mu}{2}\|\nabla
u\|_{2}^{2}+\|w\|_{1}-\alpha\|w\|_{2,1}\\\ &+\langle
y,Au-v\rangle+\frac{\beta_{1}}{2}\|Au-v\|_{2}^{2}+\langle z,\nabla
u-w\rangle+\frac{\beta_{2}}{2}\|\nabla u-w\|_{2}^{2},\end{aligned}$ (20)

where $y\in X$ and $z=(z_{x},z_{y})\in Y$ are Lagrange multipliers and
$\beta_{1},\beta_{2}$ are positive parameters. We then apply the alternating
direction method of multipliers (ADMM) to minimize (19) that consists of the
following steps per iteration $k$:

$\displaystyle u_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{u}\mathcal{L}_{\beta_{1,k},\beta_{2,k}}(u,v_{k},w_{k},y_{k},z_{k}),$
(21a) $\displaystyle v_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{v}\mathcal{L}_{\beta_{1,k},\beta_{2,k}}(u_{k+1},v,w_{k},y_{k},z_{k}),$
(21b) $\displaystyle w_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{w}\mathcal{L}_{\beta_{1,k},\beta_{2,k}}(u_{k+1},v_{k+1},w,y_{k},z_{k}),$
(21c) $\displaystyle y_{k+1}$
$\displaystyle=y_{k}+\beta_{1,k}(Au_{k+1}-v_{k+1}),$ (21d) $\displaystyle
z_{k+1}$ $\displaystyle=z_{k}+\beta_{2,k}(\nabla u_{k+1}-w_{k+1}),$ (21e)
$\displaystyle(\beta_{1,k+1},\beta_{2,k+1})$
$\displaystyle=\sigma(\beta_{1,k},\beta_{2,k}),$ (21f)

where $\sigma>1$.

######  Remark 1.

The scheme presented in (21) slightly differs from the original ADMM [6], the
latter of which has $\sigma=1$ in (21f). Having $\sigma>1$ increases the
weights of the penalty parameters $\beta_{1,k},\beta_{2,k}$ in each iteration
$k$, thus accelerating the numerical convergence speed of the proposed ADMM
algorithm. A similar technique has been used in [14, 30, 60, 61, 79].

All the subproblems (21a)-(21c) have closed-form solutions. In particular, the
first-order optimality condition for (21a) is

$\displaystyle[\beta_{1,k}A^{\top}A-(\mu+\beta_{2,k})\Delta]u^{k+1}=A^{\top}(\beta_{1,k}v_{k}-y_{k})-\nabla^{\top}(z_{k}-\beta_{2,k}w_{k}),$
(22)

where $\Delta=-\nabla^{\top}\nabla$ is the Laplacian operator. If
$\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$, then
$\beta_{1,k}A^{\top}A-(\mu+\beta_{2,k})\Delta$ is positive definite and
thereby invertible, which implies that (22) has a unique solution $u^{k+1}$.
By assuming periodic boundary condition for $u$, the operators $\Delta$ and
$A^{\top}A$ are block circulant [68], and hence (22) can be solved efficiently
by the 2D discrete Fourier transform $\mathcal{F}$. Specifically, we have the
formula

$\displaystyle
u_{k+1}=\mathcal{F}^{-1}\left(\frac{\mathcal{F}(A)^{*}\circ\mathcal{F}(\beta_{1,k}v_{k}-y_{k})-\mathcal{F}(\nabla)^{*}\circ\mathcal{F}(z_{k}-\beta_{2,k}w_{k})}{\beta_{1,k}\mathcal{F}(A)^{*}\circ\mathcal{F}(A)-(\mu+\beta_{2,k})\mathcal{F}(\Delta)}\right),$
(23)

where $\mathcal{F}^{-1}$ is the inverse discrete Fourier transform, the
superscript $*$ denotes complex conjugate, the operation $\circ$ is
componentwise multiplication, and division is componentwise. By
differentiating the objective function of (21b) and setting it to zero, we can
get a closed-form solution for $v_{k+1}$ given by

$\displaystyle
v_{k+1}=\frac{\left(\beta_{1,k}Au_{k+1}+y_{k}-\lambda\mathbbm{1}\right)+\sqrt{\left(\beta_{1,k}Au_{k+1}+y_{k}-\lambda\mathbbm{1}\right)^{2}+4\lambda\beta_{1,k}f}}{2\beta_{1,k}},$
(24)

where the square root, squaring, and division are performed componentwise.
Lastly, the $w$-subproblem (21c) can be decomposed componentwise as follows:

$\displaystyle\begin{aligned}
(w_{i,j})_{k+1}&=\operatorname*{arg\,min}_{w_{i,j}}\|w_{i,j}\|_{1}-\alpha\|w_{i,j}\|_{2}+\frac{\beta_{2,k}}{2}\left\|w_{i,j}-\left((\nabla
u_{k+1})_{i,j}+\frac{(z_{k})_{i,j}}{\beta_{2,k}}\right)\right\|_{2}^{2}\\\
&=\text{prox}\left((\nabla
u_{k+1})_{i,j}+\frac{(z_{k})_{i,j}}{\beta_{2,k}},\alpha,\frac{1}{\beta_{2,k}}\right),\end{aligned}$
(25)

where the proximal operator for $\ell_{1}-\alpha\ell_{2}$ on
$x\in\mathbb{R}^{n}$ is given by

$\displaystyle\text{prox}(x,\alpha,\beta)=\operatorname*{arg\,min}_{y}\|y\|_{1}-\alpha\|y\|_{2}+\frac{1}{2\beta}\|x-y\|_{2}^{2}.$
(26)

The proximal operator for $\ell_{1}-\alpha\ell_{2}$ has a closed form solution
summarized by Lemma 5.

######  Lemma 5 ([46]).

Given $x\in\mathbb{R}^{n}$, $\beta>0$, and $\alpha\in[0,1]$, the optimal
solution to (26) is given by one of the following cases:

  1. 1.

When $\|x\|_{\infty}>\beta$, we have

$\displaystyle x^{*}=(\|\xi\|_{2}+\alpha\beta)\frac{\xi}{\|\xi\|_{2}},$

where $\xi=\operatorname{sign}(x)\circ\max(|x|-\beta,0)$.

  2. 2.

When $(1-\alpha)\beta<\|x\|_{\infty}\leq\beta$, then $x^{*}$ is a 1-sparse
vector such that one chooses
$i\in\displaystyle\operatorname*{arg\,max}_{j}(|x_{j}|)$ and defines
$x^{*}_{i}=\left(|x_{i}|+(\alpha-1)\beta\right)\operatorname{sign}(x_{i})$ and
the remaining elements equal to 0.

  3. 3.

When $\|x\|_{\infty}\leq(1-\alpha)\beta$, then $x^{*}=0$.

In summary, we describe the ADMM scheme to solve (11) in Algorithm 2.

1 Input:  • image $f$ • blurring operator $A$ • fidelity parameter $\lambda>0$
• smoothing parameter $\mu\geq 0$ • AITV parameter $\alpha\in[0,1]$ • penalty
parameters $\beta_{1,0},\beta_{2,0}>0$ • penalty multiplier $\sigma>1$ •
relative error $\epsilon>0$ Output:$u_{k}$

2 Initialize $u_{0},w_{0},z_{0}$.

3 Set $k=0$.

4 while  _$\frac{\|u_{k}-u_{k-1}\|_{2}}{\|u_{k}\|_{2}} >\epsilon$_ do

5 $\displaystyle u_{k+1}$
$\displaystyle=\mathcal{F}^{-1}\left(\frac{\mathcal{F}(A)^{*}\circ\mathcal{F}(\beta_{1,k}v_{k}-y_{k})-\mathcal{F}(\nabla)^{*}\circ\mathcal{F}(z_{k}-\beta_{2,k}w_{k})}{\beta_{1,k}\mathcal{F}(A)^{*}\circ\mathcal{F}(A)-(\mu+\beta_{2,k})\mathcal{F}(\Delta)}\right)$
$\displaystyle v_{k+1}$
$\displaystyle=\frac{\left(\beta_{1,k}Au_{k+1}+y_{k}-\lambda\mathbbm{1}\right)+\sqrt{\left(\beta_{1,k}Au_{k+1}+y_{k}-\lambda\mathbbm{1}\right)^{2}+4\lambda\beta_{1,k}f}}{2\beta_{1,k}}$
$\displaystyle(w_{k+1})_{i,j}$ $\displaystyle=\text{prox}\left((\nabla
u_{k+1})_{i,j}+\frac{(z_{k})_{i,j}}{\beta_{2,k}},\alpha,\frac{1}{\beta_{2,k}}\right)\quad\forall(i,j)\in\Omega$
$\displaystyle y_{k+1}$ $\displaystyle=y_{k}+\beta_{1,k}(Au_{k+1}-v_{k+1})$
$\displaystyle z_{k+1}$ $\displaystyle=z_{k}+\beta_{2,k}(\nabla
u_{k+1}-w_{k+1})$ $\displaystyle(\beta_{1,k+1},\beta_{2,k+1})$
$\displaystyle=\sigma(\beta_{1,k},\beta_{2,k})$ $\displaystyle k$
$\displaystyle\coloneqq k+1$

6

7

Algorithm 2 ADMM for the AITV-Regularized Smoothing Model with Poisson
Fidelity (11)

###  3.3 Convergence Analysis

We establish the subsequential convergence of ADMM described in Algorithm 2.
The global convergence of ADMM [69] is inapplicable to our model as the
gradient operator $\nabla$ is non-surjective, which will be further
investigated in future work. For the sake of brevity, we set
$\beta=\beta_{1}=\beta_{2}$ and denote

$\displaystyle\mathcal{L}_{\beta}(u,v,w,y,z)\coloneqq\mathcal{L}_{\beta,\beta}(u,v,w,y,z).$

In addition, we introduce definitions of subdifferentials [57], which defines
a stationary point of a non-smooth objective function.

######  Definition 6.

For a proper function
$h:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\\{+\infty\\}$, define
$\text{dom}(h)\coloneqq\\{x\in\mathbb{R}^{n}:h(x)<+\infty\\}$.

  1. (a)

The regular subdifferential at $x\in\text{dom}(h)$ is given by

$\displaystyle\hat{\partial}{h}(x)\coloneqq\left\\{w:\liminf_{x^{\prime}\rightarrow
x,x^{\prime}\neq x}\frac{h(x^{\prime})-h(x)-\langle
w,x^{\prime}-x\rangle}{\|x^{\prime}-x\|}\geq 0\right\\}.$

  2. (b)

The (limiting) subdifferential at $x\in\text{dom}(h)$ is given by

$\displaystyle\partial{h}(x)\coloneqq\left\\{w:\exists\,x_{k}\rightarrow
x\text{ and }w_{k}\in\hat{\partial}{h}(x_{k})\text{ with }w_{k}\rightarrow
w\text{ and }h(x_{k})\rightarrow h(x)\right\\}.$

An important property of the limiting subdifferential is its closedness: for
any $(x_{k},v_{k})\rightarrow(x,v)$ with $v_{k}\in\partial{h}(x_{k})$, if
$h(x_{k})\rightarrow h(x)$, then $v\in\partial{h}(x)$.

######  Lemma 7.

Suppose that $\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$ and
$0\leq\alpha<1$. Let $\\{(u_{k},v_{k},w_{k},y_{k},z_{k})\\}_{k=1}^{\infty}$ be
a sequence generated by Algorithm 2. Then we have

$\displaystyle\begin{aligned}
&\mathcal{L}_{\beta_{k+1}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})-\mathcal{L}_{\beta_{k}}(u_{k},v_{k},w_{k},y_{k},z_{k})\\\
&\leq-\frac{\nu}{2}\|u_{k+1}-u_{k}\|_{2}^{2}-\frac{\beta_{0}}{2}\|v_{k+1}-v_{k}\|_{2}^{2}+\frac{1}{\sigma^{k-1}\beta_{0}}\left(\left\|y_{k+1}-y_{k}\right\|_{2}^{2}+\left\|z_{k+1}-z_{k}\right\|_{2}^{2}\right),\end{aligned}$
(28)

for some constant $\nu>0$.

###### Proof.

If $\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$, then
$\beta_{0}A^{\top}A+(\beta_{0}+\mu)\nabla^{\top}\nabla$ is positive definite,
and hence there exists $\nu>0$ such that

$\displaystyle\beta_{k}\|Au\|_{2}^{2}+(\beta_{k}+\mu)\|\nabla
u\|_{2}^{2}\geq\beta_{0}\|Au\|_{2}^{2}+(\beta_{0}+\mu)\|\nabla
u\|_{2}^{2}\geq\nu\|u\|_{2}^{2}\;\quad\forall k\in\mathbb{N},$

which implies that $\mathcal{L}_{\beta_{k}}(u,v_{k},w_{k},y_{k},z_{k})$ is
strongly convex with respect to $u$ with parameter $\nu$. Additionally,
$\mathcal{L}_{\beta_{k}}(u_{k+1},v,w_{k},y_{k},z_{k})$ is strongly convex with
respect to $v$ with parameter $\beta_{0}\leq\beta_{k}$. It follows from [5,
Theorem 5.25] that we have

$\displaystyle\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k},w_{k},y_{k},z_{k})-\mathcal{L}_{\beta_{k}}(u_{k},v_{k},w_{k},y_{k},z_{k})$
$\displaystyle\leq-\frac{\nu}{2}\|u_{k+1}-u_{k}\|_{2}^{2},$ (29)
$\displaystyle\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k},y_{k},z_{k})-\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k},w_{k},y_{k},z_{k})$
$\displaystyle\leq-\frac{\beta_{0}}{2}\|v_{k+1}-v_{k}\|_{2}^{2}.$ (30)

As $w_{k+1}$ is the optimal solution to (21c), it is straightforward to have

$\displaystyle\begin{aligned}
&\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k+1},y_{k},z_{k})-\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k},y_{k},z_{k})\leq
0.\end{aligned}$ (31)

Simple calculations by using (21d)-(21e) lead to

$\displaystyle\begin{aligned}
\mathcal{L}_{\beta_{k}}(u_{k+1},&v_{k+1},w_{k+1},y_{k+1},z_{k+1})-\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k+1},y_{k},z_{k})\\\
=&\left(\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})-\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k})\right)\\\
&+\left(\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k})-\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k+1},y_{k},z_{k})\right)\\\
=&\langle z_{k+1}-z_{k},\nabla u_{k+1}-w_{k+1}\rangle+\langle
y_{k+1}-y_{k},Au_{k+1}-v_{k+1}\rangle\\\
=&\frac{1}{\beta_{k}}\left(\|y_{k+1}-y_{k}\|_{2}^{2}+\|z_{k+1}-z_{k}\|_{2}^{2}\right).\end{aligned}$
(32)

Lastly, we have

$\displaystyle\begin{aligned}
\mathcal{L}_{\beta_{k+1}}(u_{k+1},&v_{k+1},w_{k+1},y_{k+1},z_{k+1})-\mathcal{L}_{\beta_{k}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})\\\
=&\frac{\beta_{k+1}-\beta_{k}}{2}\left(\|Au_{k+1}-v_{k+1}\|_{2}^{2}+\|\nabla
u_{k+1}-w_{k+1}\|_{2}^{2}\right)\\\
=&\frac{\beta_{k+1}-\beta_{k}}{2\beta_{k}^{2}}\left(\left\|y_{k+1}-y_{k}\right\|_{2}^{2}+\left\|z_{k+1}-z_{k}\right\|_{2}^{2}\right).\end{aligned}$
(33)

Combining (29) - (33) together with the fact that
$\beta_{k}=\sigma^{k}\beta_{0}$ for $\sigma>1$, we obtain

$\displaystyle\mathcal{L}_{\beta_{k+1}}(u_{k+1},$ $\displaystyle
v_{k+1},w_{k+1},y_{k+1},z_{k+1})-\mathcal{L}_{\beta_{k}}(u_{k},v_{k},w_{k},y_{k},z_{k})$
$\displaystyle\leq$
$\displaystyle-\frac{\nu}{2}\|u_{k+1}-u_{k}\|_{2}^{2}-\frac{\beta_{0}}{2}\|v_{k+1}-v_{k}\|_{2}^{2}+\frac{\beta_{k+1}+\beta_{k}}{2\beta_{k}^{2}}\left(\left\|y_{k+1}-y_{k}\right\|_{2}^{2}+\left\|z_{k+1}-z_{k}\right\|_{2}^{2}\right)$
$\displaystyle=$
$\displaystyle-\frac{\nu}{2}\|u_{k+1}-u_{k}\|_{2}^{2}-\frac{\beta_{0}}{2}\|v_{k+1}-v_{k}\|_{2}^{2}+\frac{\sigma+1}{2\sigma^{k}\beta_{0}}\left(\left\|y_{k+1}-y_{k}\right\|_{2}^{2}+\left\|z_{k+1}-z_{k}\right\|_{2}^{2}\right)$
$\displaystyle\leq$
$\displaystyle-\frac{\nu}{2}\|u_{k+1}-u_{k}\|_{2}^{2}-\frac{\beta_{0}}{2}\|v_{k+1}-v_{k}\|_{2}^{2}+\frac{1}{\sigma^{k-1}\beta_{0}}\left(\left\|y_{k+1}-y_{k}\right\|_{2}^{2}+\left\|z_{k+1}-z_{k}\right\|_{2}^{2}\right).$

This completes the proof. ∎

######  Lemma 8.

Suppose that $\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$ and
$0\leq\alpha<1$. Let $\\{(u_{k},v_{k},w_{k},y_{k},z_{k})\\}_{k=1}^{\infty}$ be
generated by Algorithm 2. If $\\{y_{k}\\}_{k=1}^{\infty}$ bounded, then the
sequence $\\{(u_{k},v_{k},w_{k},y_{k},z_{k})\\}_{k=1}^{\infty}$ is bounded,
$u_{k+1}-u_{k}\rightarrow 0,$ and $v_{k+1}-v_{k}\rightarrow 0$.

###### Proof.

First we show that $\\{z_{k}\\}_{k=1}^{\infty}$ is bounded. Combining (21e)
with the first-order optimality condition of (25), we have

$\displaystyle\begin{aligned}
(z_{k+1})_{i,j}=(z_{k})_{i,j}+\beta_{k}\left((\nabla
u_{k+1})_{i,j}-(w_{k+1})_{i,j}\right)&\in\partial\left(\|(w_{k+1})_{i,j}\|_{1}-\alpha\|(w_{k+1})_{i,j}\|_{2}\right)\\\
&\subseteq\partial\left(\|(w_{k+1})_{i,j}\|_{1}\right)-\alpha\partial\left(\|(w_{k+1})_{i,j}\|_{2}\right),\end{aligned}$
(34)

which implies that there exist $\xi_{1}\in\partial\|(w_{k+1})_{i,j}\|_{1}$ and
$\xi_{2}\in\partial\|(w_{k+1})_{i,j}\|_{2}$ such that
$(z_{k+1})_{i,j}=\xi_{1}-\alpha\xi_{2}$ for each $(i,j)\in\Omega$. Recall that
for $x\in\mathbb{R}^{2}$ the subgradients of the two norms are

$\displaystyle\partial\|x\|_{1}=\left\\{\xi\in\mathbb{R}^{2}:\xi_{i}=\begin{cases}\text{sign}(x_{i})&\text{
if }x_{i}\neq 0\\\ \xi_{i}\in[-1,1]&\text{ if }x_{i}=0\end{cases}\;\text{ for
}i=1,2\right\\},$ (35)
$\displaystyle\partial\|x\|_{2}=\left\\{\xi\in\mathbb{R}^{2}:\xi=\begin{cases}\frac{x}{\|x\|_{2}}&\text{
if }x\neq 0\\\ \in\\{\xi\in\mathbb{R}^{2}:\|\xi\|_{2}\leq 1\\}&\text{ if
}x=0\end{cases}\right\\}.$ (36)

Therefore, we have $\|\xi_{1}\|_{\infty}\leq 1,\|\xi_{2}\|_{\infty}\leq 1$,
and hence $\|(z_{k+1})_{i,j}\|_{\infty}\leq 1+\alpha$ (by the triangle
inequality), i.e., $\\{z_{k}\\}_{k=1}^{\infty}$ is bounded.

By the assumption $\\{(y_{k})\\}_{k=1}^{\infty}$ is bounded. There exist two
constants $C_{1},C_{2}>0$ such that $\|y_{k+1}-y_{k}\|_{2}^{2}\leq C_{1}$,
$\|z_{k+1}-z_{k}\|_{2}^{2}\leq C_{1}$, $\|y_{k}\|_{2}^{2}\leq C_{2}$, and
$\|z_{k}\|_{2}^{2}\leq C_{2}$ for all $k\in\mathbb{N}$. Hence, we have from
(LABEL:eq:lagrange_ineq) that

$\displaystyle\begin{aligned}
\mathcal{L}_{\beta_{k+1}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})\leq&\mathcal{L}_{\beta_{k}}(u_{k},v_{k},w_{k},y_{k},z_{k})-\frac{\nu}{2}\|u_{k+1}-u_{k}\|_{2}^{2}\\\
&-\frac{\beta_{0}}{2}\|v_{k+1}-v_{k}\|_{2}^{2}+\frac{2C_{1}}{\sigma^{k-1}\beta_{0}}.\end{aligned}$
(37)

A telescoping summation of (37) leads to

$\displaystyle\begin{aligned}
\mathcal{L}_{\beta_{k+1}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})\leq&\mathcal{L}_{\beta_{0}}(u_{0},v_{0},w_{0},y_{0},z_{0})+\frac{2C_{1}}{\beta_{0}}\sum_{i=0}^{k}\frac{1}{\sigma^{i-1}}\\\
&-\frac{\nu}{2}\sum_{i=0}^{k}\|u_{i+1}-u_{i}\|_{2}^{2}-\frac{\beta_{0}}{2}\sum_{i=0}^{k}\|v_{i+1}-v_{i}\|_{2}^{2}.\end{aligned}$
(38)

By completing two least-squares terms, we can rewrite
$\mathcal{L}_{\beta_{k+1}}$ as

$\displaystyle\begin{aligned}
\mathcal{L}_{\beta_{k+1}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})=&\lambda\langle
v_{k+1}-f\log v_{k+1},\mathbbm{1}\rangle+\frac{\mu}{2}\|\nabla
u_{k+1}\|_{2}^{2}+\|w_{k+1}\|_{1}-\alpha\|w_{k+1}\|_{2,1}\\\
&+\frac{\beta_{k+1}}{2}\left\|Au_{k+1}-v_{k+1}+\frac{y_{k+1}}{\beta_{k+1}}\right\|_{2}^{2}-\frac{\|y_{k+1}\|_{2}^{2}}{2\beta_{k+1}}\\\
&+\frac{\beta_{k+1}}{2}\left\|\nabla
u_{k+1}-w_{k+1}+\frac{z_{k+1}}{\beta_{k+1}}\right\|_{2}^{2}-\frac{\|z_{k+1}\|_{2}^{2}}{2\beta_{k+1}}.\end{aligned}$
(39)

Combining (38) and (39), we have

$\displaystyle\begin{aligned} \lambda\langle f-f\log
f,\mathbbm{1}\rangle+(1-\alpha)\|w_{k+1}\|_{1}-\frac{C_{2}}{\beta_{0}}\leq&\mathcal{L}_{\beta_{k+1}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})\\\
\leq&\mathcal{L}_{\beta_{0}}(u_{0},v_{0},w_{0},y_{0},z_{0})+\frac{2C_{1}}{\beta_{0}}\sum_{i=0}^{k}\frac{1}{\sigma^{i-1}}\\\
&-\frac{\nu}{2}\sum_{i=0}^{k}\|u_{i+1}-u_{i}\|_{2}^{2}-\frac{\beta_{0}}{2}\sum_{i=0}^{k}\|v_{i+1}-v_{i}\|_{2}^{2}\\\
\leq&\mathcal{L}_{\beta_{0}}(u_{0},v_{0},w_{0},y_{0},z_{0})+\frac{2C_{1}}{\beta_{0}}\sum_{i=0}^{\infty}\frac{1}{\sigma^{i-1}}.\end{aligned}$
(40)

Since $\sigma>1$, the infinite sum is finite, and hence we have $\forall
k\in\mathbb{N},$

$\displaystyle\|w_{k+1}\|_{1}\leq\frac{1}{1-\alpha}\left(\mathcal{L}_{\beta_{0}}(u_{0},v_{0},w_{0},y_{0},z_{0})-\lambda\langle
f-f\log
f,\mathbbm{1}\rangle+\frac{2C_{1}}{\beta_{0}}\sum_{i=0}^{\infty}\frac{1}{\sigma^{i-1}}+\frac{C_{2}}{\beta_{0}}\right)<\infty,$

which implies that $\\{w_{k}\\}_{k=1}^{\infty}$ is bounded. Also from (38) and
(39), we have

$\displaystyle\lambda\langle f-f\log
f,\mathbbm{1}\rangle-\frac{C_{2}}{\beta_{0}}$ $\displaystyle\leq\lambda\langle
v_{k+1}-f\log v_{k+1},\mathbbm{1}\rangle-\frac{C_{2}}{\beta_{0}}$
$\displaystyle\leq\lambda\langle v_{k+1}-f\log
v_{k+1},\mathbbm{1}\rangle-\frac{\|y_{k+1}\|_{2}^{2}}{2\beta_{k+1}}-\frac{\|z_{k+1}\|_{2}^{2}}{2\beta_{k+1}}$
$\displaystyle\leq\mathcal{L}_{\beta_{k+1}}(u_{k+1},v_{k+1},w_{k+1},y_{k+1},z_{k+1})$
$\displaystyle\leq\mathcal{L}_{\beta_{0}}(u_{0},v_{0},w_{0},y_{0},z_{0})+\frac{2C_{1}}{\beta_{0}}\sum_{i=0}^{\infty}\frac{1}{\sigma^{i-1}}<\infty.$

This shows that $\\{\langle v_{k}-f\log
v_{k},\mathbbm{1}\rangle\\}_{k=1}^{\infty}$ is bounded. By emulating the
computation in (18), it can be shown that $\\{v_{k}\\}_{k=1}^{\infty}$ is
bounded.

It suffices to prove that $\\{(Au_{k},\nabla u_{k})\\}_{k=1}^{\infty}$ is
bounded in order to prove the boundedness of $\\{u_{k}\\}_{k=1}^{\infty}$ by
Proposition 3. Using (21d), we have

$\displaystyle\|Au_{k+1}\|_{2}\leq\frac{\|y_{k+1}-y_{k}\|_{2}}{\beta_{k}}+\|v_{k+1}\|_{2}\leq\frac{\sqrt{C_{1}}}{\beta_{0}}+\|v_{k+1}\|_{2}.$

As $\\{v_{k}\\}_{k=1}^{\infty}$ is proven to be bounded, then
$\\{Au_{k}\\}_{k=1}^{\infty}$ is also bounded. We can prove $\\{\nabla
u_{k}\\}_{k=1}^{\infty}$ is bounded similarly using (21e). Altogether,
$\\{(u_{k},v_{k},w_{k},y_{k},z_{k})\\}_{k=1}^{\infty}$ is bounded.

It follows from (40) that

$\displaystyle\frac{\nu}{2}\sum_{i=0}^{k}\|u_{i+1}-u_{i}\|_{2}^{2}+\frac{\beta_{0}}{2}\sum_{i=0}^{k}\|v_{i+1}-v_{i}\|_{2}^{2}\leq$
$\displaystyle\mathcal{L}_{\beta_{0}}(u_{0},v_{0},w_{0},y_{0},z_{0})$
$\displaystyle+\frac{2C_{1}}{\beta_{0}}\sum_{i=0}^{k}\frac{1}{\sigma^{i-1}}-\lambda\langle
f-f\log f,\mathbbm{1}\rangle+\frac{C_{2}}{\beta_{0}}.$

As $k\rightarrow\infty$, we see the right-hand side is finite, which forces
the infinite summations on the left-hand side to converge, and hence we have
$u_{k+1}-u_{k}\rightarrow 0$ and $v_{k+1}-v_{k}\rightarrow 0$. ∎

######  Theorem 9.

Suppose that $\text{ker}(A)\cap\text{ker}(\nabla)=\\{0\\}$ and
$0\leq\alpha<1$. Let $\\{(u_{k},v_{k},w_{k},y_{k},z_{k})\\}_{k=1}^{\infty}$ be
generated by Algorithm 2. If $\\{y_{k}\\}_{k=1}^{\infty}$ bounded,
$\beta_{k}(v_{k+1}-v_{k})\rightarrow 0,\beta_{k}(w_{k+1}-w_{k})\rightarrow
0,y_{k+1}-y_{k}\rightarrow 0$, and $z_{k+1}-z_{k}\rightarrow 0$, then there
exists a subsequence whose limit point $(u^{*},v^{*},w^{*},y^{*},z^{*})$ is a
stationary point of (19) that satisfies the following:

$\displaystyle 0$ $\displaystyle=-\mu\Delta
u^{*}+A^{\top}y^{*}+\nabla^{\top}z^{*},$ (41a) $\displaystyle 0$
$\displaystyle=\lambda\left(\mathbbm{1}-\frac{f}{v^{*}}\right)-y^{*},$ (41b)
$\displaystyle z^{*}$
$\displaystyle\in\partial\left(\|w^{*}\|_{1}-\alpha\|w^{*}\|_{2,1}\right),$
(41c) $\displaystyle Au^{*}$ $\displaystyle=v^{*},$ (41d) $\displaystyle\nabla
u^{*}$ $\displaystyle=w^{*}.$ (41e)

###### Proof.

By Lemma 8, the sequence
$\\{(u_{k},v_{k},w_{k},y_{k},z_{k})\\}_{k=1}^{\infty}$ is bounded, so there
exists a subsequence
$\\{(u_{k_{n}},v_{k_{n}},w_{k_{n}},y_{k_{n}},z_{k_{n}})\\}_{n=1}^{\infty}$
that converges to a point $(u^{*},v^{*},w^{*},y^{*},z^{*})$. Additionally, we
have $u_{k+1}-u_{k}\rightarrow 0$ and $v_{k+1}-v_{k}\rightarrow 0$. Since
$\\{(y_{k},z_{k})\\}_{k=1}^{\infty}$ is bounded, there exists a constant $C>0$
such that $\|y_{k+1}-y_{k}\|_{2}<C$ and $\|z_{k+1}-z_{k}\|_{2}<C$ for each
$k\in\mathbb{N}$. By (21e), we have

$\displaystyle\|w_{k+1}-w_{k}\|_{2}$ $\displaystyle\leq\|w_{k+1}-\nabla
u_{k+1}\|_{2}+\|\nabla u_{k+1}-\nabla u_{k}\|_{2}+\|\nabla u_{k}-w_{k}\|_{2}$
$\displaystyle=\frac{\|z_{k+1}-z_{k}\|_{2}}{\beta_{k}}+\|\nabla u_{k+1}-\nabla
u_{k}\|_{2}+\frac{\|z_{k}-z_{k-1}\|_{2}}{\beta_{k-1}}$
$\displaystyle\leq\frac{2C}{\beta_{k-1}}+\|\nabla u_{k+1}-\nabla u_{k}\|_{2}.$

As $k\rightarrow\infty$, we have $w_{k+1}-w_{k}\rightarrow 0$. Altogether, we
can derive the following results:

$\displaystyle\lim_{n\rightarrow\infty}(u_{k_{n}+1},v_{k_{n}+1},w_{k_{n}+1})=\lim_{n\rightarrow\infty}(u_{k_{n}},v_{k_{n}},w_{k_{n}})=(u^{*},v^{*},w^{*}).$
(42)

Furthermore, the assumptions give us

$\displaystyle\lim_{n\rightarrow\infty}\beta_{k_{n}}(v_{k_{n}+1}-v_{k_{n}})$
$\displaystyle=0,$
$\displaystyle\lim_{n\rightarrow\infty}\beta_{k_{n}}(w_{k_{n}+1}-w_{k_{n}})$
$\displaystyle=0,$
$\displaystyle\lim_{n\rightarrow\infty}y_{k_{n}+1}-y_{k_{n}}$
$\displaystyle=0,$
$\displaystyle\lim_{n\rightarrow\infty}z_{k_{n}+1}-z_{k_{n}}$
$\displaystyle=0.$

By (21d)-(21e), we have

$\displaystyle\|Au^{*}-v^{*}\|_{2}$
$\displaystyle=\lim_{n\rightarrow\infty}\|Au_{k_{n}+1}-v_{k_{n}+1}\|_{2}=\lim_{n\rightarrow\infty}\frac{\|y_{k_{n}+1}-y_{k_{n}}\|_{2}}{\beta_{k_{n}}}\leq\lim_{n\rightarrow\infty}\frac{C}{\beta_{k_{n}}}=0,$
$\displaystyle\|\nabla u^{*}-w^{*}\|_{2}$
$\displaystyle=\lim_{n\rightarrow\infty}\|\nabla
u_{k_{n}+1}-w_{k_{n}+1}\|_{2}=\lim_{n\rightarrow\infty}\frac{\|z_{k+1}-z_{k}\|_{2}}{\beta_{k_{n}}}\leq\lim_{n\rightarrow\infty}\frac{C}{\beta_{k_{n}}}=0.$

Hence, we have $Au^{*}=v^{*}$ and $\nabla u^{*}=w^{*}$.

The optimality conditions at iteration $k_{n}$ are the following:

$\displaystyle-\mu\Delta
u_{k_{n}+1}+A^{\top}y_{k_{n}}+\beta_{k_{n}}A^{\top}(Au_{k_{n}+1}-v_{k_{n}})+\nabla^{\top}z_{k_{n}}+\beta_{k_{n}}\nabla^{\top}(\nabla
u_{k_{n}+1}-w_{k_{n}})=0,$ (43a)
$\displaystyle\lambda\left(\mathbbm{1}-\frac{f}{v_{k_{n}+1}}\right)-y_{k_{n}}-\beta_{k_{n}}(Au_{k_{n}+1}-v_{k_{n}+1})=0,$
(43b) $\displaystyle z_{k_{n}}+\beta_{k_{n}}(\nabla
u_{k_{n}+1}-w_{k_{n}+1})\in\partial(\|w_{k_{n}+1}\|_{1}-\alpha\|w_{k_{n}+1}\|_{2,1}).$
(43c)

Expanding (43a) by substituting in (21d)-(21e) and taking the limit, we have

$\displaystyle 0=$ $\displaystyle\lim_{n\rightarrow\infty}-\mu\Delta
u_{k_{n}+1}+A^{\top}y_{k_{n}}+\beta_{k_{n}}A^{\top}(Au_{k_{n}+1}-v_{k_{n}})+\nabla^{\top}z_{k_{n}}+\beta_{k_{n}}\nabla^{\top}(\nabla
u_{k_{n}+1}-w_{k_{n}})$ $\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}-\mu\Delta
u_{k_{n}+1}+A^{\top}y_{k_{n}}+\beta_{k_{n}}A^{\top}(Au_{k_{n}+1}-v_{k_{n}+1})+\beta_{k_{n}}A^{\top}(v_{k_{n}+1}-v_{k_{n}})+\nabla^{\top}z_{k_{n}}$
$\displaystyle+\beta_{k_{n}}\nabla^{\top}(\nabla
u_{k_{n}+1}-w_{k_{n}+1})+\beta_{k_{n}}\nabla^{\top}(w_{k_{n}+1}-w_{k_{n}})$
$\displaystyle=$ $\displaystyle\lim_{n\rightarrow\infty}-\mu\Delta
u_{k_{n}+1}+A^{\top}y_{k_{n}}+A^{\top}(y_{k_{n}+1}-y_{k_{n}})+\beta_{k_{n}}A^{\top}(v_{k_{n}+1}-v_{k_{n}})+\nabla^{\top}z_{k_{n}}$
$\displaystyle+\nabla^{\top}(z_{k_{n}+1}-z_{k_{n}})+\beta_{k_{n}}\nabla^{\top}(w_{k_{n}+1}-w_{k_{n}})$
$\displaystyle=$ $\displaystyle-\mu\Delta
u^{*}+A^{\top}y^{*}+\nabla^{\top}z^{*}.$

Substituting in (21d) into (43b) and taking the limit give us

$\displaystyle 0=$
$\displaystyle\lim_{n\rightarrow\infty}\lambda\left(\mathbbm{1}-\frac{f}{v_{k_{n}+1}}\right)-y_{k_{n}}-\beta_{k_{n}}(Au_{k_{n}+1}-v_{k_{n}+1})$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\lambda\left(\mathbbm{1}-\frac{f}{v_{k_{n}+1}}\right)-y_{k_{n}}-(y_{k_{n}+1}-y_{k_{n}})$
$\displaystyle=$
$\displaystyle\lambda\left(\mathbbm{1}-\frac{f}{v^{*}}\right)-y^{*}.$

Lastly, by substituting (21e) into (43c), we have

$\displaystyle
z_{k_{n}+1}\in\partial(\|w_{k_{n}+1}\|_{1}-\alpha\|w_{k_{n}+1}\|_{2,1}).$

By continuity, we have
$\|w_{k_{n}+1}\|_{1}-\alpha\|w_{k_{n}+1}\|_{2,1}\rightarrow\|w^{*}\|_{1}-\alpha\|w^{*}\|_{2,1}$.
Together with the fact that
$(w_{k_{n}+1},z_{k_{n}+1})\rightarrow(w^{*},z^{*})$, we have
$z^{*}\in\partial\left(\|w^{*}\|_{1}-\alpha\|w^{*}\|_{2,1}\right)$ by
closedness of the subdifferential.

Therefore, $(u^{*},v^{*},w^{*},y^{*},z^{*})$ is a stationary point. ∎

######  Remark 2.

It is true that the assumptions in Theorem 9 are rather strong, but they are
standard in the convergence analyses of other ADMM algorithms for nonconvex
problems that fail to satisfy the conditions for global convergence in [69].
For example, [34, 35, 39, 43] assumed convergence of the successive
differences of the primal variables and Lagrange multipliers. Instead, we
modify the convergence of the successive difference of the primal variables,
i.e., $\beta_{k}(v_{k+1}-v_{k})\rightarrow
0,\beta_{k}(w_{k+1}-w_{k})\rightarrow 0$. Boundedness of the Lagrange
multiplier (i.e., $\\{y_{k}\\}_{k=1}^{\infty}$) was also assumed in [44, 76],
which required a stronger assumption than ours regarding the successive
difference of the Lagrange multipliers.

##  4 Numerical Experiments

In this section, we apply the proposed method of AITV Poisson SaT/SLaT on
various grayscale and color images for image segmentation. For grayscale
images, we compare our method with the original TV SaT [16], thresholded-
Rudin-Osher-Fatemi (T-ROF) [11], and the Potts model [55] solved by either
Pock’s algorithm (Pock) [54] or Storath and Weinmann’s algorithm (Storath)
[60]. For color images, we compare with TV SLaT [10], Pock’s method [54], and
Storath’s method [60]. We can solve (9) for TV SaT/SLaT via Algorithm 2 that
utilizes the proximal operator corresponding to the $\|\cdot\|_{2,1}$ norm.
The code for T-ROF is provided by the respective
author111https://xiaohaocai.netlify.app/download/ and we can adapt it to
handle blur by using a more general data fidelity term. Pock’s method is
implemented by the lab group222Python code is available at
https://github.com/VLOGroup/pgmo-lecture/blob/master/notebooks/tv-potts.ipynb
and a translated MATLAB code is available at
https://github.com/kbui1993/MATLAB_Potts.. Storath’s method is provided by the
original author333 https://github.com/mstorath/Pottslab. Note that T-ROF,
Pock’s method, and Storath’s method are designed for images corrupted with
Gaussian noise. We apply the Anscombe transform [2] to the test images, after
which the Poisson noise becomes approximately Gaussian noise. Since Storath’s
method is not for segmentation, we perform a post-processing step of $k$-means
clustering to its piecewise-constant output. For the SLaT methods, we
parallelize the smoothing step separately for each channel.

Figure 1: The entire DRIVE dataset [59] for binary segmentation. The image
size is $584\times 565$ with background value of 200 and the pixel value for
vessels to be 255.

To quantitatively measure the segmentation performance, we use the DICE index
[24] and peak signal-to-noise ratio (PSNR). Let $S\subset\Omega$ be the
ground-truth region and $S^{\prime}\subset\Omega$ be a region obtained from
the segmentation algorithm corresponding to the ground-truth region $S$. The
DICE index is formulated by

$\displaystyle\text{DICE}=\frac{2|S\cap S^{\prime}|}{|S|+|S^{\prime}|}.$

To compare the piecewise-constant reconstruction $\tilde{f}$ according to (10)
with the original test image $f$, we compute PSNR by

$\displaystyle\text{PSNR}=20\log_{10}\frac{(M\times N)\times
P}{\sum_{i,j}(f_{i,j}-\tilde{f}_{i,j})^{2}},$

where $M\times N$ is the image size and $P=\max_{i,j}f_{i,j}$.

Poisson noise is added to the test images by the MATLAB command poissrnd. To
ease parameter tuning, we scale each test image to $[0,1]$ after its
degradation with Poisson noise and/or blur. We set $\sigma=1.25$ and
$\beta_{1,0}=\beta_{2,0}=1.0,2.0$ in Algorithm 2 for grayscale and color
images, respectively. The stopping criterion is either 300 iterations or when
the relative error of $u_{k}$ is below $\epsilon=10^{-4}$. We tune the
fidelity parameter $\lambda$ and the smoothing parameter $\mu$ for each image,
which will be specified later. For T-ROF, Pock’s method, and Storath’s method,
their parameters are manually tuned to give the best DICE indices for binary
segmentation (Section 4.1) and the PSNR values for multiphase segmentation
(Section 4.2-4.3). All experiments are performed in MATLAB R2022b on a Dell
laptop with a 1.80 GHz Intel Core i7-8565U processor and 16.0 GB RAM.

###  4.1 Grayscale, Binary Segmentation

We start with performing binary segmentation on the entire DRIVE dataset [59]
that consists of 20 images shown in Figure 1. Each image has size $584\times
565$ with modified pixel values of either 200 for the background or 255 for
the vessels. Before adding Poisson noise, we set the peak value of the image
to be $P/2$ or $P/5$, where $P=255$. Note that a lower peak value indicates
stronger noise in the image, thus more challenging for denoising. We examine
three cases: (1) $P/2$ no blur, (2) $P/5$ no blur, and (3) $P/2$ with Gaussian
blur specified by MATLAB command fspecial(’gaussian’, [10 10], 2). For the TV
SaT method, we set $\lambda=14.5,\ \mu=0.5$ for case (1), $\lambda=8.0,\
\mu=0.5$ for case (2), and $\lambda=22.5,\ \mu=0.25$ for case (3). For the
AITV SaT method, the parameters $\lambda$ and $\mu$ are set the same as TV
SaT, and we have $\alpha=0.3$ for cases (1)-(2) and $\alpha=0.8$ for case (3).

Table 1 records the DICE indices and the computational time in seconds for the
competing methods, averaged over 20 images. We observe that AITV SaT attains
the best DICE indices for all three cases with comparable computational time
to TV SaT and T-ROF, all of which are much faster than Pock and Storath. As
visually illustrated in Figure 2, AITV SaT segments more of the thinner
vessels compared to TV SaT and T-ROF in five images, thereby having the higher
average DICE indices.

Figure 2: Binary segmentation results of Figure 1 with peak $P/2$ under Gaussian blur and Poisson noise. Table 1: DICE and computational time in seconds of the binary segmentation methods averaged over 20 images in Figure 1 with standard deviations in parentheses. Bold indicates best result. |  | TV SaT | AITV SaT | T-ROF | Pock | Storath
---|---|---|---|---|---|---
$P/2$ no blur | DICE |  0.9464 (0.0091) |  0.9501 (0.0076) |  0.9463 (0.0073) |  0.8466 (0.0301) |  0.8855 (0.0181)
Time (sec.) |  4.2401 (0.3618) |  5.7342 (0.5251) |  4.9206 (1.4281) |  24.7376 (3.2454) |  19.9456 (1.8875)
$P/5$ no blur | DICE |  0.8714 (0.0134) |  0.8735 (0.0125) |  0.8570 (0.0170) |  0.6504 (0.0910) |  0.8277 (0.0191)
Time (sec.) |  4.7076 (0.6937) |  6.4027 (0.8441) |  5.4943 (0.7935) |  46.9346 (9.4969) |  21.8734 (2.7660)
$P/2$ with Gaussian Blur | DICE |  0.7244 (0.0254) |  0.7411 (0.0220) |  0.7322 (0.0251) |  0.5473 (0.0398) |  0.6944 (0.0217)
Time (sec.) |  7.4495 (1.0983) |  9.2523 (1.5959) |  11.7337 (2.1252) |  47.3911 (10.9191) |  19.9444 (2.4142)

###  4.2 Grayscale, Multiphase Segmentation

Figure 3: The entire BrainWeb dataset [3] for grayscale, multiphase segmentation. Each image is of size $104\times 87$. The pixel values are 10 (background), 48 (cerebrospinal fluid), 106 (grey matter), and 154 (white matter). Figure 4: Segmentation result of the first image of Figure 3 with peak $P/8$ under Poisson noise with no blur. From top to bottom are segmentation results for CSF, GM, and WM.  Table 2: DICE and computational time in seconds of the multiphase segmentation methods averaged over 20 images in Figure 3 with standard deviations in parentheses. Bold indicates best result. |  | TV SaT | AITV SaT | T-ROF | Pock | Storath
---|---|---|---|---|---|---
$P/8$ no blur | CSF DICE |  0.8208 (0.0270) |  0.8396 (0.0262) |  0.7398 (0.1316) |  0.4572 (0.1474) |  0.8041 (0.0244)
GM DICE |  0.8097 (0.0258) |  0.8477 (0.0157) |  0.7904 (0.0422) |  0.7507 (0.0719) |  0.7900 (0.0344)
WM DICE |  0.8449 (0.0125) |  0.8694 (0.0101) |  0.8221 (0.0132) |  0.8459 (0.0283) |  0.8138 (0.0196)
Time (sec.) |  0.2607 (0.0625) |  0.2863 (0.0470) |  0.2202 (0.0324) |  2.6139 (0.6690) |  0.3383 (0.0948)
$P/8$ with Motion Blur | CSF DICE |  0.6196 (0.0385) |  0.6260 (0.0483) |  0.6174 (0.0460) |  0.3772 (0.0561) |  0.4964 (0.0468)
GM DICE |  0.6809 (0.0304) |  0.7138 (0.0262) |  0.6528 (0.0358) |  0.6345 (0.0590) |  0.6544 (0.0399)
WM DICE |  0.7757 (0.0127) |  0.7935 (0.0110) |  0.7686 (0.0164) |  0.7529 (0.0185) |  0.7382 (0.0140)
Time (sec.) |  0.2494 (0.0443) |  0.2854 (0.0279) |  0.3647 (0.0693) |  2.5782 (0.4683) |  0.3052 (0.1399)

We examine the multiphase segmentation on the entire BrainWeb dataset [3] that
consists of 20 grayscale images as shown in Figure 3. Each image is of size
$104\times 87$ and has four regions to segment: background, cerebrospinal
fluid (CSF), grey matter (GM), and white matter (WM). The pixel values are 10
(background), 48 (CSF), 106 (GM), and 154 (WM). The maximum intensity $P=154$.
We consider two cases: (1) $P/2$ no blur and (2) $P/2$ with motion blur
specified by fspecial(’motion’, 5, 225). For the SaT methods, we have
$\mu=1.0$, $\alpha=0.6,0.7$, and $\lambda=4.0,5.0$ for case (1) and case (2),
respectively.

Across all 20 images of the BrainWeb dataset, Table 2 reports the average DICE
indices for CSF, GM, and WM and average computational times in seconds of the
segmentation methods. For both cases (1) and (2), AITV SaT attains the highest
average DICE indices for segmenting CSF, GM, and WM. AITV SaT is comparable to
TV SaT and T-ROF in terms of computational time.

Figure 4 shows the segmentation results of the first image in Figure 3 for
case (1). When segmenting CSF, the methods (TV SaT, AITV SaT, and Storath)
yield similar visual results, while Pock fails to segment roughly half of the
region. In addition, AITV SaT segments the most GM region with the least
amount of noise artifacts than the other methods. Lastly, for WM segmentation,
AITV SaT avoids the “holes” or “gaps” and segments fewer regions outside of
the ground truth, thus outperforming TV SaT and Storath. For the three
regions, T-ROF has the most noise artifacts in its segmentation results.

###  4.3 Color Segmentation

We perform color image segmentation on 10 images shown in Figure 5, which are
selected from the PASCAL VOC 2010 dataset [28]. Each image has six different
color regions. Figure 5(A) is of size $307\times 461$; Figure 5(B) is of
$500\times 367$; Figures 5(C) and 5(H)-(I) are of $500\times 375$; and Figures
5(D)-(G) and 5(J) are of $375\times 500$. Before adding Poisson noise to each
channel of each image, we set the peak value $P=10$. We choose the parameters
of the SLaT methods to be $\lambda=1.5$, $\mu=0.05$, and $\alpha=0.6$.

Figures 6-7 present the piecewise-constant approximations via (10), showing
similar segmentation results obtained by TV SLaT, AITV SLaT, and Storath.
Quantitatively in Table 3, AITV SLaT has better PSNRs than TV SLaT and Pock
for all the images and outperforms Storath for seven. Overall, the proposed
AITV SLaT has the highest PSNR on average over 10 images with the lowest
standard deviation and comparable speed as Storath.

Figure 5: Test images from the PASCAL VOC 2010 dataset [28] for color, multiphase segmentation. Each image has 6 regions. The image sizes are (A) $307\times 461$, (B) $500\times 367$, (C) $500\times 375$, (D)-(G) $375\times 500$, (H)-(I) $500\times 375$, and (J) $375\times 500$. Figure 6: Color image segmentation results of Figures 5(A)-(E). Figure 7: Color image segmentation results of Figures 5(F)-(J).  Table 3: PSNR and computation time in seconds for color image segmentation of the images in Figure 5. Bold indicates best result. |  | TV SLaT | AITV SLaT | Pock | Storath
---|---|---|---|---|---
PSNR | Figure 5(A) | 29.0337 | 30.4619 | 19.2306 | 28.7141
Figure 5(B) | 31.8673 | 31.7055 | 17.3512 | 30.9703
Figure 5(C) | 30.9605 | 33.4217 | 18.0389 | 32.2217
Figure 5(D) | 29.5265 | 32.5505 | 21.4873 | 34.7881
Figure 5(E) | 29.8903 | 31.0656 | 19.7646 | 28.3456
Figure 5(F) | 33.2308 | 34.7619 | 17.9788 | 34.4106
Figure 5(G) | 28.1136 | 30.6237 | 22.5048 | 31.2439
Figure 5(H) | 33.2682 | 33.4250 | 18.9390 | 31.9377
Figure 5(I) | 30.0659 | 31.7937 | 20.0856 | 29.3905
Figure 5(J) | 31.4164 | 34.0610 | 20.3185 | 34.3599
|  Avg. (Std.) |  30.7373 (1.7266) |  32.3870 (1.4905) |  19.5700 (1.6147) |  31.6383 (2.3686)
Time (sec.) | Figure 5(A) | 7.1992 | 8.2455 | 148.9232 | 4.9410
Figure 5(B) | 12.2238 | 11.8890 | 321.0637 | 9.9533
Figure 5(C) | 5.6775 | 7.4494 | 384.6998 | 10.1785
Figure 5(D) | 4.9511 | 7.2061 | 128.4236 | 7.0026
Figure 5(E) | 6.2380 | 7.0719 | 275.0501 | 7.9808
Figure 5(F) | 6.0509 | 7.2392 | 398.1339 | 10.0982
Figure 5(G) | 7.2821 | 7.6919 | 165.2142 | 7.9315
Figure 5(H) | 6.0269 | 9.7018 | 385.3311 | 7.9389
Figure 5(I) | 7.0375 | 8.0875 | 250.3698 | 9.7501
Figure 5(J) | 6.1873 | 10.6472 | 334.8662 | 7.2176
|  Avg. (Std.) |  6.8874 (2.0088) |  8.5230 (1.6610) |  279.2076 (102.7063) |  8.2993 (1.7031)

###  4.4 Parameter Analysis

The proposed smoothing model (11) involves the following parameters:

  * •

The fidelity parameter $\lambda$ weighs how close the approximation $Au^{*}$
is to the original image $f$. For a larger amount of noise, the value of
$\lambda$ should be chosen smaller.

  * •

The smoothing parameter $\mu$ determines how smooth the solution $u^{*}$
should be. A larger value of $\mu$ may improve denoising, but at a cost of
smearing out the edges between adjacent regions, which may be segmented
together if they have similar colors.

  * •

The sparsity parameter $\alpha\in[0,1]$ determines how sparse the gradient
vector at each pixel should be. More specifically, the closer the value
$\alpha$ is to 1, the more $\|\nabla u\|_{1}-\alpha\|\nabla u\|_{2,1}$
resembles $\|\nabla u\|_{0}$.

We perform sensitivity analysis on these parameters to understand how they
affect the segmentation quality of AITV SaT/SLaT. We consider two types of
tests in the case of $P/8$ with motion blur in Figure 3. In the first case
(Figure 8), we fix $\mu=1.0$ and vary $\alpha,\lambda$. In the second case
(Figure 9), we fix $\lambda=5.0$ and vary $\alpha,\mu$. Figure 8 reveals a
concave relationship of the DICE index of each region with respect to the
parameter $\lambda$, which implies there exists the optimal choice of
$\lambda.$ Additionally, when $\lambda$ is small, a large value for $\alpha$
can improve the DICE indices. According to Figure 9, the DICE indices of the
GM and WM regions decrease with respect to $\mu$, while the DICE index of the
CSF region is approximately constant. For $\alpha=0.8$, the DICE indices of
the GM and WM regions are the largest when $\mu\geq 1$, but the large $\alpha$
is not optimal for CSF. Hence, an intermediate value of $\alpha$, such as
$0.6$, is preferable to attain satisfactory segmentation quality for all three
regions.

Lastly, in Figure 10, we conduct sensitive analysis for the case of $P=10$ of
Figure 5. We fix $\mu=0.05$, while varying $\lambda,\alpha$ in Figure 10(A),
which indicates that the optimal value for $\alpha$ is in the range of
$0.5\leq\alpha\leq 0.7$. Then we fix $\lambda=1.5$ to examine $\mu$ and
$\alpha$ in Figure 10(B). For $0.2\leq\alpha\leq 0.7$, PSNR decreases as $\mu$
increases. Again, $\alpha=0.6$ generally yields the best PSNR.

Figure 8: Sensitivity analysis on $\lambda$ for the $P/8$ with motion blur
case of Figure 3. The parameter $\mu=1.0$ is fixed. DICE indices averaged over
10 images for each brain region are plotted with respect to $\lambda$.  Figure
9: Sensitivity analysis on $\mu$ for the $P/8$ with motion blur case of Figure
3. The parameter $\lambda=5.0$ is fixed. DICE indices averaged over 10 images
for each brain region are plotted with respect to $\mu$.  Figure 10:
Sensitivity analysis of parameters for $P=10$ case of Figure 5. (A) is the
sensitivity analysis on $\lambda$ when $\mu=0.05$ fixed; (B) is the
sensitivity analysis on $\mu$ when $\lambda=1.5$ fixed. Average PSNR is
plotted.

##  5 Conclusion and future work

In this paper, we developed the AITV Poisson SaT/SLaT framework for image
segmentation. In particular, we proposed a simplified Mumford-Shah model with
the AITV regularization and Poisson fidelity for the smoothing step. The model
was proven to have a global minimizer. Our numerical algorithm incorporated a
specific splitting scheme for ADMM and the $\ell_{1}-\alpha\ell_{2}$ proximal
operator for solving a subproblem. Convergence analysis established that the
sequence generated by ADMM has a convergent subsequence to a stationary point
of the nonconvex model. In our numerical experiments, the AITV Poisson
SaT/SLaT yielded high-quality segmentation results within seconds for various
grayscale and color images corrupted with Poisson noise and/or blur. For
future directions, we are interested in other nonconvex regularization, such
as $\ell_{1}/\ell_{2}$ on the gradient [65, 66, 72], $\ell_{p},0<p<1$, on the
gradient [32, 43, 73], and transformed total variation [33], as alternatives
to AITV. On the other hand, we can develop AITV variants of weighted TV [41]
or adaptive TV [71, 80]. Moreover, we plan to determine how to make the
sparsity parameter $\alpha$ in AITV adaptable to each image. In future work,
we will adapt other segmentation algorithms [11, 34, 35, 38, 39, 52, 78]
designed for Gaussian noise or impulsive noise to Poisson noise.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

KB performed the experiments and analysis and drafted the manuscript. All
authors contributed to the design, evaluation, discussions and production of
the manuscript.

## Funding

The work was partially supported by NSF grants DMS-1846690, DMS-1854434,
DMS-1952644, DMS-2151235, and a Qualcomm Faculty Award.

## Data Availability Statement

The images in Figure 1 are provided from the DRIVE dataset [59] at
https://drive.grand-challenge.org/DRIVE/. The images in Figure 3 are extracted
from BrainWeb [3] via the Python package “brainweb” provided at
https://github.com/casperdcl/brainweb. The images in Figure 5 are selected
from the PASCAL VOC 2010 dataset [28]. Code for AITV Poisson SaT/SLaT is
available at https://github.com/kbui1993/Official_Poisson_AITV_SaT_SLaT.

## References

  * [1] L. Ambrosio and V. M. Tortorelli, Approximation of functional depending on jumps by elliptic functional via t-convergence, Communications on Pure and Applied Mathematics, 43 (1990), pp. 999–1036.
  * [2] F. J. Anscombe, The transformation of Poisson, binomial and negative-binomial data, Biometrika, 35 (1948), pp. 246–254.
  * [3] B. Aubert-Broche, M. Griffin, G. B. Pike, A. C. Evans, and D. L. Collins, Twenty new digital brain phantoms for creation of validation image data bases, IEEE Transactions on Medical Imaging, 25 (2006), pp. 1410–1416.
  * [4] L. Bar, T. F. Chan, G. Chung, M. Jung, N. Kiryati, R. Mohieddine, N. Sochen, and L. A. Vese, Mumford and Shah model and its applications to image segmentation and image restoration, in Handbook of Mathematical Methods in Imaging, Springer, 2011, pp. 1095–1157.
  * [5] A. Beck, First-order methods in optimization, SIAM, 2017.
  * [6] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends® in Machine Learning, 3 (2011), pp. 1–122.
  * [7] K. Bui, J. Fauman, D. Kes, L. Torres Mandiola, A. Ciomaga, R. Salazar, A. L. Bertozzi, J. Gilles, D. P. Goronzy, A. I. Guttentag, et al., Segmentation of scanning tunneling microscopy images using variational methods and empirical wavelets, Pattern Analysis and Applications, 23 (2020), pp. 625–651.
  * [8] K. Bui, Y. Lou, F. Park, and J. Xin, An efficient smoothing and thresholding image segmentation framework with weighted anisotropic-isotropic total variation, arXiv preprint arXiv:2202.10115, (2022).
  * [9] K. Bui, F. Park, Y. Lou, and J. Xin, A weighted difference of anisotropic and isotropic total variation for relaxed Mumford–Shah color and multiphase image segmentation, SIAM Journal on Imaging Sciences, 14 (2021), pp. 1078–1113.
  * [10] X. Cai, R. Chan, M. Nikolova, and T. Zeng, A three-stage approach for segmenting degraded color images: Smoothing, lifting and thresholding (SLaT), Journal of Scientific Computing, 72 (2017), pp. 1313–1332.
  * [11] X. Cai, R. Chan, C.-B. Schonlieb, G. Steidl, and T. Zeng, Linkage between piecewise constant Mumford–Shah model and Rudin–Osher–Fatemi model and its virtue in image segmentation, SIAM Journal on Scientific Computing, 41 (2019), pp. B1310–B1340.
  * [12] X. Cai, R. Chan, and T. Zeng, A two-stage image segmentation method using a convex variant of the Mumford–Shah model and thresholding, SIAM Journal on Imaging Sciences, 6 (2013), pp. 368–390.
  * [13] W. Cao, J. Sun, and Z. Xu, Fast image deconvolution using closed-form thresholding formulas of $L_{q}(q=\frac{1}{2},\frac{2}{3})$ regularization, Journal of Visual Communication and Image Representation, 24 (2013), pp. 31–41.
  * [14] P. Cascarano, L. Calatroni, and E. L. Piccolomini, Efficient $\ell^{0}$ gradient-based super-resolution for simplified image segmentation, IEEE Transactions on Computational Imaging, 7 (2021), pp. 399–408.
  * [15] A. Chambolle and T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging, Journal of mathematical imaging and vision, 40 (2011), pp. 120–145.
  * [16] R. Chan, H. Yang, and T. Zeng, A two-stage image segmentation method for blurry images with Poisson or multiplicative Gamma noise, SIAM Journal on Imaging Sciences, 7 (2014), pp. 98–127.
  * [17] T. F. Chan, S. Esedoglu, and M. Nikolova, Algorithms for finding global minimizers of image segmentation and denoising models, SIAM Journal on Applied Mathematics, 66 (2006), pp. 1632–1648.
  * [18] T. F. Chan and L. A. Vese, Active contours without edges, IEEE Transactions on Image Processing, 10 (2001), pp. 266–277.
  * [19] H. Chang, Y. Lou, Y. Duan, and S. Marchesini, Total variation–based phase retrieval for Poisson noise removal, SIAM Journal on Imaging Sciences, 11 (2018), pp. 24–55.
  * [20] R. Chartrand, Exact reconstruction of sparse signals via nonconvex minimization, IEEE Signal Processing Letters, 14 (2007), pp. 707–710.
  * [21] R. Choksi, Y. G. Gennip, and A. Oberman, Anisotropic total variation regularized $L^{1}$ approximation and denoising/deblurring of 2D bar codes, Inverse Problems & Imaging, 5 (2011), pp. 591–617.
  * [22] M. R. Chowdhury, J. Qin, and Y. Lou, Non-blind and blind deconvolution under Poisson noise using fractional-order total variation, Journal of Mathematical Imaging and Vision, 62 (2020), pp. 1238–1255.
  * [23] M. R. Chowdhury, J. Zhang, J. Qin, and Y. Lou, Poisson image denoising based on fractional-order total variation, Inverse Problems & Imaging, 14 (2020).
  * [24] L. R. Dice, Measures of the amount of ecologic association between species, Ecology, 26 (1945), pp. 297–302.
  * [25] L. Ding and W. Han, $\alpha\ell_{1}-\beta\ell_{2}$ regularization for sparse recovery, Inverse Problems, 35 (2019), p. 125009.
  * [26] Y. Duan, H. Chang, W. Huang, J. Zhou, Z. Lu, and C. Wu, The ${L}_{0}$ regularized Mumford–Shah model for bias correction and segmentation of medical images, IEEE Transactions on Image Processing, 24 (2015), pp. 3927–3938.
  * [27] L. C. Evans, Partial differential equations, vol. 19, American Mathematical Soc., 2010.
  * [28] M. Everingham, L. Van Gool, C. K. Williams, J. Winn, and A. Zisserman, The PASCAL visual object classes (VOC) challenge, International Journal of Computer Vision, 88 (2009), pp. 303–308.
  * [29] H. Ge and P. Li, The Dantzig selector: recovery of signal via $\ell_{1}-\alpha\ell_{2}$ minimization, Inverse Problems, 38 (2021), p. 015006.
  * [30] S. Gu, Q. Xie, D. Meng, W. Zuo, X. Feng, and L. Zhang, Weighted nuclear norm minimization and its applications to low level vision, International Journal of Computer Vision, 121 (2017), pp. 183–208.
  * [31] W. Guo, Y. Lou, J. Qin, and M. Yan, A novel regularization based on the error function for sparse recovery, Journal of Scientific Computing, 87 (2021), pp. 1–22.
  * [32] M. Hintermüller and T. Wu, Nonconvex TVq-models in image restoration: Analysis and a trust-region regularization–based superlinearly convergent solver, SIAM Journal on Imaging Sciences, 6 (2013), pp. 1385–1415.
  * [33] L. Huo, W. Chen, H. Ge, and M. K. Ng, Stable image reconstruction using transformed total variation minimization, SIAM Journal on Imaging Sciences, 15 (2022), pp. 1104–1139.
  * [34] M. Jung, Piecewise-smooth image segmentation models with $l^{1}$ data-fidelity terms, Journal of Scientific Computing, 70 (2017), pp. 1229–1261.
  * [35] M. Jung, M. Kang, and M. Kang, Variational image segmentation models involving non-smooth data-fidelity terms, Journal of Scientific Computing, 59 (2014), pp. 277–308.
  * [36] H. Lantéri and C. Theys, Restoration of astrophysical images—the case of Poisson data with additive Gaussian noise, EURASIP Journal on Advances in Signal Processing, 2005 (2005), pp. 1–14.
  * [37] T. Le, R. Chartrand, and T. J. Asaki, A variational approach to reconstructing images corrupted by Poisson noise, Journal of Mathematical Imaging and Vision, 27 (2007), pp. 257–263.
  * [38] F. Li, M. K. Ng, T. Y. Zeng, and C. Shen, A multiphase image segmentation method based on fuzzy region competition, SIAM Journal on Imaging Sciences, 3 (2010), pp. 277–299.
  * [39] F. Li, S. Osher, J. Qin, and M. Yan, A multiphase image segmentation based on fuzzy membership functions and L1-norm fidelity, Journal of Scientific Computing, 69 (2016), pp. 82–106.
  * [40] H. Li, W. Guo, J. Liu, L. Cui, and D. Xie, Image segmentation with adaptive spatial priors from joint registration, SIAM Journal on Imaging Sciences, 15 (2022), pp. 1314–1344.
  * [41] M.-M. Li and B.-Z. Li, A novel weighted anisotropic total variational model for image applications, Signal, Image and Video Processing, 16 (2022), pp. 211–218.
  * [42] P. Li, W. Chen, H. Ge, and M. K. Ng, $\ell_{1}-\alpha\ell_{2}$ minimization methods for signal and image reconstruction with impulsive noise removal, Inverse Problems, 36 (2020), p. 055009.
  * [43] Y. Li, C. Wu, and Y. Duan, The TVp regularized Mumford-Shah model for image labeling and segmentation, IEEE Transactions on Image Processing, 29 (2020), pp. 7061–7075.
  * [44] H. Liu, K. Deng, H. Liu, and Z. Wen, An entropy-regularized ADMM for binary quadratic programming, Journal of Global Optimization, (2022), pp. 1–33.
  * [45] Y. Lou, S. Osher, and J. Xin, Computational aspects of constrained $L_{1}-L_{2}$ minimization for compressive sensing, in Modelling, Computation and Optimization in Information Systems and Management Sciences, Springer, 2015, pp. 169–180.
  * [46] Y. Lou and M. Yan, Fast L1–L2 minimization via a proximal operator, Journal of Scientific Computing, 74 (2018), pp. 767–785.
  * [47] Y. Lou, P. Yin, Q. He, and J. Xin, Computing sparse representation in a highly coherent dictionary based on difference of $L_{1}$ and $L_{2}$, Journal of Scientific Computing, 64 (2015), pp. 178–196.
  * [48] Y. Lou, T. Zeng, S. Osher, and J. Xin, A weighted difference of anisotropic and isotropic total variation model for image processing, SIAM Journal on Imaging Sciences, 8 (2015), pp. 1798–1823.
  * [49] Q.-T. Luong, Color in computer vision, in Handbook of Pattern Recognition and Computer Vision, World Scientific, 1993, pp. 311–368.
  * [50] D. B. Mumford and J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems, Communications on Pure and Applied Mathematics, (1989).
  * [51] S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics, 79 (1988), pp. 12–49.
  * [52] Z.-F. Pang, M. Geng, L. Zhang, Y. Zhou, T. Zeng, L. Zheng, N. Zhang, D. Liang, H. Zheng, Y. Dai, et al., Adaptive weighted curvature-based active contour for ultrasonic and 3T/5T MR image segmentation, Signal Processing, 205 (2023), p. 108881.
  * [53] F. Park, Y. Lou, and J. Xin, A weighted difference of anisotropic and isotropic total variation for relaxed Mumford-Shah image segmentation, in 2016 IEEE International Conference on Image Processing (ICIP), IEEE, 2016, pp. 4314–4318.
  * [54] T. Pock, A. Chambolle, D. Cremers, and H. Bischof, A convex relaxation approach for computing minimal partitions, in 2009 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2009, pp. 810–817.
  * [55] R. B. Potts, Some generalized order-disorder transformations, Mathematical Proceedings of the Cambridge Philosophical Society, 48 (1952), p. 106–109.
  * [56] Y. Rahimi, C. Wang, H. Dong, and Y. Lou, A scale-invariant approach for sparse signal recovery, SIAM Journal of Scientific Computing, 41 (2019), pp. A3649–A3672.
  * [57] R. T. Rockafellar and R. J.-B. Wets, Variational analysis, vol. 317, Springer Science & Business Media, 2009.
  * [58] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), pp. 259–268.
  * [59] J. Staal, M. D. Abràmoff, M. Niemeijer, M. A. Viergever, and B. Van Ginneken, Ridge-based vessel segmentation in color images of the retina, IEEE Transactions on Medical Imaging, 23 (2004), pp. 501–509.
  * [60] M. Storath and A. Weinmann, Fast partitioning of vector-valued images, SIAM Journal on Imaging Sciences, 7 (2014), pp. 1826–1852.
  * [61] M. Storath, A. Weinmann, and L. Demaret, Jump-sparse and sparse recovery using Potts functionals, IEEE Transactions on Signal Processing, 62 (2014), pp. 3654–3666.
  * [62] S. Tongbram, B. A. Shimray, L. S. Singh, and N. Dhanachandra, A novel image segmentation approach using fcm and whale optimization algorithm, Journal of Ambient Intelligence and Humanized Computing, (2021), pp. 1–15.
  * [63] Y. Vardi, L. A. Shepp, and L. Kaufman, A statistical model for positron emission tomography, Journal of the American Statistical Association, 80 (1985), pp. 8–20.
  * [64] L. A. Vese and T. F. Chan, A multiphase level set framework for image segmentation using the Mumford and Shah model, International Journal of Computer Vision, 50 (2002), pp. 271–293.
  * [65] C. Wang, M. Tao, C.-N. Chuah, J. Nagy, and Y. Lou, Minimizing $l_{1}$ over $l_{2}$ norms on the gradient, Inverse Problems, 38 (2022), p. 065011.
  * [66] C. Wang, M. Tao, J. G. Nagy, and Y. Lou, Limited-angle CT reconstruction via the $l_{1}/l_{2}$ minimization, SIAM Journal on Imaging Sciences, 14 (2021), pp. 749–777.
  * [67] C. Wang, M. Yan, Y. Rahimi, and Y. Lou, Accelerated schemes for the $l_{1}/l_{2}$ minimization, IEEE Transaction on Signal Processing, 68 (2020), pp. 2660–2669.
  * [68] Y. Wang, J. Yang, W. Yin, and Y. Zhang, A new alternating minimization algorithm for total variation image reconstruction, SIAM Journal on Imaging Sciences, 1 (2008), pp. 248–272.
  * [69] Y. Wang, W. Yin, and J. Zeng, Global convergence of ADMM in nonconvex nonsmooth optimization, Journal of Scientific Computing, 78 (2019), pp. 29–63.
  * [70] Y. Wen, R. H. Chan, and T. Zeng, Primal-dual algorithms for total variation based image restoration under Poisson noise, Science China Mathematics, 59 (2016), pp. 141–160.
  * [71] T. Wu, X. Gu, Y. Wang, and T. Zeng, Adaptive total variation based image segmentation with semi-proximal alternating minimization, Signal Processing, 183 (2021), p. 108017.
  * [72] T. Wu, Z. Mao, Z. Li, Y. Zeng, and T. Zeng, Efficient color image segmentation via quaternion-based $L_{1}/L_{2}$ regularization, Journal of Scientific Computing, 93 (2022), p. 9.
  * [73] T. Wu, J. Shao, X. Gu, M. K. Ng, and T. Zeng, Two-stage image segmentation based on nonconvex $\ell_{2}-\ell_{p}$ approximation and thresholding, Applied Mathematics and Computation, 403 (2021), p. 126168.
  * [74] T. Wu, Y. Zhao, Z. Mao, L. Shi, Z. Li, and Y. Zeng, Image segmentation via Fischer-Burmeister total variation and thresholding, Advances in Applied Mathematics and Mechanics, 14 (2022), pp. 960–988.
  * [75] Y. Xu, A. Narayan, H. Tran, and C. G. Webster, Analysis of the ratio of $l_{1}$ and $l_{2}$ norms in compressed sensing, Applied and Computational Harmonic Analysis, 55 (2021), pp. 486–511.
  * [76] Y. Xu, W. Yin, Z. Wen, and Y. Zhang, An alternating direction algorithm for matrix completion with nonnegative factors, Frontiers of Mathematics in China, 7 (2012), pp. 365–384.
  * [77] Z. Xu, X. Chang, F. Xu, and H. Zhang, $L_{1/2}$ regularization: A thresholding representation theory and a fast solver, IEEE Transactions on Neural Networks and Learning Systems, 23 (2012), pp. 1013–1027.
  * [78] J. Yang, Z. Guo, D. Zhang, B. Wu, and S. Du, An anisotropic diffusion system with nonlinear time-delay structure tensor for image enhancement and segmentation, Computers & Mathematics with Applications, 107 (2022), pp. 29–44.
  * [79] J. You, Y. Jiao, X. Lu, and T. Zeng, A nonconvex model with minimax concave penalty for image restoration, Journal of Scientific Computing, 78 (2019), pp. 1063–1086.
  * [80] T. Zhang, J. Chen, C. Wu, Z. He, T. Zeng, and Q. Jin, Edge adaptive hybrid regularization model for image deblurring, Inverse Problems, 38 (2022), p. 065010.
  * [81] D. Zosso, J. An, J. Stevick, N. Takaki, M. Weiss, L. S. Slaughter, H. H. Cao, P. S. Weiss, and A. L. Bertozzi, Image segmentation with dynamic artifacts detection and bias correction, Inverse Problems and Imaging, 11 (2017), pp. 577–600.