You are viewing a javascript disabled version of the site. Please enable Javascript for this site to function properly.
Go to headerGo to navigationGo to searchGo to contentsGo to footer
In content section. Select this link to jump to navigation

An Efficient Total Variation Minimization Method for Image Restoration

Article type: Research Article

Authors: Pham, Cong Thang1; ∗ | Tran, Thi Thu Thao2 | Gamard, Guilhem3

Affiliations: [1] The University of Danang – University of Science and Technology, 54 Nguyen Luong Bang, Danang, Vietnam | [2] The University of Danang – University of Economics, 71 Ngu Hanh Son, Danang, Vietnam | [3] LIP, ENS de Lyon, 46 allée d’Italie, 69364 Lyon, France

Correspondence: [∗] Corresponding author.

Keywords: total variation, image restoration, mixed Poisson–Gaussian noise, convex optimization, split-Bregman method

DOI: 10.15388/20-INFOR407

Journal: Informatica, vol. 31, no. 3, pp. 539-560, 2020

Received April 2019
|
Accepted February 2020
|
Published: 2020
Get PDF open access

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.

1Introduction

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 u∗. 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:

(1)
f=P(u)+W,
where f is observed image, u is the unknown image, P(u) means that the image u is corrupted by Poisson noise, and W∼N(0,σ2) 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 L2 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):

(2)
u∗=arg minu∈S(Ω)∫Ω|∇u|+λ12∫Ω(u−f)2+λ2∫Ω(u−flogu),
where f is the observed image, Ω⊂R2 is a bounded domain, and S(Ω) is the set of positive functions from Ω to R; finally, λ1,λ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:

(3)
u∗=arg minu∈S(Ω)E(u),E(u)=∫Ωα(x)|∇u(x)|dx+λ12∫Ω(u(x)−f(x))2dx+λ2∫Ω(u(x)−f(x)logu(x))dx,
where f is the observed image, λ1 and λ2 are positive regularization parameters, S(Ω)={u∈BV(Ω):u>0} is closed and convex, with BV(Ω) being the space of functions Ω→R with bounded variation; and finally α(x) is a continuous function in S(Ω).

The function α(x) 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 α(x) is chosen as follows:

α(x)=11+l·|v(x)|2,
where l is a threshold value and v(x)=|∇Gσ(x)∗f|, in which ∗ denotes the convolution with Gσ(x)=12πσ2exp(−x22σ2), 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 E(·) is convex, which enables us to use larger time-step parameters during gradient descent when solving (3). We introduce the influence function α(x), 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.

2Preliminaries

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 P(u|f). Bayes’s rule gives:

(4)
P(u|f)=P(f|u)P(u)P(f).
The probability density function of the observed image f corrupted by Gaussian noise PN (respectively, by Poisson noise PP) is:
PN(f|u)=1σ2πexp(−(u−f)22σ2),PP(f|u)=ufexp(−u)f!,
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:
Pmixed(f|u)=ufexp(−u)f!12πσexp(−(u−f)22σ2).
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 M×N, and let I={1,…,M}×{1,…,N} denote the domain of f. For i in I, we write fi the pixel of f at position i (and similarly ui the pixel of u at position i). Then,
Pmixed(f|u)=∏i∈I(ui)fie(−ui)fi!sexp(−(ui−fi)22σ2)
with s=(2πσ)−1. Maximizing Pmixed is equivalent to minimizing −logPmixed, so let us compute the quantity −log(Pmixed(f|u)):
(5)
∑i∈Iui−filog(ui)+log(fi!)+y(ui−fi)2,
for some constant y>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(fi!) altogether.

Now we assume that P(u) follows a Gibbs prior (Le et al., 2007):

(6)
P(u)=1zexp(−∫|∇u|),
where z is a normalization factor. We need to make a couple of comments here. First, u is not a function R2→R, but rather a discrete array of pixels; thus the integral in that expression is going to be translated to a sum, while ∇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 P(u) 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 P(u|f). By equation (4), this amounts to minimize the expression −log(P(f|u))−log(P(u)), so we can plug in equations (5) and (6) to get:

(7)
u∗=arg minu∑i∈I1z|∇ui|+y(ui−fi)2+(ui−filog(ui)),
and we can view this expression as a discrete approximation of the functional E(·) defined as:
(8)
E(u)=∫Ω|∇u|dx+λ12∫Ω(u−f)2dx+λ2∫Ω(u−flogu)dx,
with λ1=2yz and λ2=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, |∇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 λ1 and λ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).

3Existence 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 E(·) is a convex functional; then, we show that E(·) has a lower bound. These two facts together imply the existence and uniqueness of the minimizer of E(·).

Theorem 1.

The functional u↦E(u), where E is defined in (3), is strictly convex.

Proof.

Let us set: h(u)=λ12(u−f)2+λ2(u−flogu). The first and the second order derivative of h are:

h′(u)=λ1u2−u(λ1f−λ2)−λ2fu
and
h″(u)=λ1u2+λ2fu2.
Since f is a positive, and u∈S(Ω), we have: h″(u)>0, i.e. h(u) is strictly convex. Moreover, the TV regularization is convex, thence E(u) is also strictly convex.  □

Theorem 2.

Let f∈S(Ω)∩L∞(Ω), then the problem (3) has an exactly one solution u∈BV(Ω) and satisfying:

infΩf⩽u⩽supΩf.

Proof.

Let us denote that a=inf(f),b=sup(f), and

Edata(u)=λ12∫Ω(u−f)2dx+λ2∫Ω(u−flogu)dx.

Fixing x∈Ω and denoting the data fidelity term with h on R+, where

g(t)=λ12(t−f(x))2+λ2(t−f(x)logt).

Easily, we have that the first order derivative of g satisfies:

g′(t)=(t−f(x))(λ1+λ2t).

The function g decreases if t∈(0,f(x)) and increases if t∈(f(x),+∞). Therefore, for every V⩾f(x), we have

g(inf(t,V))⩽g(t).

Hence, if V=b, we have

Edata(inf(u,V))⩽Edata(u).

Furthermore, from Kornprobst et al. (1999), we have: ∫Ω|∇inf(u,b)|⩽∫Ω|∇u|. Hence, E(inf(u,b))⩽E(u). In the same way, we have: E(sup(u,a))⩽E(u), where a=inf(f). Thence, we can assume a⩽un⩽b, the sequence {un} is bounded in L1(Ω).

Since {un} is a minimizing sequence, we know that E(un) is bounded. Hence, also the regularization term ∫Ω|∇u| is bounded and {un} is bounded in BV(Ω).

There exists u∗∈BV(Ω) such that up to a subsequence, we have that un converges weakly to u∗∈BV(Ω) and un converges strongly to u∗∈L1(Ω). We have S(Ω) is closed and convex. Using 0<a⩽u∗⩽b, the lower semicontinuity of the total variation and Fatou’s lemma, we get that u∗ is a minimizer of the problem (3).  □

4Numerical Method

4.1Discretization Scheme

Our scheme allows to perform both deblurring and denoising simultaneously. To do so, we need to compute:

(9)
u∗=arg minu∈S(Ω)E(u),E(u)=∫Ωα(x)|∇u|dx+λ12∫Ω(Ku−f)2dx+λ2∫Ω(Ku−flogKu)dx,
where K is a blurring operator (convex), f is the observed image, S(Ω) is the set of positive functions defined over Ω with bounded total variation, and λ1,λ2 are positive regularization parameters. This functional u↦E(u) 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 R2→R. Therefore we have to choose a discretization scheme for numerical computations. If u is a image, we write uj,k for the pixel at coordinates (j,k) in u. We define the following quantities:

∇1uj,k=uj+1,k−uj−1,k,∇2uj,k=uj,k+1−uj,k−1,∇uj,k=(∇1uj,k,∇2uj,k),|∇uj,k|=(∇1uj,k)2+(∇2uj,k)2+ε2,
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
δE(u)δuj,k=div(∇uj,k|∇uj,k|)−λ1KT(Kuj,k−fj,k)−λ2(K−fj,kuj,k).
The operator divergence div(∇u|∇u|) is defined by
(∇11u)(∇2u)2−2(∇1u)(∇2u)(∇12u)+(∇22u)(∇1u)2((∇1u)2+(∇1u)2+ε2)3/2,
where
∇11uj,k=∇1(∇1uj,k)=uj+1,k−2uj,k+uj−1,k,∇22uj,k=∇2(∇2uj,k)=uj,k+1−2uj,k+uj,k−1,∇12uj,k=∇1(∇2uj,k)=uj+1,k+1+uj−1,k−1−uj+1,k−1−uj−1,k+1.
Thus, for a small parameter δt>0, a solution of the minimization problem (2) may be computed by
u(t+1)−u(t)δt=div(α(x)(∇u(t)|∇u(t)|))−λ1KT(Ku−f)−λ2(K−fu).
When the time-step parameter δ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.2Proposed Algorithm

First, let us first review the split-Bregman algorithm (Goldstein and Osher, 2009). Suppose that we have a scalar γ and two convex functionals Ψ(·) and G(·); and that we need to solve the following constrained optimization problem:

(10)
findarg minu,d‖d‖1+γ2G(u),s.t.d=Ψ(u).

We convert (10) into an unconstrained problem:

(11)
findarg minu,d‖d‖1+γ2G(u)+ρ2‖d−Ψ(u)−b‖22,
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):
u(k+1)=arg minuγ2G(u)+ρ2‖d(k)−Ψ(u)−b(k)‖22,d(k+1)=arg mind‖d‖1+ρ2‖d−Ψ(u(k+1))−b(k)‖22,b(k+1)=b(k)+Ψ(u(k+1))−d(k+1).
Now we return to the problem (9). We define
G(u)=λ12(Ku−f)2+λ2(Ku−flogKu)andΨ(u)=α∇u.
We set ν=Ku; then, based on equation (11), the split-Bregman problem for (9) is defined as:
(12)
arg minu,d(‖d‖1+γ2G(ν)+ρ12‖ν−Ku−c‖22+ρ22∑i=1,2‖di−α∇iu−bi‖22),
where the parameters ρ1, ρ2 and γ are positive, d=(d1,d2), b=(b1,b2) and ∇u=(∇1u,∇2u).

The split-Bregman method for solving (12) is described as follows:

u(k+1)=arg minuρ12‖ν(k)−Ku−c(k)‖22+ρ22∑i=1,2‖di(k)−α∇iu−bi(k)‖22,ν(k+1)=arg minνγ2G(ν)+ρ12‖ν−Ku(k+1)−c(k)‖22,di(k+1)=arg mind‖di‖1+ρ22‖di−α∇iu(k+1)−bi(k)‖22,c(k+1)=c(k)+Ku(k+1)−ν(k+1),bi(k+1)=bi(k)+α∇iu(k+1)−di(k+1).

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

Subproblem 1. The u subproblem is a quadratic optimization problem, whose optimality condition reads:

(13)
(ρ1KT·K+ρ2α∑i=1,2∇iT∇i)u(k+1)=ρ1KT(ν(k)−c(k))+ρ2∑i=1,2∇iT(di(k)−bi(k)),
under considering periodic boundary conditions. Note that left-hand-side matrix in (13) includes a Laplacian matrix (∇1T∇1+∇2T∇2=−Δ) 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:
(14)
u=F−1(F(r)ρ1F(KT)·F(K)−ρ2F(α)·F(Δ)),
where
r=ρ1KT(ν(k)−b1(k))+ρ2∑i=1,2∇iT(di(k)−bi(k)),
F and F−1 are the forward and inverse Fourier transform operators.

Subproblem 2. The optimality condition for the ν subproblem is given by

γ2(λ1(ν−f)+λ2(1−fν))+ρ1(ν−Ku(k+1)−c(k))=0.
This equation canbe rewritten as:
(γ2λ1+ρ1)(ν(k+1))2−(γ2λ1f−λ2γ2+ρ1(Ku(k+1)+c(k)))ν(k+1)−γ2λ2f=0.
The positive solution is given by
(15)
ν(k+1)=S(k)+(S(k))2+γλ2fγλ1+2ρ1,
where
Sk=λ1γf−λ2γ+2ρ1(Ku(k+1)+bν(k))2(γλ1+2ρ1).

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 d1 and d2 separately:

(16)
di(k+1)=shrink(α∇iu(k+1)+bi(k),1ρ2).
The problem (16) is solved as:
(17)
di(k+1)=α∇iu(k+1)+bi(k)|α∇iu(k+1)+bi(k)|·max(|α∇iu(k+1)+bi(k)|−1ρ2,0).

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 ς:

‖u(k)−u(k−1)‖2‖u(k)‖2<ς,
where ς is a small positive parameter. For our experiments, we set tolerance ς=0.0005 and N=500.

Algorithm 1

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

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

5Numerical Simulations

5.1Implementation 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-L1 (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: λ1=0.4, λ2=0.6. We set ρ1=1, ρ2=1. The parameter σ in α(x) is set to 1. The threshold value l in the function α(x) and the parameters γ are chosen to keep the poise between noise removal and detail preservation capabilities.

The test images11 are 8-bit gray scale standard images of size 256×256 shown in Fig. 1.

Fig. 1

Original images.

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

PSNR=10log10(MNImax2‖u∗−u‖22),
where u,u∗ are, respectively, the original image and the reconstructed (or noisy) image, Imax 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

SSIM(u,u∗)=(2μuμu∗+c1)(2σu,u∗+c2)(μu2+μu∗2+c1)(σu2+σu∗2+c2),
where μu, μu∗ are the means of u, u∗, respectively; σu,σu∗, their standard deviations; σu,u∗, the covariance of two images u and u∗; c1=(K1L)2; c2=(K2L)2; L is the dynamic range of the pixel values (255 for 8-bit grayscale images); and finally K1≪1, K2≪1 are small constants.

5.2Numerical Results and Discussion

5.2.1Image 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 Imax and by Gaussian noise with standard deviation σG to the test images. In Figs. 2, 4 and 5, we give the results for denoising the corrupted images for different noise levels Imax and σG=10.

Fig. 2

Recovered results for the test images. (a) Noisy image with Imax=120, σG=10, (b) TV L1, (c) Bregman, (d) MS-MPG, (e) Our proposed.

Recovered results for the test images. (a) Noisy image with Imax=120, σG=10, (b) TV L1, (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-L1 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.

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-L1 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 Imax=60, σG=10, (b) TV L1, (c) Bregman, (d) MS-MPG, (e) Ours.

Recovered results for the test images. (a) Noisy image with Imax=60, σG=10, (b) TV L1, (c) Bregman, (d) MS-MPG, (e) Ours.
Fig. 5

Recovered results for the test images. (a) Noisy image with Imax=60, σG=10, (b) TV-L1, (c) Bregman, (d) MS-MPG, (e) Ours.

Recovered results for the test images. (a) Noisy image with Imax=60, σG=10, (b) TV-L1, (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 L1 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.

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 L1 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 L1 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.

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 L1 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 L1. 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 Imax and σG=10.

ImageMethodCPU time (s)
Imax=120Imax=60
TV L14.34495.6730
ClockBregman0.94600.8212
MS-MPG4.14654.8734
Ours1.09451.1081
TV L15.62297.4171
CocoBregman1.02650.8414
MS-MPG4.08445.0879
Ours1.12391.2251
TV L14.30966.4129
LampBregman0.92250.9473
MS-MPG4.18104.8758
Ours0.94311.1266
Table 2

PSNR values and SSIM measures for noisy images and recovered images with Imax=120.

ImagePSNRMSSIM
NoisyTV L1BregmanMS-MPGOursNoisyTV L1BregmanMS-MPGOurs
Imax=120, σ=10
Jetplane18.941622.720324.119024.784825.32510.40450.70610.75140.75110.7748
Lake19.641321.367522.590622.997224.47980.52350.63600.68120.70690.7603
Aerial17.447118.955019.584019.305119.88060.55820.50830.58080.57110.7130
Clock18.385224.604025.794524.884426.12010.29970.83390.88220.77960.8970
Car19.138521.469422.155922.879324.06200.48480.61060.65420.68040.7256
Coco16.911920.424220.421520.342620.65390.27550.85510.87980.82960.8950
Lamp17.877024.280824.359424.106224.63390.24460.85220.88910.78890.8985
Poulina18.838125.256725.720325.978126.06530.32500.76480.79340.79820.8074
Spine21.000425.256124.685525.534926.10100.61800.79250.77630.79670.8206
Head21.778724.356726.234826.906127.09790.63240.80330.80430.82730.8400
Average18.996022.869123.566623.771924.44200.43660.73630.76930.75300.8132
Imax=120, σ=15
Jetplane16.715022.203323.491523.691824.14150.33200.67610.72480.69590.7320
Lake17.257420.821522.082722.226023.04420.43840.60210.67320.67090.7040
Aerial15.800618.767119.279519.106019.47060.46220.45940.57400.51390.6472
Clock16.461924.216525.364524.237125.77400.24400.81050.86010.81860.8805
Car16.858920.951221.773522.126922.76080.40150.58090.64020.63380.6727
Coco15.419320.339820.410920.133220.54880.21810.82650.85990.77410.8789
Lamp16.046123.897223.909023.516924.30630.19640.82250.86950.72100.8799
Poulina16.662724.919525.270925.275325.41420.24520.73460.76590.74910.7739
Spine18.558223.730124.401524.312224.92720.53780.74180.76890.75210.7794
Head19.351224.54925.419925.835625.98930.55880.75670.78360.78540.7991
Average16.913122.439523.140423.046123.63770.36340.70110.75200.71150.7748
Table 3

PSNR values and SSIM measures for noisy images and recovered images with with Imax=60.

ImagePSNRMSSIM
NoisyTV L1BregmanMS-MPGOursNoisyTV L1BregmanMS-MPGOurs
Imax=60, σ=10
Jetplane14.092921.451522.511622.370522.80570.25700.63960.64820.67300.6854
Lake14.719020.148020.988520.733521.55860.34880.55670.59870.59450.6325
Aerial13.909118.712218.938618.692919.28980.34650.40360.52960.38010.5635
Clock13.994123.755424.760724.316625.06820.18660.77590.78650.79310.8439
Car14.239320.339020.899320.898821.49200.31240.54170.57230.57090.5864
Coco13.437319.960919.908220.045920.26650.15730.79690.78150.82180.8535
Lamp13.623523.311823.287023.489223.61010.14660.78980.78230.80670.8568
Poulina14.169224.242924.876824.870424.92720.18040.69010.71690.72520.7316
Spine16.091022.74923.398123.326623.50110.45970.68210.72860.71530.7308
Head16.971823.66724.299124.278024.47630.49700.70590.74940.72840.7550
Average14.524721.833822.386822.302222.69960.28920.65820.68940.68090.7240
Imax=60, σ=15
ImageNoisyTV L1BregmanMS-MPGOursNoisyTV L1BregmanMS-MPGOurs
Jetplane11.431420.560421.031721.272721.37290.19110.58830.52080.62470.6319
Lake12.045019.367619.778919.910220.09110.24410.50530.52090.55260.5545
Aerial11.621618.402118.900118.563219.14820.24250.34350.42160.35180.4362
Clock11.450622.991423.673723.425024.31870.13650.72970.63870.72980.8123
Car11.535419.603119.870520.108120.24980.21630.48980.47760.52540.5311
Coco11.147719.669419.180919.758019.83150.11490.74320.60100.76280.8227
Lamp11.118222.673422.100522.814522.95510.10140.73410.61510.74300.8263
Poulina11.592723.490423.839823.880824.00400.12570.63530.61060.68180.6960
Spine13.455120.808521.968222.012922.11220.38440.61150.64930.65480.6658
Head14.344222.479922.495422.510522.96980.43700.64450.68900.68530.6991
Average11.974221.004621.284021.425621.70530.21940.60250.57450.63120.6676

5.2.2Image 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×9, and standard deviation of 1. After the blurring operation, we corrupt the images by Possion noise Imax=120 and σ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 L1 (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 L1, (c) Bregman, (d) MS-MPG, (e) Our proposed.

Recovered results for the test images. (a) Blurring and Noisy image, (b) TV L1, (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 L1 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.

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 L1 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 L1, (c) Bregman, (d) MS-MPG, (e) Our proposed.

Recovered results for the test images. (a) Blurring and Noisy image, (b) TV L1, (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 L1 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.

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 L1 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 Imax=120, σ=15.

ImagePSNRMSSIM
NoisyTV L1BregmanMS-MPGOursNoisyTV L1BregmanMS-MPGOurs
Imax=120, σ=15
Jetplane14.952218.707918.7201218.042019.00290.22820.63840.66000.58830.6860
Lake16.053519.641919.610019.647220.26750.28760.55060.54490.56540.6090
Aerial15.370118.332518.854918.749518.99210.31070.46470.49020.49600.5030
Clock16.189123.234823.589822.589323.71330.17610.77580.81960.65750.8313
Car15.590519.620219.660019.377420.14860.24080.53750.54860.52860.5863
Coco15.382920.147920.176219.002520.35720.14100.80820.84680.75240.8608
Lamp15.947723.400523.631521.665723.76350.12960.80010.85970.69190.8679
Poulina15.347519.551819.680120.270520.43920.17100.68710.70990.69230.7205
Spine16.047619.190718.679718.884719.35440.38650.56940.58320.60510.6286
Head14.981216.788816.699116.604418.34810.45900.63520.65620.67110.7061
Average15.586219.861719.930119.483320.43870.25310.64670.67190.62490.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.

6Conclusion

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.

Acknowledgements

The authors would like to thank professor S.D. Dvoenko and professor A.V. Kopylov, Tula State University, Tula, Russia, for their advice and comments.

References

1 

Aubert, G., Aujol, J. ((2008) ). A variational approach to remove multiplicative noise. SIAM Journal on Applied Mathematics, 68: (4), 925–946.

2 

Barcelos, C.A.Z., Chen, Y. ((2000) ). Heat flows and related minimization problem in image restoration. Computers & Mathematics with Applications, 39: (5–6), 81–97.

3 

Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J. ((2010) ). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning, 3: , 1–122.

4 

Calatroni, L., De Los Reyes, J.C., Schönlieb, C.B. ((2017) ). Infimal convolution of data discrepancies for mixed noise removal. SIAM Journal on Imaging Sciences, 10: (3), 1196–1233.

5 

Chambolle, A. ((2004) ). An algorithm for total variation minimization and applications. Journal of Mathematical Imaging and Vision, 20: , 89–97.

6 

Chambolle, A., Caselles, V., Novaga, M., Cremers, D., Pock, T. ((2010) ). An introduction to total variation for image analysis. In: Theoretical Foundations and Numerical Methods for Sparse Recovery, Radon Series on Computational and Applied Mathematics, Vol. 9: , pp. 263–340.

7 

Chan, T.F., Shen, J. ((2005) ). Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods. Society for Industrial and Applied Mathematics. 400 pp.

8 

Chouzenoux, E., Jezierska, A., Pesquet, J.C., Talbot, H. ((2015) ). A convex approach for image restoration with exact Poisson–Gaussian likelihood. SIAM Journal on Imaging Sciences, 8: (4), 2662–2682.

9 

De Los Reyes, J.C., Schönlieb, C.B. ((2013) ). Image denoising: learning the noise model via nonsmooth PDE-constrained optimization. Inverse Problems & Imaging, 7: (4), 1183–1214.

10 

Dong, Y., Zeng, T. ((2013) ). A convex variational model for restoring blurred images with multiplicative noise. SIAM Journal on Imaging Sciences, 6: (3), 1598–1625.

11 

Foi, A., Trimeche, M., Katkovnik, V., Egiazarian, K. ((2008) ). Practical Poissonian–Gaussian noise modeling and fitting for single-image raw-data. IEEE Transactions on Image Processing, 17: (10), 1737–1754.

12 

Goldstein, T., Osher, S. ((2009) ). The split Bregman method for L1-regularized problems. SIAM Journal on Imaging Sciences, 2: (2), 323–343.

13 

Jezierska, A., Chaux, C., Pesquet, J.C., Talbot, H. ((2011) ). An EM approach for Poisson–Gaussian noise modeling In: 19th European Signal Processing Conference (EUSIPCO), Barcelona, Spain, pp. 2244–2248.

14 

Kornprobst, P., Deriche, R., Aubert, G. ((1999) ). Image sequence analysis via partial differential equations. Journal of Mathematical Imaging and Vision, 11: , 5–26.

15 

Lanza, A., Morigi, S., Sgallari, F., Wen, Y.W. ((2014) ). Image restoration with Poisson–Gaussian mixed noise. Computer Methods in Biomechanics and Biomedical Engineering: Imaging and Visualization, 2: (1), 12–24.

16 

Le Montagner, Y., Angelini, E.D., Olivo-Marin, J.C. ((2014) ). An unbiased risk estimator for image denoising in the presence of mixed Poisson–Gaussian noise. IEEE Transactions on Image Processing, 23: (3), 1255–1268.

17 

Le, T., Chartrand, R., Asaki, T.J. ((2007) ). A variational approach to reconstructing images corrupted by Poisson noise. Journal of Mathematical Imaging and Vision, 27: , 257–263.

18 

Li, J., Shen, Z., Yin, R., Zhang, X. ((2015) ). A reweighted l2 method for image restoration with Poisson and mixed Poisson–Gaussian noise. Inverse Problems & Imaging, 9: (3), 875–894.

19 

Li, J., Luisier, F., Blu, T. ((2018) ). PURE-LET image deconvolution. IEEE Transactions on Image Processing, 27: (1), 92–105.

20 

Liu, L., Chen, L. Chen, C.L.P., Tang, Y.Y., Man Pun, C. ((2017) ). Weighted joint sparse representation for removing mixed noise in image. IEEE Transactions on Cybernetics, 47: (3), 600–611.

21 

Luisier, F., Blu, T., Unser, M. ((2011) ). Image denoising in mixed Poisson–Gaussian noise. IEEE Transactions on Image Processing, 20: (3), 696–708.

22 

Makitalo, M., Foi, A. ((2013) ). Optimal inversion of the generalized Anscombe transformation for Poisson–Gaussian noise. IEEE Transactions on Image Processing, 22: (1), 91–103.

23 

Marnissi, Y., Zheng, Y., Pesquei, J.C. ((2016) ). Fast variational Bayesian signal recovery in the presence of Poisson–Gaussian noise. In: IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3964–3968.

24 

Micchelli, C.A., Shen, L., Xu, Y. ((2011) ). Proximity algorithms for image models: denoising. Inverse Problems, 27: (4), 045009.

25 

Pham, C.T., Kopylov, A. ((2015) ). Multi-quadratic dynamic programming procedure of edge-preserving denoising for medical images. ISPRS-International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XL–5/W6: , 101–106.

26 

Pham, C.T., Gamard, G., Kopylov, A., Tran, T.T.T. ((2018) ). An algorithm for image restoration with mixed noise using total variation regularization. Turkish Journal of Electrical Engineering and Computer Sciences, 26: (6), 2831–2845.

27 

Pham, C.T., Tran, T.T.T., Phan, T.D.K., Dinh, V.S., Pham, M.T., Nguyen M.H. ((2019) ). An adaptive algorithm for restoring image corrupted by mixed noise. Journal Cybernetics and Physics, 8: (2), 73–82.

28 

Rudin, L., Osher, S., Fatemi, E. ((1992) ). Nonlinear total variation-based noise removal algorithms. Physica D, 60: , 259–268.

29 

Sciacchitano, F., Dong, Y., Zeng, T. ((2015) ). Variational approach for restoring blurred images with Cauchy noise. SIAM Journal on Imaging Sciences, 8: (3), 1894–1922.

30 

Tran, T.T.T., Pham, C.T., Kopylov, A.V., Nguyen, V.N. ((2019) ). An adaptive variational model for medical images restoration. ISPRS – International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XLII–2/W12: , 219–224.

31 

Wang, Z., Bovik, A.C. ((2006) ). Modern Image Quality Assessment, Synthesis Lectures on Image, Video, and Multimedia Processing. Morgan and Claypool Publishers. 156 pp.

32 

Wang, Y., Yang, J., Yin, W., Zhang, Y. ((2008) ). A new alternating minimization algorithm for total variation image reconstruction. SIAM Journal on Imaging Sciences, 1: (3), 248–272.

33 

Yan, M. ((2013) ). Restoration of images corrupted by impulse noise and mixed Gaussian impulse noise using blind inpainting. SIAM Journal on Imaging Sciences, 6: (3), 1227–1245.

34 

Zhang, B., Fadili, M.J., Starck, J.L., Olivo-Marin, J.C. ((2007) ). Multiscale variance-stabilizing transform for mixed-Poisson–Gaussian processes and its applications in bioimaging. In: IEEE International Conference on Image Processing (ICIP), Vol. 6, 233–236.

Show more