An Efficient Total Variation Minimization Method for Image Restoration

Abstract

In this paper, we present an effective algorithm for solving the Poisson–Gaussian total variation model. The existence and uniqueness of solution for the mixed Poisson–Gaussian model are proved. Due to the strict convexity of the model, the split-Bregman method is employed to solve the minimization problem. Experimental results show the effectiveness of the proposed method for mixed Poisson–Gaussion noise removal. Comparison with other existing and well-known methods is provided as well.

1 Introduction

Image acquisition is an ubiquitous technology, found for example in photography, medical imagery, astronomy, etc. Nevertheless, in almost all situations, the image-capturing devices are imperfect: some unwanted noise is added to the signal. Therefore, the obtained images are post-processed by numerical algorithms before being delivered to the users; those algorithms have to solve the image restoration problem.

In the image restoration problem, an original image u is corrupted by some random noise η, resulting in a noisy image f. Our task is to reconstruct u, knowing both f and the distribution of η. Of course, there is in general no way to find the exact image u; image restoration algorithms rather yield a good approximation of u, usually noted ${\mathit{u}}^{\ast }$ . To do so, they exploit a priori knowledge on the image u.

Various distributions have been considered for the noise, e.g. Gaussian (Rudin et al., 1992; Pham and Kopylov, 2015), Poisson (Chan and Shen, 2005; Le et al., 2007), Cauchy (Sciacchitano et al., 2015), as well as some mixed noise models, e.g. mixed Gaussian-Impulse noise (Yan, 2013), mixed Gaussian–Salt and Pepper noise (Liu et al., 2017), mixed Poisson–Gaussian (Calatroni et al., 2017; Pham et al., 2018; Tran et al., 2019).

A growing interest in Poisson–Gaussian probabilistic models has recently arisen (Chouzenoux et al., 2015). The mixture of Poisson and Gaussian noise occurs in several practical setups (e.g. microscopy, astronomy), where the sensors used to capture images have two sources of noise: a signal-dependent source which comes from the way light intensity is measured; and a signal-independent source which is simply thermal and electronic noise. Gaussian noise is just additive, so it cannot properly approximate the Poisson–Gaussian distributions observed in practice, which are strongly signal-dependent.

In general, the mixed Poisson–Gaussian noise model can be expressed as follows: $$\begin{array}{c}\mathit{f}=\mathcal{P}\mathrm{(}\mathit{u}\mathrm{)}+\mathit{W}\mathrm{,}\end{array}$$(1) where f is observed image, u is the unknown image, $\mathcal{P}\mathrm{(}\mathit{u}\mathrm{)}$ means that the image u is corrupted by Poisson noise, and $\mathit{W}\sim \mathcal{N}\mathrm{(}0\mathrm{,}{\mathit{\sigma }}^2\mathrm{)}$ is a Gaussian noise with zero mean and variance σ.

Recently, several approaches have been devoted to the mixed Poisson–Gaussian noise model (Foi et al., 2008; Jezierska et al., 2011; Lanza et al., 2014; Le Montagner et al., 2014). Many algorithms for denoising images corrupted by mixed Poisson–Gaussian noise have been investigated using approximations based on variance stabilization transforms (Zhang et al., 2007; Makitalo and Foi, 2013) or PURE-LET based approaches (Luisier et al., 2011; Li et al., 2018). Variational models based on the Bayesian framework have been also proposed for removing and denoising and deconvolution of mixed Poisson–Gaussian noise (Calatroni et al., 2017). This framework is perhaps a popular approach to mixed Poisson–Gaussian noise model. Authors in De Los Reyes and Schönlieb (2013) proposed a nonsmooth PDE-constrained optimization approach for the determination of the correct noise model in total variation image denoising. Authors in Lanza et al. (2014) focused on the maximum a posteriori approach to derive a variational formulation composed of the total variation (TV) regularization term and two fidelities. A weighted squared ${\mathit{L}}_2$ norm noise approximation was proposed for mixed Poisson–Gaussian noise in Li et al. (2015), or an efficient primal-dual algorithm was also proposed in Chouzenoux et al. (2015) by investigating the properties of the Poisson–Gaussian negative log-likelihood as a convex Lipschitz differentiable function. Recently, authors in Marnissi et al. (2016) proposed a variational Bayesian method for Poisson–Gaussian noise, using an exact Poisson–Gaussian likelihood. Similarily, authors in Calatroni et al. (2017) proposed a variational approach which includes an infimal convolution combination of standard data delities classically associated to one single-noise distribution, and a TV regularization as regularizing energy. Generally, image restoration by variational models based on TV can be a good solution to the mixed Poisson–Gaussian noise removal with the following formula (Calatroni et al., 2017; Pham et al., 2019): $$\begin{array}{c}{\mathit{u}}^{\ast }=\underset{\mathit{u}\in \mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}}{arg\; min}{{\displaystyle \int }}_{\mathrm{\Omega }}|\nabla \mathit{u}|+{\displaystyle \frac{{\mathit{\lambda }}_1}{2}}{{\displaystyle \int }}_{\mathrm{\Omega }}{\mathrm{(}\mathit{u}-\mathit{f}\mathrm{)}}^2+{\mathit{\lambda }}_2{{\displaystyle \int }}_{\mathrm{\Omega }}\mathrm{(}\mathit{u}-\mathit{f}log\mathit{u}\mathrm{)}\mathrm{,}\end{array}$$(2) where f is the observed image, $\mathrm{\Omega }\subset {\mathbb{R}}^2$ is a bounded domain, and $\mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ is the set of positive functions from Ω to $\mathbb{R}$ ; finally, ${\mathit{\lambda }}_1\mathrm{,}{\mathit{\lambda }}_2$ are positive regularization parameters (see Chan and Shen, 2005, for details on this method).

However, in some cases, intermediate solutions of (2) obtained during the execution of algorithms may contain pixels with negative values. To avoid this problem, authors in Pham et al. (2018) proposed a modified scheme of gradient descent (MSGD) that guarantee positive values for each pixel in the image domain.

In this work, we focus on the model (2) and consider the following model: $$\begin{array}{rl}{\mathit{u}}^{\ast }& =\underset{\mathit{u}\in \mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}}{arg\; min}\mathit{E}\mathrm{(}\mathit{u}\mathrm{)}\mathrm{,}\\ \mathit{E}\mathrm{(}\mathit{u}\mathrm{)}& ={{\displaystyle \int }}_{\mathrm{\Omega }}\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}|\nabla \mathit{u}\mathrm{(}\mathit{x}\mathrm{)}|\mathit{d}\mathit{x}+{\displaystyle \frac{{\mathit{\lambda }}_1}{2}}{{\displaystyle \int }}_{\mathrm{\Omega }}{\mathrm{(}\mathit{u}\mathrm{(}\mathit{x}\mathrm{)}-\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}\mathrm{)}}^2\mathit{d}\mathit{x}\\ & +{\mathit{\lambda }}_2{{\displaystyle \int }}_{\mathrm{\Omega }}\mathrm{(}\mathit{u}\mathrm{(}\mathit{x}\mathrm{)}-\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}log\mathit{u}\mathrm{(}\mathit{x}\mathrm{)}\mathrm{)}\mathit{d}\mathit{x}\mathrm{,}\end{array}$$(3) where f is the observed image, ${\mathit{\lambda }}_1$ and ${\mathit{\lambda }}_2$ are positive regularization parameters, $\mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}=\{\mathit{u}\in \text{BV}\mathrm{(}\mathrm{\Omega }\mathrm{)}:\mathit{u}\mathrm{>}0\}$ is closed and convex, with $\text{BV}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ being the space of functions $\mathrm{\Omega }\to \mathbb{R}$ with bounded variation; and finally $\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}$ is a continuous function in $\mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ .

The function $\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}$ is used to control the intensity of the diffusion, which is an edge indicator for spatially adaptive image restoration (Barcelos and Chen, 2000). Typically, the function $\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}$ is chosen as follows: $$\begin{array}{c}\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}={\displaystyle \frac{1}{1+\mathit{l}·|\mathit{v}\mathrm{(}\mathit{x}\mathrm{)}{|}^2}}\mathrm{,}\end{array}$$ where l is a threshold value and $\mathit{v}\mathrm{(}\mathit{x}\mathrm{)}=|\nabla {\mathit{G}}_{\mathit{\sigma }}\mathrm{(}\mathit{x}\mathrm{)}\ast \mathit{f}|$ , in which ∗ denotes the convolution with ${\mathit{G}}_{\mathit{\sigma }}\mathrm{(}\mathit{x}\mathrm{)}=\frac{1}{2\mathit{\pi }{\mathit{\sigma }}^2}exp\mathrm{(}-\frac{{\mathit{x}}^2}{2{\mathit{\sigma }}^2}\mathrm{)}$ , i.e. the Gaussian filter with standard deviation σ.

The main contributions of this paper are the following. We give an elementary proof of the existence and uniqueness of model (3). Moreover, we check that the functional $\mathit{E}\mathrm{(}·\mathrm{)}$ is convex, which enables us to use larger time-step parameters during gradient descent when solving (3). We introduce the influence function $\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}$ , which acts as an edge-detection function, to get the model (3) in order to improve the ability of edge preservation and to control the speed of smoothing. In addition, we propose a new method to solve (3) that perceptibly improves the quality of the denoised images. By changing the time-step parameter, users can either get faster denoising with comparable results to previous methods, or better quality denoising with comparable running times. Our method is a technical improvement over the split-Bregman algorithm. We report experimental results for the aforementioned method, for various parameters in the noise distribution. The quality of denoising is measured with the SSIM and PSNR metrics. If we tune the time-step parameter to get similar quality result as the original split-Bregman method, we get faster running times.

The rest of the paper is organized as follows. In Section 2, we describe the Poisson–Gaussian model and introduce the notation used in this work. In Section 3, we prove the existence and uniqueness of the solution. In Section 4, using the split-Bregman algorithm, we present the proposed optimization framework. Next, in Section 5, we show some numerical results of our proposed method and we compare them with the results obtained with other existing methods. Finally, some conclusions are drawn in Section 6.

2 Preliminaries

We recall the principle behind equation (2). Note that the contents of this section are not a rigorous proof; we simply provide a bit of context around the equation, why it was considered in the first place, and one possible reason for its practical efficiency. We also state our assumptions on both the initial image and the noise along the way.

Our goal is to recover the original image u, knowing the noisy image f. Our strategy is to find the image u which maximizes the conditional probability $\mathit{P}\mathrm{(}\mathit{u}|\mathit{f}\mathrm{)}$ . Bayes’s rule gives: $$\begin{array}{c}\mathit{P}\mathrm{(}\mathit{u}|\mathit{f}\mathrm{)}={\displaystyle \frac{\mathit{P}\mathrm{(}\mathit{f}|\mathit{u}\mathrm{)}\mathit{P}\mathrm{(}\mathit{u}\mathrm{)}}{\mathit{P}\mathrm{(}\mathit{f}\mathrm{)}}}.\end{array}$$(4) The probability density function of the observed image f corrupted by Gaussian noise ${\mathit{P}}_{\mathcal{N}}$ (respectively, by Poisson noise ${\mathit{P}}_{\mathcal{P}}$ ) is: $$\begin{array}{c}{\mathit{P}}_{\mathcal{N}}\mathrm{(}\mathit{f}|\mathit{u}\mathrm{)}={\displaystyle \frac{1}{\mathit{\sigma }\sqrt{2\mathit{\pi }}}}exp\mathrm{(}-{\displaystyle \frac{{\mathrm{(}\mathit{u}-\mathit{f}\mathrm{)}}^2}{2{\mathit{\sigma }}^2}}\mathrm{)}\mathrm{,}{\mathit{P}}_{\mathcal{P}}\mathrm{(}\mathit{f}|\mathit{u}\mathrm{)}={\displaystyle \frac{{\mathit{u}}^{\mathit{f}}exp\mathrm{(}-\mathit{u}\mathrm{)}}{\mathit{f}!}}\mathrm{,}\end{array}$$ where σ is the variance of the Gaussian noise. As we explained in the introduction, the two sources of noise are independent of each other, so the distribution of the mixed noise may be expressed as: $$\begin{array}{c}{\mathit{P}}_{\mathrm{mixed}}\mathrm{(}\mathit{f}|\mathit{u}\mathrm{)}={\displaystyle \frac{{\mathit{u}}^{\mathit{f}}exp\mathrm{(}-\mathit{u}\mathrm{)}}{\mathit{f}!}}{\displaystyle \frac{1}{\sqrt{2\mathit{\pi }}\mathit{\sigma }}}exp\mathrm{(}-{\displaystyle \frac{{\mathrm{(}\mathit{u}-\mathit{f}\mathrm{)}}^2}{2{\mathit{\sigma }}^2}}\mathrm{)}.\end{array}$$ We assume that the values of the pixels in the original image are independent, and that the noise is also independent on each pixel. (However, we do not assume that the noise and the original image are independent of each other.) Suppose that f has size $\mathit{M}\times \mathit{N}$ , and let $\mathit{I}=\{1\mathrm{,}\dots \mathrm{,}\mathit{M}\}\times \{1\mathrm{,}\dots \mathrm{,}\mathit{N}\}$ denote the domain of f. For i in I, we write ${\mathit{f}}_{\mathit{i}}$ the pixel of f at position i (and similarly ${\mathit{u}}_{\mathit{i}}$ the pixel of u at position i). Then, $$\begin{array}{c}{\mathit{P}}_{\mathrm{mixed}}\mathrm{(}\mathit{f}|\mathit{u}\mathrm{)}={\displaystyle \prod }_{\mathit{i}\in \mathit{I}}{\displaystyle \frac{{\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{)}}^{{\mathit{f}}_{\mathit{i}}}{\mathit{e}}^{\mathrm{(}-{\mathit{u}}_{\mathit{i}}\mathrm{)}}}{{\mathit{f}}_{\mathit{i}}!}}\mathit{s}exp\mathrm{(}-{\displaystyle \frac{{\mathrm{(}{\mathit{u}}_{\mathit{i}}-{\mathit{f}}_{\mathit{i}}\mathrm{)}}^2}{2{\mathit{\sigma }}^2}}\mathrm{)}\end{array}$$ with $\mathit{s}={\mathrm{(}\sqrt{2\mathit{\pi }}\mathit{\sigma }\mathrm{)}}^{-1}$ . Maximizing ${\mathit{P}}_{\mathrm{mixed}}$ is equivalent to minimizing $-log{\mathit{P}}_{\mathrm{mixed}}$ , so let us compute the quantity $-log\mathrm{(}{\mathit{P}}_{\mathrm{mixed}}\mathrm{(}\mathit{f}|\mathit{u}\mathrm{)}\mathrm{)}$ : $$\begin{array}{c}{\displaystyle \sum }_{\mathit{i}\in \mathit{I}}{\mathit{u}}_{\mathit{i}}-{\mathit{f}}_{\mathit{i}}log\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{)}+log\mathrm{(}{\mathit{f}}_{\mathit{i}}!\mathrm{)}+\mathit{y}{\mathrm{(}{\mathit{u}}_{\mathit{i}}-{\mathit{f}}_{\mathit{i}}\mathrm{)}}^2\mathrm{,}\end{array}$$(5) for some constant $\mathit{y}\mathrm{>}0$ . In the above equation, u varies but f is constant. Since our goal is to minimize the whole expression, we can ignore the term $log\mathrm{(}{\mathit{f}}_{\mathit{i}}!\mathrm{)}$ altogether.

Now we assume that $\mathit{P}\mathrm{(}\mathit{u}\mathrm{)}$ follows a Gibbs prior (Le et al., 2007): $$\begin{array}{c}\mathit{P}\mathrm{(}\mathit{u}\mathrm{)}={\displaystyle \frac{1}{\mathit{z}}}exp\mathrm{(}-{\displaystyle \int }|\nabla \mathit{u}|\mathrm{)}\mathrm{,}\end{array}$$(6) where z is a normalization factor. We need to make a couple of comments here. First, u is not a function ${\mathbb{R}}^2\to \mathbb{R}$ , but rather a discrete array of pixels; thus the integral in that expression is going to be translated to a sum, while $\nabla \mathit{u}$ will be translated as a linear approximation. Second, this assumption appears to contradict the previous one, that the pixels of the original image are independent of one another. However, the assumption on $\mathit{P}\mathrm{(}\mathit{u}\mathrm{)}$ is local: each pixel depends (weakly) on the neighbouring pixels only, so we do not lose much by assuming independence. This turns out to yield good results in practice (Chan and Shen, 2005).

We now have all the ingredients to maximize $\mathit{P}\mathrm{(}\mathit{u}|\mathit{f}\mathrm{)}$ . By equation (4), this amounts to minimize the expression $-log\mathrm{(}\mathit{P}\mathrm{(}\mathit{f}|\mathit{u}\mathrm{)}\mathrm{)}-log\mathrm{(}\mathit{P}\mathrm{(}\mathit{u}\mathrm{)}\mathrm{)}$ , so we can plug in equations (5) and (6) to get: $$\begin{array}{c}{\mathit{u}}^{\ast }=\underset{\mathit{u}}{arg\; min}{\displaystyle \sum }_{\mathit{i}\in \mathit{I}}{\displaystyle \frac{1}{\mathit{z}}}|\nabla {\mathit{u}}_{\mathit{i}}|+\mathit{y}{\mathrm{(}{\mathit{u}}_{\mathit{i}}-{\mathit{f}}_{\mathit{i}}\mathrm{)}}^2+\mathrm{(}{\mathit{u}}_{\mathit{i}}-{\mathit{f}}_{\mathit{i}}log\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{)}\mathrm{)}\mathrm{,}\end{array}$$(7) and we can view this expression as a discrete approximation of the functional $\mathit{E}\mathrm{(}·\mathrm{)}$ defined as: $$\begin{array}{c}\mathit{E}\mathrm{(}\mathit{u}\mathrm{)}={{\displaystyle \int }}_{\mathrm{\Omega }}|\nabla \mathit{u}|\mathit{d}\mathit{x}+{\displaystyle \frac{{\mathit{\lambda }}_1}{2}}{{\displaystyle \int }}_{\mathrm{\Omega }}{\mathrm{(}\mathit{u}-\mathit{f}\mathrm{)}}^2\mathit{d}\mathit{x}+{\mathit{\lambda }}_2{{\displaystyle \int }}_{\mathrm{\Omega }}\mathrm{(}\mathit{u}-\mathit{f}log\mathit{u}\mathrm{)}\mathit{d}\mathit{x}\mathrm{,}\end{array}$$(8) with ${\mathit{\lambda }}_1=2\mathit{y}\mathit{z}$ and ${\mathit{\lambda }}_2=\mathit{z}$ . (We multiplied by z, which is positive and constant, so the minimum is the same.) Intuitively, the last two terms are data fidelity terms, which ensure that the restored image u is not “too far” from the original image u (taking the distribution of the noise into account). By contrast, $|\nabla \mathit{u}|$ is a smoothness term, which guarantees that the reconstructed image is not too irregular (this is where our a priori knowledge on the original picture lies). The parameters ${\mathit{\lambda }}_1$ and ${\mathit{\lambda }}_2$ will have to be determined experimentally later on.

In the following sections, we introduce some theoretical results about the existence and uniqueness result for solution of (3).

3 Existence and Unicity of the Solution

Motivated by Aubert and Aujol (2008), Dong and Zeng (2013), we have the following existence and uniqueness results for the optimization problem (3). We prove that (3) has an unique solution in two steps: first, we show that $\mathit{E}\mathrm{(}·\mathrm{)}$ is a convex functional; then, we show that $\mathit{E}\mathrm{(}·\mathrm{)}$ has a lower bound. These two facts together imply the existence and uniqueness of the minimizer of $\mathit{E}\mathrm{(}·\mathrm{)}$ .

Theorem 1.

The functional $\mathit{u}\mapsto \mathit{E}\mathrm{(}\mathit{u}\mathrm{)}$ , where E is defined in (3), is strictly convex.

Proof.

Let us set: $\mathit{h}\mathrm{(}\mathit{u}\mathrm{)}=\frac{{\mathit{\lambda }}_1}{2}{\mathrm{(}\mathit{u}-\mathit{f}\mathrm{)}}^2+{\mathit{\lambda }}_2\mathrm{(}\mathit{u}-\mathit{f}log\mathit{u}\mathrm{)}$ . The first and the second order derivative of h are: $$\begin{array}{c}{\mathit{h}}^{\prime }\mathrm{(}\mathit{u}\mathrm{)}={\displaystyle \frac{{\mathit{\lambda }}_1{\mathit{u}}^2-\mathit{u}\mathrm{(}{\mathit{\lambda }}_1\mathit{f}-{\mathit{\lambda }}_2\mathrm{)}-{\mathit{\lambda }}_2\mathit{f}}{\mathit{u}}}\end{array}$$ and $$\begin{array}{c}{\mathit{h}}^{″}\mathrm{(}\mathit{u}\mathrm{)}={\displaystyle \frac{{\mathit{\lambda }}_1{\mathit{u}}^2+{\mathit{\lambda }}_2\mathit{f}}{{\mathit{u}}^2}}.\end{array}$$ Since f is a positive, and $\mathit{u}\in \mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ , we have: ${\mathit{h}}^{″}\mathrm{(}\mathit{u}\mathrm{)}\mathrm{>}0$ , i.e. $\mathit{h}\mathrm{(}\mathit{u}\mathrm{)}$ is strictly convex. Moreover, the TV regularization is convex, thence $\mathit{E}\mathrm{(}\mathit{u}\mathrm{)}$ is also strictly convex.  □

Theorem 2.

Let $\mathit{f}\in \mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}\cap {\mathit{L}}^{\infty }\mathrm{(}\mathrm{\Omega }\mathrm{)}$ , then the problem (3) has an exactly one solution $\mathit{u}\in \mathit{B}\mathit{V}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ and satisfying: $$\begin{array}{c}\underset{\mathrm{\Omega }}{inf}\mathit{f}⩽\mathit{u}⩽\underset{\mathrm{\Omega }}{sup}\mathit{f}.\end{array}$$

Proof.

Let us denote that $\mathit{a}=inf\mathrm{(}\mathit{f}\mathrm{)}\mathrm{,}\mathit{b}=sup\mathrm{(}\mathit{f}\mathrm{)}$ , and $$\begin{array}{c}{\mathit{E}}_{\mathrm{data}}\mathrm{(}\mathit{u}\mathrm{)}={\displaystyle \frac{{\mathit{\lambda }}_1}{2}}{{\displaystyle \int }}_{\mathrm{\Omega }}{\mathrm{(}\mathit{u}-\mathit{f}\mathrm{)}}^2\mathit{d}\mathit{x}+{\mathit{\lambda }}_2{{\displaystyle \int }}_{\mathrm{\Omega }}\mathrm{(}\mathit{u}-\mathit{f}log\mathit{u}\mathrm{)}\mathit{d}\mathit{x}.\end{array}$$

Fixing $\mathit{x}\in \mathrm{\Omega }$ and denoting the data fidelity term with h on ${\mathbb{R}}^{+}$ , where $$\begin{array}{c}\mathit{g}\mathrm{(}\mathit{t}\mathrm{)}={\displaystyle \frac{{\mathit{\lambda }}_1}{2}}{\mathrm{(}\mathit{t}-\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}\mathrm{)}}^2+{\mathit{\lambda }}_2\mathrm{(}\mathit{t}-\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}log\mathit{t}\mathrm{)}.\end{array}$$

Easily, we have that the first order derivative of g satisfies: $$\begin{array}{c}{\mathit{g}}^{\prime }\mathrm{(}\mathit{t}\mathrm{)}=\mathrm{(}\mathit{t}-\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}\mathrm{)}\mathrm{(}{\mathit{\lambda }}_1+{\displaystyle \frac{{\mathit{\lambda }}_2}{\mathit{t}}}\mathrm{)}.\end{array}$$

The function g decreases if $\mathit{t}\in \mathrm{(}0\mathrm{,}\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}\mathrm{)}$ and increases if $\mathit{t}\in \mathrm{(}\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}\mathrm{,}+\infty \mathrm{)}$ . Therefore, for every $\mathit{V}⩾\mathit{f}\mathrm{(}\mathit{x}\mathrm{)}$ , we have $$\begin{array}{c}\mathit{g}\mathrm{(}inf\mathrm{(}\mathit{t}\mathrm{,}\mathit{V}\mathrm{)}\mathrm{)}⩽\mathit{g}\mathrm{(}\mathit{t}\mathrm{)}.\end{array}$$

Hence, if $\mathit{V}=\mathit{b}$ , we have $$\begin{array}{c}{\mathit{E}}_{\mathrm{data}}\mathrm{(}inf\mathrm{(}\mathit{u}\mathrm{,}\mathit{V}\mathrm{)}\mathrm{)}⩽{\mathit{E}}_{\mathrm{data}}\mathrm{(}\mathit{u}\mathrm{)}.\end{array}$$

Furthermore, from Kornprobst et al. (1999), we have: ${\int }_{\mathrm{\Omega }}|\nabla inf\mathrm{(}\mathit{u}\mathrm{,}\mathit{b}\mathrm{)}|⩽{\int }_{\mathrm{\Omega }}|\nabla \mathit{u}|$ . Hence, $\mathit{E}\mathrm{(}inf\mathrm{(}\mathit{u}\mathrm{,}\mathit{b}\mathrm{)}\mathrm{)}⩽\mathit{E}\mathrm{(}\mathit{u}\mathrm{)}$ . In the same way, we have: $\mathit{E}\mathrm{(}sup\mathrm{(}\mathit{u}\mathrm{,}\mathit{a}\mathrm{)}\mathrm{)}⩽\mathit{E}\mathrm{(}\mathit{u}\mathrm{)}$ , where $\mathit{a}=inf\mathrm{(}\mathit{f}\mathrm{)}$ . Thence, we can assume $\mathit{a}⩽{\mathit{u}}_{\mathit{n}}⩽\mathit{b}$ , the sequence $\{{\mathit{u}}_{\mathit{n}}\}$ is bounded in ${\mathit{L}}^1\mathrm{(}\mathrm{\Omega }\mathrm{)}$ .

Since $\{{\mathit{u}}_{\mathit{n}}\}$ is a minimizing sequence, we know that $\mathit{E}\mathrm{(}{\mathit{u}}_{\mathit{n}}\mathrm{)}$ is bounded. Hence, also the regularization term ${\int }_{\mathrm{\Omega }}|\nabla \mathit{u}|$ is bounded and $\{{\mathit{u}}_{\mathit{n}}\}$ is bounded in $\mathit{B}\mathit{V}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ .

There exists ${\mathit{u}}^{\ast }\in \mathit{B}\mathit{V}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ such that up to a subsequence, we have that ${\mathit{u}}_{\mathit{n}}$ converges weakly to ${\mathit{u}}^{\ast }\in \mathit{B}\mathit{V}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ and ${\mathit{u}}_{\mathit{n}}$ converges strongly to ${\mathit{u}}^{\ast }\in {\mathit{L}}^1\mathrm{(}\mathrm{\Omega }\mathrm{)}$ . We have $\mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ is closed and convex. Using $0\mathrm{<}\mathit{a}⩽{\mathit{u}}^{\ast }⩽\mathit{b}$ , the lower semicontinuity of the total variation and Fatou’s lemma, we get that ${\mathit{u}}^{\ast }$ is a minimizer of the problem (3).  □

4 Numerical Method

4.1 Discretization Scheme

Our scheme allows to perform both deblurring and denoising simultaneously. To do so, we need to compute: $$\begin{array}{rl}& {\mathit{u}}^{\ast }=\underset{\mathit{u}\in \mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}}{arg\; min}\mathit{E}\mathrm{(}\mathit{u}\mathrm{)}\mathrm{,}\\ & \mathit{E}\mathrm{(}\mathit{u}\mathrm{)}={{\displaystyle \int }}_{\mathrm{\Omega }}\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}|\nabla \mathit{u}|\mathit{d}\mathit{x}+{\displaystyle \frac{{\mathit{\lambda }}_1}{2}}{{\displaystyle \int }}_{\mathrm{\Omega }}{\mathrm{(}\mathit{K}\mathit{u}-\mathit{f}\mathrm{)}}^2\mathit{d}\mathit{x}+{\mathit{\lambda }}_2{{\displaystyle \int }}_{\mathrm{\Omega }}\mathrm{(}\mathit{K}\mathit{u}-\mathit{f}log\mathit{K}\mathit{u}\mathrm{)}\mathit{d}\mathit{x}\mathrm{,}\end{array}$$(9) where K is a blurring operator (convex), f is the observed image, $\mathit{S}\mathrm{(}\mathrm{\Omega }\mathrm{)}$ is the set of positive functions defined over Ω with bounded total variation, and ${\mathit{\lambda }}_1\mathrm{,}{\mathit{\lambda }}_2$ are positive regularization parameters. This functional $\mathit{u}\mapsto \mathit{E}\mathrm{(}\mathit{u}\mathrm{)}$ is still strictly convex, because K is assumed to be convex.

The images we are handling are discrete, i.e. matrices of pixel values rather than functions from ${\mathbb{R}}^2\to \mathbb{R}$ . Therefore we have to choose a discretization scheme for numerical computations. If u is a image, we write ${\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}$ for the pixel at coordinates $\mathrm{(}\mathit{j}\mathrm{,}\mathit{k}\mathrm{)}$ in u. We define the following quantities: $$\begin{array}{rl}& {\nabla }_1{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}={\mathit{u}}_{\mathit{j}+1\mathrm{,}\mathit{k}}-{\mathit{u}}_{\mathit{j}-1\mathrm{,}\mathit{k}}\mathrm{,}{\nabla }_2{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}={\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}+1}-{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}-1}\mathrm{,}\\ & \nabla {\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}=\mathrm{(}{\nabla }_1{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{,}{\nabla }_2{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{)}\mathrm{,}|\nabla {\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}|=\sqrt{{\mathrm{(}{\nabla }_1{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{)}}^2+{\mathrm{(}{\nabla }_2{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{)}}^2+{\mathit{\epsilon }}^2}\mathrm{,}\end{array}$$ where ε is a small positive number, added to avoid divisions by 0 in the implementation of the algorithms. Finding a minimum for the problem (2) can be achieved via the steepest gradient descent method $$\begin{array}{c}{\displaystyle \frac{\mathit{\delta }\mathit{E}\mathrm{(}\mathit{u}\mathrm{)}}{\mathit{\delta }{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}}}=\text{div}\mathrm{\left(}{\displaystyle \frac{\nabla {\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}}{|\nabla {\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}|}}\mathrm{\right)}-{\mathit{\lambda }}_1{\mathit{K}}^{\mathit{T}}\mathrm{(}\mathit{K}{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}-{\mathit{f}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{)}-{\mathit{\lambda }}_2\mathrm{(}\mathit{K}-{\displaystyle \frac{{\mathit{f}}_{\mathit{j}\mathrm{,}\mathit{k}}}{{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}}}\mathrm{)}.\end{array}$$ The operator divergence $\text{div}\mathrm{\left(}\frac{\nabla \mathit{u}}{|\nabla \mathit{u}|}\mathrm{\right)}$ is defined by $$\begin{array}{c}{\displaystyle \frac{\mathrm{(}{\nabla }_{11}\mathit{u}\mathrm{)}{\mathrm{(}{\nabla }_2\mathit{u}\mathrm{)}}^2-2\mathrm{(}{\nabla }_1\mathit{u}\mathrm{)}\mathrm{(}{\nabla }_2\mathit{u}\mathrm{)}\mathrm{(}{\nabla }_{12}\mathit{u}\mathrm{)}+\mathrm{(}{\nabla }_{22}\mathit{u}\mathrm{)}{\mathrm{(}{\nabla }_1\mathit{u}\mathrm{)}}^2}{{\mathrm{(}{\mathrm{(}{\nabla }_1\mathit{u}\mathrm{)}}^2+{\mathrm{(}{\nabla }_1\mathit{u}\mathrm{)}}^2+{\mathit{\epsilon }}^2\mathrm{)}}^{3\mathrm{/}2}}}\mathrm{,}\end{array}$$ where $$\begin{array}{rl}& {\nabla }_{11}{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}={\nabla }_1\mathrm{(}{\nabla }_1{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{)}={\mathit{u}}_{\mathit{j}+1\mathrm{,}\mathit{k}}-2{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}+{\mathit{u}}_{\mathit{j}-1\mathrm{,}\mathit{k}}\mathrm{,}\\ & {\nabla }_{22}{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}={\nabla }_2\mathrm{(}{\nabla }_2{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{)}={\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}+1}-2{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}+{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}-1}\mathrm{,}\\ & {\nabla }_{12}{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}={\nabla }_1\mathrm{(}{\nabla }_2{\mathit{u}}_{\mathit{j}\mathrm{,}\mathit{k}}\mathrm{)}={\mathit{u}}_{\mathit{j}+1\mathrm{,}\mathit{k}+1}+{\mathit{u}}_{\mathit{j}-1\mathrm{,}\mathit{k}-1}-{\mathit{u}}_{\mathit{j}+1\mathrm{,}\mathit{k}-1}-{\mathit{u}}_{\mathit{j}-1\mathrm{,}\mathit{k}+1}.\end{array}$$ Thus, for a small parameter $\mathit{\delta }\mathit{t}\mathrm{>}0$ , a solution of the minimization problem (2) may be computed by $$\begin{array}{c}{\displaystyle \frac{{\mathit{u}}^{\mathrm{(}\mathit{t}+1\mathrm{)}}-{\mathit{u}}^{\mathrm{(}\mathit{t}\mathrm{)}}}{\mathit{\delta }\mathit{t}}}=\text{div}\mathrm{(}\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}\mathrm{(}{\displaystyle \frac{\nabla {\mathit{u}}^{\mathrm{(}\mathit{t}\mathrm{)}}}{|\nabla {\mathit{u}}^{\mathrm{(}\mathit{t}\mathrm{)}}|}}\mathrm{\left)}\mathrm{\right)}-{\mathit{\lambda }}_1{\mathit{K}}^{\mathit{T}}\mathrm{(}\mathit{K}\mathit{u}-\mathit{f}\mathrm{)}-{\mathit{\lambda }}_2\mathrm{(}\mathit{K}-{\displaystyle \frac{\mathit{f}}{\mathit{u}}}\mathrm{)}.\end{array}$$ When the time-step parameter $\mathit{\delta }\mathit{t}$ becomes small, the convergence speed becomes so slow that larger images are proceeded with poor efficiency. There are many methods (Chambolle, 2004; Micchelli et al., 2011; Boyd et al., 2010) which can be used for the minimization problem in (2). In this paper, we extend the split-Bregman algorithm (Goldstein and Osher, 2009) to solve the minimization problem.

4.2 Proposed Algorithm

First, let us first review the split-Bregman algorithm (Goldstein and Osher, 2009). Suppose that we have a scalar γ and two convex functionals $\mathrm{\Psi }\mathrm{(}·\mathrm{)}$ and $\mathit{G}\mathrm{(}·\mathrm{)}$ ; and that we need to solve the following constrained optimization problem: $$\begin{array}{rl}& \text{find}\underset{\mathit{u}\mathrm{,}\mathit{d}}{arg\; min}\Vert \mathit{d}{\Vert }_1+{\displaystyle \frac{\mathit{\gamma }}{2}}\mathit{G}\mathrm{(}\mathit{u}\mathrm{)}\mathrm{,}\\ & \text{s.t.}\mathit{d}=\mathrm{\Psi }\mathrm{(}\mathit{u}\mathrm{)}.\end{array}$$(10)

We convert (10) into an unconstrained problem: $$\begin{array}{c}\text{find}\underset{\mathit{u}\mathrm{,}\mathit{d}}{arg\; min}\Vert \mathit{d}{\Vert }_1+{\displaystyle \frac{\mathit{\gamma }}{2}}\mathit{G}\mathrm{(}\mathit{u}\mathrm{)}+{\displaystyle \frac{\mathit{\rho }}{2}}\Vert \mathit{d}-\mathrm{\Psi }\mathrm{(}\mathit{u}\mathrm{)}-\mathit{b}{\Vert }_2^2\mathrm{,}\end{array}$$(11) where ρ is a penalty parameter (a positive constant) and b is a variable related to the split-Bregman iteration algorithm (to be explicited later). The solution to problem (11) can be approximated by the split-Bregman Algorithm (Goldstein and Osher, 2009): $$\begin{array}{rl}& {\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}=\underset{\mathit{u}}{arg\; min}{\displaystyle \frac{\mathit{\gamma }}{2}}\mathit{G}\mathrm{(}\mathit{u}\mathrm{)}+{\displaystyle \frac{\mathit{\rho }}{2}}{\Vert {\mathit{d}}^{\mathrm{(}\mathit{k}\mathrm{)}}-\mathrm{\Psi }\mathrm{(}\mathit{u}\mathrm{)}-{\mathit{b}}^{\mathrm{(}\mathit{k}\mathrm{)}}\Vert }_2^2\mathrm{,}\\ & {\mathit{d}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}=\underset{\mathit{d}}{arg\; min}\Vert \mathit{d}{\Vert }_1+{\displaystyle \frac{\mathit{\rho }}{2}}{\Vert \mathit{d}-\mathrm{\Psi }\mathrm{\left(}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}\mathrm{\right)}-{\mathit{b}}^{\mathrm{(}\mathit{k}\mathrm{)}}\Vert }_2^2\mathrm{,}\\ & {\mathit{b}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}={\mathit{b}}^{\mathrm{(}\mathit{k}\mathrm{)}}+\mathrm{\Psi }\mathrm{\left(}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}\mathrm{\right)}-{\mathit{d}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}.\end{array}$$ Now we return to the problem (9). We define $$\begin{array}{c}\mathit{G}\mathrm{(}\mathit{u}\mathrm{)}={\displaystyle \frac{{\mathit{\lambda }}_1}{2}}{\mathrm{(}\mathit{K}\mathit{u}-\mathit{f}\mathrm{)}}^2+{\mathit{\lambda }}_2\mathrm{(}\mathit{K}\mathit{u}-\mathit{f}log\mathit{K}\mathit{u}\mathrm{)}\text{and}\mathrm{\Psi }\mathrm{(}\mathit{u}\mathrm{)}=\mathit{\alpha }\nabla \mathit{u}.\end{array}$$ We set $\mathit{\nu }=\mathit{K}\mathit{u}$ ; then, based on equation (11), the split-Bregman problem for (9) is defined as: $$\begin{array}{c}\underset{\mathit{u}\mathrm{,}\mathit{d}}{arg\; min}\mathrm{(}\Vert \mathit{d}{\Vert }_1+{\displaystyle \frac{\mathit{\gamma }}{2}}\mathit{G}\mathrm{(}\mathit{\nu }\mathrm{)}+{\displaystyle \frac{{\mathit{\rho }}_1}{2}}\Vert \mathit{\nu }-\mathit{K}\mathit{u}-\mathit{c}{\Vert }_2^2+{\displaystyle \frac{{\mathit{\rho }}_2}{2}}{\displaystyle \sum }_{\mathit{i}=1\mathrm{,}2}\Vert {\mathit{d}}_{\mathit{i}}-\mathit{\alpha }{\nabla }_{\mathit{i}}\mathit{u}-{\mathit{b}}_{\mathit{i}}{\Vert }_2^2\mathrm{)}\mathrm{,}\end{array}$$(12) where the parameters ${\mathit{\rho }}_1$ , ${\mathit{\rho }}_2$ and γ are positive, $\mathit{d}=\mathrm{(}{\mathit{d}}_1\mathrm{,}{\mathit{d}}_2\mathrm{)}$ , $\mathit{b}=\mathrm{(}{\mathit{b}}_1\mathrm{,}{\mathit{b}}_2\mathrm{)}$ and $\nabla \mathit{u}=\mathrm{(}{\nabla }_1\mathit{u}\mathrm{,}{\nabla }_2\mathit{u}\mathrm{)}$ .

The split-Bregman method for solving (12) is described as follows: $$\begin{array}{rl}& {\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}=\underset{\mathit{u}}{arg\; min}{\displaystyle \frac{{\mathit{\rho }}_1}{2}}{\Vert {\mathit{\nu }}^{\mathrm{(}\mathit{k}\mathrm{)}}-\mathit{K}\mathit{u}-{\mathit{c}}^{\mathrm{(}\mathit{k}\mathrm{)}}\Vert }_2^2+{\displaystyle \frac{{\mathit{\rho }}_2}{2}}{\displaystyle \sum }_{\mathit{i}=1\mathrm{,}2}{\Vert {\mathit{d}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}-\mathit{\alpha }{\nabla }_{\mathit{i}}\mathit{u}-{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}\Vert }_2^2\mathrm{,}\\ & {\mathit{\nu }}^{\mathrm{(}\mathit{k}+1\mathrm{)}}=\underset{\mathit{\nu }}{arg\; min}{\displaystyle \frac{\mathit{\gamma }}{2}}\mathit{G}\mathrm{(}\mathit{\nu }\mathrm{)}+{\displaystyle \frac{{\mathit{\rho }}_1}{2}}{\Vert \mathit{\nu }-\mathit{K}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}-{\mathit{c}}^{\mathrm{(}\mathit{k}\mathrm{)}}\Vert }_2^2\mathrm{,}\\ & {\mathit{d}}_{\mathit{i}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}=\underset{\mathit{d}}{arg\; min}\Vert {\mathit{d}}_{\mathit{i}}{\Vert }_1+{\displaystyle \frac{{\mathit{\rho }}_2}{2}}{\Vert {\mathit{d}}_{\mathit{i}}-\mathit{\alpha }{\nabla }_{\mathit{i}}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}-{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}\Vert }_2^2\mathrm{,}\\ & {\mathit{c}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}={\mathit{c}}^{\mathrm{(}\mathit{k}\mathrm{)}}+\mathit{K}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}-{\mathit{\nu }}^{\mathrm{(}\mathit{k}+1\mathrm{)}}\mathrm{,}\\ & {\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}={\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}+\mathit{\alpha }{\nabla }_{\mathit{i}}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}-{\mathit{d}}_{\mathit{i}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}.\end{array}$$

There are three subproblems to solve: u, ν and d.

Subproblem 1. The u subproblem is a quadratic optimization problem, whose optimality condition reads: $$\begin{array}{rl}& \mathrm{(}{\mathit{\rho }}_1{\mathit{K}}^{\mathit{T}}·\mathit{K}+{\mathit{\rho }}_2\mathit{\alpha }{\displaystyle \sum }_{\mathit{i}=1\mathrm{,}2}{\nabla }_{\mathit{i}}^{\mathit{T}}{\nabla }_{\mathit{i}}\mathrm{)}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}\\ & ={\mathit{\rho }}_1{\mathit{K}}^{\mathit{T}}\mathrm{(}{\mathit{\nu }}^{\mathrm{(}\mathit{k}\mathrm{)}}-{\mathit{c}}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{)}+{\mathit{\rho }}_2{\displaystyle \sum }_{\mathit{i}=1\mathrm{,}2}{\nabla }_{\mathit{i}}^{\mathit{T}}\mathrm{(}{\mathit{d}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}-{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{)}\mathrm{,}\end{array}$$(13) under considering periodic boundary conditions. Note that left-hand-side matrix in (13) includes a Laplacian matrix ( ${\nabla }_1^{\mathit{T}}{\nabla }_1+{\nabla }_2^{\mathit{T}}{\nabla }_2=-\mathrm{\Delta }$ ) and is strictly diagonally dominant. Following Wang et al. (2008), equation (13) can be solved efficiently with one fast Fourier transform (FFT) operation and one inverse FFT operation as: $$\begin{array}{c}\mathit{u}={\mathcal{F}}^{-1}\mathrm{\left(}{\displaystyle \frac{\mathcal{F}\mathrm{(}\mathit{r}\mathrm{)}}{{\mathit{\rho }}_1\mathcal{F}\mathrm{(}{\mathit{K}}^{\mathit{T}}\mathrm{)}·\mathcal{F}\mathrm{(}\mathit{K}\mathrm{)}-{\mathit{\rho }}_2\mathcal{F}\mathrm{(}\mathit{\alpha }\mathrm{)}·\mathcal{F}\mathrm{(}\mathrm{\Delta }\mathrm{)}}}\mathrm{\right)}\mathrm{,}\end{array}$$(14) where $$\begin{array}{c}\mathit{r}={\mathit{\rho }}_1{\mathit{K}}^{\mathit{T}}\mathrm{(}{\mathit{\nu }}^{\mathrm{(}\mathit{k}\mathrm{)}}-{\mathit{b}}_1^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{)}+{\mathit{\rho }}_2{\displaystyle \sum }_{\mathit{i}=1\mathrm{,}2}{\nabla }_{\mathit{i}}^{\mathit{T}}\mathrm{(}{\mathit{d}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}-{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{)}\mathrm{,}\end{array}$$ $\mathcal{F}$ and ${\mathcal{F}}^{-1}$ are the forward and inverse Fourier transform operators.

Subproblem 2. The optimality condition for the ν subproblem is given by $$\begin{array}{c}{\displaystyle \frac{\mathit{\gamma }}{2}}\mathrm{(}{\mathit{\lambda }}_1\mathrm{(}\mathit{\nu }-\mathit{f}\mathrm{)}+{\mathit{\lambda }}_2\mathrm{(}1-{\displaystyle \frac{\mathit{f}}{\mathit{\nu }}}\mathrm{\left)}\mathrm{\right)}+{\mathit{\rho }}_1\mathrm{(}\mathit{\nu }-\mathit{K}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}-{\mathit{c}}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{)}=0.\end{array}$$ This equation canbe rewritten as: $$\begin{array}{rl}& \mathrm{(}{\displaystyle \frac{\mathit{\gamma }}{2}}{\mathit{\lambda }}_1+{\mathit{\rho }}_1\mathrm{)}{\mathrm{\left(}{\mathit{\nu }}^{\mathrm{(}\mathit{k}+1\mathrm{)}}\mathrm{\right)}}^2-\mathrm{(}{\displaystyle \frac{\mathit{\gamma }}{2}}{\mathit{\lambda }}_1\mathit{f}-{\mathit{\lambda }}_2{\displaystyle \frac{\mathit{\gamma }}{2}}+{\mathit{\rho }}_1\mathrm{(}\mathit{K}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}+{\mathit{c}}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{\left)}\mathrm{\right)}{\mathit{\nu }}^{\mathrm{(}\mathit{k}+1\mathrm{)}}-{\displaystyle \frac{\mathit{\gamma }}{2}}{\mathit{\lambda }}_2\mathit{f}\\ & =0.\end{array}$$ The positive solution is given by $$\begin{array}{c}{\mathit{\nu }}^{\mathrm{(}\mathit{k}+1\mathrm{)}}={\mathit{S}}^{\mathrm{(}\mathit{k}\mathrm{)}}+\sqrt{{\mathrm{\left(}{\mathit{S}}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{\right)}}^2+{\displaystyle \frac{\mathit{\gamma }{\mathit{\lambda }}_2\mathit{f}}{\mathit{\gamma }{\mathit{\lambda }}_1+2{\mathit{\rho }}_1}}}\mathrm{,}\end{array}$$(15) where $$\begin{array}{c}{\mathit{S}}^{\mathit{k}}={\displaystyle \frac{{\mathit{\lambda }}_1\mathit{\gamma }\mathit{f}-{\mathit{\lambda }}_2\mathit{\gamma }+2{\mathit{\rho }}_1\mathrm{(}\mathit{K}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}+{\mathit{b}}_{\mathit{\nu }}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{)}}{2\mathrm{(}\mathit{\gamma }{\mathit{\lambda }}_1+2{\mathit{\rho }}_1\mathrm{)}}}.\end{array}$$

Subproblem 3. The solution of the d subproblem can readily be obtained by applying the soft thresholding operator (see Micchelli et al., 2011). We can use shrinkage operators to compute the optimal values of ${\mathit{d}}_1$ and ${\mathit{d}}_2$ separately: $$\begin{array}{c}{\mathit{d}}_{\mathit{i}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}=\text{shrink}\mathrm{(}\mathit{\alpha }{\nabla }_{\mathit{i}}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}+{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}\mathrm{,}{\displaystyle \frac{1}{{\mathit{\rho }}_2}}\mathrm{)}.\end{array}$$(16) The problem (16) is solved as: $$\begin{array}{c}{\mathit{d}}_{\mathit{i}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}={\displaystyle \frac{\mathit{\alpha }{\nabla }_{\mathit{i}}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}+{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}}{|\mathit{\alpha }{\nabla }_{\mathit{i}}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}+{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}|}}·max\mathrm{\left(}\right|\mathit{\alpha }{\nabla }_{\mathit{i}}{\mathit{u}}^{\mathrm{(}\mathit{k}+1\mathrm{)}}+{\mathit{b}}_{\mathit{i}}^{\mathrm{(}\mathit{k}\mathrm{)}}|-{\displaystyle \frac{1}{{\mathit{\rho }}_2}}\mathrm{,}0\mathrm{)}.\end{array}$$(17)

The algorithm. The complete method is summarized in Algorithm 1. We need a stopping criterion for the iteration; we end the loop if the maximum number of allowed outer iterations N has been carried out (to guarantee an upper bound on running time) or the following condition is satisfied for some prescribed tolerance ς: $$\begin{array}{c}{\displaystyle \frac{\Vert {\mathit{u}}^{\mathrm{(}\mathit{k}\mathrm{)}}-{\mathit{u}}^{\mathrm{(}\mathit{k}-1\mathrm{)}}{\Vert }_2}{\Vert {\mathit{u}}^{\mathrm{(}\mathit{k}\mathrm{)}}{\Vert }_2}}\mathrm{<}\mathit{\varsigma }\mathrm{,}\end{array}$$ where ς is a small positive parameter. For our experiments, we set tolerance $\mathit{\varsigma }=0.0005$ and $\mathit{N}=500$ .

Algorithm 1

Adaptive split-Bregman algorithm for solving the model (9).

5 Numerical Simulations

5.1 Implementation Issues

In this section, we show some numerical reconstructions obtained applying our proposed method for mixed Poisson–Gaussian noise. We compare our reconstructions with other images obtained other well known methods, such as TV- ${\mathit{L}}_1$ (Chambolle et al., 2010), the Modified scheme for Mixed Poisson–Gaussian model (MS-MPG) (Pham et al., 2018) and the Bregman method (Goldstein and Osher, 2009). All of the compared methods perform image denoising with their optimal parameters. For a fair comparison, the regularization parameters of both MS-MPG and our proposed are the same: ${\mathit{\lambda }}_1=0.4$ , ${\mathit{\lambda }}_2=0.6$ . We set ${\mathit{\rho }}_1=1$ , ${\mathit{\rho }}_2=1$ . The parameter σ in $\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}$ is set to 1. The threshold value l in the function $\mathit{\alpha }\mathrm{(}\mathit{x}\mathrm{)}$ and the parameters γ are chosen to keep the poise between noise removal and detail preservation capabilities.

The test images1 are 8-bit gray scale standard images of size $256\times 256$ shown in Fig. 1.

Fig. 1

Original images.

All the experiments were run on a machine with Core i7-CPU 2 GHz, SDRAM 4 GB-DDR III 2 Ghz, Windows 10 (64 bit), and implemented in MATLAB. To compare the efficiency of algorithms, we use the Peak Signal-to-Noise Ratio (PSNR) and the Structure Similarity Index (SSIM) (Wang and Bovik, 2006).

The first metric, PSNR (db), is defined by $$\begin{array}{c}\mathit{PSNR}=10{log}_{10}\mathrm{\left(}{\displaystyle \frac{\mathit{\text{MN}}{\mathit{I}}_{max}^2}{\Vert {\mathit{u}}^{\ast }-\mathit{u}{\Vert }_2^2}}\mathrm{\right)}\mathrm{,}\end{array}$$ where $\mathit{u}\mathrm{,}{\mathit{u}}^{\ast }$ are, respectively, the original image and the reconstructed (or noisy) image, ${\mathit{I}}_{max}$ is the maximum intensity of the original image, M and N are the number of image pixels in rows and columns.

The second metric, SSIM, is defined by $$\begin{array}{c}\mathit{SSIM}\mathrm{(}\mathit{u}\mathrm{,}{\mathit{u}}^{\ast }\mathrm{)}={\displaystyle \frac{\mathrm{(}2{\mathit{\mu }}_{\mathit{u}}{\mathit{\mu }}_{{\mathit{u}}^{\ast }}+{\mathit{c}}_1\mathrm{)}\mathrm{(}2{\mathit{\sigma }}_{\mathit{u}\mathrm{,}{\mathit{u}}^{\ast }}+{\mathit{c}}_2\mathrm{)}}{\mathrm{(}{\mathit{\mu }}_{\mathit{u}}^2+{\mathit{\mu }}_{{\mathit{u}}^{\ast }}^2+{\mathit{c}}_1\mathrm{)}\mathrm{(}{\mathit{\sigma }}_{\mathit{u}}^2+{\mathit{\sigma }}_{{\mathit{u}}^{\ast }}^2+{\mathit{c}}_2\mathrm{)}}}\mathrm{,}\end{array}$$ where ${\mathit{\mu }}_{\mathit{u}}$ , ${\mathit{\mu }}_{{\mathit{u}}^{\ast }}$ are the means of u, ${\mathit{u}}^{\ast }$ , respectively; ${\mathit{\sigma }}_{\mathit{u}}\mathrm{,}{\mathit{\sigma }}_{{\mathit{u}}^{\ast }}$ , their standard deviations; ${\mathit{\sigma }}_{\mathit{u}\mathrm{,}{\mathit{u}}^{\ast }}$ , the covariance of two images u and ${\mathit{u}}^{\ast }$ ; ${\mathit{c}}_1={\mathrm{(}{\mathit{K}}_1\mathit{L}\mathrm{)}}^2$ ; ${\mathit{c}}_2={\mathrm{(}{\mathit{K}}_2\mathit{L}\mathrm{)}}^2$ ; L is the dynamic range of the pixel values (255 for 8-bit grayscale images); and finally ${\mathit{K}}_1\ll 1$ , ${\mathit{K}}_2\ll 1$ are small constants.

5.2 Numerical Results and Discussion

5.2.1 Image Denoising

Our method can perform image deblurring and denoising simultaneously. In this section, we perform only image denoising. Noisy observations are generated by Poisson noise with some peak ${\mathit{I}}_{max}$ and by Gaussian noise with standard deviation ${\mathit{\sigma }}_{\mathrm{G}}$ to the test images. In Figs. 2, 4 and 5, we give the results for denoising the corrupted images for different noise levels ${\mathit{I}}_{max}$ and ${\mathit{\sigma }}_{\mathrm{G}}=10$ .

Fig. 2

Recovered results for the test images. (a) Noisy image with ${\mathit{I}}_{max}=120$ , ${\mathit{\sigma }}_{\mathrm{G}}=10$ , (b) TV ${\mathit{L}}_1$ , (c) Bregman, (d) MS-MPG, (e) Our proposed.

Fig. 3

The zoomed-in part of the recovered images in Fig. 2. (a) First column: details of original images; (b) Second column: details of observed images; (c) Third column: details of restored images by TV- ${\mathit{L}}_1$ method; (d) Fourth column: details of restored images by Bregman method; (e) Fifth column: details of restored images by MS-MPG method; (f) Sixth column: details of restored images by our proposed method.

Fig. 4

Recovered results for the test images. (a) Noisy image with ${\mathit{I}}_{max}=60$ , ${\mathit{\sigma }}_{\mathrm{G}}=10$ , (b) TV ${\mathit{L}}_1$ , (c) Bregman, (d) MS-MPG, (e) Ours.

Fig. 5

Recovered results for the test images. (a) Noisy image with ${\mathit{I}}_{max}=60$ , ${\mathit{\sigma }}_{\mathrm{G}}=10$ , (b) TV- ${\mathit{L}}_1$ , (c) Bregman, (d) MS-MPG, (e) Ours.

Fig. 6

The zoomed-in part of the recovered images in Fig. 4. (a) Details of original images; (b) details of observed images; (c) details of restored images by TV ${\mathit{L}}_1$ method; (d) details of restored images by Bregman method; (e) details of restored images by MS-MPG method; (f) details of restored images by our proposed method.

Fig. 7

The zoomed-in part of the recovered images in Fig. 5. (a) Details of original images; (b) details of observed images; (c) details of restored images by TV ${\mathit{L}}_1$ method; (d) details of restored images by Bregman method; (e) details of restored images by MS-MPG method; (f) details of restored images by our proposed method.

For a better visual comparison, we have enlarged some details of the restored images in Figs. 3, 6 and 7 (we include the original images in the first column). It can be seen that our method gives even better visual quality than other methods. Table 1 shows the computation time in second(s) of the compared methods for Fig. 2. We see from Table 1 that the computation time of the restored images by the proposed method and the Bregman method is about the same. However, the computational time required by the proposed method is less than that required by the MS-MPG and TV ${\mathit{L}}_1$ . The comparison metrics PSNR, SSIM are also computed using various noise levels and shown in Table 2 and Table 3. The best values among all the methods are shown in bold. We give the values of the PSNR and SSIM for the noisy and recovered images. The results shown in Tables 1, 2 and 3 prove that the proposed method is convergent and gets higher PSNR and SSIM values than others.

Table 1

Execution time for different denoising methods (in seconds) with noise level ${\mathit{I}}_{max}$ and ${\mathit{\sigma }}_{\mathrm{G}}=10$ .

Image	Method	CPU time (s)
		${\mathit{I}}_{max}=120$	${\mathit{I}}_{max}=60$
			TV ${\mathit{L}}_1$	4.3449	5.6730
Clock	Bregman	0.9460	0.8212
	MS-MPG	4.1465	4.8734
	Ours	1.0945	1.1081
	TV ${\mathit{L}}_1$	5.6229	7.4171
Coco	Bregman	1.0265	0.8414
	MS-MPG	4.0844	5.0879
	Ours	1.1239	1.2251
	TV ${\mathit{L}}_1$	4.3096	6.4129
Lamp	Bregman	0.9225	0.9473
	MS-MPG	4.1810	4.8758
	Ours	0.9431	1.1266

Table 2

PSNR values and SSIM measures for noisy images and recovered images with ${\mathit{I}}_{max}=120$ .

Image	PSNR					MSSIM
	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours
	${\mathit{I}}_{max}=120$ , $\mathit{\sigma }=10$
Jetplane	18.9416	22.7203	24.1190	24.7848	25.3251	0.4045	0.7061	0.7514	0.7511	0.7748
Lake	19.6413	21.3675	22.5906	22.9972	24.4798	0.5235	0.6360	0.6812	0.7069	0.7603
Aerial	17.4471	18.9550	19.5840	19.3051	19.8806	0.5582	0.5083	0.5808	0.5711	0.7130
Clock	18.3852	24.6040	25.7945	24.8844	26.1201	0.2997	0.8339	0.8822	0.7796	0.8970
Car	19.1385	21.4694	22.1559	22.8793	24.0620	0.4848	0.6106	0.6542	0.6804	0.7256
Coco	16.9119	20.4242	20.4215	20.3426	20.6539	0.2755	0.8551	0.8798	0.8296	0.8950
Lamp	17.8770	24.2808	24.3594	24.1062	24.6339	0.2446	0.8522	0.8891	0.7889	0.8985
Poulina	18.8381	25.2567	25.7203	25.9781	26.0653	0.3250	0.7648	0.7934	0.7982	0.8074
Spine	21.0004	25.2561	24.6855	25.5349	26.1010	0.6180	0.7925	0.7763	0.7967	0.8206
Head	21.7787	24.3567	26.2348	26.9061	27.0979	0.6324	0.8033	0.8043	0.8273	0.8400
Average	18.9960	22.8691	23.5666	23.7719	24.4420	0.4366	0.7363	0.7693	0.7530	0.8132
${\mathit{I}}_{max}=120$ , $\mathit{\sigma }=15$
Jetplane	16.7150	22.2033	23.4915	23.6918	24.1415	0.3320	0.6761	0.7248	0.6959	0.7320
Lake	17.2574	20.8215	22.0827	22.2260	23.0442	0.4384	0.6021	0.6732	0.6709	0.7040
Aerial	15.8006	18.7671	19.2795	19.1060	19.4706	0.4622	0.4594	0.5740	0.5139	0.6472
Clock	16.4619	24.2165	25.3645	24.2371	25.7740	0.2440	0.8105	0.8601	0.8186	0.8805
Car	16.8589	20.9512	21.7735	22.1269	22.7608	0.4015	0.5809	0.6402	0.6338	0.6727
Coco	15.4193	20.3398	20.4109	20.1332	20.5488	0.2181	0.8265	0.8599	0.7741	0.8789
Lamp	16.0461	23.8972	23.9090	23.5169	24.3063	0.1964	0.8225	0.8695	0.7210	0.8799
Poulina	16.6627	24.9195	25.2709	25.2753	25.4142	0.2452	0.7346	0.7659	0.7491	0.7739
Spine	18.5582	23.7301	24.4015	24.3122	24.9272	0.5378	0.7418	0.7689	0.7521	0.7794
Head	19.3512	24.549	25.4199	25.8356	25.9893	0.5588	0.7567	0.7836	0.7854	0.7991
Average	16.9131	22.4395	23.1404	23.0461	23.6377	0.3634	0.7011	0.7520	0.7115	0.7748

Table 3

PSNR values and SSIM measures for noisy images and recovered images with with ${\mathit{I}}_{max}=60$ .

Image	PSNR					MSSIM
	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours
	${\mathit{I}}_{max}=60$ , $\mathit{\sigma }=10$
Jetplane	14.0929	21.4515	22.5116	22.3705	22.8057	0.2570	0.6396	0.6482	0.6730	0.6854
Lake	14.7190	20.1480	20.9885	20.7335	21.5586	0.3488	0.5567	0.5987	0.5945	0.6325
Aerial	13.9091	18.7122	18.9386	18.6929	19.2898	0.3465	0.4036	0.5296	0.3801	0.5635
Clock	13.9941	23.7554	24.7607	24.3166	25.0682	0.1866	0.7759	0.7865	0.7931	0.8439
Car	14.2393	20.3390	20.8993	20.8988	21.4920	0.3124	0.5417	0.5723	0.5709	0.5864
Coco	13.4373	19.9609	19.9082	20.0459	20.2665	0.1573	0.7969	0.7815	0.8218	0.8535
Lamp	13.6235	23.3118	23.2870	23.4892	23.6101	0.1466	0.7898	0.7823	0.8067	0.8568
Poulina	14.1692	24.2429	24.8768	24.8704	24.9272	0.1804	0.6901	0.7169	0.7252	0.7316
Spine	16.0910	22.749	23.3981	23.3266	23.5011	0.4597	0.6821	0.7286	0.7153	0.7308
Head	16.9718	23.667	24.2991	24.2780	24.4763	0.4970	0.7059	0.7494	0.7284	0.7550
Average	14.5247	21.8338	22.3868	22.3022	22.6996	0.2892	0.6582	0.6894	0.6809	0.7240
${\mathit{I}}_{max}=60$ , $\mathit{\sigma }=15$
Image	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours
Jetplane	11.4314	20.5604	21.0317	21.2727	21.3729	0.1911	0.5883	0.5208	0.6247	0.6319
Lake	12.0450	19.3676	19.7789	19.9102	20.0911	0.2441	0.5053	0.5209	0.5526	0.5545
Aerial	11.6216	18.4021	18.9001	18.5632	19.1482	0.2425	0.3435	0.4216	0.3518	0.4362
Clock	11.4506	22.9914	23.6737	23.4250	24.3187	0.1365	0.7297	0.6387	0.7298	0.8123
Car	11.5354	19.6031	19.8705	20.1081	20.2498	0.2163	0.4898	0.4776	0.5254	0.5311
Coco	11.1477	19.6694	19.1809	19.7580	19.8315	0.1149	0.7432	0.6010	0.7628	0.8227
Lamp	11.1182	22.6734	22.1005	22.8145	22.9551	0.1014	0.7341	0.6151	0.7430	0.8263
Poulina	11.5927	23.4904	23.8398	23.8808	24.0040	0.1257	0.6353	0.6106	0.6818	0.6960
Spine	13.4551	20.8085	21.9682	22.0129	22.1122	0.3844	0.6115	0.6493	0.6548	0.6658
Head	14.3442	22.4799	22.4954	22.5105	22.9698	0.4370	0.6445	0.6890	0.6853	0.6991
Average	11.9742	21.0046	21.2840	21.4256	21.7053	0.2194	0.6025	0.5745	0.6312	0.6676

5.2.2 Image Deblurring and Denoising

In this section, we perform image denoising and delurring simultaneously. In our simulation, we use the Gaussian blur with a window size $9\times 9$ , and standard deviation of 1. After the blurring operation, we corrupt the images by Possion noise ${\mathit{I}}_{max}=120$ and ${\mathit{\sigma }}_{\mathrm{G}}=15$ . As in the previous experiment, we compare our results with those obtained by employing the Bregman method, the MS-MPG and the TV ${\mathit{L}}_1$ (see recovered results in Figs. 8, 10, and their zoom-in part in Figs. 9, 11).

Fig. 8

Recovered results for the test images. (a) Blurring and Noisy image, (b) TV ${\mathit{L}}_1$ , (c) Bregman, (d) MS-MPG, (e) Our proposed.

Fig. 9

The zoomed-in part of the recovered images in Fig. 8. (a) Details of original images; (b) details of observed images; (c) details of restored images by TV ${\mathit{L}}_1$ method; (d) details of restored images by Bregman method; (e) details of restored images by MS-MPG method; (f) details of restored images by our proposed method.

Fig. 10

Recovered results for the test images. (a) Blurring and Noisy image, (b) TV ${\mathit{L}}_1$ , (c) Bregman, (d) MS-MPG, (e) Our proposed.

Fig. 11

The zoomed-in part of the recovered images in Fig. 10. (a) Details of original images; (b) details of observed images; (c) details of restored images by TV ${\mathit{L}}_1$ method; (d) details of restored images by Bregman method; (e) details of restored images by MS-MPG method; (f) details of restored images by our proposed method.

Table 4

PSNR values and SSIM measures for noisy and blurring images and recovered images with with ${\mathit{I}}_{max}=120$ , $\mathit{\sigma }=15$ .

Image	PSNR					MSSIM
	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours	Noisy	TV ${\mathit{L}}_1$	Bregman	MS-MPG	Ours
	${\mathit{I}}_{max}=120$ , $\mathit{\sigma }=15$
Jetplane	14.9522	18.7079	18.72012	18.0420	19.0029	0.2282	0.6384	0.6600	0.5883	0.6860
Lake	16.0535	19.6419	19.6100	19.6472	20.2675	0.2876	0.5506	0.5449	0.5654	0.6090
Aerial	15.3701	18.3325	18.8549	18.7495	18.9921	0.3107	0.4647	0.4902	0.4960	0.5030
Clock	16.1891	23.2348	23.5898	22.5893	23.7133	0.1761	0.7758	0.8196	0.6575	0.8313
Car	15.5905	19.6202	19.6600	19.3774	20.1486	0.2408	0.5375	0.5486	0.5286	0.5863
Coco	15.3829	20.1479	20.1762	19.0025	20.3572	0.1410	0.8082	0.8468	0.7524	0.8608
Lamp	15.9477	23.4005	23.6315	21.6657	23.7635	0.1296	0.8001	0.8597	0.6919	0.8679
Poulina	15.3475	19.5518	19.6801	20.2705	20.4392	0.1710	0.6871	0.7099	0.6923	0.7205
Spine	16.0476	19.1907	18.6797	18.8847	19.3544	0.3865	0.5694	0.5832	0.6051	0.6286
Head	14.9812	16.7888	16.6991	16.6044	18.3481	0.4590	0.6352	0.6562	0.6711	0.7061
Average	15.5862	19.8617	19.9301	19.4833	20.4387	0.2531	0.6467	0.6719	0.6249	0.6999

In Table 4, we give the values of the PSNR and SSIM for different images and different variational methods. The best values among all the methods are shown in bold. Comparing the values of the PSNR and SSIM, we can clearly see that our method outperforms the others even in presence of blur.

6 Conclusion

In this paper, we have studied a fast total variation minimization method for image restoration. We propose an adaptive model for mixed Poisson–Gaussion noise removal. It is proved that the adaptive model is strictly convex. Then, we have employed split Bregman method to solve the proposed minimization problem. Our experimental results have shown that the quality of restored images by the proposed method are competitive with those restored by the existing total variation restoration methods. The most important contribution is that the proposed algorithm is very efficient.