Primal-dual method to smoothing TV-based model for image denoising

Abstract

The total variation-based Rudin–Osher–Fatemi model is an effective and popular prior model in the image processing problem. Different to frequently using the splitting scheme to directly solve this model, we propose the primal dual method to solve the smoothing total variation-based Rudin–Osher–Fatemi model and give some convergence analysis of proposed method. Numerical implements show that our proposed model and method can efficiently improve the numerical results compared with the Rudin–Osher–Fatemi model.

Keywords

Smoothing total variation primal dual method image denoising convex conjugation

Introduction

The degradation of an image is usually unavoidable during its acquisition, and perhaps the most fundamental image restoration problem is denoising. It forms a significant preliminary step in many machine vision tasks such as object detection and recognition. A major concern in designing image denoising models is to preserve important image features such as edges while removing noise. Variational PDE based models^1–3 have been widely used over the past decades for image denoising with edge preservation. A classical variational denoising method is the total variation (TV)-based minimizing process proposed by Rudin–Osher–Fatemi (ROF),⁴ where they used the properties that TV can reduce oscillations and regularize the geometry of level sets without penalizing discontinuities. Formally, the ROF can be written as

min_{u} \frac{λ}{2} \int Ω (u - f) 2 dx + | u | TV

(1)

where the regularization weight

λ > 0

should be tuned to match the noise level contaminating f and

| u | TV : = \int Ω | \nabla u | dx

can be extended to the space

BV (Ω)

that contains functions with discontinuities. Here, Ω denotes the image domains.

However, from the view of numerical computation, solving (1) is a difficult task to find a good method since its corresponding Euler–Lagrangian equation includes the term $\frac{\nabla u}{| \nabla u |}$ , which is degenerated when $| \nabla u | = 0$ . Especially, this fact usually happens in the flat areas of image. A lot of methods have been proposed to overcome this drawback. The first class of techniques is to slightly perturb the TV-norm functional to get

min_{u} \underset{E (u, α)}{\underset{︸}{\frac{λ}{2} \int Ω (u - f) 2 dx + \int Ω \sqrt{| \nabla u | 2 + α} dx}}

(2)

where α is a small positive parameter, then to solve the associated Euler–Lagrange equation,

g (u) : = λ (u - f) - div (\frac{\nabla u}{\sqrt{| \nabla u | 2 + α}}) = 0

(3)

with the homogeneous Neumann boundary condition. A number of methods have been proposed to solve the equation (3). Rudin et al.⁴ used a time marching scheme to reach the steady state of the parabolic equation

u t = g (u)

with the initial condition f. It has been proved to be quite efficient for regularizing images without smoothing the boundaries of the objects. However, the convergence is usually slow due to the Courant–Friedrichs–Lewy (CFL) condition, especially in flat regions where

| \nabla u | \approx 0

. To completely get rid of CFL conditions, Vogel and Oman⁵ proposed a lagged diffusion fixed point iteration scheme (FP) to solve the stationary Euler–Lagrange directly, where the

(n + 1)

-th iteration is obtained by solving the sparse linear equation

λ (u n + 1 - f) - div (\frac{\nabla u n + 1}{\sqrt{| \nabla u n | 2 + α}}) = 0 \Rightarrow u n + 1 = (λ I + L (u n)) - 1 λ f

by defining

L (u n) : = - \frac{\nabla}{\sqrt{| \nabla u n | 2 + α}}

. They showed that the fixed point iteration converged to the minimizer of the functional from the ROF model. However, inverse matrices still need to be computed in this method. There are some other methods based on using the dual problem of the original problem (1) such as semi-gradient and the primal dual method (PDM)^6–8 or the splitting scheme to overcome the nonsmoothing operator norm

| \nabla \cdot | 1

such as the alternating direction method of multipliers (ADMM).^9–12

This article focuses on the problem (2) by assuming $0 \leq α \leq 1$ . It is easy to find that $\int Ω \sqrt{| \nabla u | 2 + α} dx$ as $α \to 0$ tends to the original TV norm and $\int Ω \sqrt{| \nabla u | 2 + α} dx$ as $α \to 1$ tends to the surface area $\int Ω \sqrt{| \nabla u | 2 + 1} dx$ . Furthermore, using the Taylor approximation $\sqrt{1 + x} \approx 1 + x / 2$ , we have

\int Ω \sqrt{| \nabla u | 2 + α} dx \approx \int Ω (\sqrt{α} + | \nabla u | 2 / 2) dx

(4)

Then the Euler–Lagrange equation of problem (2) approximates to Laplace equation in flat region of image where $| \nabla u | \approx 0$ . In this case, the problem (2) can reduce some circumstances such as staircasing effect. Actually the problem (2) was introduced in literature,^13,14 where they considered the image as embedded maps and minimal surfaces. Then they solved it in the sense of Laplace–Beltrami operator equation. Different from these aforementioned methods, we here propose the PDM to solve the problem (2). In this method, each iteration updates both a primal and a dual variable. It is thus able to avoid some of the difficulties that arise when working only on the primal or dual side.^4,8 Therefore, we first employ the Legrndre–Fenchel transform to transform the problem (2) to a saddle point problem. Then, we can solve it based on the PDM. Based on the triviality following the work in Chambolle and Pock,⁸ the convergence of the proposed PDM can be kept. Our numerical results demonstrate that the algorithm instances achieve the solution of difficult convex optimization problems in reasonable time and without parameter tuning of algorithm.

The contents of the article are organized as follows. In “PDM” section, we give some preliminaries of the PDM and use it to solve our proposed model. We arrange some numerical implementations in “Numerical implementations” section. Some concluding remarks are presented in “Conclusions” section.

PDM

PDM was first presented by Arrow et al.¹⁵ and named as the primal-dual hybrid gradient (PDHG) method in Arrow et al.¹⁵ and Zhu et al.¹⁶ This method enjoys to solve the primal variable and the dual variable alternatively and then widely applied to compute the saddle point of min–max problems.^17–19 A gradient descent method is employed to find the primal and the dual variables alternatively in Zhu et al.²⁰ and Esser et al.²¹ A predictor-corrector scheme is used in the alternating direction iterations for finding the dual variable in Chen and Teboulle.²² Recently, Chambolle and Pock⁸ presented an unified form of PDM. They demonstrated that, with a properly specified step-size policy and an averaging scheme, this method can achieve the $O (1 / N)$ rate of convergence, where $N$ is total number of iterations. In this section, we will apply the algorithm in the problem (2) to seek the saddle point of our min–max problem. We therefore first give a brief introduction of the method here and then propose it to solve the problem (2).

PDM

Let X and Y be two finite-dimensional real vector spaces equipped with the Euclidian inner product $〈 \cdot, \cdot 〉$ and the norm $‖ \cdot ‖ = 〈 \cdot, \cdot 〉 \frac{1}{2}$ . Define the operator norm as

‖ A ‖ = max {‖ A x ‖, x \in X with ‖ x ‖ \leq 1}

where

A : X \to Y

Assume that functions $F : X \to ℝ \cup {\infty}$ and $G : Y \to ℝ \cup {\infty}$ are proper convex function, we consider the following problem:

min_{x} F (x) + G (Ax)

(5)

where A is a matrix operator. Based on the definition of convex conjugation as

H * (s) = sup x {〈 s, t 〉 - H (t)}

for a function

H : X \to ℝ \cup {\infty}

, the dual problem of the problem (5) can be written as

max_{y} - F * (y) - G * (- A T y)

(6)

Similar to the work in Chambolle and Pock,⁸ the primal problem (5) and its dual problem (6) can be connected by a generic saddle-point problem as

min_{x} max_{y} F (x, y) : = 〈 Ax, y 〉 + F (x) - G * (y), \forall x \in X, y \in Y

(7)

Furthermore, the above fact holds if and only if $A T y * \in \partial F (x *)$ or also $Ax * \in \partial G (y *)$ . Thus, the problem (5) can be solved by the primal dual algorithm⁸

{x k + 1 = (I + τ \partial F) - 1 (x k - τ A T y k), (8 a) \bar{x} k = 2 x k + 1 - x k, (8 b) y k + 1 = (I + σ \partial G *) - 1 (y k - σ A \bar{x} k), (8 c)

where

τ > 0

and

σ > 0

are step sizes of the primal and dual variables, respectively.

Obviously, the strategy (8a)–(8c) can be regarded as using the gradient descent method to solve the primal problem

min_{x} 〈 Ax, y 〉 + F (x)

and using the gradient ascent method to solve the dual problem

max_{y} 〈 Ax, y 〉 - G * (y)

with the modified step (8b), where the step length are chosen as τ and σ corresponding to the variables x and y, respectively.

Theorem 2.1 Assume that the min–max problem (7) has a saddle point $(x *, y *)$ and $τ σ ‖ A ‖ 2 < 1$ , then the sequence ${(x n, y n)}$ generated by the strategy (8a)–(8c) converges to $(x *, y *)$ by choosing some suitable original values $(x 0, y 0) \in X \times Y$ and $\bar{x} = x 0$ .

For the proof, see Theorem 1 in Chambolle and Pock.⁸ Obviously, the PDM is matrix inversion-free if the functions F(x) and $G * (y)$ are separable for the variables x and y, respectively. So this method can be used to solve large scale problem such as machine learning and medical image processing problem.^23,24

The problem (2) based on the PDM

In this section, we consider solving the problem (2) based on the PDM. For simplicity, we assume that the image region Ω is squared. The problem (2) is computed numerically by projecting an image of $N = n \times n$ pixels. Set related Euclidean spaces $S = R n \times n$ and $T = S \times S$ . The usual scalar product can be denoted as $〈 v, w 〉 S : = \sum_{i = 1}^{n} \sum_{j = 1}^{n} v i, j w i, j$ for $v, w \in S$ . If $u \in S$ , we use $\nabla u = (\nabla_{x}^{+} u, \nabla_{y}^{+} u) \in T$ to denote the first-order forward difference operators with and use $div (p, q) = \nabla_{x}^{-} p + \nabla_{y}^{-} q$ to denote the first-order backward difference operators with

We also use the following notational brevities to represent $ℓ 1$ and $ℓ \infty$ norms as:

| \nabla u i, j | 1 = \sqrt{\nabla x u_{i, j}^{2} + \nabla y u_{i, j}^{2}} and | q | \infty = max_{1 \leq i, j \leq n} {| q i, j |}

for

| q i, j | = \sqrt{(q_{i, j}^{x}) 2 + (q_{i, j}^{y}) 2}

with

q \in T

. Based on these notations, we denote the discrete equivalent of (2) as

min_{u} {\sum_{i = 1}^{n} \sum_{j = 1}^{n} (\frac{λ}{2} (u i, j - f i, j) 2 + \sqrt{| \nabla u i, j | 2 + α})}

(8)

where

| \nabla u i, j | 2 = (\nabla_{x}^{+} u i, j) 2 + (\nabla_{y}^{+} u i, j) 2

. To solve the problem (5), we set X = S, Y = T, and

F (u) = \sum_{i = 1}^{n} \sum_{j = 1}^{n} \frac{λ}{2} (u i, j - f i, j) 2 and G (Au) = \sum_{i = 1}^{n} \sum_{j = 1}^{n} \sqrt{| \nabla u i, j | 2 + α}

where

G (\cdot) = | \cdot | 1

and

A = \nabla

. Using the property of the convex conjugate of a proper convex function,²⁵ as

g (t) = \sqrt{α + t 2} \Leftrightarrow g * * (t) = sup | s | \leq 1 {〈 t, s 〉 + \sqrt{α (1 - s 2)}}

then the problem (9) can be written to the min–max problem

min_{u} max_{| p | \infty \leq 1} {\sum_{i = 1}^{n} \sum_{j = 1}^{n} (\frac{λ}{2} (u i, j - f i, j) 2 - 〈 u i, j, div p i, j 〉 + \sqrt{α (1 - | p i, j | 2)})}

(9)

where

p i, j = (p_{i, j}^{x}, p_{i, j}^{y}) \in T

. Since the subjective function in (10) is proper convex, we can transfer the order of min and max. Thus, we can use the primal dual scheme (8) to solve it.

For the original variation u: In the problem (10), ignoring the unrelated term to u, we can obtain

min_{u} {\sum_{i = 1}^{n} \sum_{j = 1}^{n} (\frac{λ}{2} (u i, j - f i, j) 2 - 〈 u i, j, div p i, j 〉)}

whose optimality condition is

λ (u i, j - f i, j) - div p i, j = 0

for $1 \leq i, j \leq n$ . Then, we can use gradient ascent method to solve it as

u_{i, j}^{k} - u_{i, j}^{k + 1} = τ (λ (u_{i, j}^{k + 1} - f i, j) - div p i, j) \Rightarrow u_{i, j}^{k + 1} = \frac{u_{i, j}^{k} + λ τ f i, j + τ div p i, j}{1 + λ τ}

(10)

For the dual variation p: In the problem (10), ignoring the unrelated term to p and introducing an indicator function

χ K (p i, j) = {0, if p i, j \in K + \infty if p i, j \neg\in K

where

K : = {| p i, j | \infty \leq 1}

, we then have

which optimality condition is

\nabla u i, j \sqrt{1 - | p i, j | 2} - \sqrt{α} p i, j + \partial χ K (p i, j) = 0

for

1 \leq i, j \leq n

. Note that

\partial χ K (p i, j) = 0

since indicator function is a constant function. Here partial denotes the subgradient defined by

\partial ℎ (\bar{y}) : = {v | ℎ (\bar{x}) - ℎ (y) \geq (v, \bar{x} - \bar{y})}

at the point

\bar{y}

for a function ℎ. Then we use the projection gradient method as

p_{i, j}^{k + 1} = \frac{p_{i, j}^{k} + σ (\nabla u i, j \sqrt{1 - | p_{i, j}^{k} | 2} - \sqrt{α} p_{i, j}^{k})}{max {1, | p_{i, j}^{k} + σ (\nabla u i, j \sqrt{1 - | p_{i, j}^{k} | 2} - \sqrt{α} p_{i, j}^{k}) |}}

(11)

to solve it.

Based on two above subproblems, we summarize the PDM to solve the problem (9) as follows.

Algorithm 2.1. PDM to the problem (9)

Set original values $u 0 = f$ and $p 0 = 0$ . Let $k : = 0$ .

Repeat

{u_{i, j}^{k + 1} = \frac{u_{i, j}^{k} + λ τ f i, j + τ div p_{i, j}^{k}}{1 + λ τ} {\bar{u}}_{i, j}^{k} = 2 u_{i, j}^{k + 1} - u_{i, j}^{k} p_{i, j}^{k + 1} = \frac{p_{i, j}^{k} + σ (\nabla {\bar{u}}_{i, j}^{k} \sqrt{1 - | p_{i, j}^{k} | 2} - \sqrt{α} p_{i, j}^{k})}{max {1, | p_{i, j}^{k} + σ (\nabla {\bar{u}}_{i, j}^{k} \sqrt{1 - | p_{i, j}^{k} | 2} - \sqrt{α} p_{i, j}^{k}) |}}

(12)

Until Convergence.

Remark 2.1 Using the gradient method to compute the variables u and p in (11) and (12) actually corresponds to the original primal dual scheme (8). For example, when we use the gradient ascent method to obtain the original variable $u i, j$ in (11), it is also the optimization condition of the problem

Thus, we easily find to the formulation (8a) as

u_{i, j}^{k + 1} = (I + τ \nabla F) - 1 (u i, j - (u_{i, j}^{k} + τ div p_{i, j}^{k}))

(13)

For the dual variable $p i, j$ in (12), we also obtain the similar formulation as (8c) by setting $G * (p i, j) : = \sqrt{α (1 - | p i, j | 2)} + χ K (p i, j)$ .

Remark 2.2 To the problem (10), it is easy to obtain that it has an unique saddle $(u, p)$ since the subjective function in (9) is strict convex and coercive. For the PDM, it is very important to compute the operator norm of ∇ to keep its convergence. Actually, it is not difficult to find that $‖ \nabla ‖ 2 2 \leq 8$ .

Based on the above remarks, we have the following assertion for Algorithm 2.1.

Theorem 2.2 Assume that $τ σ < 1 / 8$ in Algorithm 2.1, then the sequence ${(u k, p k)}$ converges to the saddle point $(u *, p *)$ . Furthermore, $u *$ is the solution of the problem (9).

Numerical implementations

In numerical implementations, we will terminate the iteration if the stop conditions is satisfied in term of the $ℓ 2$ relative error $max {‖ u k + 1 - u k ‖ / ‖ u k ‖ 2, ‖ p k + 1 - p k ‖ 2 / ‖ p k ‖ 2} < 10 - 5$ or the iteration arriving to 500. To standardize discussions, we begin by normalizing pixel values of the inputting image $\bar{f}$ to [0,255]. We do this by using the linear-stretch formula as $f = 255 \times (\bar{f} - min (\bar{f})) | / (max (\bar{f}) - min (\bar{f}))$ , where max and min, respectively, represent maximum and minimum of $\bar{f}$ . The simulations are preformed in Matlab R20014a on a PC with an Intel Core i5 M520 at 2.40 GHz and 4 GB of memory. Furthermore, to obtain noise images, we use the matlab function “randn” to add the white Gaussian noise to the original clean image.

For numerical comparisons, we consider to compare our proposed PDM with the numerical methods such as the time marching method (TM)⁴ and the lagged diffusion fixed point method (FPM) in Vogel and Oman⁵ to solve the smoothing TV-based ROF model.¹⁵ Actually, our PDM degenerates to the PDM used in Chambolle and Pock⁸ to solve the ROF model while α = 0. Apart from some default settings, like the maximum number of iterations, the related parameters in the codes are chosen by trials and errors to give the best results of the respective methods. Here, we choose three different images with the size 256 × 256 as testing images shown in Figure 1. To measure more quantitatively, the similarity between the recovered image and the original (supposedly noise-free) image, we use both the signal-to-noise ratio (SNR) and the structural similarity index (SSIM), the latter being well known to better reflect perceived visual quality.

Figure 1.

Original images in our implementations. (a) Cameraman (b) Lena (c) Synthetic image.

Examples 3.1 Here we consider to use the PDM, FDM, and TM to solve the smoothing TV-based model (2) by setting $α = 0.01$ and adding white Gaussian noise to Figure 1(a) with three different deviations as 10, 15, and 20 ( Table 1 ). From convergence curves in Figure 2 , it can be noticed that all methods can efficiently suppress noise. However, the energy curve generated by the TM need much more iteration number to attain steady state due to the effect by the CFL condition. Simultaneously, we also notice that the quality of restored images are reduced and the regularization parameter λ in the model (2) needs to be decreased while increasing the noise in the image. That is because that the capability of reconstructing becomes harder and the regularization term $\int Ω \sqrt{| \nabla u | 2 + α} dx$ to be dominance. From the averages of CPU time, our proposed method is obviously superior to other two methods, especially for the large scale image such in image denoising problem in the medical image. Simultaneously, we notice that our method will give the best SNR and SSIM when choosing a suitable stopping condition by comparing the convergence curves in the second and third rows of Figure 2 .

Figure 2.

Numerical curves of $SNR$ , $SSIM$ , and the energy of the problem (2). Top-row denotes the original images by adding to different white Gaussian noise.

Table 1.

The related results in Example 3.1.

(σ, λ)	(10, 0.08)			(15, 0.175)			(20, 0.032)
Method	SNR	SSIM	Time	SNR	SSIM	Time	SNR	SSIM	Time
TM	19.1310	0.5056	64.7656	17.3887	0.4571	90.1563	15.8731	0.4037	78.9375
PDM	20.5998	0.5519	50.2969	17.8345	0.4559	50.0781	14.1136	0.3080	42.7813
FPM	20.5990	0.5523	32.0313	17.7935	0.4542	49.2656	14.0446	0.3081	78.8594

Here (σ, λ) denotes the level of noise and regularization parameter in the problem (2).

Examples 3.2 We consider the effect while choosing the different parameter α. Here the Lena image as the testing image is added to different noise levels with different deviations shown in Table 2 . As we can see from the data, it is easy to find that the smoothing TV-based model can improve the restoration quality then the results generated by using suitable α. Note that it is of using the PDM to solve the original ROF model when setting α = 0. Furthermore, the regularization parameter λ is need to be reduced when we increase the smoothing parameter α. This is because that we need a higher λ to penalty

\int Ω \sqrt{| \nabla u | 2 + α} dx

while decreasing α.

Table 2.

The related results in Example 3.2.

Noise	σ=10			σ=20			σ=30
Parameter	SNR	SSIM	λ	SNR	SSIM	λ	SNR	SSIM	λ
$α = 0$	19.1983	0.7457	0.18	15.4873	0.6320	0.07	13.3730	0.5478	0.039
$α = 0.001$	19.2249	0.7492	0.18	15.4639	0.6413	0.07	13.3582	0.5594	0.039
$α = 0.01$	19.2005	0.7448	0.18	15.4999	0.6320	0.07	13.3768	0.5495	0.038
$α = 0.1$	19.2091	0.7469	0.17	15.5119	0.6306	0.07	13.3844	0.5522	0.036
$α = 1$	19.2206	0.7474	0.15	15.5275	0.6338	0.06	13.4102	0.5573	0.030
$α = 1.5$	19.2249	0.7492	0.14	15.4639	0.6413	0.05	13.3582	0.5594	0.026

Examples 3.3 Since the original ROF model (1) usually generates the staircase effect in the smoothing region while restoring degraded image. Simultaneously, from the introduction, we can find the smoothing TV-based ROF model (2) approximately own the diffusion ability with isotropy, that is to say $\int Ω \sqrt{| \nabla u | 2 + α} dx \approx \int Ω (1 + | \nabla u | 2 / 2) dx$ , when $| \nabla u | \approx 0$ . So we expect that the model (2) can efficiently reduce staircase effect and keep image edges when choosing suitable the smoothing parameter α. This fact has been observed from Example 3.2 by choosing α = 1. Here we mainly focus on reducing the staircase effect by choosing suitable α, and numerical results are arranged in Figure 3 . To clear comparisons, we also plot the slices of the 60th in the restored image. Obviously, the smoothing TV-based model (2) can reduce the staircase effect compared with the ROF model (2) and keep image edges, this phenomenon is especially obvious in peak regions and the corner regions in Figure 3(c) and ( d ).

Figure 3.

Numerical results by using the models (1) and (2) to solve the noisy images. To clear comparison, we only show the upper-left part of restoration images with size $128 \times 128$ and plot the slice of the 64th column. First to third columns are added to the white Gaussian noise with three different deviations, such as 10, 20, and 30.

Conclusions

In this article, we have introduced the PDM based on the work in Chambolle and Pock⁸ to solve the smoothing TV-based ROF model (2). Numerical implementations show that our method is superior to the classic methods such as the time marching method⁴ and the FPM⁵ on the stability and high efficiency, especially to restore the large-scale image. Furthermore, the proposed model can effectively reduce the staircase effect by choosing suitable smoothing parameter α.

Footnotes

Acknowledgments

The author is grateful to the anonymous referee for his/her careful checking of details and helpful comments that improved this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: NSF of China (No. 11401170) and Foundation of Henan Educational Committee of China (No.14A110018).

References

Aubert G and Kornprobst P. Mathematical problem in image processing: partial differential equations and the calculus of variations. New York, NY: Springer, 2008.

Chan

Shen

. Image processing and analysis-variational, PDE, wavelet, and stochastic methods, Philadelphia: SIAM, 2005.

Scherzer

. Handbook of mathematical methods in imaging, New York: Springer, 2015.

Rudin

Osher

Fatemi

. Nonlinear total variation based noise removal algorithms. Physica D 1992; 60: 259–268.

Vogel

Oman

. Iterative methods for total variation denoising. SIAM J Sci Comput 1996; 17: 227–238.

Chan

Golub

Mulet

. A nonlinear primal-dual method for total variation-based image restoration. SIAM J Sci Comput 1999; 20: 1964–1977.

Chambolle

. An algorithm for total variation minimization and applications. J Math Imag Vis 2004; 20: 89–97.

Chambolle

Pock

. A first-order primal-dual algorithm for convex problems with applications to imaging. J Math Imag Vis 2011; 40: 120–145.

Boyd

Parikh

Chu

. Distributed optimization and statistical learning via the alternating direction method of multiplier. Found Trends Mach Learn 2011; 3: 1–122.

10.

Tai

. Augmented Lagrangian method, dual methods, and split Bregman iteration for ROF, vectorial TV, and high order models. SIAM J Imag Sci 2010; 3: 300–339.

11.

Setzer S. Split Bregman algorithm, Douglas-Rachford splitting and frame shrinkage. In: Proceedings of the second international conference on scale space and variational methods in computer vision, 2009, Voss, Norway, vol. 5567, pp.464–476.

12.

Goldstein

Osher

. The split Bregman method for L1 regularized problems. SIAM J Imag Sci 2009; 2: 323–343.

13.

Malladi

Sethian

. Image processing: Flows under min/max curvature and mean curvature. Graph Models Image Process 1996; 58: 127–141.

14.

Fallah

Ford

. On mean curvature diffusion in nonlinear image filtering. Pattern Recog Lett 1998; 19: 433–437.

15.

Arrow

Hurwicz

Uzawa

. Studies in linear and non-linear programming, Stanford, CA: Stanford University Press, 1958.

16.

Zhu

Wright

Chan

. Duality-based algorithms for total-variation-regularized image restoration. Comput Opt Appl 2010; 47: 377–400.

17.

Combettes

Pesquet

. A proximal decomposition method for solving convex variational inverse problems. Inverse Prob 2008; 24: 065014.

18.

Dupe

Fadili

Starck

. A proximal iteration for deconvolving Poisson noisy images using sparse representations. IEEE Trans Image Process 2009; 18: 310–321.

19.

Eckstein

Bertsekas

. On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators. Math Prog A 1992; 55: 293–318.

20.

Zhu

Wright

Chan

. Duality-based algorithms for total-variation-regularized image restoration. Comput Opt Appl 2010; 47: 377–400.

21.

Esser

Zhang

Chan

. A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science. SIAM J Imag Sci 2010; 3: 1015–1046.

22.

Chen

Teboulle

. A proximal-based decomposition method for convex minimization problems. Math Prog A 1994; 64: 81–101.

23.

Sidky

Jøgensen

Pan

. Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm. Phys Med Biol 2012; 57: 3065–3091.

24.

Zhang Y and Xiao L. Stochastic primal-dual coordinate method for regularized empirical risk minimization. In: The 32nd international conference on machine learning, Lille, France, 2015.

25.

Ekeland I and Turbull T. Infinite dimensional optimization and convexity. Chicago, IL: The University of Chicago Press, 1983.