An Efficient Variational Method for Image Restoration

Abstract

Image restoration is one of the most fundamental issues in imaging science. Total variation regularization is widely used in image restoration problems for its capability to preserve edges. In this paper, we consider a constrained minimization problem with double total variation regularization terms. To solve this problem, we employ the split Bregman iteration method and the Chambolle’s algorithm. The convergence property of the algorithm is established. The numerical results demonstrate the effectiveness of the proposed method in terms of peak signal-to-noise ratio (PSNR) and the structure similarity index (SSIM).

1. Introduction

Image restoration is a fundamental problem in the literature of image processing. It plays an important role in various areas such as remote sensing, astronomy, medical imaging, and microscopy [1, 2]. Many image restoration tasks can be posed as linear inverse problems of the form: 𝑔=𝐻𝑓+𝑛, (1) where 𝐻∈ℝ^𝑛2×𝑛2 is a blurring matrix constructed from a discretized point spread function (PSF), 𝑓∈ℝ^𝑛2 is an original 𝑛×𝑛 gray-scale image, 𝑔∈ℝ^𝑛2 is a degraded observation, and 𝑛∈ℝ^𝑛2 is additive noise. We remark that the matrix 𝐻 has special structure that can be exploited in computations, when the special boundary conditions such as periodic and Dirichlet boundary conditions are imposed [3, 4]. In this work, the PSF is assumed to be known. In fact, if the PSF is unknown, there are a variety of means of techniques available for estimating it [5, 6].

Image restoration problems are frequently ill conditioned; thus, the straightforward solution of (1) typically does not yield a meaningful approximation [7, 8]. In order to avoid this difficulty, one typical method is to replace the linear system (1) by a nearby system that is less sensitive to the error 𝑛 in 𝑔 and considers the computed solutions of the latter system. This replacement is known as regularization.

One of the most popular regularization approaches is Tikhonov regularization [9] which seeks to minimize a penalized least squares problem of the form: min_𝑓𝑔−𝐻𝑓²₂+𝛼𝑅𝑓²₂, (2) where the first term is the data fidelity of the solution 𝑓, and the regularization term ‖𝑅𝑓‖²₂ restricts smoothness of the solution. The positive regularization parameter 𝛼 plays the role of balancing the tradeoff between the fidelity and noise sensitivity. The regularization operator 𝑅 is a carefully chosen matrix, often the identity matrix or a discrete approximation of the first or second order derivative operator.

In application to image processing, however, the Tikhonov regularization can produce poor solutions (with overly smoothed edges) when the desired solution comes from an image with edges; that is, it overly penalizes discontinuities in the solutions [10]. In this regard, Rudin et al. [11] proposed a total variation (TV) regularization which has the ability to preserve edges well and remove noise at the same time. The resulting model (commonly referred to as the Rudin-Osher-Fatemi (ROF) model) has been proven to be successful in a wide range of applications in image processing [12]. We should note that there are many other edge-preserving restoration techniques in the literature, such as the anisotropic diffusion methods of Perona and Malik for denoising [13], morphological wavelets, and dyadic tree-based edge-preserving method proposed by Xiang and Ramadge [14]. In this paper, we focus on the the minimization problem with TV regularization min_𝑓𝑔−𝐻𝑓²₂+𝛼𝑓_TV, (3) where ‖⋅‖_TV denotes the discrete TV norm. To define the discrete TV norm, we first introduce the discrete gradient ∇𝑓 [15]: ∇𝑓_𝑖,𝑗=∇𝑓^𝑥_𝑖,𝑗,∇𝑓^𝑦_𝑖,𝑗 (4) with ∇𝑓^𝑥_𝑖,𝑗=𝑓_𝑖+1,𝑗−𝑓_𝑖,𝑗,if𝑖<𝑛,0,if𝑖=𝑛,∇𝑓^𝑦_𝑖,𝑗=𝑓_𝑖,𝑗+1−𝑓_𝑖,𝑗,if𝑗<𝑛,0,if𝑗=𝑛, (5) for 𝑖,𝑗=1,…,𝑛, and 𝑓_𝑖,𝑗 represents the value of pixel (𝑖,𝑗) in the image (the (𝑗−1)𝑛+𝑖th entry of the vector 𝑓). Then the discrete TV norm of 𝑓 is defined as follows: 𝑓_TV∶=_{1≤𝑖,𝑗≤𝑛}∇𝑓_𝑖,𝑗=_{1≤𝑖,𝑗≤𝑛}∇𝑓^𝑥_𝑖,𝑗²+∇𝑓^𝑦_𝑖,𝑗², (6) where |𝑦|=𝑦²₁+𝑦²₂ for every 𝑦=(𝑦₁,𝑦₂)∈ℝ².

In recent years, a great many of algorithms have been developed for total variation based image restoration and proved to be effective for reducing blur and noise while preserving edges. In the original TV regularization paper [11], the authors proposed a time marching scheme to solve the associated Euler-Lagrange equation of (3). The drawback of this method is very slow due to stability constraints. Later, Vogel and Oman [16] proposed a lagged diffusivity fixed point method to solve the same Euler-Lagrange equation of (3). They proved that this method had a global convergent property and was asymptotically faster than the explicit time marching scheme [17]. In [18], Chan et al. applied the Newton’s method to solve the nonlinear primal-dual system of the system (3). Chambolle [15] considered a dual formulation of the TV denoising problem and proposed a semi-implicit gradient descent algorithm to solve the resulting constrained optimization problem. This method is globally convergent with a suitable step size. Recently, Wang et al. [19] proposed a fast total variation deconvolution method which uses splitting technique and constructs an iterative procedure of alternately solving a pair of easy subproblems associated with an increasing sequence of penalty parameter values. In [20], Goldstein and Osher proposed the novel split Bregman iterative algorithm to deal with the artificial constraints; their method has several advantages such as fast convergence rate and stability.

More recently, Huang et al. [2] proposed a minimization problem of the form min_𝑓,𝑢𝒥𝑓,𝑢≡min_𝑓,𝑢𝐻𝑓−𝑔²₂+𝛼₁𝑓−𝑢²₂+𝛼₂‖𝑢‖_TV, (7) where 𝛼₁ and 𝛼₂ are positive regularization parameters. The authors employed an alternating minimization algorithm to solve the system (7). The numerical results on image restoration show the efficiency of their method. The idea of this method is similar to the one proposed in [19]. Both of them use the penalty method by introducing an auxiliary variable.

In [2], the minimization problem (7) can be solved by two steps: one is the deblurring step which is employed by the Tikhonov regularization and the other step is denoising. Although the noise can be removed to a certain extent in the second step, we may lose some details before the denoising step because of the Tikhonov regularization used in the first step (deblurring), as we know that Tikhonov regularization penalizes edges.

In [21], Chavent and Kunisch considered a total bounded variation regularization minimization problem given by min_𝑓𝐻𝑓−𝑔²₂+𝛼𝑓²₂+𝛽𝑓_TV, (8) where 𝛼,𝛽 are both positive parameters. The authors proved that the solution of system (8) is unique. In [22], Hinrermüller and Kunisch applied the semismooth Newton method to solve the Fenchel predual of (8). Numerical results for image denoising and zooming/resizing showed the efficiency of their approach. In [23], Liu and Huang introduced an extended split Bregman iteration to solve the minimization problem (8). Numerical simulations illustrated the excellent reconstruction performance of their method.

Note that the unconstrained problem (3) is equivalent to the following constrained minimization problem: min_𝑓,𝑢𝒥𝑓,𝑢≡min_𝑓𝐻𝑓−𝑔²₂+𝛼₂𝑓_TV+𝛼₃‖𝑢‖_TVs.t.𝑓=𝑢, (9) where 𝛼₂,𝛼₃>0. Then using penalty method, we obtain proposed minimization problem as follows: min_𝑓,𝑢𝒥𝑓,𝑢≡min_𝑓𝐻𝑓−𝑔²₂+𝛼₁𝑓−𝑢²₂+𝛼₂𝑓_TV+𝛼₃‖𝑢‖_TV, (10) where 𝛼₁, 𝛼₂, and 𝛼₃ are positive regularization parameters. We note that the problem (9) is the same with the problem (7) if we set the parameter 𝛼₂=0. The minimization problem (10) can be rewritten as follows: min_𝑓,𝑢𝒥𝑓,𝑢≡min_𝑢min_𝑓𝐻𝑓−𝑔²₂+𝛼₁𝑓−𝑢²₂+𝛼₂𝑓_TV+𝛼₃‖𝑢‖_TV. (11) The method for solving the minimization problem (11) will be discussed in Section 2. Our numerical results will show that the proposed method yields state of the art results both in terms of SSIM and PSNR.

This paper is outlined as follows. In the next section, we first give a brief introduction of the split Bregman method, and then we propose an iterative algorithm for solving (10). The convergence property of the proposed method is given in Section 3. In Section 4, we present numerical experiments to show the efficiency of the proposed method. Finally, the concluding remarks can be found in Section 5.

2. Alternating Minimization Iterative Scheme

In this section, we derive an algorithm to solve the minimization problem (11). Before we discuss the alternating iterative algorithm for solving (11), we would like to give a brief introduction of the split Bregman iteration [20].

2.1. Split Bregman Iteration

2.1.1. Bregman Iteration

Bregman iteration is a concept that originated in functional analysis for finding minimizer of convex functionals [24]. Osher et al. [25] first applied the Bregman iteration to the ROF model for denoising problem in image processing. The basic idea of the Bregman iteration is to transform a constrained optimization problem to an unconstrained problem. The objective functional in the transformed unconstrained problem is defined by means of the Bregman distance of a convex functional. Suppose the unconstrained problem is formulated as min_𝑓𝐽𝑓+𝜆2𝐹𝑓,𝑔, (12) where 𝐽(𝑓) is a convex function and 𝐹(𝑓,𝑔) is convex and differentiable. The Bregman distance of a convex function 𝐽(𝑓) at the point 𝑣 is defined as the following (nonnegative) quantity: 𝐷^𝑝_𝐽𝑓,𝑣≡𝐽𝑓−𝐽(𝑣)−⟨𝑝,𝑓−𝑣⟩, (13) where 𝑝∈𝜕𝐽; that is, 𝑝 is one of the subgradients of 𝐽 at 𝑣. Then the authors in [25] employed the Bregman iteration method to solve the unconstrained problem (11). Assume 𝐹(𝑓,𝑔)=‖𝐻𝑓−𝑔‖², and then the Bregman iteration method is to alternatively iterate the following scheme: 𝑓^𝑘+1=argmin_𝑓𝐷^𝑝𝑘_𝐽𝑓,𝑓^𝑘+𝜆2𝐻𝑓−𝑔²₂=argmin_𝑓𝐽𝑓−𝑝^𝑘,𝑓−𝑓^𝑘+𝜆2𝐻𝑓−𝑔²₂,𝑝^𝑘+1=𝑝^𝑘−𝜆𝐻^𝑇𝐻𝑓^𝑘+1−𝑔. (14) As shown in [25, 26], when 𝐻 is linear, the iteration (13) can be reformulated into the simplified method 𝑓^𝑘+1=argmin_𝑢𝐽𝑓+𝜆2𝐻𝑓−𝑏^𝑘2₂,𝑏^𝑘+1=𝑏^𝑘+𝑔−𝐻𝑓^𝑘+1. (15)(16) This Bregman iteration technique has mainly two advantages over tradition penalty function/continuation methods. One is that it converges very quickly when applied to certain types of objective functions, especially for problems where 𝐽(𝑓) contains an 𝑙₁-regularization term. The other advantage is that the value of 𝜆 in (14) remains constant. See [20, 25, 26] for further details.

2.1.2. Split Bregman Iteration

In [20], the authors considered the problem min_𝑢,𝑑‖𝑑‖₁+𝐸𝑓suchthat𝑑=𝜙𝑓, (17) where 𝐸(𝑓) and 𝜙(𝑓) are convex and differentiable. To solve the problem (16), they convert the constrained problem into an unconstrained problem: min_𝑓,𝑑‖𝑑‖₁𝑓+𝜆+𝐸2𝑑−𝜙(𝑓)²₂. (18)

Let 𝐽(𝑓,𝑑)=‖𝑑‖₁+𝐸(𝑓) and 𝐹(𝑓,𝑑)=‖𝑑−𝜙(𝑓)‖²₂, and the aforementioned Bregman iteration (14) can be similarly applied to (17). Then, they obtained the following elegant two-phase iterative algorithm (split Bregman iteration scheme).

Split Bregman iteration: 𝑓^𝑘+1,𝑑^𝑘+1=min_𝑓,𝑑‖𝑑‖₁𝑓+𝜆+𝐸2𝑑−𝜙(𝑓)−𝑏^𝑘2₂,𝑏^𝑘+1=𝑏^𝑘+𝜙𝑓^𝑘+1−𝑑^𝑘+1. (19)(20) The split Bregman iteration has stable convergence property, and it is extremely fast and very simple to program. For more details on split Bregman and its applications, see [20, 23, 27, 28].

2.2. Proposed Alternating Iteration Scheme

In this section, we propose an alternating minimization algorithm to solve the problem (10). Given an initial 𝑢⁰, we get 𝒫_ℎ𝑢^𝑘−1∶=𝑓^𝑘=argmin_𝑓𝐻𝑓−𝑔²₂+𝛼₁𝑓−𝑢^𝑘−12₂+𝛼₂𝑓_TV,𝒫_𝑡𝑣𝑓^𝑘∶=𝑢^𝑘=argmin_𝑢𝛼₁𝑢−𝑓^𝑘2₂+𝛼₃‖𝑢‖_TV, (21) for 𝑖=1,2,…. As a matter of convenience, we express the relationship between 𝑢^𝑘 and 𝑢^𝑘−1 as below: 𝑢^𝑖=𝒫_tv𝒫_ℎ𝑢^𝑖−1, (22) we denote 𝑢^𝑘=𝒯(𝑢^𝑘−1) for simplicity, where 𝒯(⋅)=𝒫_tv(𝒫_ℎ(⋅)).

Considering 𝑓, the minimization problem (10) is reduced to a minimization problem with respect to 𝑓min_𝑓𝐻𝑓−𝑔²₂+𝛼₁𝑓−𝑢²₂+𝛼₂𝑓_TV. (23) Although there are many methods to solve (22), we focus on the split Bregman iteration here. Let 𝐸(𝑓)=(1/𝛼₂)‖𝐻𝑓−𝑔‖²₂+(𝛼₁/𝛼₂)‖𝑓−𝑢‖²,𝐽(𝑓,𝑑)=‖𝑑‖₁+𝐸(𝑓), and 𝐹(𝑓,𝑑)=‖𝑑−∇𝑓‖². According to (18) and (19), we have 𝑓^𝑖+1,𝑑^𝑖+1=argmin_𝑓,𝑑‖𝑑‖₁+1𝛼₂𝐻𝑓−𝑔²₂+𝛼₁𝛼₂𝑓−𝑢²₂+𝜆2𝑑−∇𝑓−𝑏^𝑖2₂,𝑏^𝑖+1=𝑏^𝑖+∇𝑓^𝑖+1−𝑑^𝑖+1. (24) Clearly, the minimization with respect to 𝑓^𝑖+1 and 𝑑^𝑖+1 in (23) is decoupled, and thus they can be solved separately. We perform the subminimization problem (23) efficiently by iteratively minimizing the following subproblems with respect to 𝑓 and 𝑑 separately: 𝑓^𝑖+1=argmin_𝑓𝐻𝑓−𝑔²₂+𝛼₁𝑓−𝑢²₂+𝛼₂𝜆2𝑑^𝑖−∇𝑓−𝑏^𝑖2₂,𝑑^𝑖+1=argmin_𝑑‖𝑑‖₁+𝜆2𝑑−∇𝑓^𝑖+1−𝑏^𝑖2₂. (25)(26) The minimizer 𝑓^𝑖+1 is given by normal equations: 𝐻^𝑇𝐻+𝛼₂𝜆2∇^𝑇∇+𝛼₁𝐼𝑓^𝑖+1=𝐻^𝑇𝑔+𝛼₁𝑢+𝛼₂𝜆2∇^𝑇𝑑^𝑖−𝑏^𝑖, (27) where 𝐼 is the identity matrix and ∇^𝑇=−div represents the adjoint of ∇. When an appropriate boundary condition is given, the normal equation (26) can be solved by fast algorithms. In this paper, we impose the periodic boundary condition, and then the matrix 𝐻^𝑇𝐻,∇^𝑇∇ are all block circulant [4, 8], which can be diagonalized by the two-dimensional discrete Fourier transform ℱ [19]. By applying the convolution theorem, we obtain𝑓^𝑖+1=ℱ⁻¹2ℱ(𝐻)^∗∘ℱ𝑔+2𝛼₁ℱ(𝑢)+𝛼₂𝜆ℱ∇_𝑥^∗∘ℱ𝑣^𝑖_𝑥+𝛼₂𝜆ℱ∇_𝑦^∗∘ℱ𝑣^𝑖_𝑦2ℱ(𝐻)^∗∘ℱ(𝐻)+2𝛼₁𝐼+𝛼₂𝜆ℱ∇_𝑥^∗∘ℱ∇_𝑥+𝛼₂𝜆ℱ∇_𝑦^∗∘ℱ∇_𝑦, (28)where “*” denotes complex conjugacy, “∘” denotes component-wise multiplication, 𝑣^𝑖=𝑑^𝑖−𝑏^𝑖=(𝑑^𝑖_𝑥−𝑏^𝑖_𝑥,𝑑^𝑖_𝑦−𝑏^𝑖_𝑦)=(𝑣^𝑖_𝑥,𝑣^𝑖_𝑦), ∇=(∇_𝑥,∇_𝑦), and the division is component-wise as well.

The minimizer 𝑑^𝑖+1 of (25) can be determined by the following shrinkage formula: 𝑑^𝑘+1=max∇𝑓^𝑖+1+𝑏^𝑖∇𝑓,0^𝑖+1+𝑏^𝑖∇𝑓^𝑖+1+𝑏^𝑖,0withtheconvention0⋅0=0. (29)

To sum up, we get Algorithm 1 for solving the subproblem (22).

(1) initialization: 𝜆,𝛼1,𝛼2,𝑢,𝑓0=0,𝑑0=𝑏0=0,𝑖=0; (2) iteration:    compute 𝑓𝑖+1 by using formula (27)    𝑑𝑖+1=max∇𝑓𝑖+1+𝑏𝑖∇𝑓,0𝑖+1+𝑏𝑖∇𝑓𝑖+1+𝑏𝑖, with the convention 00⋅0=0    𝑏𝑖+1=𝑏𝑖+∇𝑓𝑖+1−𝑑𝑖+1    stop or set 𝑖=𝑖+1

Algorithm 1 

The split Bregman scheme for subproblem (22).

Considering 𝑢, the outer minimization problem can be interpreted as the TV minimization scheme to denoise the recovered 𝑓 generated by the previous step. The minimum problem is as follows: 𝒫_tv𝑓=argmin_𝑢𝛼₁𝑢−𝑓²+𝛼₃‖𝑢‖_TV. (30) There are several efficient methods for solving this problem, such as the primal-dual method [18, 29], the lagged diffusivity fixed point iteration proposed by Vogel and Oman [16], the semismooth Newton method [30], and Chambolle’s dual algorithm [15]. For the simplicity of Chambolle’s dual algorithm, we adopt it here for TV denoising problem (29). The idea given by Chambolle is to replace the optimization of the image 𝑢 by the optimization of a vector field 𝑝 that is related to 𝑢 by 𝑢=𝑓−(𝛼₃/2𝛼₁)div𝑝. For a noisy image 𝑓, the vector field is the one that minimizes min_𝑝𝛼𝑓−³2𝛼₁div𝑝²,s.t.𝑝_𝑖,𝑗≤1,∀𝑖,𝑗=1,…,𝑛, (31) where 𝑝_𝑖,𝑗=𝑝^𝑥_𝑖,𝑗𝑝^𝑦_𝑖,𝑗 (32) is the dual variable at the (𝑖,𝑗)th pixel location, 𝑝 is the concatenation of all 𝑝_𝑖,𝑗, and the discrete divergence of 𝑝 is given by div𝑝_𝑖,𝑗≡𝑝^𝑥_𝑖,𝑗−𝑝^𝑥_{𝑖−1,𝑗}+𝑝^𝑦_𝑖,𝑗−𝑝^𝑦_{𝑖,𝑗−1} (33) with 𝑝^𝑥_0,𝑗=𝑝^𝑦_𝑖,0=0. The vector div𝑝 is the concatenation of all (div𝑝)_𝑖,𝑗. For simplicity, we denote 𝛽=𝛼₃/2𝛼₁. The iterative scheme proposed by Chambolle for computing the optimal solution 𝑝 is as follows: 𝑝^𝑙+1_𝑖,𝑗=𝑝^𝑙_𝑖,𝑗+𝜏∇(div𝑝^𝑙−𝑓/𝛽)_𝑖,𝑗1+𝜏∇(div𝑝^𝑙−𝑓/𝛽)_𝑖,𝑗,∀1≤𝑖,𝑗≤𝑛, (34) where 𝜏 is the step size and 𝑝^𝑙_𝑖,𝑗 is the 𝑙th iterate of the iterative method for minimizer; see [15] for more details. After the minimizer 𝑝^∗ of the constrained optimization problem in (31) is determined, the denoised image can be computed by 𝑢^∗=𝑓−𝛽div𝑝^∗. (35)

In summary, we obtain Algorithm 2 by using alternating minimization scheme to solve the minimization problem (10).

(1) initialization: 𝑢0,𝑝0𝑖,𝑗𝐰0,𝐮0,𝛼1,𝛼2,𝜆,𝑘=0; (2) iteration    compute 𝑓𝑘 using Algorithm 1 for fixed 𝑢𝑘−1    compute 𝑢𝑘 according to Chambolle method (31) and    (34) for fixed 𝑓𝑘    stop or set 𝑘=𝑘+1

Algorithm 2 

Alternating iteration scheme for solving (10).

3. Convergence Analysis

In this section, we make use of a theorem proposed in [31] to give the convergence property of the proposed algorithm. The theorem is given by the following.

Theorem 1 (see [31]). Let 𝒯∶ℝ^𝑛2→ℝ^𝑛2 be a 𝛽-averaged nonexpansive operator and the set of fixed points of 𝒯 be nonempty. Then for any 𝑥⁰, the sequence 𝑥^𝑘=𝑇^𝑘𝑥⁰ converges weakly to a fixed point in ℝ^𝑛2.

Definition 2. An operator 𝒯 is called nonexpansive if, for all 𝑥 and 𝑦∈ℝ^𝑛2, 𝒯𝑥−𝒯𝑦≤𝑥−𝑦. (36) Given a nonexpansive operator 𝒩, let 𝒯=(1−𝛼)ℐ+𝛼𝑁; for some 𝛼∈(0,1), the operator 𝒯 is said to be 𝛼-average.

Definition 3. An operator 𝒯 is called 𝜈-inverse strongly monotone (ism) if there is 𝜈>0, such that 𝒯𝑥−𝒯𝑦,𝑥−𝑦≥𝜈𝒯𝑥−𝒯𝑦². (37) Let ℰ=ℐ−𝒯 be complement of the operator 𝒯, and then we can easily get the following identity: 𝑥−𝑦²−𝒯𝑥−𝒯𝑦²=2⟨ℰ𝑥−ℰ𝑦,𝑥−𝑦⟩−ℰ𝑥−ℰ𝑦². (38) An operator ℱ is called firmly nonexpansive if it is 1-ism.

Lemma 4 (see [31]). An operator 𝒩 is nonexpansive if and only if its complement ℰ=ℐ−𝒩 is 1/2-ism. If ℰ is 𝜈-ism and 𝛾>0, then the operator 𝛾ℰ is (𝜈/𝛾)-ism.

Lemma 5. An operator 𝒫 is 𝛽-average nonexpansive if and only if its complement ℰ=ℐ−𝒫 is (1/2𝛽)-ism.

Proof. Firstly, suppose 𝒫 is 𝛽-average; from Definition 2, there exists a nonexpansive 𝒩 such that 𝒫=(1−𝛽)ℐ+𝛽𝒩, and then ℰ=ℐ−𝒫=𝛽(ℐ−𝒩). Since 𝒩 is nonexpansive, from Lemma 4, we have that ℐ−𝒩 is 1/2-ism and ℰ=𝛽(ℐ−𝒩) is (1/2𝛽)-ism.
Now we assume that ℰ is (1/2𝛽)-ism, we write 𝒫=(1−𝛽)𝐼+𝛽𝒩, for 𝒩=ℐ−(1/𝛽)ℰ, from Lemma 4, we have that ℐ−𝒩 is (1/2)-ism and 𝒩 is nonexpansive, and then from Definition 2, the operator 𝒫 is 𝛽-average.

Lemma 6 (see [32]). Let be convex and semicontinuous and 𝛼>0. Suppose 𝑥 is defined as follows: 𝑥=argmin_𝑥𝑦−𝑥²+𝛼𝜑(𝑥). (39) Define 𝒮 such that 𝑥=𝒮(𝑦) for every 𝑦. Then 𝒮 and ℐ−𝒮 are firmly nonexpansive.

Theorem 7. Let 𝛼₁ and 𝛼₂ be positive numbers. Suppose 𝑓 is defined as follows: 𝑓=argmin_𝑓𝐻𝑓−𝑔²+𝛼₁𝑓−𝑢²+𝛼₂𝑓_𝑇𝑉. (40) Define 𝒫_𝑛 such that 𝑓=𝒫_𝑛(𝑢) for every 𝑢. Then 𝒫_𝑛 is 1/2-averaged nonexpansive.

Proof. Let 𝜙(𝑓)=‖𝑓‖_TV, for every 𝑢, the minimum in (22) is achieved at a unique point 𝒫_𝑛 which is characterized by the inclusion 𝐻^𝑇𝑔+𝛼₁𝑢−𝐻^𝑇𝐻+𝛼₁𝐼𝒫_𝑛(𝑢)∈𝛼₂2𝜕𝜙𝑓 (41) from the property of subdifferential of 𝜙(𝑓), ∀𝑢,𝑤∈ℝ^𝑛2, and we have the following inequalities: 2𝛼₂𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢),𝐻^𝑇𝑔+𝛼₁𝑢−𝐻^𝑇𝐻+𝛼₁𝐼𝒫_𝑛(𝑢𝒫)+𝜙_𝑛(𝑢)≤𝜙𝒫_𝑛(𝑤),2𝛼₂𝒫_𝑛(𝑢)−𝒫_𝑛(𝑤),𝐻^𝑇𝑔+𝛼₁𝑤−𝐻^𝑇𝐻+𝛼₁𝐼𝒫_𝑛(𝑤𝒫)+𝜙_𝑛(𝑤)≤𝜙𝒫_𝑛(𝑢). (42)(43) Adding these two inequalities, we obtain 𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢)^𝑇1𝛼₁𝐻^𝑇𝐻+𝐼𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢≤𝒫_𝑛)(𝑤)−𝒫_𝑛(𝑢)^𝑇(𝑤−𝑢). (44) It is obvious that we have the following inequality: 𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢)²≤𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢)^𝑇×1𝛼₁𝐻^𝑇𝐻+𝐼𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢). (45) From (43) and (44), we obtain 𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢)²≤𝒫_𝑛(𝑤)−𝒫_𝑛(𝑢)^𝑇(𝑤−𝑢). (46) Then from Definition 3, the operator 𝒫_𝑛 is firmly nonexpansive. Analogously we can easily obtain that the operator ℐ−𝒫_𝑛 is also firmly nonexpansive. Therefore, it follows from Lemma 5 that the operator 𝒫_𝑛 and ℐ−𝒫_𝑛 is 1/2-averaged nonexpansive.

Corollary 8. The operator 𝒯=𝒫_𝑡𝑣(𝒫_ℎ) is 3/4-average nonexpansive.

Proof. From Lemmas 5 and 6, and Theorem 7, we know that 𝒫_tv and 𝒫_ℎ are both 1/2-averaged nonexpansive operators, and there exist nonexpansive operators 𝒩_ℎ and 𝒩_tv such that 𝒫_ℎ=12ℐ+𝒩_ℎ,𝒫_tv=12ℐ+𝒩_tv. (47) Thus we have 𝒯=𝒫_tv𝒫_ℎ=14ℐ+𝒩_ℎ+𝒩_tv+𝒩_tv𝒩_ℎ. (48) Set 𝒩=(1/3)(𝒩_ℎ+𝒩_tv+𝒩_tv𝒩_ℎ), then for any 𝑥 and 𝑦, =1𝒩𝑥−𝒩𝑦3𝒩_tv𝑥−𝑦+𝒩_ℎ𝑥−𝑦+𝒩_tv𝒩_ℎ≤1𝑥−𝑦3𝒩_tv𝑥−𝒩_tv𝑦+𝒩_ℎ𝑥−𝒩_ℎ𝑦+𝒩_tv𝒩_ℎ(𝑥)−𝒩_tv𝒩_ℎ𝑦≤132𝑥−𝑦+𝒩_ℎ(𝑥)−𝒩_ℎ𝑦≤133𝑥−𝑦=𝑥−𝑦. (49) Consequently, 𝒩 is nonexpansive, and we rewrite 𝒫_ℎ as 3𝒯=1−43ℐ+4𝒩. (50) It follows from Definition 2 that the operator 𝒯=𝒫_𝑡𝑣(𝒫_ℎ) is 3/4-averaged.

According to Theorem 1 and Corollary 8, we conclude that for any initial guess 𝑢⁰∈ℝ^𝑛2, {𝑢^𝑘} generated by (20) converges to the minimizer of 𝒥 in (10).

4. Numerical Experiments

In this section, we present several numerical experiments to illustrate the behavior of the proposed method for image restoration problems. The quality of the restoration results by different methods is compared quantitatively by using the PSNR and SSIM. PSNR is an engineering term for the ratio between the maximum possible power of a signal and the power of corrupting noise that affects the fidelity of its representation. The higher PSNR value, the higher image quality. The SSIM is a well-known quality metric used to measure the similarity between two images. The method is developed by Wang et al. [33] and is based on three specific statistical measures that are much closer to how the human eye perceives differences between images. The higher SSIM value, the better restoration. We also use blur signal-to-noise ratio (BSNR) to describe how much noise is added in the blurry image. Suppose 𝑓,𝑔,𝑓, and 𝑛 are the original image, the blurred and noisy image, the restored image, and the noise, respectively. The PSNR, SSIM, and BSNR are defined as follows [2, 34]: 𝑛PSNR=10log²Max²_𝐼𝑓−𝑓², (51) where Max_𝐼 is the maximum possible pixel value of the image (i.e., when the pixels are represented by using 8 bits per sample, this is 255): SSIM=𝑙𝑓,𝑓𝑐𝑓,𝑓𝑠𝑓,𝑓, (52) where 𝑙𝑓,𝑓=2𝜇_𝑓𝜇𝑓+𝐶₁𝜇²_𝑓+𝜇²𝑓+𝐶₁,𝑐𝑓,𝑓=2𝜎_𝑓𝜎𝑓+𝐶₂𝜎²_𝑓+𝜎²𝑓+𝐶₂,𝑠𝑓,𝑓=2𝜎_𝑓𝑓+𝐶₃𝜎_𝑓𝜎𝑓+𝐶₃. (53) The first term in (52) is the luminance comparison function which measures the closeness of the two images’ mean luminance (𝜇_𝑓 and 𝜇𝑓). The unique maximum of this factor equals 1 if and only if 𝜇_𝑓𝑓=𝜇. The second term 𝑐(𝑓,𝑓) is the contrast comparison function which measures the closeness of the contrast of the two images. Here the contrast is measured by the standard deviation 𝜎_𝑓 and 𝜎𝑓. This term achieves maximum value 1 if and only if 𝜎_𝑓𝑓=𝜎. The third term is the structure comparison function which measures the correlation coefficient between the two images 𝑓 and 𝑓. Note that 𝜎_𝑓𝑓 is the covariance between 𝑓 and 𝑓. The positive values of the SSIM index are in [0,1]. A value of 0 means no correlation between images, and 1 means that 𝑓=𝑓. The positive constants 𝐶₁,𝐶₂, and 𝐶₃ can be thought of as stabilizing constants for near-zero denominator values. In the following experiments, we will also use SSIM index map to reveals areas of high/low similarity between two images, the whiter SSIM index map, the closer between the two images. We refer the reader to see [33, 34] for further details on SSIM and SSIM index map.

The BSNR is given by BSNR=20log₁₀𝑔‖𝑛‖. (54)

In the following experiments, we compare our proposed method (we call the proposed method FNDTV later) with FastTV [2]. For FastTV, based on the suggestions in [2], we fixed its parameters 𝛼₁=0.003 for BSNR = 40 dB and 𝛼₁=0.006 for BSNR = 30 dB, and we determine the best value of 𝛼₂ such as the restored images with best performance. For our method, we also determine the best value of the regularization parameters to give the best performance.

The stopping criterion of the proposed method is that the relative difference between the successive iteration of the restored image should satisfy the following inequality: 𝑓^𝑘+1−𝑓^𝑘𝑓^𝑘+1≤tol, (55) where 𝑓^𝑘 is the restored image at the 𝑘th iteration of the proposed method. We set tol=1×10⁻⁴ in all tests for both methods.

Four test images, “Cameraman,” “Bridge,” “Lena,” and “Resolution Chart,” which are commonly used in the literature, are shown in Figure 1. We test several kinds of blurring kernels including average, motion, gaussian, and out-of-focus. These different blurring kernels can be generated by the Matlab built-in function 𝑓𝑠𝑝𝑒𝑐𝑖𝑎𝑙. In all tests, we add the Gaussian white noise of different BSNR to the blurred images.

(a)

(b)

(c)

(d)

Figure 1 

Original images. (a) “Cameraman,” (b) “Resolution Chart,” (c) “Bridge,” (d) “Lena.”

All experiments are carried out on Windows XP 32-bit and Matlab v7.10 running on a desktop equipped with an Intel Core2 Duo CPU 2.93 GHz and 2 GB of RAM.

4.1. Average Blur Example

In this example, we consider the well-known “Cameraman” image (256×256) which is shown in Figure 1(a). The image is blurred by a 3×3 box average kernel and contaminated by BSNR = 40 dB Gaussian noise. The blurred and noisy image is shown in Figure 2(a).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

Results of different methods when restoring blurred and noisy image “Cameraman” degraded by average 9×9 uniform blur and a noise with BSNR = 40 dB: (a) blurred and noisy image “Cameraman”; (b) restored image by FastTV; (c) restored image by FNDTV; (d) SSIM index map of the corrupted image; (e) SSIM index map of the recovered image by FastTV; (f) SSIM index map of the recovered image by FNDTV.

Figures 2(b)–2(d) show the restored images by FNDTV and FastTV. We can see that the visual quality of reconstruction by FNDTV is slightly better than the outcome of FastTV. From Table 1, it is not difficult to see that both of PSNR and SSIM of the recovered image by FNDTV are higher than the one obtained by FastTV.

Table 1 

Numerical results for the experiments in terms of PSNR (dB) and SSIM.

We also show the SSIM index maps of the restored images recovered by the two methods in Figures 2(e)–2(f), and the map can deliver more information about the quality degradation of the restored images. In contrast, the SSIM map of the restored image by the proposed method is slightly whiter than the SSIM map by FastTV.

4.2. Motion Blur Example

Motion blur is considered in this example. The observed image is the so-called “Resolution Chart” [11], it is degraded by a motion bur with length 7 and contaminated by a white Gaussian noise with BSNR = 30 dB. The degraded image is shown in Figure 3(a). The restored images by FastTV and FNDTV are presented in Figures 3(b) and 3(c). We also report the PSNR and SSIM values by these methods in Table 1. We see that both the PSNR and SSIM values of the restored image by the FNDTV method are higher than FastTV. In addition, the SSIM map obtained by the proposed method is slightly whiter than the map by FastTV.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Results of different methods when restoring blurred and noisy image “Resolution Chart” degraded by motion blur with length 7 and a noise with BSNR = 30 dB: (a) blurred and noisy image “Resolution Chart”; (b) restored image by FastTV; (c) restored image by FNDTV; (d) SSIM index map of the corrupted image; (e) SSIM index map of the recovered image by FastTV; (f) SSIM index map of the recovered image by FNDTV.

4.3. Gaussian Blur Example

The original image “Bridge” is shown in Figure 1(c). The blurred and noisy image is degraded by a Gaussian blur with radius 3 and standard deviation 𝜎=3 and then contaminated by Gaussian noise with BSNR = 40 dB. Figure 4(a) shows the blurred and noisy observation. The restored images obtained by FastTV and FNDTV are shown in Figures 4(b) and 4(c). The numerical results of two different methods in terms of PSNR and SSIM are given in Table 1.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Results of different methods when restoring blurred and noisy image “Bridge” degraded by Gaussian blur with radius 3 and standard deviation 𝜎 = 3 and a noise with BSNR = 40 dB: (a) blurred and noisy image “Bridge”; (b) restored image by FastTV; (c) restored image by FNDTV; (d) SSIM index map of the corrupted image; (e) SSIM index map of the recovered image by FastTV; (f) SSIM index map of the recovered image by FNDTV.

From Table 1, it is not difficult to see that the PSNR and SSIM of the restored image by FNDTV are higher than the one obtained by FastTV.

4.4. Out-of-Focus Blur Example

This example consists in restoring the image “Lena” degraded by an out-of-focus blur with radius 3 and contaminated by BSNR = 30 white Gaussian noise. The image “Lena” is a good test image because it has a nice mixture of detail, flat regions, shading area, and texture and has been widely used in the literature to test image restoration algorithms [19]. Figure 5(a) shows the blurred and noisy image. The restored results by both methods are shown in Figures 5(b)–5(d), and Table 1 lists the PSNR and SSIM values. From the table, we observe that both the PSNR and SSIM values of the restored image by the proposed method are better than what obtained by FastTV.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

Results of different methods when restoring blurred and noisy image “Lena” degraded by out-of-focus blur with radius 3 and a noise with BSNR = 30 dB: (a) blurred and noisy image “Lena”; (b) restored image by FastTV; (c) restored image by FNDTV; (d) SSIM index map of the corrupted image; (e) SSIM index map of the recovered image by FastTV; (f) SSIM index map of the recovered image by FNDTV.

5. Conclusion

In this paper, we have presented a new efficient algorithm for image restoration based on total variation regularization. We give a convergence proof for the algorithm, and the numerical results show that the proposed method is competitive with the state of the art method FastTV. In addition, an important feature is that the proposed method can suppress noise very well while it can preserve details of the restored image. We will consider extending the proposed method for color or other multichannel image restoration in the future.

Conflict of Interests

All of the coauthors do not have a direct financial relation with the trademarks mentioned in our paper that might lead to a conflict of interests for any of the coauthors.

Acknowledgments

The work of Jun Liu, Ting-Zhu Huang, and Si Wang is supported by 973 Program (2013CB329404), NSFC (61170311), Chinese Universities Specialized Research Fund for the Doctoral Program (20110185110020), Sichuan Province Sci., and Tech. Research Project (2012GZX0080). The work of Xiao-Guang Lv is supported by Nature science foundation of Jiangsu Province (BK20131209).