# A Fast High-Order Total Variation Minimization Method for Multiplicative Noise Removal

Authors: Xiao-Guang Lv; Jiang Le; Jin Huang; Liu Jun
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Multiplicative noise removal problem has  received considerable attention in recent years. The total variation regularization method for the solution of the noise removal problem can preserve edges well but has the sometimes undesirable staircase effect. In this paper, we propose a fast high-order total variation minimization method to restore multiplicative noisy images. The proposed method is able to preserve edges and at the same time avoid the staircase effect in the smooth regions. An alternating minimization algorithm is employed to solve the proposed high-order total variation minimization problem. We discuss the convergence of the alternating minimization algorithm. Some numerical results show that the proposed method gives restored images of higher quality than some existing multiplicative noise removal methods.


---

## Body

## 1. Introduction

Image denoising problem has been widely studied in the areas of image processing. The problem includes additive noise removal and multiplicative noise removal. Most of the literature deals with the additive noise model [1–4]. Given a noisy image g=u+v,  where u is the original image and v is the noise, the denoising problem is to recover u from the observed image g. The additive noise model has been extensively studied. We refer the reader to [5–11] and references therein for a review of various methods.In this paper, we consider the problem of seeking the original image from a multiplicative noisy image. In the multiplicative noise model, the recorded imageg is the multiplication of the original image u and the noise v:
(1)g=uv.
This problem arises in ultrasound imaging, synthetic aperture radar (SAR) and sonar (SAS), laser imaging, and magnetic field inhomogeneity in MRI [12, 13]. In this work, we concentrate on the assumption that v follows a Gamma distribution, which commonly occurs in SAR. With respect to a given resolution cell of the imaging device, the complex SAR image is computed after the SAR system receives the coherent sum of reflected monochromatic microwaves. The complex SAR image at one point is usually modeled by a complex zero-mean circular Gaussian density (i.e., the real and imaginary parts are independent Gaussian variable with a common variance). The intensity image is defined as the square of the magnitude of the corresponding complex SAR image;  its variables follow exponentially distribution independently [14, 15]. Since the complex observation is acquired in highly coherent system, the intensity image occurs to have peculiar granular appearance (known as speckle). The conventional SAR system reduces the speckle effect of the intensity image by the multilooking process (L-look, in the case of L-looks), that is, averaging the observed intensity images from slightly different angles of the same resolution cell [16].Due to the coherent nature of these image acquisition processes, nearly all the information of the original image may vanish when it is distorted by the multiplicative noise. In general, the standard additive denoising method, so prevalent in image processing, is inadequate. Therefore, it is necessary to devise efficient and reliable algorithms for recovering the true image from the observed multiplicative noisy image. Until the past decade, a few variational approaches devoted to multiplicative noise removal have been proposed. The early variational approach for multiplicative noise removal is the one by Rudin et al. [17]. According to the statistical properties of the multiplicative noise v, the recovery of the image u was based on solving the following constrained optimization problem [17]:
(2)minu∫Ω‍|Du|dxdysubjectto:∫Ω‍gudxdy=1,subjectto:∫Ω‍(gu-1)2dxdy=σ2,
where ∫Ω‍|Du|dxdy is denoted by the seminorm on boundary variation (BV) space that coincides with  ∫Ω‍ux2+uy2dxdy  when u is smooth. The two constraints state that the mean of the noise is equal to 1, and the variance is equal to  σ2. In (2), only basic statistical properties, the mean, and the variance of the noise v are considered, which somehow limits the restored results.By using a maximum a posteriori (MAP) estimator, Aubert and Aujol [18] proposed a function whose minimizer corresponds to the denoised image to be recovered. This function is
(3)minu∫Ω‍(logu+gu)dxdy+α∫Ω‍|Du|dxdy,
where the total variation of u is utilized as the regularization term, and α is the regularization parameter which controls the trade-off between a good fit of  g and a smoothness requirement due to the total variation regularization. Although the functional they proposed is not convex, they still proved the existence of a minimizer and showed the capability of their model on some numerical examples.As a result of the drawback of the objective function (3) that is not convex for all u, the obtained solution is likely not the global optimal solution of (3). To overcome this problem, Huang et al. [19] proposed and studied a strictly convex objective function for multiplicative noise removal problem. In [19], the authors introduced an auxiliary variable z=logu in (3). With the new variable, the proposed unconstrained total variation denoising problem is described as follows:
(4)minz,w∫Ω‍(z+ge-z)dxdy+α1∥z-w∥L22+α2∫Ω|Dw|dxdy,
where α1 and α2 are positive regularization parameters. The parameter α1 measures the trade-off between an image obtained by a maximum likelihood estimation from the first term and a total variation denoised image w. The parameter α2 measures the amount of regularization to a denoised image w. The main advantage of the proposed method is that the total variation norm can be used in the noise removal process in an efficient manner. Therefore the method proposed in [19] has the ability to preserve edges very well in the denoised image.Recently, Bioucas-Dias and Figueiredo [20] proposed an efficient multiplicative noise removal method by using variable splitting and constrained optimization. After converting the multiplicative noise model into an additive one by taking logarithms, they used variable splitting to obtain an equivalent constrained problem and then solved this optimization problem by using the augmented Lagrangian method. A set of experiments has shown that the proposed method, which they named MIDAL (multiplicative image denoising by augmented Lagrangian), yields state-of-the-art results both in terms of speed and denoising performance. In [21], Steidl and Teuber considered a variational restoration model consisting of the I-divergence as data fitting term and the total variation seminorm or nonlocal means as regularizer for removing multiplicative noise. They proposed to compute the minimizers of their restoration functionals by using Douglas-Rachford splitting techniques. For a particular splitting, they presented a semi-implicit scheme to solve the involved nonlinear systems of equations and proved its Q-linear convergence.It is known that the total variation denoising method is a PDE-based technique that preserves edges well but has the sometimes undesirable staircase effect, namely, the transformation of smooth regions (ramps) into piecewise constant regions (stairs) [22–26]. Therefore, the approaches involving the traditional total variation often cause staircase effect since they favor solutions that are piecewise constant. Staircase solutions fail to satisfy the evaluation of visual quality and they can develop false edges that do not exist in the true image. To attenuate the staircase effect, there is a growing interest in the literature for replacing the total variation norm by the high-order total variation norm. The motivation behind such attempt is to restore potentially a wider of images, which comprise more than merely piecewise constant regions. The majority of the high-order norms involve second-order differential operators because piecewise-vanishing second-order derivatives lead to piecewise-linear solutions that better fit smooth intensity changes; see [27] for more details. There are two classes of second-order regularization methods for image restoration problems. The first class employs a second-order regularizer in a standalone way. For example, in [28], the authors considered a fourth-order partial differential equations (PDE) model for noise removal and employed the dual algorithm of Chambolle for solving the high-order problems. In [29], Hu and Jacob applied higher degree total variation (HDTV) regularization for image recovery. The second class is to combine the TV norm with a second-order regularizer. For example, a technique in [30, 31] combining the TV filter with a fourth-order PDE filter was proposed to preserve edges and at the same time avoid the staircase effect in smooth regions for noise removal. In [32], Papafitsoros and Schönlieb considered a high-order model involving convex functions of the TV and the TV of the first derivatives for image restoration problems and used the split Bregman method to solve numerically the corresponding discretized problem.In this paper, we present a fast high-order total variation minimization method for multiplicative noise removal problem. This technique substantially reduces the staircase effect, while preserving sharp jump discontinuities (edges). An alternating minimization algorithm is employed to solve the proposed high-order total variation minimization model. We discuss the convergence of the proposed alternating minimization algorithm. Some numerical results show that the proposed method gives restored images of higher quality than some existing multiplicative noise removal methods.The organization of this paper is outlined as follows. In the next section, we present the high-order total variation minimization method for multiplicative noise removal problem. In Section3, an alternating minimization algorithm is employed to find the minimizer of the proposed minimization problem. We analyze the convergence of the proposed alternating minimization algorithm in Section 4. Some numerical experiments are given to illustrate the performance of the proposed algorithm in Section 5. Concluding remarks are given in the last section.


## 2. A High-Order Total Variation Multiplicative Denoising Model

In this section, our aim is to propose a strictly convex objective function for restoring images distorted by the multiplicative noise. The most important is that we apply a high-order total variation to the objective function to achieve sharp edges and staircase reduction efficiently. We introduce our multiplicative denoising model from the statistical perspective using the Bayesian formula.Letg, u, and v denote samples of instances of some random variables  F,  U, and V. Moreover, we assume that the random variable v on each pixel is mutually independent and identically distributed (i.i.d). Moreover, the random variable v in each pixel follows Gamma distribution; that is, its probability density function is
(5)PV(v;L,M)=MLvL-1Γ(L)e-Mvforv>0,
where  Γ(·) is the usual Gamma-function, and L and M denote the shape scale and inverse parameters in the Gamma distribution, respectively. We note that the mean of a gamma-distributed variable is L/M and the variance is L/M2. As the multiplicative noise, in this work we assume that the mean of  v equals 1, and then we have L=M.According to the posteriori estimation, the restored imageu can be determined by
(6)u=argmaxuPU∣G(u∣g).
From the Bayes rule, we have PU∣G(u∣g)=PG∣U(g∣u)PU(u)/PG(g). It then follows that
(7)u=argmaxuPG∣U(g∣u)PU(u)PG(g).
Using the proposition in [18]: PV(g/u)(1/u)=PG∣U(g∣u), we obtain
(8)PG∣U(g∣u)=LLgL-1uLΓ(L)e-Lg/u.
Taking the logarithm transformation into account, we assume that the image prior PU(u) is given by the following:
(9)PU(u)=PU∣W(u∣w)PW(w)
with
(10)PU∣W(u∣w)∝e-α1∥logu-w∥22,(11)PW(w)=e-α2∥w∥HTV,
where α1 and α2 are two positive scalars and  ∥w∥HTV is a high-order total variation of  w. Here we assume that the difference between logu and w follows a Gaussian distribution, and w follows a high-order total variation prior. Thus, maximizing PU∣G(u∣g) amounts to minimizing the log-likelihood:
(12)-log(PU∣G(u∣g))=-log(PG∣U(g∣u))-log(PU∣W(u∣w))-log(PW(w))+log(PG(g)).
Since  PG(g) is constant, (6) can be rewritten as the following optimization problem:
(13)u=argminu(-log(PG∣U(g∣u))-log(PU∣W(u∣w))=argminu-log(PW(w))).
The previous computation leads us to propose the following high-order total variation model for restoring images corrupted by the Gamma noise by introducing a new variable z=logu:
(14)minz,w𝒥(z,w)≐minz,w∫Ω‍(z+ge-z)dxdy+α1∥z-w∥L22+α2∫Ω‍|D2w|dxdy,
where α1 and α2 are positive regularization parameters. In the above model, D2w  denotes the Hessian of w, and |D2w|=wxx2+wxy2+wyx2+wyy2 is just the Frobenius norm of the Hessian D2w. The parameter α1 measures the trade-off between an image obtained by a maximum likelihood estimation from the first term and a high-order total variation denoised image w. The parameter α2 measures the amount of regularization to a denoised image w. The main advantage of the proposed method is that the high-order total variation norm can be used in the noise removal process in an efficient manner. Therefore, the proposed method has the ability to preserve edges very well and reduce staircase effect in the denoised image.Moreover, it is obvious that whenu includes an edge, z also contains an edge at the same point; that is, the logarithm transformation preserves image edges; see [19] for more details. Based on this idea, we can view z as an image in the logarithm domain. We can apply a high-order total variation regularization to z to recover its edges and reduce staircase effect. It is interesting to note that z in the proposed model can be positive, zero, or negative, and the corresponding  u=ez is still positive, although u in the objective function (3) should be positive. Especially, one can find that when z in the proposed model is a large negative value, u is just close to zero.


## 3. The Alternating Minimization Algorithm

In this section, an alternating minimization algorithm is employed to solve the proposed high-order total variation minimization problem efficiently. To solve (14), we need to consider the following two minimization subproblems.(i)
Forw fixed, determine the solution of
(15)minz∫Ω‍(z+ge-z)dxdy+α1∥z-w∥L22.(ii)
Forz fixed, determine the solution of
(16)minwα1∥z-w∥L22+α2∫Ω‍|D2w|dxdy.In the discrete setting, minimizing the first problem amounts to solve the following optimization problem:
(17)minz∑1n2‍(zi+ge-zi)+α1∥z-w∥22.
Straightforward computations can show that the minimizer of the above problem is the solution of the following nonlinear systems:
(18)1-ge-zi+2α1(zi-wi)=0fori=1,2,…,n2.
There are n2 decoupled nonlinear equations to be solved. The second derivative with respect to z of the first term in (17) is equal to ge-z, which is always greater than zero. Here we assume that g>0 in the multiplicative noise model (1). Therefore, the first term of the objective function is strictly convex for all z. Hence, we know that the objective function (17) is also strictly convex. Therefore, the corresponding nonlinear equation (18) has a unique solution and the solution can be determined by using the Newton method very efficiently.Subproblem (16) is a high-order total variation regularization process for image denoising. Some effective algorithms for the minimizer of the high-order total variation problem have been proposed in the literature. For example, a lot of PDE-based methods are widely used to preserve edges and reduce the staircase effect for the high-order total variation problem; see [6, 22, 23, 30] for more details. In particular, inspired by the work from Chambolle [5], a dual algorithm and its convergence analysis for minimization of the high-order model are presented in [24, 28]. The algorithm overcomes the numerical difficulties related to the nondifferentiability of the high-order model. Some numerical results show that the convergence of the dual algorithm is faster than the gradient descent time-marching algorithm.In the present paper, we employ the dual algorithm proposed by [24, 28] for the high-order total variation denoising problem (16). The discrete version of  problem (16) is represented as
(19)minwα1α2∥z-w∥22+∥w∥HTV.
In the dual algorithm, we need to solve the following minimization problem:
(20)minp∥z-βdiv2p∥,subjectto|pj,k|≤1,∀1≤j,k≤n,
where β=α2/2α1 and p=(p11p12p21p22) with phk∈Rn×n for h,k=1,2. This problem can be solved by a semi-implicit fixed point iteration:
(21)pj,k(i+1)=pj,k(i)-τAj,k(i)1+τ|Aj,k(i)|,
where p0=0, Aj,k(i)=(∇2(div2p(i)-z/β))j,k, and τ≤1/64 is the step size. The operator ∇2 is defined as ∇2w=(wxxwxywyxwyy). Here,  div2pj,k(i)=Dxx((p11)j,k(i))+Dyx-((p12)j,k(i))+Dxy+((p21)j,k(i))+Dyy((p22)j,k(i)); see [24, 28] for more details. After the minimizer p* of the constrained optimization problem in (20) is determined, the denoised image w can be computed by w=z-βdiv2p*. In addition, it is shown that if τ≤1/64, then div2p(i) is convergent and p(i) with the dual method is the solution of the minimization problem (20) as i→∞.Starting from an initial guessw(0), the alternating minimization method computes a sequence of iterates
(22)z(1),w(1),z(2),w(2),…,z(k),w(k),…
such that
(23)z(k)≐𝒮h(w(k-1))=argminz∑1n2‍(zi+ge-zi)+α1∥z-w(k-1)∥22,w(k)≐𝒮p(z(k))=argminwα1α2∥z(k)-w∥22+∥w∥HTV
for k=1,2,….We are now in a position to describe the alternating minimization algorithm for the multiplicative noise removal problem.Algorithm 1.
The alternating minimization algorithm for multiplicative noise removal:
Choose an initial guessw(0), and k=1.(1)
Forw(k-1) fixed, employ the Newton iterative method to compute
(24)z(k)=argminz∑1n2‍(zi+ge-zi)+α1∥z-w(k-1)∥22.(2)
Forz(k) fixed, use a fixed point iteration to compute
(25)p*=argminp∥z(k)-βdiv2p∥and obtain(26)w(k)=z(k)-βdiv2p*.(3)
Check the stopping criteria. If a stopping criterion is satisfied, then exit with an approximate restored imageu(k)=ez(k); otherwise, let  k=k+1 and go to step (1).


## 4. Convergence Analysis

The main aim of this section is to analyze the convergence of the proposed alternating minimization algorithm. We first remark that𝒥(z,w) in (14) is strictly convex, as the sum of a convex function and of two strictly convex functions;  see [33] for  more details. Hence, it suffices to show that there exists a unique minimizer of the objective function 𝒥(z,w). We denote by 𝒥* the minimum value of 𝒥(z,w). In the following  we show that our alternating minimization algorithm gives asymptotically the solution of the discrete problem associated with (14). We have the following theorem.Theorem 2.
For any initial guessw(0)∈ℝn2, (z(k),w(k)) generated by Algorithm 1 converges to a minimizer of the objective function 𝒥(z,w) as  k→∞.Proof.
It follows from the alternating iterative process in Algorithm1 that
(27)𝒥(z(k),w(k))≥𝒥(z(k+1),w(k))≥𝒥(z(k+1),w(k+1)).
It is obvious that the sequence𝒥(z(k),w(k)) is nonincreasing. We note that it is bounded from below by the minimum value 𝒥*. Let θ1(z)=∑1n2‍(zi+ge-zi), θ2(w)=∥w∥HTV and DΦ(z,w)=α1∥z-w∥22. Since θ1 and θ2 are convex, Φ is convex and differentiable and DΦ is the associated Bregman distance, we know from [34, 35] that 𝒥(z,w)=θ1(z)+θ2(w)+DΦ(z,w) has the five-point property. Hence, 𝒥(z(k),w(k)) thus converges in  ℝ to 𝒥*. So we have
(28)limk→∞𝒥(z(k),w(k))=𝒥*.
We note that if the objective function  𝒥(z,w) is coercive, the sequence (z(k),w(k)) is bounded since the sequence 𝒥(z(k),w(k)) converges.
Now we need to show that the objective function𝒥(f,u,v) is coercive. Let Dxx, Dxy, Dyx, and Dyy be the second order difference matrices in the xx, xy, yx, and yy direction, respectively, and
(29)S=(DxxDxyDyxDyy).
It is not difficult to obtain the lower bound of the discrete high-order total variation as follows:
(30)∥w∥HTV=∑1≤j,k≤n‍((wxx)j,k)2+((wxy)j,k)2+((wyx)j,k)2+((wyy)j,k)2≥12∑1≤j,k≤n‍|(wxx)j,k|+|(wxy)j,k|+|(wyx)j,k|+|(wyy)j,k|=12∥Sw∥1.
Let  𝒥1(z,w)=α1∥z-w∥22 and 𝒥2(z,w)=(α2/2)∥Sw∥1. Denote
(31)Ω={(z,w)𝒥1(z,w)+𝒥2(z,w)=0}.
If  (z,w)∈Ω, we obtain  z=w. As ∥(z,w)∥2→∞, it is easy to show that
(32)𝒥(z,w)=∑in2‍(zi+ge-zi)⟶∞.
On the other hand, we obtain
(33)∑in2‍(zi+ge-zi)≥∑in2‍(loggi+1)
since it is a strictly convex function. By using the above inequality, we have
(34)𝒥(z,w)≥∑in2‍(loggi+1)+𝒥1(z,w)+𝒥2(z,w).
Hence, if (z,w)∉Ω with ∥(z,w)∥2→∞, it is not difficult to obtain that 𝒥1(z,w)+𝒥2(z,w)→∞. Thus, we get that 𝒥(z,w) also tends to infinity. It follows from the definition of coercive that 𝒥(z,w) is coercive.
Since the sequence(z(k),w(k)) is bounded, we can thus extract a convergent subsequence {(z(ki),w(ki))} from {(z(k),w(k))} such that limi→∞{(z(ki),w(ki))}=(z^,w^). Then we obtain m=𝒥(z^,w^). Moreover, we have, for all  ki∈N and all  z,
(35)J(z,w(ki))≥𝒥(z(ki+1),z(ki)),
and for all ki∈ℕ and all w,
(36)𝒥(z(ki+1),w)≥𝒥(z(ki+1),w(ki+1)).
Let us denote byz~ a cluster point of  {z(ki+1)}. We may immediately obtain from (27) that
(37)𝒥*=𝒥(z^,w^)=𝒥(z~,w^).
Thus, we have that ∑in2‍(z^i+ge-z^i)+α1∥z^-w^∥22=∑in2‍(z~i+ge-z~i)+α1∥z~-w^∥22. We conclude that z~=z^ since z~ and z^ are the minimizers of the strictly convex function 𝒮h(w^) defined in (23). Hence, limki→∞z(ki+1)=z^. Following a similar analysis, we can show that limki→∞w(ki+1)=w^. By passing to the limit in (35) and (36), we have
(38)𝒥(z,w^)≥𝒥(z^,w^),∀z,𝒥(z^,w)≥𝒥(z^,w^),∀w.
We know that (38) can be rewritten as follows:
(39)𝒥(z^,w^)=minz𝒥(z,w^),𝒥(z^,w^)=minw𝒥(z^,w).
From the definition of  𝒥(z,w), (39) is equivalent to the following:
(40)0∈ψ(z^,w^),0∈ϕ(z^,w^),
where ψ(z,w)=1-ge-z+2α1(z-w) and ϕ(z,w)=2α1(z-w)+α2∂(∥w∥HTV). The subdifferential of  𝒥(z,w) at (z^,w^) is given by
(41)∂𝒥(z^,w^)=(ψ(z^,w^)ϕ(z^,w^)).
Therefore, we have
(42)(00)∈∂𝒥(z^,w^),
which implies that (z^,w^) is the minimizer of 𝒥(z,w). Hence the whole sequence 𝒥(z(k),w(k)) converges towards 𝒥*, the unique minimum of  𝒥(z,w). We deduce that the sequence (z(k),w(k)) converges to (z^,w^), the minimizer of  𝒥(z,w)   as k tends to infinity.


## 5. Numerical Experiments

In this section we provide numerical results from multiplicative noise removal problem to demonstrate the performance of the proposed regularization method. All computations of the present paper were carried out in Matlab 7.10. The results were obtained by running the Matlab codes on an Intel Core i3 CPU (2.27 GHz, 2.27 GHz) computer with RAM of 2048 M. The initial guess is chosen to be the noisy image in all tests.The quality of the restoration results with different methods is compared quantitatively by using the Peak-Signal-to-Noise Ratio (PSNR), the relative error (ReErr) and Structural SIMilarity index (SSIM). They are defined as follows:(43)PSNR=20log10(255n∥u(k)-utrue∥2),ReErr=∥u(k)-utrue∥2∥utrue∥2,
where utrue and u(k) are the ideal image and the restored image of order n respectively, and
(44)SSIM=(2μuμu~+C1)(2σuu~+C2)(μu2+μu~2+C1)(σu2+σu~2+C2),
where μu and μu~ are averages of u and u~, respectively. σu and σu~ are the variance of u and u~, respectively. σuu~ is the covariance of  u and u~. The positive constants C1 and C2 can be thought of as stabilizing constants for near-zero denominator values. In general, a high PSNR-value or a small ReErr-value indicates that the restoration is more accurate. The SSIM is a well-known quality metric used to measure the similarity between two images. This method developed by Wang et al.  [36] is based on three specific statistical measures that are much closer to how the human eye perceives differences between images. The whiter SSIM map, the closer two images are. Therefore, we use the SSIM map to reveal areas of high or low similarity between two images in this work.We note that a straightforward high-order total variation regularization scheme for multiplicative noise removal is to solve the following optimization problem:(45)minu∫Ω‍(logu+gu)dxdy+α∫Ω‍|D2u|dxdy,
where α is a positive regularization parameter which measures the trade-off between a good fit and a regularized solution. The Euler-Lagrange equation for minimization problem (45) amounts to solve (omitting the boundary conditions):
(46)u-gαu2+(uxx|D2u|)xx+(uxy|D2u|)yx+(uyx|D2u|)xy+(uyy|D2u|)yy=0.
Hence, the time-marching (TM) algorithm solving the high-order total variation model is as follows:
(47)ut=g-uαu2-(uxx|D2u|+ϵ)xx-(uxy|D2u|+ϵ)yx-(uyx|D2u|+ϵ)xy-(uyy|D2u|+ϵ)yy,
where the parameter  ϵ=10-5 is introduced to avoid divisions by zero.In the first test, we compare our method with the straightforward high-order total variation regularization scheme. In the time-marching algorithm for solving the straightforward high-order total variation regularization problem, we set the step size as 0.1 in order to obtain a stable iterative procedure. The algorithm is stopped when the maximum number of iterations is reached. In addition, we carry out many experiments with differentα-values in the model (45), the ones with the best results are presented in this work. For the proposed high-order total variation method (HighTV method), we terminate the iterations for the method and accept u(k) as the computed approximation of the ideal image as soon as the maximum number of allowed outer iterations has been carried out or the relative differences between consecutive iterates u(1),u(2),u(3),… satisfy
(48)∥u(k+1)-u(k)∥2∥u(k+1)∥2<10-4.
In the proposed HighTV regularization method, there are two regularization parameters. We know that the quality of the restored image is highly depended on the regularization parameters. Similar to [19], we fix α1=19 in the proposed method in order to reduce the computational time in search for good regularization parameters. We determine the best value of α2 such that the relative error of such a restored image with respect to such an ideal image is the smallest.In the first experiment, we consider the “Cameraman” image of size256×256. The original image in Figure 1(a) is distorted by the Gamma noise with  L=10 and L=20. The noisy images are shown in Figures 1(b) and 1(c). It is seen from Figures 1(b) and 1(c) that the smaller L is, the more noisy the observed images are. When L=10, the relative error (ReErr) between the noisy image and the original image is 0.3150. When L=20, the ReErr between the noisy image and the original image is 0.2234.The related images in Test 1. (a) Original image. (b) Noisy image withL=10. (c) Noisy image with  L=20.
(a)
(b)
(c)In Figures2(a)–2(d), we show the restoration results of the two different methods for L=10 and L=20. According to the figures, we see that the proposed algorithm can provide better visual quality of restored images than those by the time-marching algorithm. In Figures 3(a) and 3(b), we compare the speed of convergence of the two algorithms for L=10 and L=20. In the figures, the x-axis is the number of iterations and the y-axis is the the relative error between the restored image and the original image. We see that the speed of convergence of the proposed algorithm is faster and the relative error of the proposed algorithm is smaller. In Table 1, we compare the results of the computational time and the the relative error when the best restored images are obtained. From Table 1, we observe that the computational time required by our method is lesser than that of the time-marching algorithm.Table 1
Restoration results with TM and High TV methods for Test 1.
Method
Parameter
CPU time
ReErr
L
=
10
TM
λ
=
240
12.1623
0.1150
HighTV
α
2
=
0.005
5.9285
0.1039
L
=
20
TM
α
=
800
12.0281
0.0921
HighTV
α
2
=
0.003
4.6362
0.0871Restored results for Test 1. (a) Restored image by the TM algorithm forL=10. (b) Restored image by the high-order method for L=10. (c) Restored image by the TM algorithm for  L=20. (d) Restored image by the high-order method for  L=20.
(a)
(b)
(c)
(d)Comparison of the TM algorithm and the proposed method for Test 1. (a)L=10. (b)  L=20.
(a)
(b)In the following, we compare the proposed HighTV method with the one proposed in [18] (AA method), the one proposed in [19] (HNW method) and the one proposed in [20] (MIDAL method). For the AA method, we use the time-marching algorithm with ϵ=10-4 to solve the minimization models as proposed in [18]. Similarly, we set the step size as 0.1 in order to obtain a stable iterative procedure. The algorithm is stopped when the maximum number of iterations is reached. In addition, we carry out many experiments with different α-values in the model (3), the ones  with the best results are presented in this work. For the HNW and MIDAL methods, we use the same stopping rule as the HighTV method.For the MIDAL method, the authors have verified experimentally that setting two parameter be equal yields good results. Therefore, we only selectα by searching for the value leading to the smallest relative error of such a restored image with respect to the ideal image. In the HNW method, there are two regularization parameters. In [19], the authors fix α1=19  in all tests in order to reduce the computational time in search for good regularization parameters. We use the best value of  α2 such that the relative error of such a restored image with respect to such an ideal image is the smallest. In this way, we have a fair comparison since we compare about the best restorations for four different methods.The aim of the second test is to give evidence of the effectiveness of employing the propose algorithm over the other three algorithms for multiplicative noise removal problems. Especially, we show that the proposed method can reduce the staircase effect. In this test, we consider the well-known “Lena” image of size256×256. The original image in Figure 4(a) is distorted by the Gamma noise with L=3 and L=13. The noisy images are shown in Figures 4(b) and 4(c). When L=3, the PSNR and the relative error (ReErr) between the noisy image and the original image are 11.95 dB and 0.5768. When L=13, the PSNR and the ReErr between the noisy image and the original image are 18.37 dB and 0.2765.The related images in Test 2. (a) Original image. (b) Noisy image withL=3. (c) Noisy image with  L=13.
(a)
(b)
(c)The restored images by the AA method, the HNW method, the MIDAL method and the proposed method are shown in Figures5(a)–5(h). From these figures, compared with the AA method, the HNW method and the MIDAL method, the proposed approach yields better results in image restoration since it avoids the staircase effect of the general total variation methods while at same time preserving edges as well as the total variation methods. In Figure 6, we have enlarged some details of the eight restored images. It is clear that the models involving the general total variation transform smooth regions into piecewise constant regions. As it is seen in the zoomed parts, the proposed method outperforms another three methods since the proposed model process smooth regions better than another three models. For the comparison of the performance quantitatively and the computational efficiency, in Table 2, we give their restoration results in SNRs and ReErrs. We observe from Table 2 that both the SNR and ReErr values of the restored images by the proposed method are better than those by the AA method, HNW method and the MIDAL method.Table 2
Restoration results with four different methods for Test 2 and Test 3.
Image
Method
Parameter
PSNR
ReErr
Lena (L=3)
AA
λ
=
15
22.96
0.1625
HNW
α
2
=
0.02
23.39
0.1543
MIDAL
α
=
2
23.85
0.1451
HighTV
α
2
=
0.012
24.11
0.1425
Lena (L=13)
AA
λ
=
300
26.97
0.1022
HNW
α
2
=
0.008
27.14
0.0992
MIDAL
α
=
4
27.31
0.0979
HighTV
α
2
=
0.004
27.62
0.0951
Barche (L=6)
AA
λ
=
80
23.22
0.1271
HNW
α
2
=
0.012
23.69
0.1223
MIDAL
α
=
3
23.70
0.1221
HighTV
α
2
=
0.008
23.88
0.1187
Barche (L=16)
AA
λ
=
15
25.14
0.1024
HNW
α
2
=
0.007
25.67
0.0973
MIDAL
α
=
4.5
25.66
0.0974
HighTV
α
2
=
0.0035
26.01
0.0939Restored results by different methods for Test 2. (a) Restored image by the AA method forL=3. (b) Restored image by the HNW method for  L=3. (c) Restored image by the MIDAL method for L=3. (d) Restored image by our method for  L=3. (e) Restored image by the AA method for  L=13. (f) Restored image by the HNW method for  L=13. (g) Restored image by the MIDAL method for L=13. (h) Restored image by our method for L=13.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)A small portion of Lena image is emphasized to show the difference of the four methods for Test 2. (a) Original image. (b) Noisy image withL=3. (c) Noisy image with L=13. (d) Restored image by the AA method for L=3. (e) Restored image by the HNW method for  L=3. (f) Restored image by the MIDAL method for L=3. (g) Restored image by our method for L=3. (h) Restored image by the AA method for L=13. (i) Restored image by the HNW method for L=13. (j) Restored image by the MIDAL method for L=13. (k) Restored image by our method for L=13.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)In the third test, the “barche” image is considered. The original image in Figure7(a) is distorted by the Gamma noise with  L=6 and L=16, see Figures 7(b) and 7(c). When L=6, the PSNR and the ReErr between the noisy image and the original image are 13.20 dB and 0.4090. When L=16, the PSNR and the ReErr between the noisy image and the original image are 17.44 dB and 0.2514.The related images in Test 3. (a) Original image. (b) Noisy image withL=6. (c) Noisy image with  L=16.
(a)
(b)
(c)The information of restored images by the AA method, the HNW method, the MIDAL method and the proposed method are displayed in Figures8(a)–8(h).  From the visual quality of restored images, the proposed regularization method is quite competitive with the other three algorithms. In addition, in Figures 9(a)–9(h) we show the SSIM maps of the restored images. It is not difficult to observe that the SSIM maps by the proposed method are whiter than the AA method, the HNW method and the MIDAL method, which means that the proposed method behaves slightly better. The restoration results of these four different methods in SNRs, ReErrs are also presented in Table 2. It is shown from Table 2 that both the SNR and ReErr values of the restored images by the proposed method are better than those by the AA method, the HNW method and the MIDAL method.Restored results by different methods for Test 3. (a) Restored image by the AA method forL=6. (b) Restored image by the HNW method for  L=6. (c) Restored image by the MIDAL method for L=6. (d) Restored image by our method for L=6. (e) Restored image by the AA method for L=16. (f) Restored image by the HNW method for L=16. (g) Restored image by the MIDAL method for L=16. (h) Restored image by our method for L=16.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Comparison of the SSIM maps for Test 3. (a) The SSIM map by the AA method forL=6. (b) The SSIM map by the HNW method for L=6. (c) The SSIM map by the MIDAL method for L=6. (d) The SSIM map by our method for L=6. (e) The SSIM map by the AA method for L=16. (f) The SSIM map by the HNW method for L=16. (g) The SSIM map by the MIDAL method for L=16. (h) The SSIM map by our method for L=16.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 6. Conclusions

In this paper, we investigate a fast and efficient way to restore blurred and noisy images with a high-order total variation minimization technique. An alternating minimization algorithm is employed to solve the high-order total variation model. The edges in the restored image can be preserved quite well and the staircase effect is reduced effectively in the proposed algorithm. Our experimental results show that the performance of the proposed method is superior to that of some existing total variation restoration methods.



---

*Source: 834035-2013-04-10.xml*