# Image Restoration using Nonlocal Regularized Variational Model with Spatially Adapted Regularization Parameter

Authors: Chen Pan; Helin Feng
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Restoration of blurred images corrupted with impulse noise has attracted increasing attention due to its great potential toward practical applications. To guarantee the noisy image deblurring performance, total variation (TV), one of the most popular regularizers, has been extensively exploited in the literature. This local regularizer is known for its capability of preserving edges and penalizing oscillations but not sharp discontinuities. However, TV regularized variational models commonly suffer from undesirable staircase-like artifacts in homogeneous regions. The reason behind this phenomenon is that TV regularizer favors solutions that are piecewise constant in numerical experiments. To overcome the model-dependent deficiency, we proposed to replace the commonly used local regularizer by its nonlocal extension, i.e., nonlocal total variation (NLTV) regularizer. By taking advantage of the high degree of self-similarity within images, the proposed nonlocal variational model is able to remove both blurring and impulse noise while preserving meaningful image details. To guarantee solution stability and efficiency, the resulting nonsmooth image restoration problem was effectively dealt with using an alternating optimization algorithm. In addition, to further enhance the image quality, a local variance estimator-based automatic method was introduced to calculate the spatially adapted regularization parameters during image restoration. Extensive experiments have demonstrated that the proposed method could achieve superior imaging performance compared with several state-of-the-art image restoration methods.


---

## Body

## 1. Introduction

### 1.1. Background and Related Work

Real-world images are often degraded by both blurring and additive/multiplicative noise during image acquisition and analog-to-digital conversion. Over the past decades, restoration of such images has attracted increasing attention in the fields of image processing, such as astronomical imaging [1], remote sensing [2], biomedical imaging [3], and many others [4].Image restoration is a typical ill-conditioned inverse problem, which requires regularization techniques to guarantee a unique and stable solution. In the current literature, the classical regularization for image restoration is a quadratic regularization term or Tikhonov regularization [5]. However, the restoration results often suffer from oversmoothing of important signal features (e.g., edges and textures in images) thus compromising the image visual quality. To overcome this limitation, total variation (TV) regularizer, first proposed by Rudin et al. [6] for Gaussian noise removal, has become one of the most popular regularization tools due to its ability to preserve edges and penalize oscillations but not sharp discontinuities [7, 8]. In addition, TV regularized variational models have also been extensively used for removal of impulse [9], Poisson [10], Rician [11], and Speckle (Gamma) [12] noises.In the case of images degraded by both blurring and additive Gaussian noise, TV regularized L2 variational models based on theMaximum A Posteriori (MAP) framework have been developed to enhance image quality [13]. In particular, the models were proposed by combining the squared L2-norm data-fidelity term and TV regularizer. However, these models often tend to cause staircase-like artifacts due to the nature of TV in favoring piecewise constant solution. To overcome the model-dependent deficiency, many extensions of TV have been reported in the current literature. Second-order total variation (TV2) regularizer [14, 15] was proposed to replace the TV regularizer to suppress the undesirable staircase-like artifacts. Many efforts have been made to combine both TV and TV2 regularizers to further improve image quality [16, 17]. Motivated by the observation that gradients of natural images follow a hyper-Laplacian distribution [18], the nonconvex hybrid TV regularizer has been developed based on the combination of nonconvex TV and TV2 regularizers [19]. From a statistical point of view, the nonconvex regularizer is well capable of edge-preserving image restoration. The second-order total generalized variation (TGV) regularizer [20–22] has also been proposed to remove the staircase-like artifacts yielded by TV in recovered images. All of these traditional regularizers are developed based on local image operators, which are capable of preserving edges and reducing noise very well but may not perform well in preserving fine texture because texture is nonlocal in nature [23]. Inspired by the concepts of nonlocal means (NLM) [24] and graph Laplacian [25], Gilboa and Osher [26] defined a variational framework based on the nonlocal operators. The corresponding nonlocal total variation (NLTV) regularizer [27, 28] has gained great success in image restoration when images are degraded by both blurring and additive Gaussian noise. In particular, the NLTV regularizer is capable of taking advantage of the high degree of self-similarity within images. It is able to preserve the meaningful image details while removing both blurring and impulse noise in practice.In this paper, we mainly focus on the variational models for restoring blurred images with impulse noise. Thus, the abovementioned variational models with squared L2 data-fidelity term cannot be used directly to deal with the image restoration problems addressed in this paper. To restore blurred images with impulse noise, TV regularized L1 variational (TV-L1) models have attracted considerable attention during recent years [29–31]. In particular, both L1-norm data-fidelity term and TV regularizer were combined in these variational models. Besides, in [32, 33], the authors successfully extended the TV-L1 model to deblurring multichannel images corrupted by impulse noise. On the other hand, to recover the image with pixel values in the same range, a simple but useful box constraint was incorporated into the TV-L1 model to further enhance image quality [34, 35]. Hintermüller and Rincon-Camacho [36] developed a spatially adapted regularization parameter choice rule in the TV-L1 model where small details were desired to be preserved. However, these proposed models still tend to cause staircase-like artifacts since TV favors a piecewise constant solution. More recently, Liu et al. [37] proposed a spatially adapted variational model by replacing TV with a hybrid TV (TV1,2), which was based on the combination of both TV and TV2 regularizers.


### 1.2. Motivation and Contributions

It is well known that nonlocal self-similarity property has been widely exploited in the fields of both low- and high-level computer vision [38–40]. This property is based on the observation that image patches in natural images tend to redundantly repeat themselves many times [41]. Therefore, if images are degraded by blurring and impulse noise, the variational models based on local regularizers (e.g., TV and TV1,2) may not perform well in preserving fine textures in restored images. To guarantee high image quality, it is necessary to incorporate the nonlocal regularizer into the variational image restoration model. Inspired by the previous works in [27, 28, 42, 43], we tend to propose a NLTV regularized L1 variational model by replacing squared L2-norm data-fidelity term with L1-norm data-fidelity term to handle blurred images with impulse noise.On the other hand, the regularization parameter plays an important role in TV regularized variational models. Since images are always comprised of image features of different scales, the constant parameter can lead to significant degradation of the image quality [44]. Thus, to achieve high image quality, the regularization parameter should be selected adaptively according to local image features. As discussed in [36, 44–46], the authors proposed to develop a local variance estimator-based automated method for selection of spatially adapted regularization parameter. Similar to these previous works, we will introduce a spatially adapted regularization parameter selection scheme for our proposed nonlocal variational model. To the best of our knowledge, no research has been conducted on NLTV regularized L1 variational model with spatially adapted regularization parameter selection thus far.Due to the nondifferentiability of NLTV regularizer and nonsmoothness of L1-norm data-fidelity term, it is difficult to solve our proposed nonlocal variational model through the commonly used optimization algorithms. In the current literature, the alternating optimization algorithm [29] has been widely used for image restoration problems and achieved satisfactory results on optimization performance. Motivated by these previous works, to get a numerically stable solution, we will develop an alternating optimization algorithm to solve the proposed image restoration model. In particular, the original nonsmooth optimization problem will be decomposed into simpler subproblems by introducing several intermediate variables. Each of these subproblems has a closed-form solution or could be efficiently dealt with using a simple numerical method.In conclusion, the main ideas and contributions of this paper, given the current state-of-the-art research work, can be summarized as follows:(1)
A spatially adapted nonlocal variational model is proposed to restore blurred images with impulse noise. This model considers the image nonlocal self-similarity property, and the spatially adapted regularization parameter is selected automatically to enhance image quality.(2)
An alternating optimization algorithm is developed to optimize the proposed nonlocal variational image restoration model, in which the optimization of the image restoration model is decomposed into two subproblems. Each subproblem has a closed-form solution or can be efficiently solved using existing optimization algorithms.(3)
The experimental results on several grayscale images have shown the superior performance of our proposed method in terms of both quantitative and visual quality evaluations. It may well provide a new idea for restoration of blurred images with impulse noise.The main benefit of our proposed method is that it takes full advantage of the image nonlocal self-similarity property and spatially adapted regularization parameters. Therefore, the proposed method is capable of restoring degraded images while preserving edges and fine details in practice.


### 1.3. Organization

The remainder of this paper is organized into several sections. In Section2, the image degradation model is introduced, and following this, the MAP-based image restoration framework is presented. The NLTV regularized L1 variational model for restoration of noisy blurred images is specifically developed in Section 3. In Section 4, an alternating optimization algorithm is developed to effectively solve the resulting nonlocal regularized image restoration problem. To further improve image quality, a local variance estimator-based automatic method for selection of spatially adapted regularization parameter is introduced in Section 5. Section 6 presents numerous experiments using several grayscale images, which verify the effectiveness of the proposed method. Finally, we conclude this paper by summarizing our contributions and discussing the future work in Section 7.


## 1.1. Background and Related Work

Real-world images are often degraded by both blurring and additive/multiplicative noise during image acquisition and analog-to-digital conversion. Over the past decades, restoration of such images has attracted increasing attention in the fields of image processing, such as astronomical imaging [1], remote sensing [2], biomedical imaging [3], and many others [4].Image restoration is a typical ill-conditioned inverse problem, which requires regularization techniques to guarantee a unique and stable solution. In the current literature, the classical regularization for image restoration is a quadratic regularization term or Tikhonov regularization [5]. However, the restoration results often suffer from oversmoothing of important signal features (e.g., edges and textures in images) thus compromising the image visual quality. To overcome this limitation, total variation (TV) regularizer, first proposed by Rudin et al. [6] for Gaussian noise removal, has become one of the most popular regularization tools due to its ability to preserve edges and penalize oscillations but not sharp discontinuities [7, 8]. In addition, TV regularized variational models have also been extensively used for removal of impulse [9], Poisson [10], Rician [11], and Speckle (Gamma) [12] noises.In the case of images degraded by both blurring and additive Gaussian noise, TV regularized L2 variational models based on theMaximum A Posteriori (MAP) framework have been developed to enhance image quality [13]. In particular, the models were proposed by combining the squared L2-norm data-fidelity term and TV regularizer. However, these models often tend to cause staircase-like artifacts due to the nature of TV in favoring piecewise constant solution. To overcome the model-dependent deficiency, many extensions of TV have been reported in the current literature. Second-order total variation (TV2) regularizer [14, 15] was proposed to replace the TV regularizer to suppress the undesirable staircase-like artifacts. Many efforts have been made to combine both TV and TV2 regularizers to further improve image quality [16, 17]. Motivated by the observation that gradients of natural images follow a hyper-Laplacian distribution [18], the nonconvex hybrid TV regularizer has been developed based on the combination of nonconvex TV and TV2 regularizers [19]. From a statistical point of view, the nonconvex regularizer is well capable of edge-preserving image restoration. The second-order total generalized variation (TGV) regularizer [20–22] has also been proposed to remove the staircase-like artifacts yielded by TV in recovered images. All of these traditional regularizers are developed based on local image operators, which are capable of preserving edges and reducing noise very well but may not perform well in preserving fine texture because texture is nonlocal in nature [23]. Inspired by the concepts of nonlocal means (NLM) [24] and graph Laplacian [25], Gilboa and Osher [26] defined a variational framework based on the nonlocal operators. The corresponding nonlocal total variation (NLTV) regularizer [27, 28] has gained great success in image restoration when images are degraded by both blurring and additive Gaussian noise. In particular, the NLTV regularizer is capable of taking advantage of the high degree of self-similarity within images. It is able to preserve the meaningful image details while removing both blurring and impulse noise in practice.In this paper, we mainly focus on the variational models for restoring blurred images with impulse noise. Thus, the abovementioned variational models with squared L2 data-fidelity term cannot be used directly to deal with the image restoration problems addressed in this paper. To restore blurred images with impulse noise, TV regularized L1 variational (TV-L1) models have attracted considerable attention during recent years [29–31]. In particular, both L1-norm data-fidelity term and TV regularizer were combined in these variational models. Besides, in [32, 33], the authors successfully extended the TV-L1 model to deblurring multichannel images corrupted by impulse noise. On the other hand, to recover the image with pixel values in the same range, a simple but useful box constraint was incorporated into the TV-L1 model to further enhance image quality [34, 35]. Hintermüller and Rincon-Camacho [36] developed a spatially adapted regularization parameter choice rule in the TV-L1 model where small details were desired to be preserved. However, these proposed models still tend to cause staircase-like artifacts since TV favors a piecewise constant solution. More recently, Liu et al. [37] proposed a spatially adapted variational model by replacing TV with a hybrid TV (TV1,2), which was based on the combination of both TV and TV2 regularizers.


## 1.2. Motivation and Contributions

It is well known that nonlocal self-similarity property has been widely exploited in the fields of both low- and high-level computer vision [38–40]. This property is based on the observation that image patches in natural images tend to redundantly repeat themselves many times [41]. Therefore, if images are degraded by blurring and impulse noise, the variational models based on local regularizers (e.g., TV and TV1,2) may not perform well in preserving fine textures in restored images. To guarantee high image quality, it is necessary to incorporate the nonlocal regularizer into the variational image restoration model. Inspired by the previous works in [27, 28, 42, 43], we tend to propose a NLTV regularized L1 variational model by replacing squared L2-norm data-fidelity term with L1-norm data-fidelity term to handle blurred images with impulse noise.On the other hand, the regularization parameter plays an important role in TV regularized variational models. Since images are always comprised of image features of different scales, the constant parameter can lead to significant degradation of the image quality [44]. Thus, to achieve high image quality, the regularization parameter should be selected adaptively according to local image features. As discussed in [36, 44–46], the authors proposed to develop a local variance estimator-based automated method for selection of spatially adapted regularization parameter. Similar to these previous works, we will introduce a spatially adapted regularization parameter selection scheme for our proposed nonlocal variational model. To the best of our knowledge, no research has been conducted on NLTV regularized L1 variational model with spatially adapted regularization parameter selection thus far.Due to the nondifferentiability of NLTV regularizer and nonsmoothness of L1-norm data-fidelity term, it is difficult to solve our proposed nonlocal variational model through the commonly used optimization algorithms. In the current literature, the alternating optimization algorithm [29] has been widely used for image restoration problems and achieved satisfactory results on optimization performance. Motivated by these previous works, to get a numerically stable solution, we will develop an alternating optimization algorithm to solve the proposed image restoration model. In particular, the original nonsmooth optimization problem will be decomposed into simpler subproblems by introducing several intermediate variables. Each of these subproblems has a closed-form solution or could be efficiently dealt with using a simple numerical method.In conclusion, the main ideas and contributions of this paper, given the current state-of-the-art research work, can be summarized as follows:(1)
A spatially adapted nonlocal variational model is proposed to restore blurred images with impulse noise. This model considers the image nonlocal self-similarity property, and the spatially adapted regularization parameter is selected automatically to enhance image quality.(2)
An alternating optimization algorithm is developed to optimize the proposed nonlocal variational image restoration model, in which the optimization of the image restoration model is decomposed into two subproblems. Each subproblem has a closed-form solution or can be efficiently solved using existing optimization algorithms.(3)
The experimental results on several grayscale images have shown the superior performance of our proposed method in terms of both quantitative and visual quality evaluations. It may well provide a new idea for restoration of blurred images with impulse noise.The main benefit of our proposed method is that it takes full advantage of the image nonlocal self-similarity property and spatially adapted regularization parameters. Therefore, the proposed method is capable of restoring degraded images while preserving edges and fine details in practice.


## 1.3. Organization

The remainder of this paper is organized into several sections. In Section2, the image degradation model is introduced, and following this, the MAP-based image restoration framework is presented. The NLTV regularized L1 variational model for restoration of noisy blurred images is specifically developed in Section 3. In Section 4, an alternating optimization algorithm is developed to effectively solve the resulting nonlocal regularized image restoration problem. To further improve image quality, a local variance estimator-based automatic method for selection of spatially adapted regularization parameter is introduced in Section 5. Section 6 presents numerous experiments using several grayscale images, which verify the effectiveness of the proposed method. Finally, we conclude this paper by summarizing our contributions and discussing the future work in Section 7.


## 2. Problem Formulation

### 2.1. Image Degradation Model

Images are often degraded by blurring and different kinds of noise. In this paper, we will concentrate mainly on the impulse noise, such as salt-and-pepper noise and random-valued impulse noise. We denote byΩ⊂ℝ2 the image domain, which is often considered to be a bounded and piecewise smooth open subset. Let f:Ω⟶ℝ be a real function defined on Ω. The corrupted image g:Ω⟶ℝ is then given by(1)g=Nℋf,where N denotes the image degradation by impulse noise. For the sake of simplicity, we assume that the linear and continuous blurring operator ℋ∈ℒL2Ω is known beforehand. Denote the dynamic range of ℋf to be fmin,fmax, i.e., fmin≤ℋfx≤fmax,∀x∈Ω, the degradation model N is defined as follows [47].(i)
Salt-and-pepper noise: the gray level ofg at pixel location x∈Ω is(2)gx=fmin,with probabilityp1/2,fmax,with probabilityp1/2,ℋfx,with probability1−p1,wherep1 determines the level of the salt-and-pepper noise. The gray values of degraded pixels change to either the minimum value fmin or the maximum value fmax.(ii)
Random-valued impulse noise: the gray level ofg at pixel location x∈Ω is(3)gx=f^x,with probabilityp2,ℋfx,with probability1−p2,where0≤p2≤1 is the level of the random-valued impulse noise. The gray values of f^x are identically and uniformly distributed random numbers in fmin,fmax.Recently, many variational image restoration models have been developed by combining the nonsmooth L1 data-fidelity term with different regularizers. In this paper, we will restrict our attention to the local TV regularizer and its nonlocal extension. The corresponding MAP-based image restoration model will be proposed in the next subsection.


### 2.2. MAP-Based Image Restoration Model

The MAP estimation theory, which inherently includes prior constraints in the form of prior probability density functions, has attracted increasing attention in the field of image processing, computer vision, and medical imaging. For a given degraded imageg, the a-posteriori likelihood function pf|g via Bayesian modeling is given by pf|g=pg|fpf/pg. To determine an approximation to the unknown image f, the MAP estimator can be used to minimize the negative log-likelihood function, i.e.,(4)f∗=minf−logpfg=minf−logpgf−logpf,where the additive term logpg is neglected due to its independence to f. In practice, the probability density pg|f should be assumed suitably according to the model of the noise in the given image g. In the case of impulse noise, the negative log-likelihood function −logpg|f is commonly defined as a L1 data-fidelity term, i.e., −logpg|f=∫Ωℋf−gdx. To stabilize the image restoration, it is necessary to incorporate the additional prior information via the a-priori probability density function (PDF) pf into the restoration process. The most frequently used PDFs are Gibbs functions in current literature which satisfy(5)pf∝exp−ςℛf,where ς is a positive scale parameter and ℛf denotes a regularization functional. If the Gibbs a-priori PDF equation (5) is incorporated into the MAP framework equation (4), minimizing the negative log-likelihood function −logpf|g is equivalent to solving the following minimization problem:(6)f∗=minfλ∫Ωℋf−gdx+ℛf,where λ>0 is the regularization parameter which controls the tradeoff between the data-fidelity and regularization terms. In [29, 33], the regularization term ℛf is selected as a TV prior, i.e., ℛf=∫Ω∇fdx. However, the results often tend to cause staircase-like artifacts since TV regularizer favors a piecewise constant solution in bounded variation (BV) space. The undesirable solutions could introduce false edges that do not exist within the true image and then fail to guarantee high visual quality. To improve the image restoration capability, Liu et al. [37] proposed to replace TV by TV1,2, i.e., ℛf=∫Ωξ∇f+1−ξ∇2fdx, where ξ∈0,1 denotes the gradient dependent weight. Essentially, both TV and TV1,2 are local regularizers, which have constrained the further performance improvement in image restoration. There is a significant potential to incorporate the nonlocal regularizer into the image restoration framework equation (6). In particular, the nonlocal regularizer is well capable of taking full advantage of the pattern redundancy and self-similarity within an image [24]. Thus, it can generate superior restoration performance compared with local regularizers. In next section, we will propose a nonlocal regularized L1 variational model for restoration of blurred images corrupted by impulse noise.


## 2.1. Image Degradation Model

Images are often degraded by blurring and different kinds of noise. In this paper, we will concentrate mainly on the impulse noise, such as salt-and-pepper noise and random-valued impulse noise. We denote byΩ⊂ℝ2 the image domain, which is often considered to be a bounded and piecewise smooth open subset. Let f:Ω⟶ℝ be a real function defined on Ω. The corrupted image g:Ω⟶ℝ is then given by(1)g=Nℋf,where N denotes the image degradation by impulse noise. For the sake of simplicity, we assume that the linear and continuous blurring operator ℋ∈ℒL2Ω is known beforehand. Denote the dynamic range of ℋf to be fmin,fmax, i.e., fmin≤ℋfx≤fmax,∀x∈Ω, the degradation model N is defined as follows [47].(i)
Salt-and-pepper noise: the gray level ofg at pixel location x∈Ω is(2)gx=fmin,with probabilityp1/2,fmax,with probabilityp1/2,ℋfx,with probability1−p1,wherep1 determines the level of the salt-and-pepper noise. The gray values of degraded pixels change to either the minimum value fmin or the maximum value fmax.(ii)
Random-valued impulse noise: the gray level ofg at pixel location x∈Ω is(3)gx=f^x,with probabilityp2,ℋfx,with probability1−p2,where0≤p2≤1 is the level of the random-valued impulse noise. The gray values of f^x are identically and uniformly distributed random numbers in fmin,fmax.Recently, many variational image restoration models have been developed by combining the nonsmooth L1 data-fidelity term with different regularizers. In this paper, we will restrict our attention to the local TV regularizer and its nonlocal extension. The corresponding MAP-based image restoration model will be proposed in the next subsection.


## 2.2. MAP-Based Image Restoration Model

The MAP estimation theory, which inherently includes prior constraints in the form of prior probability density functions, has attracted increasing attention in the field of image processing, computer vision, and medical imaging. For a given degraded imageg, the a-posteriori likelihood function pf|g via Bayesian modeling is given by pf|g=pg|fpf/pg. To determine an approximation to the unknown image f, the MAP estimator can be used to minimize the negative log-likelihood function, i.e.,(4)f∗=minf−logpfg=minf−logpgf−logpf,where the additive term logpg is neglected due to its independence to f. In practice, the probability density pg|f should be assumed suitably according to the model of the noise in the given image g. In the case of impulse noise, the negative log-likelihood function −logpg|f is commonly defined as a L1 data-fidelity term, i.e., −logpg|f=∫Ωℋf−gdx. To stabilize the image restoration, it is necessary to incorporate the additional prior information via the a-priori probability density function (PDF) pf into the restoration process. The most frequently used PDFs are Gibbs functions in current literature which satisfy(5)pf∝exp−ςℛf,where ς is a positive scale parameter and ℛf denotes a regularization functional. If the Gibbs a-priori PDF equation (5) is incorporated into the MAP framework equation (4), minimizing the negative log-likelihood function −logpf|g is equivalent to solving the following minimization problem:(6)f∗=minfλ∫Ωℋf−gdx+ℛf,where λ>0 is the regularization parameter which controls the tradeoff between the data-fidelity and regularization terms. In [29, 33], the regularization term ℛf is selected as a TV prior, i.e., ℛf=∫Ω∇fdx. However, the results often tend to cause staircase-like artifacts since TV regularizer favors a piecewise constant solution in bounded variation (BV) space. The undesirable solutions could introduce false edges that do not exist within the true image and then fail to guarantee high visual quality. To improve the image restoration capability, Liu et al. [37] proposed to replace TV by TV1,2, i.e., ℛf=∫Ωξ∇f+1−ξ∇2fdx, where ξ∈0,1 denotes the gradient dependent weight. Essentially, both TV and TV1,2 are local regularizers, which have constrained the further performance improvement in image restoration. There is a significant potential to incorporate the nonlocal regularizer into the image restoration framework equation (6). In particular, the nonlocal regularizer is well capable of taking full advantage of the pattern redundancy and self-similarity within an image [24]. Thus, it can generate superior restoration performance compared with local regularizers. In next section, we will propose a nonlocal regularized L1 variational model for restoration of blurred images corrupted by impulse noise.


## 3. NLTV Regularized L1 Image Restoration Model

### 3.1. Nonlocal Total Variation Functional

The neighborhood filter NLM, first proposed by Buades et al. [24] for Gaussian noise removal, has attracted increasing attention in recent years. The underlying idea of NLM is to make full use of the nonlocal similar patches within an image. To further enhance image quality, Gilboa and Osher [26] have extended TV to NLTV by combining the concepts of NLM and graph Laplacian. In the following, we will introduce some definitions and notations regarding the NLTV functional [26, 27, 48] that will be used throughout this paper. Let Ω⊂ℝ2 denote the image domain. Given a reference image u, the similarity between two pixels x,y∈Ω can be measured using a weight function ωx,y, which satisfies(7)ωx,y=exp−∫ΩNGσ˜zux+z−uy+z2dzh2,where ΩN=−RN−1/2,RN−1/2×−RN−1/2,RN−1/2. of size RN×RN denotes the neighborhood window, Gσ˜ is the Gaussian kernel with standard deviation σ˜, and h is a positive smoothing factor to control the degree of filtering. Let fx be a real function Ω⟶ℝ. Given the symmetric weight function ωu from a reference image u, the nonlocal gradient ∇NLfx,° at point x∈Ω is defined as follows:(8)∇NLfx,x+s≔fx+s−fxωux,x+s,where s∈ΩS with ΩS=−RS−1/2,RS−1/2×−RS−1/2,RS−1/2 of size RS×RS denoting the search region. The magnitude of ∇NLf at point x∈Ω is then obtained:(9)∇NLfx≔∫ΩSfx+s−fx2ωux,x+sds.LetdivNL denote the nonlocal divergence operator, which is the adjoint of the nonlocal gradient, i.e., divNL=∇NLT. The divergence for a nonlocal vector q∈Ω×Ω⟶ℝ can be defined as follows:(10)divNLqx≔∫ΩSqx,x+s−qx+s,xωx,x+sds,which satisfies the conjugation with the nonlocal gradient, i.e., ∇NLf,q=f,−divNLq. Thus, the nonlocal Laplacian operator ΔNLfx and NLTV regularizer JNLTV can now be, respectively, defined by(11)ΔNLfx≔12divNL∇NLfx=∫ΩSfx+s−fxωx,x+sds,and(12)JNLTV≔∫Ω∇NLfdx=∫Ω∫ΩSfx+s−fx2ωx,x+sdsdx.NLTV regularizer has achieved significant success due to its weighted nonlocal gradient. In conventional local regularizers (e.g., TV [6], TV2 [14, 15], TGV [20], and TV1,2 [17, 37]), the gradient is calculated using only a few neighbors of the target pixel; whereas all pixels are considered as the neighbors in NLTV regularizer. More detailedly, weights in local regularizers are applied equally for the gradients of all pixels within an image, whereas the weights in NLTV regularizer are spatially adaptive according to the similarities between the neighborhoods of the targeting and neighboring pixels [24].


### 3.2. Nonlocal Regularized L1 Variational Model

Due to the significant advantage of NLTV regularizer, it is necessary to incorporate NLTV into the MAP-based image restoration model equation (6). In this case, we assume that the original image f follows a Gibbs prior of the form:(13)pf∝exp−ςJNLTV.Using the MAP estimation theory introduced in Section2.2, the minimization problem equation (6) can be rewritten as follows:(14)f∗=minfλ∫Ωℋf−gdx+∫Ω∇NLfdx,where λ>0 is a predefined regularization parameter. This variational image restoration model (termed NLTV-L1 for short) consists of two terms: the first term, called the L1 data-fidelity term, denotes a measure of the difference between the restored data and the observed version. The second term is the NLTV regularizer which can stabilize the restoration result. Due to nonsmooth natures of these two terms, it is intractable to efficiently solve the unconstrained minimization problem equation (14) via commonly used numerical methods. To achieve a robust and efficient solution, an alternating optimization algorithm will be developed to solve the NLTV-L1 model.


## 3.1. Nonlocal Total Variation Functional

The neighborhood filter NLM, first proposed by Buades et al. [24] for Gaussian noise removal, has attracted increasing attention in recent years. The underlying idea of NLM is to make full use of the nonlocal similar patches within an image. To further enhance image quality, Gilboa and Osher [26] have extended TV to NLTV by combining the concepts of NLM and graph Laplacian. In the following, we will introduce some definitions and notations regarding the NLTV functional [26, 27, 48] that will be used throughout this paper. Let Ω⊂ℝ2 denote the image domain. Given a reference image u, the similarity between two pixels x,y∈Ω can be measured using a weight function ωx,y, which satisfies(7)ωx,y=exp−∫ΩNGσ˜zux+z−uy+z2dzh2,where ΩN=−RN−1/2,RN−1/2×−RN−1/2,RN−1/2. of size RN×RN denotes the neighborhood window, Gσ˜ is the Gaussian kernel with standard deviation σ˜, and h is a positive smoothing factor to control the degree of filtering. Let fx be a real function Ω⟶ℝ. Given the symmetric weight function ωu from a reference image u, the nonlocal gradient ∇NLfx,° at point x∈Ω is defined as follows:(8)∇NLfx,x+s≔fx+s−fxωux,x+s,where s∈ΩS with ΩS=−RS−1/2,RS−1/2×−RS−1/2,RS−1/2 of size RS×RS denoting the search region. The magnitude of ∇NLf at point x∈Ω is then obtained:(9)∇NLfx≔∫ΩSfx+s−fx2ωux,x+sds.LetdivNL denote the nonlocal divergence operator, which is the adjoint of the nonlocal gradient, i.e., divNL=∇NLT. The divergence for a nonlocal vector q∈Ω×Ω⟶ℝ can be defined as follows:(10)divNLqx≔∫ΩSqx,x+s−qx+s,xωx,x+sds,which satisfies the conjugation with the nonlocal gradient, i.e., ∇NLf,q=f,−divNLq. Thus, the nonlocal Laplacian operator ΔNLfx and NLTV regularizer JNLTV can now be, respectively, defined by(11)ΔNLfx≔12divNL∇NLfx=∫ΩSfx+s−fxωx,x+sds,and(12)JNLTV≔∫Ω∇NLfdx=∫Ω∫ΩSfx+s−fx2ωx,x+sdsdx.NLTV regularizer has achieved significant success due to its weighted nonlocal gradient. In conventional local regularizers (e.g., TV [6], TV2 [14, 15], TGV [20], and TV1,2 [17, 37]), the gradient is calculated using only a few neighbors of the target pixel; whereas all pixels are considered as the neighbors in NLTV regularizer. More detailedly, weights in local regularizers are applied equally for the gradients of all pixels within an image, whereas the weights in NLTV regularizer are spatially adaptive according to the similarities between the neighborhoods of the targeting and neighboring pixels [24].


## 3.2. Nonlocal Regularized L1 Variational Model

Due to the significant advantage of NLTV regularizer, it is necessary to incorporate NLTV into the MAP-based image restoration model equation (6). In this case, we assume that the original image f follows a Gibbs prior of the form:(13)pf∝exp−ςJNLTV.Using the MAP estimation theory introduced in Section2.2, the minimization problem equation (6) can be rewritten as follows:(14)f∗=minfλ∫Ωℋf−gdx+∫Ω∇NLfdx,where λ>0 is a predefined regularization parameter. This variational image restoration model (termed NLTV-L1 for short) consists of two terms: the first term, called the L1 data-fidelity term, denotes a measure of the difference between the restored data and the observed version. The second term is the NLTV regularizer which can stabilize the restoration result. Due to nonsmooth natures of these two terms, it is intractable to efficiently solve the unconstrained minimization problem equation (14) via commonly used numerical methods. To achieve a robust and efficient solution, an alternating optimization algorithm will be developed to solve the NLTV-L1 model.


## 4. Numerical Optimization Algorithm

In this section, we propose to develop an alternating optimization algorithm for solving the NLTV-L1 model equation (14). The main advantage of this algorithm is that the two nonsmooth terms in equation (14) can be effectively decoupled. In each step of this algorithm, the corresponding subproblem can be efficiently solved through existing optimization algorithms. In particular, we will deal with the L1-related minimization problem in the first step via a soft shrinkage operator. To further enhance the image quality, the second step related with NLTV regularized L2 deconvolution will be implemented using the preconditioned Bregmanized operator splitting (PBOS) proposed by Zhang et al. [27].


### 4.1. Overview of PBOS

We tend to give a brief overview of PBOS used in this paper. Consider the following unconstrained optimization problem:(15)f¯∗=minf¯12∫Ωℋf¯−g¯2dx+βℛf¯,where β>0 is a predefined regularization parameter and ℛf¯ denotes a general convex regularizer, such as TV [6], TV2 [14, 15], TGV [20], TV1,2 [17, 37]), and NLTV [27]. By combining the Bregman iterative scheme [49] and forward-background operator splitting technique [50], the Bregmanized operator splitting (BOS) method is introduced to effectively solve equation (15) as follows:(16)v¯k+1=f¯k−δℋ⊤ℋf¯k−w¯k,f¯k+1=minf¯12δ∫Ωf¯−v¯k+12dx+βℛf¯,w¯k+1=w¯k−ℋf¯k+1−g¯,for a positive parameter 0<δ<1/ℋ⊤ℋ2. The superscript ⊤ denotes the transpose (conjugate transpose) operator for real (complex) matrices or vectors. In this paper, we mainly focus on the NLTV regularizer, i.e., ℛf¯=∫Ω∇NLf¯dx. However, the image restoration performance is highly dependent on the updating of weight function ω in equation (7). In the case of image denoising, ω can be estimated accurately because the noisy image still contains the most image similarity information. In contrast, it is difficult to achieve the satisfactory ω for image deblurring due to loss of structural similarities in blurred image. To overcome this limitation, Lou et al. [51] and Zhang et al. [27] proposed to estimate ω from the preconditioned image generated by Tikhonov’s regularization method. Although the preconditioned image would suffer from noise amplification, the main geometric structure preserved will be beneficial for the accurate estimation of ω since the updating scheme equation (7) is insensitive to noise. Thus, the preconditioned Bregmanized operator splitting (PBOS) for solving equation (15) can be defined as follows:(17)v¯k+1=f¯k−δℋ⊤ℋℋ⊤+ε−1ℋf¯k−w¯k,f¯k+1=minf¯12δ∫Ωf¯−v¯k+12dx+βℛf¯,w¯k+1=w¯k−ℋf¯k+1−g¯,with 0<δ<1/ℋ⊤ℋℋ⊤+ϵ−1ℋ2, where ϵ>0 is a small regularization parameter to stabilize the solution of v¯-subproblem in equation (17).


### 4.2. Alternating Optimization Algorithm for the NLTV-L1 Model

We introduce one auxiliary variableu=ℋf−g in equation (14), which can be written as an equivalent unconstrained optimization problem as follows:(18)minf,uλ∫Ωudx+α2∫Ωu−ℋf−g2dx+∫Ω∇NLfdx,where α>0 is a penalty parameter. It could be easily proved that, as α⟶+∞, the solution of equation (18) is convergent to the solution of equation (14). However, when α is not so large, the approximation in equations (18) to (14) is inaccurate in practice. We note that equation (18) can be rewritten as follows:(19)minf,uλ∫Ωudx+α2∫Ωu−ℋf−g2dx+∫Ω∇NLfdx=minfminuλ∫Ωudx+α2∫Ωu−ℋf−g2dx+∫Ω∇NLfdx.In this work, the solution of equation (18) can be obtained by applying the alternating optimization algorithm. Staring from an initial guess u0 and f0, we compute the sequence(20)u1,f1,u2,f2,…,uk,fk,uk+1,fk+1,…,via the following formulas:(21)uk+1=minuλ∫Ωudx+α2∫Ωu−ℋfk−g2dx.(22)fk+1=minfα2∫Ωℋf−uk+1+g2dx+∫Ω∇NLfdx.Now, we alternatively solve equation (18) with respect to u and f until the solution converges to an optimal solution.


#### 4.2.1. Solving theu-Subproblem Equation (21)

The first step of the alternating optimization algorithm is to perform the L1 minimization problem. For fixedfk, the solution uk+1 of u-subproblem in equation (21) is equivalent to solving the following minimization problem:(23)uk+1=minuν∫Ωudx+∫Ωu−vk2dx,where ν=2λ/α and vk=ℋfk−g. The exact minimizer of ϕu=u−v2+νu can be obtained by a soft shrinkage operator [52, 53] as follows(24)u=ψv=vvmaxv−ν2,0.Therefore, the solution of the optimization problem equation (21) can be obtained using uk+1=ψℋfk−g, i.e.,(25)uk+1=ℋfk−gℋfk−gmaxℋfk−g−λα,0.


#### 4.2.2. Solving thef-Subproblem Equation (22)

The last step is to apply a nonlocal variational model to recover the image generated by the previous step. For fixeduk+1, the solution fk+1 of f-subproblem in equation (22) is equivalent to optimizing the following NLTV regularized L2 minimization model:(26)fk+1=minf12∫Ωℋf−uk+1+g2dx+1α∫Ω∇NLfdx.As introduced in Section4.1, PBOS can be directly adopted for NLTV regularized L2 minimization problem. We first replace uk+1+g by g¯k+1, the solution of f-subproblem equation (26) can be achieved by considering the following alternating iteration scheme:(27)vi+1=fk,i−δℋ⊤ℋℋ⊤+ϵ−1ℋfk,i−wi,fk,i+1=minf12δ∫Ωf−vi+12dx+1α∫Ω∇NLfdx,wi+1=wi−ℋfk,i+1−g¯k+1,for i=0,1,…,Imax. The variable f satisfies the following updating scheme, i.e., fk,0=fk and fk+1=fk,Imax. Both the first and third steps in equation (27) can be updated easily, but it is difficult to effectively update the second step due to the nonsmooth nature of NLTV regularizer. In the current literature, the split Bregman method [54] has been successfully used in NLTV-based multiplicative noise removal [48]. To apply this method, we first replace ∇NLf by d and add a penalty function term. The solution of f-subproblem in equation (27) can be achieved by considering the following split Bregman iteration equation:(28)dj+1,fk,i,j+1=mind,f12∫Ωf−vi+12dx+η2∫Ωd−∇NLf−bj2dx+μ∫Ωddx,with μ=δ/α and fk,i,0=fk,i. Here, η>0 is a penalty parameter and the variable b is updated using bj+1=bj+∇NLfk,i,j+1−dj+1 for j=0,1,…,Jmax. The variable fk,i+1=fk,i,Jmax can be obtained for updating the third step in equation (27).Thed-subproblem in equation (28) is in essence a convex L1 minimization problem:(29)dj+1=mindη2∫Ωd−∇NLfk,i,j−bj2dx+μ∫Ωddx,which can be explicitly computed using the following soft shrinkage operator [52, 54], i.e.,(30)dj+1=∇NLfk,i,j+bj∇NLfk,i,j+bjmax∇NLfk,i,j+bj−μη,0.This shrinkage operator is extremely fast and requires only a few operations per element ofdj+1. Given the fixed vi+1, dj+1, and bj, solution fk,i,j+1 of the f-subproblem in equation (28) is equivalent to solving the following minimization probem:(31)fk,i,j+1=minf12∫Ωf−vi+12dx+η2∫Ωdj+1−∇NLf−bj2dx.We can easily obtain the optimal solutionfk,i,j+1 because the f-subproblem is both smooth and differentiable. In particular, the solution fk,i,j+1 consists in solving the following system of linear equation:(32)fk,i,j+1−vi+1−ηdivNL∇NLfk,i,j+1+bj−dj+1=0,which provides(33)fk,i,j+1=I−2ηΔNL−1vi+1+ηdivNLbj−dj+1,where I denotes the identity element. Since the nonlocal Laplacian ΔNL is negative semidefinite, the operator I−2ηΔNL is diagonally dominant with nonlocal weight ω. To ensure optimal efficiency, the widely used GaussSeidel method can be adopted to achieve an approximate solution of equation (33). This method is capable of guaranteeing a fast convergence. The implementation of PBOS for equation (26) is described by the following pseudocode shown in Algorithm 1.Algorithm 1: PBOS for NLTV regularized L2 minimization problem (26).
(1)
Input: variables ℋ,g¯k+1 and parameters ϵ,α,δ,η.(2)
Initialize: variables fk,0,w0,b0.(3)
fori=0 to Imaxdo(4)
updatevi+1: vi+1=fk,i−δℋ⊤ℋℋ⊤+ϵ−1ℋfk,i−wi.(5)
update weight functionω using formula (7).(6)
fk,i,0←fk,i.(7)
whilea stopping criterion is not satisfieddo(8)
updatedj+1: dj+1=∇NLfk,i,j+bj/∇NLfk,i,j+bjmax∇NLfk,i,j+bj−δ/αη,0.(9)
updatefk,i,j+1: solve fk,i,j+1=I−2ηΔNL−1vi+1+ηdivNLbj−dj+1 using the Gauss-Seidel method.(10)
updatebj+1: bj+1=bj+∇NLfk,i,j+1−dj+1.(11)
end while(12)
fk,i+1←fk,i,Jmax.(13)
updatewi+1: wi+1=wi−ℋfk,i+1−uk+1+g.(14)
end for(15)
fk+1←fk,Imax.


### 4.3. Optimization Procedure

As described in Section4.2, an alternating optimization algorithm was proposed to decompose the original NLTV-L1 minimization problem equation (14) into two simpler subproblems, i.e., u-subproblem equation (21) and f-subproblem equation (22). Each subproblem has a closed-form solution or could be efficiently handled using existing optimization algorithms. In particular, the u-subproblem was solved by a soft shrinkage operator. Solution of the f-subproblem was achieved using PBOS. In conclusion, Algorithm 2 shows the pseudocode of our alternating optimization algorithm for the NLTV-L1 minimization problem equation (14).Algorithm 2: Alternating optimization algorithm for NLTV-L1 model equation (14).
(1)
Input: variables ℋ,g and parameters λ,α.(2)
Initialize: variables u0,f0.(3)
fork=0 to Kmaxdo//Stage 1: u-subproblem(4)
updateuk+1: uk+1=ℋfk−g/ℋfk−gmaxℋfk−g−λ/α,0.//Stage 2: f-subproblem(5)
Solve the following optimization problem by PBOS shown in Algorithm 1fk+1=minf1/2∫Ωℋf−uk+1+g2dx+1/α∫Ω∇NLfdx(6)
end for(7)
Output: Recovered image f.


## 4.1. Overview of PBOS

We tend to give a brief overview of PBOS used in this paper. Consider the following unconstrained optimization problem:(15)f¯∗=minf¯12∫Ωℋf¯−g¯2dx+βℛf¯,where β>0 is a predefined regularization parameter and ℛf¯ denotes a general convex regularizer, such as TV [6], TV2 [14, 15], TGV [20], TV1,2 [17, 37]), and NLTV [27]. By combining the Bregman iterative scheme [49] and forward-background operator splitting technique [50], the Bregmanized operator splitting (BOS) method is introduced to effectively solve equation (15) as follows:(16)v¯k+1=f¯k−δℋ⊤ℋf¯k−w¯k,f¯k+1=minf¯12δ∫Ωf¯−v¯k+12dx+βℛf¯,w¯k+1=w¯k−ℋf¯k+1−g¯,for a positive parameter 0<δ<1/ℋ⊤ℋ2. The superscript ⊤ denotes the transpose (conjugate transpose) operator for real (complex) matrices or vectors. In this paper, we mainly focus on the NLTV regularizer, i.e., ℛf¯=∫Ω∇NLf¯dx. However, the image restoration performance is highly dependent on the updating of weight function ω in equation (7). In the case of image denoising, ω can be estimated accurately because the noisy image still contains the most image similarity information. In contrast, it is difficult to achieve the satisfactory ω for image deblurring due to loss of structural similarities in blurred image. To overcome this limitation, Lou et al. [51] and Zhang et al. [27] proposed to estimate ω from the preconditioned image generated by Tikhonov’s regularization method. Although the preconditioned image would suffer from noise amplification, the main geometric structure preserved will be beneficial for the accurate estimation of ω since the updating scheme equation (7) is insensitive to noise. Thus, the preconditioned Bregmanized operator splitting (PBOS) for solving equation (15) can be defined as follows:(17)v¯k+1=f¯k−δℋ⊤ℋℋ⊤+ε−1ℋf¯k−w¯k,f¯k+1=minf¯12δ∫Ωf¯−v¯k+12dx+βℛf¯,w¯k+1=w¯k−ℋf¯k+1−g¯,with 0<δ<1/ℋ⊤ℋℋ⊤+ϵ−1ℋ2, where ϵ>0 is a small regularization parameter to stabilize the solution of v¯-subproblem in equation (17).


## 4.2. Alternating Optimization Algorithm for the NLTV-L1 Model

We introduce one auxiliary variableu=ℋf−g in equation (14), which can be written as an equivalent unconstrained optimization problem as follows:(18)minf,uλ∫Ωudx+α2∫Ωu−ℋf−g2dx+∫Ω∇NLfdx,where α>0 is a penalty parameter. It could be easily proved that, as α⟶+∞, the solution of equation (18) is convergent to the solution of equation (14). However, when α is not so large, the approximation in equations (18) to (14) is inaccurate in practice. We note that equation (18) can be rewritten as follows:(19)minf,uλ∫Ωudx+α2∫Ωu−ℋf−g2dx+∫Ω∇NLfdx=minfminuλ∫Ωudx+α2∫Ωu−ℋf−g2dx+∫Ω∇NLfdx.In this work, the solution of equation (18) can be obtained by applying the alternating optimization algorithm. Staring from an initial guess u0 and f0, we compute the sequence(20)u1,f1,u2,f2,…,uk,fk,uk+1,fk+1,…,via the following formulas:(21)uk+1=minuλ∫Ωudx+α2∫Ωu−ℋfk−g2dx.(22)fk+1=minfα2∫Ωℋf−uk+1+g2dx+∫Ω∇NLfdx.Now, we alternatively solve equation (18) with respect to u and f until the solution converges to an optimal solution.


### 4.2.1. Solving theu-Subproblem Equation (21)

The first step of the alternating optimization algorithm is to perform the L1 minimization problem. For fixedfk, the solution uk+1 of u-subproblem in equation (21) is equivalent to solving the following minimization problem:(23)uk+1=minuν∫Ωudx+∫Ωu−vk2dx,where ν=2λ/α and vk=ℋfk−g. The exact minimizer of ϕu=u−v2+νu can be obtained by a soft shrinkage operator [52, 53] as follows(24)u=ψv=vvmaxv−ν2,0.Therefore, the solution of the optimization problem equation (21) can be obtained using uk+1=ψℋfk−g, i.e.,(25)uk+1=ℋfk−gℋfk−gmaxℋfk−g−λα,0.


### 4.2.2. Solving thef-Subproblem Equation (22)

The last step is to apply a nonlocal variational model to recover the image generated by the previous step. For fixeduk+1, the solution fk+1 of f-subproblem in equation (22) is equivalent to optimizing the following NLTV regularized L2 minimization model:(26)fk+1=minf12∫Ωℋf−uk+1+g2dx+1α∫Ω∇NLfdx.As introduced in Section4.1, PBOS can be directly adopted for NLTV regularized L2 minimization problem. We first replace uk+1+g by g¯k+1, the solution of f-subproblem equation (26) can be achieved by considering the following alternating iteration scheme:(27)vi+1=fk,i−δℋ⊤ℋℋ⊤+ϵ−1ℋfk,i−wi,fk,i+1=minf12δ∫Ωf−vi+12dx+1α∫Ω∇NLfdx,wi+1=wi−ℋfk,i+1−g¯k+1,for i=0,1,…,Imax. The variable f satisfies the following updating scheme, i.e., fk,0=fk and fk+1=fk,Imax. Both the first and third steps in equation (27) can be updated easily, but it is difficult to effectively update the second step due to the nonsmooth nature of NLTV regularizer. In the current literature, the split Bregman method [54] has been successfully used in NLTV-based multiplicative noise removal [48]. To apply this method, we first replace ∇NLf by d and add a penalty function term. The solution of f-subproblem in equation (27) can be achieved by considering the following split Bregman iteration equation:(28)dj+1,fk,i,j+1=mind,f12∫Ωf−vi+12dx+η2∫Ωd−∇NLf−bj2dx+μ∫Ωddx,with μ=δ/α and fk,i,0=fk,i. Here, η>0 is a penalty parameter and the variable b is updated using bj+1=bj+∇NLfk,i,j+1−dj+1 for j=0,1,…,Jmax. The variable fk,i+1=fk,i,Jmax can be obtained for updating the third step in equation (27).Thed-subproblem in equation (28) is in essence a convex L1 minimization problem:(29)dj+1=mindη2∫Ωd−∇NLfk,i,j−bj2dx+μ∫Ωddx,which can be explicitly computed using the following soft shrinkage operator [52, 54], i.e.,(30)dj+1=∇NLfk,i,j+bj∇NLfk,i,j+bjmax∇NLfk,i,j+bj−μη,0.This shrinkage operator is extremely fast and requires only a few operations per element ofdj+1. Given the fixed vi+1, dj+1, and bj, solution fk,i,j+1 of the f-subproblem in equation (28) is equivalent to solving the following minimization probem:(31)fk,i,j+1=minf12∫Ωf−vi+12dx+η2∫Ωdj+1−∇NLf−bj2dx.We can easily obtain the optimal solutionfk,i,j+1 because the f-subproblem is both smooth and differentiable. In particular, the solution fk,i,j+1 consists in solving the following system of linear equation:(32)fk,i,j+1−vi+1−ηdivNL∇NLfk,i,j+1+bj−dj+1=0,which provides(33)fk,i,j+1=I−2ηΔNL−1vi+1+ηdivNLbj−dj+1,where I denotes the identity element. Since the nonlocal Laplacian ΔNL is negative semidefinite, the operator I−2ηΔNL is diagonally dominant with nonlocal weight ω. To ensure optimal efficiency, the widely used GaussSeidel method can be adopted to achieve an approximate solution of equation (33). This method is capable of guaranteeing a fast convergence. The implementation of PBOS for equation (26) is described by the following pseudocode shown in Algorithm 1.Algorithm 1: PBOS for NLTV regularized L2 minimization problem (26).
(1)
Input: variables ℋ,g¯k+1 and parameters ϵ,α,δ,η.(2)
Initialize: variables fk,0,w0,b0.(3)
fori=0 to Imaxdo(4)
updatevi+1: vi+1=fk,i−δℋ⊤ℋℋ⊤+ϵ−1ℋfk,i−wi.(5)
update weight functionω using formula (7).(6)
fk,i,0←fk,i.(7)
whilea stopping criterion is not satisfieddo(8)
updatedj+1: dj+1=∇NLfk,i,j+bj/∇NLfk,i,j+bjmax∇NLfk,i,j+bj−δ/αη,0.(9)
updatefk,i,j+1: solve fk,i,j+1=I−2ηΔNL−1vi+1+ηdivNLbj−dj+1 using the Gauss-Seidel method.(10)
updatebj+1: bj+1=bj+∇NLfk,i,j+1−dj+1.(11)
end while(12)
fk,i+1←fk,i,Jmax.(13)
updatewi+1: wi+1=wi−ℋfk,i+1−uk+1+g.(14)
end for(15)
fk+1←fk,Imax.


## 4.2.1. Solving theu-Subproblem Equation (21)

The first step of the alternating optimization algorithm is to perform the L1 minimization problem. For fixedfk, the solution uk+1 of u-subproblem in equation (21) is equivalent to solving the following minimization problem:(23)uk+1=minuν∫Ωudx+∫Ωu−vk2dx,where ν=2λ/α and vk=ℋfk−g. The exact minimizer of ϕu=u−v2+νu can be obtained by a soft shrinkage operator [52, 53] as follows(24)u=ψv=vvmaxv−ν2,0.Therefore, the solution of the optimization problem equation (21) can be obtained using uk+1=ψℋfk−g, i.e.,(25)uk+1=ℋfk−gℋfk−gmaxℋfk−g−λα,0.


## 4.2.2. Solving thef-Subproblem Equation (22)

The last step is to apply a nonlocal variational model to recover the image generated by the previous step. For fixeduk+1, the solution fk+1 of f-subproblem in equation (22) is equivalent to optimizing the following NLTV regularized L2 minimization model:(26)fk+1=minf12∫Ωℋf−uk+1+g2dx+1α∫Ω∇NLfdx.As introduced in Section4.1, PBOS can be directly adopted for NLTV regularized L2 minimization problem. We first replace uk+1+g by g¯k+1, the solution of f-subproblem equation (26) can be achieved by considering the following alternating iteration scheme:(27)vi+1=fk,i−δℋ⊤ℋℋ⊤+ϵ−1ℋfk,i−wi,fk,i+1=minf12δ∫Ωf−vi+12dx+1α∫Ω∇NLfdx,wi+1=wi−ℋfk,i+1−g¯k+1,for i=0,1,…,Imax. The variable f satisfies the following updating scheme, i.e., fk,0=fk and fk+1=fk,Imax. Both the first and third steps in equation (27) can be updated easily, but it is difficult to effectively update the second step due to the nonsmooth nature of NLTV regularizer. In the current literature, the split Bregman method [54] has been successfully used in NLTV-based multiplicative noise removal [48]. To apply this method, we first replace ∇NLf by d and add a penalty function term. The solution of f-subproblem in equation (27) can be achieved by considering the following split Bregman iteration equation:(28)dj+1,fk,i,j+1=mind,f12∫Ωf−vi+12dx+η2∫Ωd−∇NLf−bj2dx+μ∫Ωddx,with μ=δ/α and fk,i,0=fk,i. Here, η>0 is a penalty parameter and the variable b is updated using bj+1=bj+∇NLfk,i,j+1−dj+1 for j=0,1,…,Jmax. The variable fk,i+1=fk,i,Jmax can be obtained for updating the third step in equation (27).Thed-subproblem in equation (28) is in essence a convex L1 minimization problem:(29)dj+1=mindη2∫Ωd−∇NLfk,i,j−bj2dx+μ∫Ωddx,which can be explicitly computed using the following soft shrinkage operator [52, 54], i.e.,(30)dj+1=∇NLfk,i,j+bj∇NLfk,i,j+bjmax∇NLfk,i,j+bj−μη,0.This shrinkage operator is extremely fast and requires only a few operations per element ofdj+1. Given the fixed vi+1, dj+1, and bj, solution fk,i,j+1 of the f-subproblem in equation (28) is equivalent to solving the following minimization probem:(31)fk,i,j+1=minf12∫Ωf−vi+12dx+η2∫Ωdj+1−∇NLf−bj2dx.We can easily obtain the optimal solutionfk,i,j+1 because the f-subproblem is both smooth and differentiable. In particular, the solution fk,i,j+1 consists in solving the following system of linear equation:(32)fk,i,j+1−vi+1−ηdivNL∇NLfk,i,j+1+bj−dj+1=0,which provides(33)fk,i,j+1=I−2ηΔNL−1vi+1+ηdivNLbj−dj+1,where I denotes the identity element. Since the nonlocal Laplacian ΔNL is negative semidefinite, the operator I−2ηΔNL is diagonally dominant with nonlocal weight ω. To ensure optimal efficiency, the widely used GaussSeidel method can be adopted to achieve an approximate solution of equation (33). This method is capable of guaranteeing a fast convergence. The implementation of PBOS for equation (26) is described by the following pseudocode shown in Algorithm 1.Algorithm 1: PBOS for NLTV regularized L2 minimization problem (26).
(1)
Input: variables ℋ,g¯k+1 and parameters ϵ,α,δ,η.(2)
Initialize: variables fk,0,w0,b0.(3)
fori=0 to Imaxdo(4)
updatevi+1: vi+1=fk,i−δℋ⊤ℋℋ⊤+ϵ−1ℋfk,i−wi.(5)
update weight functionω using formula (7).(6)
fk,i,0←fk,i.(7)
whilea stopping criterion is not satisfieddo(8)
updatedj+1: dj+1=∇NLfk,i,j+bj/∇NLfk,i,j+bjmax∇NLfk,i,j+bj−δ/αη,0.(9)
updatefk,i,j+1: solve fk,i,j+1=I−2ηΔNL−1vi+1+ηdivNLbj−dj+1 using the Gauss-Seidel method.(10)
updatebj+1: bj+1=bj+∇NLfk,i,j+1−dj+1.(11)
end while(12)
fk,i+1←fk,i,Jmax.(13)
updatewi+1: wi+1=wi−ℋfk,i+1−uk+1+g.(14)
end for(15)
fk+1←fk,Imax.


## 4.3. Optimization Procedure

As described in Section4.2, an alternating optimization algorithm was proposed to decompose the original NLTV-L1 minimization problem equation (14) into two simpler subproblems, i.e., u-subproblem equation (21) and f-subproblem equation (22). Each subproblem has a closed-form solution or could be efficiently handled using existing optimization algorithms. In particular, the u-subproblem was solved by a soft shrinkage operator. Solution of the f-subproblem was achieved using PBOS. In conclusion, Algorithm 2 shows the pseudocode of our alternating optimization algorithm for the NLTV-L1 minimization problem equation (14).Algorithm 2: Alternating optimization algorithm for NLTV-L1 model equation (14).
(1)
Input: variables ℋ,g and parameters λ,α.(2)
Initialize: variables u0,f0.(3)
fork=0 to Kmaxdo//Stage 1: u-subproblem(4)
updateuk+1: uk+1=ℋfk−g/ℋfk−gmaxℋfk−g−λ/α,0.//Stage 2: f-subproblem(5)
Solve the following optimization problem by PBOS shown in Algorithm 1fk+1=minf1/2∫Ωℋf−uk+1+g2dx+1/α∫Ω∇NLfdx(6)
end for(7)
Output: Recovered image f.


## 5. Automatic Selection of Spatially Adapted Regularization Parameter

It is well known that the regularization parameterλ plays an important role in image restoration. In particular, larger λ easily leads to little noise reduction; whereas smaller λ will result in oversmoothing of texture regions. Since images are always comprised of image features of different scales, the constant parameter λ will degrade the image visual quality [44]. Therefore, to achieve satisfactory restoration performance, the regularization parameter should be selected adaptively according to the local image features. Roughly speaking, the smaller parameter enables effective noise reduction in homogeneous regions (with large-scale features); whereas larger one performs well in fine detail preservation in texture regions (with small-scale features). In [36, 44–46], many efforts have been made to take full advantage of the local image features and achieved remarkable restoration results. To the best of our knowledge, no research has been conducted on NLTV-L1 model with spatially adapted regularization parameter selection thus far. In this section, we will investigate how to adaptively select the regularization parameter for our proposed NLTV-L1 model in terms of both salt-and-pepper and random-valued impulse noise.Motivated by the previous works in [36, 44–46], we employ the similar automated method in spatially adapted regularization parameter selection for our NLTV-L1 model equation (14). Before the selection step is performed, we define a local window centered at pixel x∈Ω as follows:(34)Ωxρ=y:x−y∞≤ρ2,and assume that ϖx,y is a local window filter, i.e.,(35)ϖx,y=1Ωxρ,ifx−y∞≤ρ2,ε0,else,,with ∫Ω∫Ωϖx,ydxdy=1 and 0<ε0≤min1,1/Ωxρ. In this paper, the local expected absolute value estimator is defined as follows:(36)Sfx=∫Ωϖx,yℋf−gydy,for x∈Ω. Based on the formula in (36), we can obtain the following NLTV regularized minimization problem with local constraints:(37)minf∫Ω∇NLfdxs.t.Sf≤πa.e.inΩ,where “a.e.” stands for “almost everywhere,” and π∈ℝ denotes the expected absolute value that depends on the type of noise [36]. It should be noted that the constrained minimization problem (37) is related to an unconstrained optimization model with spatially adapted regularization parameters λ∈L2Ω as follows:(38)f∗=minf∫Ωλℋf−gdx+∫Ω∇NLfdx.For fixed regularization parameterλ, the minimization problem (38) can be solved via the alternating optimization algorithm given in Algorithm 2. Similar to the previous works in [36, 44–46], we propose the updating scheme for the regularization parameter λ as follows:(39)λ˜t+1x=κminλ˜tx+τmaxSftx−πx,0,L.(40)λt+1x=∫Ωϖx,yλ˜t+1ydy,where τ>0 is a step size, κ=2 is a constant parameter adopted in this paper, L is a large positive parameter to ensure that λ˜t+1 stays bounded, and ft is the current estimate of the latent image f. More details about the determination of the parameter τ and expected absolute value π in equation (39) will be discussed, respectively, in Sections 5.1 and 5.2.


### 5.1. Selection ofλ in the Case of Salt-and-Pepper Noise

In the case of salt-and-pepper noise, we are able to obtain the expected absolute valueπx=p1/2 for any x∈Ω, where p1 is the noise level defined in equation (2). To keep the updated λ˜t+1 at the same scale as λ˜t, the parameter τ is set as τ=τt=2λ˜t∞/p1 as done in [36, 44]. Therefore, for given recovered image fk and intermediate variable λ˜t, we can define the following scheme to update the regularization parameter λt+1 as follows:(41)λ˜t+1x=κminλ˜tx+τtmaxSftx−πx,0,L,(42)λt+1x=∫Ωϖx,yλ˜t+1ydy,(43)τt+1=2λ˜t+1∞/p1,where ft is the current estimate of the latent image f which is obtained by implementing the Algorithm 2 with fixed λ=λt. According to the iterative update of λ in equations (41)–(43), we can obtain the NLTV-L1 model with spatially adapted regularization parameters for restoring blurred images with salt-and-pepper noise given by Algorithm 3.Algorithm 3: NLTV-L1 model with spatially adapted regularization parameter selection.
(1)
Input: variables ℋ,g and parameters p1,p2,κ,L.(2)
Initialize: variable f0 and parameters λ0,τ0,π0.(3)
whilea stopping criterion is not satisfieddo(4) Solve the following optimization problem by Algorithm 2ft=minf∫Ωλtℋf−gdx+∫Ω∇NLfdx(5) Updateλt+1 based on the previous ftSftx=∫Ωϖx,yℋft−gydy,λ˜t+1x=κminλ˜tx+τtmaxSftx−πtx,0,L,λt+1x=∫Ωϖx,yλ˜t+1ydy,whereτt=2λ˜t∞/p1,for salt−and−pepper noise,2λ˜t∞/p2,for random−valued impulse noise.andπtx=p1/2,for salt−and−pepper noise,p2∫Ωϖx,yℋft2−ℋft+1/2ydy,for random−valued impulse noise.(6)
end while(7)
Output: Recovered image f.


### 5.2. Selection ofλ in the Case of Random-Valued Impulse Noise

Analogous to the salt-and-pepper noise,τ is set as τ=τt=2λ˜t∞/p2 to keep the updated λ˜t+1 at the same scale as λ˜t. For a given λ˜t associated with a recovered image ft, the current expected absolute value can be calculated as follows (44)πtx=p2∫Ωϖx,yℋft2−ℋft+12ydy,where πx∈p2/4,p2/2 for any x∈Ω. This spatially adapted expected absolute value is significantly more complex than the constant value πx=p1/2 in equation (41). Then, we propose the updating scheme of λ as follows:(45)λ˜t+1x=κminλ˜tx+τtmaxSftx−πtx,0,L,λt+1x=∫Ωϖx,yλ˜t+1ydy,τt+1=2λ˜t+1∞p2.For this type of noise, the recovered imageft+1 can be obtained by implementing the Algorithm 2 with fixed λ=λt+1. As a conclusion of this subsection, the NLTV-L1 model with spatially adapted regularization parameters for restoring blurred images with random-valued impulse noise can also be found in Algorithm 3. For more details on automatic selection of spatially adapted regularization parameters in terms of both salt-and-pepper noise and random-valued impulse noise, we refer the interested reader to see [36] and references therein. Owing to the proposed NLTV-L1 model and automatic selection technique, the new image restoration version developed in this paper is called spatially adapted NLTV-L1 algorithm (termed SA-NLTV-L1 for short).


## 5.1. Selection ofλ in the Case of Salt-and-Pepper Noise

In the case of salt-and-pepper noise, we are able to obtain the expected absolute valueπx=p1/2 for any x∈Ω, where p1 is the noise level defined in equation (2). To keep the updated λ˜t+1 at the same scale as λ˜t, the parameter τ is set as τ=τt=2λ˜t∞/p1 as done in [36, 44]. Therefore, for given recovered image fk and intermediate variable λ˜t, we can define the following scheme to update the regularization parameter λt+1 as follows:(41)λ˜t+1x=κminλ˜tx+τtmaxSftx−πx,0,L,(42)λt+1x=∫Ωϖx,yλ˜t+1ydy,(43)τt+1=2λ˜t+1∞/p1,where ft is the current estimate of the latent image f which is obtained by implementing the Algorithm 2 with fixed λ=λt. According to the iterative update of λ in equations (41)–(43), we can obtain the NLTV-L1 model with spatially adapted regularization parameters for restoring blurred images with salt-and-pepper noise given by Algorithm 3.Algorithm 3: NLTV-L1 model with spatially adapted regularization parameter selection.
(1)
Input: variables ℋ,g and parameters p1,p2,κ,L.(2)
Initialize: variable f0 and parameters λ0,τ0,π0.(3)
whilea stopping criterion is not satisfieddo(4) Solve the following optimization problem by Algorithm 2ft=minf∫Ωλtℋf−gdx+∫Ω∇NLfdx(5) Updateλt+1 based on the previous ftSftx=∫Ωϖx,yℋft−gydy,λ˜t+1x=κminλ˜tx+τtmaxSftx−πtx,0,L,λt+1x=∫Ωϖx,yλ˜t+1ydy,whereτt=2λ˜t∞/p1,for salt−and−pepper noise,2λ˜t∞/p2,for random−valued impulse noise.andπtx=p1/2,for salt−and−pepper noise,p2∫Ωϖx,yℋft2−ℋft+1/2ydy,for random−valued impulse noise.(6)
end while(7)
Output: Recovered image f.


## 5.2. Selection ofλ in the Case of Random-Valued Impulse Noise

Analogous to the salt-and-pepper noise,τ is set as τ=τt=2λ˜t∞/p2 to keep the updated λ˜t+1 at the same scale as λ˜t. For a given λ˜t associated with a recovered image ft, the current expected absolute value can be calculated as follows (44)πtx=p2∫Ωϖx,yℋft2−ℋft+12ydy,where πx∈p2/4,p2/2 for any x∈Ω. This spatially adapted expected absolute value is significantly more complex than the constant value πx=p1/2 in equation (41). Then, we propose the updating scheme of λ as follows:(45)λ˜t+1x=κminλ˜tx+τtmaxSftx−πtx,0,L,λt+1x=∫Ωϖx,yλ˜t+1ydy,τt+1=2λ˜t+1∞p2.For this type of noise, the recovered imageft+1 can be obtained by implementing the Algorithm 2 with fixed λ=λt+1. As a conclusion of this subsection, the NLTV-L1 model with spatially adapted regularization parameters for restoring blurred images with random-valued impulse noise can also be found in Algorithm 3. For more details on automatic selection of spatially adapted regularization parameters in terms of both salt-and-pepper noise and random-valued impulse noise, we refer the interested reader to see [36] and references therein. Owing to the proposed NLTV-L1 model and automatic selection technique, the new image restoration version developed in this paper is called spatially adapted NLTV-L1 algorithm (termed SA-NLTV-L1 for short).


## 6. Experimental Results and Discussion

In this section, numerical results are presented to illustrate the superior performance of our proposed SA-NLTV-L1 for restoring blurred images with impulse noise. We compare the proposed method with several state-of-the-art image restoration techniques on different numerical experiments. All image restoration methods are implemented using Matlab (The MathWorks, Natick, Inc., MA) on a machine with 3.10 GHz Intel Core i5-2500 CPU and 4 GB RAM. Eight different grayscale images (with size256×256) used throughout this paper are displayed in Figure 1. Using the Matlab built-in function fspecial, both Disk and Gaussian blur kernels with size 15×15, visually illustrated in Figure 2, are generated in our experiments. In particular, these two blur kernels are fspecial (‘disk’, 7) and fspecial (‘gaussian’, [15, 15], 3), respectively. To objectively evaluate the restoration performance, the quality of restored images by different methods is compared quantitatively using three popular assessment methods, i.e., peak-to-noise ratio (PSNR), structural similarity (SSIM) index [55], and feature similarity (FSIM) index [56].Figure 1
Eight different grayscale images of size256×256. From top-left to bottom-right: (a) Barbara. (b) Boats. (c) Cameraman. (d) Lena. (e) Man. (f) Hill. (g) Peppers. (h) House.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 2
Two different blur kernels with size15×15. From left to right: (a) Disk blur kernel fspecial (‘disk’, 7), and (b) Gaussian blur kernel fspecial (‘gaussian’, [15, 15], 3).
(a)(b)


### 6.1. Comparison with Other Image Restoration Techniques

The proposed method will be compared with five image restoration techniques as follows:(1)
TV-L1:TV Regularized L1 Variational Model [29]. TV-L1 was proposed by combining the L1 data-fidelity term with TV regularizer. In particular, this model is typically capable of preserving edge features without amplifying the impulse noise. In this paper, TV-L1 was solved using an alternating minimization algorithm.(2)
NLTV-L1:NLTV Regularized L1 Variational Model. This model takes full advantage of the natural redundancy of patterns within an image. Thus it is capable of suppressing the staircase-like artifacts often yielded by TV-L1 in recovered image. NLTV-L1 was solved using the alternating optimization algorithm given in Algorithm 2.(3)
SA-TV-L1:Spatially Adapted TV Regularized L1 Variational Model [36]. This method is an improved version of the TV-L1 by introducing an automated method for selection of spatially adapted regularization parameter. For the sake of comparison, the corresponding TV-L1 model was also solved using the alternating minimization algorithm proposed in [29].(4)
SA-TV1,2-L1: Spatially Adapted TV1,2Regularized L1 Variational Model [37]. Compared with SA-TV-L1, Liu et al. [37] proposed to replace TV by TV1,2 to further improve the restoration performance. In this paper, the alternating direction method of multipliers (ADMM) was adopted to efficiently solve the resulting optimization problem.The first two variational methods (i.e., TV-L1 and NLTV-L1) were implemented by using only constant regularization parameter. In contrast, SA-TV-L1 and SA-TV1,2-L1 take full advantage of the spatially adapted regularization parameter. However, these three variational methods were proposed based on local regularizers, which could constrain the further performance improvement in image restoration. In contrast, the proposed SA-NLTV-L1 takes into consideration both nonlocal regularizer and spatially adapted regularization parameters. These two properties are capable of enhancing the performance of our image restoration method.


### 6.2. Optimal Parameter Selection

Nonlocal variational image restoration techniques are dependent on several parameters, i.e.,h, RS, and RN. In the current literature, many efforts have been made to investigate the influence of these parameters on image restoration [27, 57, 58]. To the best of our knowledge, no research has been conducted on the NLTV-L1 model for restoring blurred images with impulse noise. It is generally thought that the patch size RN×RN is less important (common choices are between 5 and 9) [58]. Without loss of generality, we set the patch size to be 7×7 throughout the rest of this paper. To keep a good balance between computational complexity and restoration performance, we only select the 10 best and 4 nearest neighbors for each pixel in the searching window ΩS [27]. The searching window size RS×RS is empirically set as 11×11 as conducted in [27]. The smoothing parameter h is also an important parameter that we will study more carefully. To estimate the optimal parameter h, manual experiments are implemented on the eight different grayscale images, shown in Figure 1. In our experiments, the images are corrupted by two different blur kernels and different levels p (i.e., p1 and p2) of impulse noise ranged from 10% to 70% with 10% in step.In Figure3 and Figure 4, the optimal value of h for different images was determined using the image quality measures in terms of PSNR and SSIM. It can be observed that the optimum h gradually becomes larger as noise level p increases. The relationship between h and p can be expressed as a linear approximation formula, i.e., h=α1p+α0. The fitting coefficients α1 and α0 can be estimated empirically from Figures 3 and 4. To illustrate, the approximation results are summarized in Tables 1 and 2. It can be observed that most linear approximation formulas have the similar coefficient α1 which is more important than α0 in practice. For the sake of simplicity, we, respectively, obtained the mean version for salt-and-pepper and random-valued impulse noise by averaging the formulas in Tables 1 and 2, i.e.,(46)h=3.6658p1+0.0156andh=3.6307p2+0.0645,where p1 and p2 denote the noise levels of salt-and-pepper and random-valued impulse noise, respectively. The other parameter values in our numerical experiments were empirically set as ϵ=1.0×10−1, α=2.0×103, δ=1, η=δ/4α. λ was adaptively calculated using the automatic parameter selection technique introduced in Section 5. The maximal iterations Imax and Kmax in Algorithms 1 and 2 were set as 5 for all examples since the algorithms could attain a steady state. We determined the iterative scheme in Algorithm 3 by the following stopping criterion:(47)ft+1−ft2ft+12≤1.0×10−3,or the maximal number of iterations Tmax=10 was reached. Here, ft denotes the recovered image at the tth iteration of the proposed method. All other image restoration methods were implemented with their optimum default parameter setting in numerical experiments. We refer the interested reader to see [29, 36, 37] for more details on the optimal setting of parameters.Figure 3
Determining the optimal relationship between smoothing parameterh and noise level p by means of PSNR under different degradation conditions: (a) Type I blur and salt-and-pepper noise. (b) Type II blur and salt-and-pepper noise. (c) Type I blur and random-valued impulse noise and (d) Type II blur and random-valued impulse noise. The noise level p ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)Figure 4
Determining the optimal relationship between smoothing parameterh and noise level p (both p1 and p2) by means of SSIM under different degradation conditions: (a) Type I blur and salt-and-pepper noise, (b) Type II blur and salt-and-pepper noise, (c) Type I blur and random-valued impulse noise and (d) Type II blur and random-valued impulse noise. The noise level p ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)Table 1
Optimal values ofh for four different images (Barbara, Boat, Cameraman, and Lena) in the case of salt-and-pepper noise. Both PSNR and SSIM were adopted to evaluate the image restoration results.
hType I blurType II blurPSNRSSIMPSNRSSIMBarbarah=4.29p1−0.23h=3.88p1−0.14h=3.42p1+0.15h=3.71p1+0.01Boath=3.37p1+0.15h=3.24p1+0.22h=4.76p1−0.39h=3.57p1+0.19Cameramanh=3.82p1−0.07h=3.35p1+0.18h=4.04p1−0.11h=3.84p1−0.08Lenah=3.51p1+0.10h=3.73p1+0.02h=2.49p1+0.31h=3.73p1−0.06Table 2
Optimal values ofh for four different images (Barbara, Boat, Cameraman and Lena) in the case of random-valued impulse noise. Both PSNR and SSIM were adopted to evaluate the image restoration results.
hType I blurType II blurPSNRSSIMPSNRSSIMBarbarah=3.85p2−0.01h=4.49p2−0.15h=3.46p2+0.10h=3.51p2+0.04Boath=4.05p2−0.14h=3.55p2+0.13h=3.68p2+0.11h=3.11p2+0.22Cameramanh=2.99p2+0.29h=3.27p2+0.06h=3.70p2+0.07h=4.47p2−0.18Lenah=3.57p2+0.01h=3.50p2+0.16h=3.69p2+0.17h=3.19p2+0.16


### 6.3. Quantitative Comparison of Restoration Performance

This section is devoted to compare SA-NLTV-L1 with some current state-of-the-art image restoration techniques. To evaluate the stability of our restoration method, four different images were corrupted by different Gaussian blur kernels and salt-and-pepper/random-valued impulse noise with differentp ranged from 10% to 50% with 10% in step. To evaluate the image restoration performance, both PSNR and SSIM metrics were used simultaneously. Figures 5 and 6 depict the quantitative results. Essentially, both NLTV-L1 and SA-NLTV-L1 are nonlocal variational methods which take full advantage of the natural redundancy of patterns within an image. The nonlocal property is reliably capable of noise suppression and detail preservation during image restoration. Thus both NLTV-L1 and SA-NLTV-L1 can yield the most satisfactory restoration results compared with other local variational methods (i.e., TV-L1, SA-TV-L1, and SA-TV1,2-L1). Among these local methods, SA-TV1,2-L1 outperforms other local methods under consideration. In particular, TV1,2 can suppress the staircase-like artifacts generated by TV in recovered images. Compared with TV-L1, the image quality can be further improved using the spatially adapted regularization parameter for SA-TV-L1. The superior performance of our SA-NLTV-L1 also benefits from the spatially adapted regularization parameter. As shown in Figures 5 and 6, the advantage of SA-NLTV-L1 over other competing methods increases with higher noise level p in terms of both PSNR and SSIM.Figure 5
Quantitative evaluation (in terms of PSNR) of various restoration methods for different images under different degradation conditions: (a) Barbara (type I blur and salt-and-pepper noise), (b) Boat (type I blur and random-valued impulse noise), (c) Cameraman (type II blur and salt-and-pepper noise) and (d) Lena (type II blur and random-valued impulse noise). The noise level (bothp1 and p2) ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)Figure 6
Quantitative evaluation (in terms of SSIM) of various restoration methods for different images under different degradation conditions: (a) Barbara (type I blur and salt-and-pepper noise), (b) Boat (type I blur and random-valued impulse noise), (c) Cameraman (type II blur and salt-and-pepper noise) and (d) Lena (type II blur and random-valued impulse noise). The noise level (bothp1 and p2) ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)


### 6.4. Qualitative Visual Quality Assessment

#### 6.4.1. Influence of Nonlocal Regularizer

Qualitative visual quality comparison plays a crucial role in evaluating the performance of an image restoration method. In particular, the recovered image should contain important geometrical structures, few (or ideally) no visible artifacts. In this subsection, we will investigate the influence of nonlocal regularizer on image restoration performance. Both “Barbara” and “Boat” were selected to evaluate the visual quality of restoration results. As shown in Figures7 and 8, these two images were corrupted by type II blur and different levels of salt-and-pepper/random-valued impulse noise ranged from 20% to 50% with 10% in step. The smoothing parameter h for both NLTV-L1 and SA-NLTV-L1 was determined by the linear approximation formula introduced in Section 6.2. Some of the visual results in Figures 7 and 8 show that both TV-L1 and SA-TV-L1 suffer from some disadvantages such as undesirable staircase-like artifacts and blocky-like structures. This is mainly due to the fact that TV favors piecewise constant solutions. SA-TV1,2-L1 could further enhance the image quality owing to the TV1,2 regularizer and spatially adapted regularization parameter. However, the recovered results easily produce blurred-edge effects which significantly degrade the visual quality. In contrast, NLTV-L1 could achieve the most satisfactory visual quality in terms of noise/blur removal and detail preserving. The advantage of NLTV-L1 was further confirmed by the local magnification views of regions of interest in recovered images shown in Figures 9 and 10. NLTV-L1 performed well on these two images, and produced more similar geometric structures to those of the original version. The superior performance of NLTV-L1 benefits from the self-similarity properties within an image. The qualitative analysis in this subsection could be confirmed by the quantitative results shown in Figures 5 and 6.Figure 7
First column: image “Barbara” corrupted by type II blur and different levelsp1 of salt-and-pepper noise ranged from 20% to 50%. (a) Degraded. From second to fifth column: restoration results yielded by (b) TV-L1, (c) SA-TV-L1, (d) SA-TV1,2-L1 and (e) NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)Figure 8
First column: image “Boat” corrupted by type II blur and different levelsp2 of random-valued impulse noise ranged from 20% to 50%. (a) Degraded. From second to fifth column: restoration results yielded by (b) TV-L1, (c) SA-TV-L1, (d) SA-TV1,2-L1, and (e) NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)Figure 9
Local magnification views of the restoration results for image “Barbara” corrupted by type II blur and salt-and-pepper noise withp1=50%. From top-left to bottom-right: (a) original image, (b) degraded image and recovered images generated by (c) TV-L1, (d) SA-TV-L1, (e) SA-TV1,2-L1, and (f) NLTV-L1, respectively.
(a)(b)(c)(d)(e)(f)Figure 10
Local magnification views of the restoration results for image “Boat” corrupted by type II blur and random-valued impulse noise withp2=50%. From top-left to bottom-right: (a) original image, (b) degraded image and recovered images generated by (c) TV-L1, (d) SA-TV-L1, (e) SA-TV1,2-L1 and (f) NLTV-L1, respectively.
(a)(b)(c)(d)(e)(f)


#### 6.4.2. Influence of Spatially Adapted Regularization Parameter Selection

As shown in Figures5 and 6, SA-NLTV-L1 generates better quantitative results than NLTV-L1 in terms of PSNR and SSIM under different image degradation conditions. In this subsection, we will investigate the influence of spatially adapted regularization parameter selection on qualitative visual appearance. In Figure 11, we first assess the image restoration performance for image “Cameraman” corrupted by type I blur and salt-and-pepper noise with p1=50%. In the case of small type I blur, both NLTV-L1 and SA-NLTV-L1 were capable of maintaining a good balance between noise/blur removal and detail preserving. To facilitate the qualitative visual evaluation, we also presented the absolute difference images between original and degraded/recovered images shown in the bottom panel of Figure 11. To achieve high-quality image restoration, the corresponding absolute difference images should contain as little structural details as possible. As observed in Figure 11, the geometrical structures were more noticeable in absolute difference images yielded by NLTV-L1. In contrast, our SA-NLTV-L1 generated smaller absolute differences due to the spatially adapted regularization parameter shown in Figure 12. As can be observed, the values of the adaptive parameter λ are illustrated in a gray scale. Light gray regions refer to larger λ, whereas dark gray belongs to regions where λ is smaller. From the figure, we find that λ Figure 13 becomes large in the regions of the camera and tripod which contain significative structural information. Note that larger regularization parameters can perform well in fine detail preservation, whereas smaller ones in homogeneous regions demonstrate the results of artifact reduction.Figure 11
Example restoration results for image “Cameraman” corrupted by type I blur and salt-and-pepper noise withp1=50%. Top: degraded image and recovered versions yielded by NLTV-L1 and SA-NLTV-L1, respectively. Bottom: absolute difference images of the degraded and recovered versions shown in top figures. As can be observed, SA-NLTV-L1 shows better performance in both noise/blur removal and detail preserving. (a) Degraded image. (b) NLTV-L1. (c) SA-NLTV-L1.Figure 12
The spatially adapted regularization parametersλ for SA-NLTV-L1 in the case of image “Cameraman” corrupted by type I blur and random-valued impulse noise with p1=50%. Larger λ tends to perform well in preservation of significative structures, whereas smaller ones can demonstrate the results of undesirable artifacts suppression.Figure 13
First row: image “Lena” corrupted by type II blur and different levelsp2 of random-valued impulse noise ranged from (a–d) 20% to 50%. From second to third row: restoration results yielded by (e–h) NLTV-L1 and (i–l) SA-NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)In the case of large type II blur and random-valued impulse noise withp2 ranged from 20% to 50%, we will assess the restoration performance for image “Lena” under different degradation conditions. The restoration results and their associated magnified views are displayed in Figures 13 and 14, respectively. Similar effects as those obtained in the preceding experiments can be found in this example. Both NLTV-L1 and SA-NLTV-L1 were capable of strong impulse noise removal. However, due to the large type II blur, restoration results generated by NLTV-L1 suffered from the ripple artifacts in homogeneous regions (e.g., face and bare shoulder). In contrast, SA-NLTV-L1 could effectively remove the ripple artifacts and preserved the meaningful geometrical structures. The improvements of the recovered images are due to the spatially adapted regularization parameter λ shown in Figure 15. It can be observed that λ is large in the regions of hair and is small in the homogeneous regions. Thus, our SA-NLTV-L1 could establish a proper balance between detail preservation and artifact reduction and generate better visual appearance. In conclusion, the superior performance of the proposed method benefits from the nonlocal regularizer and spatially adapted regularization parameter selection.Figure 14
Local magnification views of the degraded and recovered images show in. As can be observed, SA-NLTV-L1 yields better recovered geometrical structures compared to NLTV-L1 in the local regions of interest. (a–d)20% to 50%. (e–h) NLTV-L1. (i–l) SA-NLTV-L1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 15
The spatially adapted regularization parameterλ for SA-NLTV-L1 in the case of image “Lena” corrupted by type II blur and different levels p2 of random-valued impulse noise ranged from 20% to 50%. Larger λ tends to perform well in preservation of significative structures; whereas smaller ones can demonstrate the results of artifact reduction. (a)20%. (b) 30%. (c) 40%. (d) 50%.
(a)(b)(c)(d)


## 6.1. Comparison with Other Image Restoration Techniques

The proposed method will be compared with five image restoration techniques as follows:(1)
TV-L1:TV Regularized L1 Variational Model [29]. TV-L1 was proposed by combining the L1 data-fidelity term with TV regularizer. In particular, this model is typically capable of preserving edge features without amplifying the impulse noise. In this paper, TV-L1 was solved using an alternating minimization algorithm.(2)
NLTV-L1:NLTV Regularized L1 Variational Model. This model takes full advantage of the natural redundancy of patterns within an image. Thus it is capable of suppressing the staircase-like artifacts often yielded by TV-L1 in recovered image. NLTV-L1 was solved using the alternating optimization algorithm given in Algorithm 2.(3)
SA-TV-L1:Spatially Adapted TV Regularized L1 Variational Model [36]. This method is an improved version of the TV-L1 by introducing an automated method for selection of spatially adapted regularization parameter. For the sake of comparison, the corresponding TV-L1 model was also solved using the alternating minimization algorithm proposed in [29].(4)
SA-TV1,2-L1: Spatially Adapted TV1,2Regularized L1 Variational Model [37]. Compared with SA-TV-L1, Liu et al. [37] proposed to replace TV by TV1,2 to further improve the restoration performance. In this paper, the alternating direction method of multipliers (ADMM) was adopted to efficiently solve the resulting optimization problem.The first two variational methods (i.e., TV-L1 and NLTV-L1) were implemented by using only constant regularization parameter. In contrast, SA-TV-L1 and SA-TV1,2-L1 take full advantage of the spatially adapted regularization parameter. However, these three variational methods were proposed based on local regularizers, which could constrain the further performance improvement in image restoration. In contrast, the proposed SA-NLTV-L1 takes into consideration both nonlocal regularizer and spatially adapted regularization parameters. These two properties are capable of enhancing the performance of our image restoration method.


## 6.2. Optimal Parameter Selection

Nonlocal variational image restoration techniques are dependent on several parameters, i.e.,h, RS, and RN. In the current literature, many efforts have been made to investigate the influence of these parameters on image restoration [27, 57, 58]. To the best of our knowledge, no research has been conducted on the NLTV-L1 model for restoring blurred images with impulse noise. It is generally thought that the patch size RN×RN is less important (common choices are between 5 and 9) [58]. Without loss of generality, we set the patch size to be 7×7 throughout the rest of this paper. To keep a good balance between computational complexity and restoration performance, we only select the 10 best and 4 nearest neighbors for each pixel in the searching window ΩS [27]. The searching window size RS×RS is empirically set as 11×11 as conducted in [27]. The smoothing parameter h is also an important parameter that we will study more carefully. To estimate the optimal parameter h, manual experiments are implemented on the eight different grayscale images, shown in Figure 1. In our experiments, the images are corrupted by two different blur kernels and different levels p (i.e., p1 and p2) of impulse noise ranged from 10% to 70% with 10% in step.In Figure3 and Figure 4, the optimal value of h for different images was determined using the image quality measures in terms of PSNR and SSIM. It can be observed that the optimum h gradually becomes larger as noise level p increases. The relationship between h and p can be expressed as a linear approximation formula, i.e., h=α1p+α0. The fitting coefficients α1 and α0 can be estimated empirically from Figures 3 and 4. To illustrate, the approximation results are summarized in Tables 1 and 2. It can be observed that most linear approximation formulas have the similar coefficient α1 which is more important than α0 in practice. For the sake of simplicity, we, respectively, obtained the mean version for salt-and-pepper and random-valued impulse noise by averaging the formulas in Tables 1 and 2, i.e.,(46)h=3.6658p1+0.0156andh=3.6307p2+0.0645,where p1 and p2 denote the noise levels of salt-and-pepper and random-valued impulse noise, respectively. The other parameter values in our numerical experiments were empirically set as ϵ=1.0×10−1, α=2.0×103, δ=1, η=δ/4α. λ was adaptively calculated using the automatic parameter selection technique introduced in Section 5. The maximal iterations Imax and Kmax in Algorithms 1 and 2 were set as 5 for all examples since the algorithms could attain a steady state. We determined the iterative scheme in Algorithm 3 by the following stopping criterion:(47)ft+1−ft2ft+12≤1.0×10−3,or the maximal number of iterations Tmax=10 was reached. Here, ft denotes the recovered image at the tth iteration of the proposed method. All other image restoration methods were implemented with their optimum default parameter setting in numerical experiments. We refer the interested reader to see [29, 36, 37] for more details on the optimal setting of parameters.Figure 3
Determining the optimal relationship between smoothing parameterh and noise level p by means of PSNR under different degradation conditions: (a) Type I blur and salt-and-pepper noise. (b) Type II blur and salt-and-pepper noise. (c) Type I blur and random-valued impulse noise and (d) Type II blur and random-valued impulse noise. The noise level p ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)Figure 4
Determining the optimal relationship between smoothing parameterh and noise level p (both p1 and p2) by means of SSIM under different degradation conditions: (a) Type I blur and salt-and-pepper noise, (b) Type II blur and salt-and-pepper noise, (c) Type I blur and random-valued impulse noise and (d) Type II blur and random-valued impulse noise. The noise level p ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)Table 1
Optimal values ofh for four different images (Barbara, Boat, Cameraman, and Lena) in the case of salt-and-pepper noise. Both PSNR and SSIM were adopted to evaluate the image restoration results.
hType I blurType II blurPSNRSSIMPSNRSSIMBarbarah=4.29p1−0.23h=3.88p1−0.14h=3.42p1+0.15h=3.71p1+0.01Boath=3.37p1+0.15h=3.24p1+0.22h=4.76p1−0.39h=3.57p1+0.19Cameramanh=3.82p1−0.07h=3.35p1+0.18h=4.04p1−0.11h=3.84p1−0.08Lenah=3.51p1+0.10h=3.73p1+0.02h=2.49p1+0.31h=3.73p1−0.06Table 2
Optimal values ofh for four different images (Barbara, Boat, Cameraman and Lena) in the case of random-valued impulse noise. Both PSNR and SSIM were adopted to evaluate the image restoration results.
hType I blurType II blurPSNRSSIMPSNRSSIMBarbarah=3.85p2−0.01h=4.49p2−0.15h=3.46p2+0.10h=3.51p2+0.04Boath=4.05p2−0.14h=3.55p2+0.13h=3.68p2+0.11h=3.11p2+0.22Cameramanh=2.99p2+0.29h=3.27p2+0.06h=3.70p2+0.07h=4.47p2−0.18Lenah=3.57p2+0.01h=3.50p2+0.16h=3.69p2+0.17h=3.19p2+0.16


## 6.3. Quantitative Comparison of Restoration Performance

This section is devoted to compare SA-NLTV-L1 with some current state-of-the-art image restoration techniques. To evaluate the stability of our restoration method, four different images were corrupted by different Gaussian blur kernels and salt-and-pepper/random-valued impulse noise with differentp ranged from 10% to 50% with 10% in step. To evaluate the image restoration performance, both PSNR and SSIM metrics were used simultaneously. Figures 5 and 6 depict the quantitative results. Essentially, both NLTV-L1 and SA-NLTV-L1 are nonlocal variational methods which take full advantage of the natural redundancy of patterns within an image. The nonlocal property is reliably capable of noise suppression and detail preservation during image restoration. Thus both NLTV-L1 and SA-NLTV-L1 can yield the most satisfactory restoration results compared with other local variational methods (i.e., TV-L1, SA-TV-L1, and SA-TV1,2-L1). Among these local methods, SA-TV1,2-L1 outperforms other local methods under consideration. In particular, TV1,2 can suppress the staircase-like artifacts generated by TV in recovered images. Compared with TV-L1, the image quality can be further improved using the spatially adapted regularization parameter for SA-TV-L1. The superior performance of our SA-NLTV-L1 also benefits from the spatially adapted regularization parameter. As shown in Figures 5 and 6, the advantage of SA-NLTV-L1 over other competing methods increases with higher noise level p in terms of both PSNR and SSIM.Figure 5
Quantitative evaluation (in terms of PSNR) of various restoration methods for different images under different degradation conditions: (a) Barbara (type I blur and salt-and-pepper noise), (b) Boat (type I blur and random-valued impulse noise), (c) Cameraman (type II blur and salt-and-pepper noise) and (d) Lena (type II blur and random-valued impulse noise). The noise level (bothp1 and p2) ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)Figure 6
Quantitative evaluation (in terms of SSIM) of various restoration methods for different images under different degradation conditions: (a) Barbara (type I blur and salt-and-pepper noise), (b) Boat (type I blur and random-valued impulse noise), (c) Cameraman (type II blur and salt-and-pepper noise) and (d) Lena (type II blur and random-valued impulse noise). The noise level (bothp1 and p2) ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)


## 6.4. Qualitative Visual Quality Assessment

### 6.4.1. Influence of Nonlocal Regularizer

Qualitative visual quality comparison plays a crucial role in evaluating the performance of an image restoration method. In particular, the recovered image should contain important geometrical structures, few (or ideally) no visible artifacts. In this subsection, we will investigate the influence of nonlocal regularizer on image restoration performance. Both “Barbara” and “Boat” were selected to evaluate the visual quality of restoration results. As shown in Figures7 and 8, these two images were corrupted by type II blur and different levels of salt-and-pepper/random-valued impulse noise ranged from 20% to 50% with 10% in step. The smoothing parameter h for both NLTV-L1 and SA-NLTV-L1 was determined by the linear approximation formula introduced in Section 6.2. Some of the visual results in Figures 7 and 8 show that both TV-L1 and SA-TV-L1 suffer from some disadvantages such as undesirable staircase-like artifacts and blocky-like structures. This is mainly due to the fact that TV favors piecewise constant solutions. SA-TV1,2-L1 could further enhance the image quality owing to the TV1,2 regularizer and spatially adapted regularization parameter. However, the recovered results easily produce blurred-edge effects which significantly degrade the visual quality. In contrast, NLTV-L1 could achieve the most satisfactory visual quality in terms of noise/blur removal and detail preserving. The advantage of NLTV-L1 was further confirmed by the local magnification views of regions of interest in recovered images shown in Figures 9 and 10. NLTV-L1 performed well on these two images, and produced more similar geometric structures to those of the original version. The superior performance of NLTV-L1 benefits from the self-similarity properties within an image. The qualitative analysis in this subsection could be confirmed by the quantitative results shown in Figures 5 and 6.Figure 7
First column: image “Barbara” corrupted by type II blur and different levelsp1 of salt-and-pepper noise ranged from 20% to 50%. (a) Degraded. From second to fifth column: restoration results yielded by (b) TV-L1, (c) SA-TV-L1, (d) SA-TV1,2-L1 and (e) NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)Figure 8
First column: image “Boat” corrupted by type II blur and different levelsp2 of random-valued impulse noise ranged from 20% to 50%. (a) Degraded. From second to fifth column: restoration results yielded by (b) TV-L1, (c) SA-TV-L1, (d) SA-TV1,2-L1, and (e) NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)Figure 9
Local magnification views of the restoration results for image “Barbara” corrupted by type II blur and salt-and-pepper noise withp1=50%. From top-left to bottom-right: (a) original image, (b) degraded image and recovered images generated by (c) TV-L1, (d) SA-TV-L1, (e) SA-TV1,2-L1, and (f) NLTV-L1, respectively.
(a)(b)(c)(d)(e)(f)Figure 10
Local magnification views of the restoration results for image “Boat” corrupted by type II blur and random-valued impulse noise withp2=50%. From top-left to bottom-right: (a) original image, (b) degraded image and recovered images generated by (c) TV-L1, (d) SA-TV-L1, (e) SA-TV1,2-L1 and (f) NLTV-L1, respectively.
(a)(b)(c)(d)(e)(f)


### 6.4.2. Influence of Spatially Adapted Regularization Parameter Selection

As shown in Figures5 and 6, SA-NLTV-L1 generates better quantitative results than NLTV-L1 in terms of PSNR and SSIM under different image degradation conditions. In this subsection, we will investigate the influence of spatially adapted regularization parameter selection on qualitative visual appearance. In Figure 11, we first assess the image restoration performance for image “Cameraman” corrupted by type I blur and salt-and-pepper noise with p1=50%. In the case of small type I blur, both NLTV-L1 and SA-NLTV-L1 were capable of maintaining a good balance between noise/blur removal and detail preserving. To facilitate the qualitative visual evaluation, we also presented the absolute difference images between original and degraded/recovered images shown in the bottom panel of Figure 11. To achieve high-quality image restoration, the corresponding absolute difference images should contain as little structural details as possible. As observed in Figure 11, the geometrical structures were more noticeable in absolute difference images yielded by NLTV-L1. In contrast, our SA-NLTV-L1 generated smaller absolute differences due to the spatially adapted regularization parameter shown in Figure 12. As can be observed, the values of the adaptive parameter λ are illustrated in a gray scale. Light gray regions refer to larger λ, whereas dark gray belongs to regions where λ is smaller. From the figure, we find that λ Figure 13 becomes large in the regions of the camera and tripod which contain significative structural information. Note that larger regularization parameters can perform well in fine detail preservation, whereas smaller ones in homogeneous regions demonstrate the results of artifact reduction.Figure 11
Example restoration results for image “Cameraman” corrupted by type I blur and salt-and-pepper noise withp1=50%. Top: degraded image and recovered versions yielded by NLTV-L1 and SA-NLTV-L1, respectively. Bottom: absolute difference images of the degraded and recovered versions shown in top figures. As can be observed, SA-NLTV-L1 shows better performance in both noise/blur removal and detail preserving. (a) Degraded image. (b) NLTV-L1. (c) SA-NLTV-L1.Figure 12
The spatially adapted regularization parametersλ for SA-NLTV-L1 in the case of image “Cameraman” corrupted by type I blur and random-valued impulse noise with p1=50%. Larger λ tends to perform well in preservation of significative structures, whereas smaller ones can demonstrate the results of undesirable artifacts suppression.Figure 13
First row: image “Lena” corrupted by type II blur and different levelsp2 of random-valued impulse noise ranged from (a–d) 20% to 50%. From second to third row: restoration results yielded by (e–h) NLTV-L1 and (i–l) SA-NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)In the case of large type II blur and random-valued impulse noise withp2 ranged from 20% to 50%, we will assess the restoration performance for image “Lena” under different degradation conditions. The restoration results and their associated magnified views are displayed in Figures 13 and 14, respectively. Similar effects as those obtained in the preceding experiments can be found in this example. Both NLTV-L1 and SA-NLTV-L1 were capable of strong impulse noise removal. However, due to the large type II blur, restoration results generated by NLTV-L1 suffered from the ripple artifacts in homogeneous regions (e.g., face and bare shoulder). In contrast, SA-NLTV-L1 could effectively remove the ripple artifacts and preserved the meaningful geometrical structures. The improvements of the recovered images are due to the spatially adapted regularization parameter λ shown in Figure 15. It can be observed that λ is large in the regions of hair and is small in the homogeneous regions. Thus, our SA-NLTV-L1 could establish a proper balance between detail preservation and artifact reduction and generate better visual appearance. In conclusion, the superior performance of the proposed method benefits from the nonlocal regularizer and spatially adapted regularization parameter selection.Figure 14
Local magnification views of the degraded and recovered images show in. As can be observed, SA-NLTV-L1 yields better recovered geometrical structures compared to NLTV-L1 in the local regions of interest. (a–d)20% to 50%. (e–h) NLTV-L1. (i–l) SA-NLTV-L1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 15
The spatially adapted regularization parameterλ for SA-NLTV-L1 in the case of image “Lena” corrupted by type II blur and different levels p2 of random-valued impulse noise ranged from 20% to 50%. Larger λ tends to perform well in preservation of significative structures; whereas smaller ones can demonstrate the results of artifact reduction. (a)20%. (b) 30%. (c) 40%. (d) 50%.
(a)(b)(c)(d)


## 6.4.1. Influence of Nonlocal Regularizer

Qualitative visual quality comparison plays a crucial role in evaluating the performance of an image restoration method. In particular, the recovered image should contain important geometrical structures, few (or ideally) no visible artifacts. In this subsection, we will investigate the influence of nonlocal regularizer on image restoration performance. Both “Barbara” and “Boat” were selected to evaluate the visual quality of restoration results. As shown in Figures7 and 8, these two images were corrupted by type II blur and different levels of salt-and-pepper/random-valued impulse noise ranged from 20% to 50% with 10% in step. The smoothing parameter h for both NLTV-L1 and SA-NLTV-L1 was determined by the linear approximation formula introduced in Section 6.2. Some of the visual results in Figures 7 and 8 show that both TV-L1 and SA-TV-L1 suffer from some disadvantages such as undesirable staircase-like artifacts and blocky-like structures. This is mainly due to the fact that TV favors piecewise constant solutions. SA-TV1,2-L1 could further enhance the image quality owing to the TV1,2 regularizer and spatially adapted regularization parameter. However, the recovered results easily produce blurred-edge effects which significantly degrade the visual quality. In contrast, NLTV-L1 could achieve the most satisfactory visual quality in terms of noise/blur removal and detail preserving. The advantage of NLTV-L1 was further confirmed by the local magnification views of regions of interest in recovered images shown in Figures 9 and 10. NLTV-L1 performed well on these two images, and produced more similar geometric structures to those of the original version. The superior performance of NLTV-L1 benefits from the self-similarity properties within an image. The qualitative analysis in this subsection could be confirmed by the quantitative results shown in Figures 5 and 6.Figure 7
First column: image “Barbara” corrupted by type II blur and different levelsp1 of salt-and-pepper noise ranged from 20% to 50%. (a) Degraded. From second to fifth column: restoration results yielded by (b) TV-L1, (c) SA-TV-L1, (d) SA-TV1,2-L1 and (e) NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)Figure 8
First column: image “Boat” corrupted by type II blur and different levelsp2 of random-valued impulse noise ranged from 20% to 50%. (a) Degraded. From second to fifth column: restoration results yielded by (b) TV-L1, (c) SA-TV-L1, (d) SA-TV1,2-L1, and (e) NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)Figure 9
Local magnification views of the restoration results for image “Barbara” corrupted by type II blur and salt-and-pepper noise withp1=50%. From top-left to bottom-right: (a) original image, (b) degraded image and recovered images generated by (c) TV-L1, (d) SA-TV-L1, (e) SA-TV1,2-L1, and (f) NLTV-L1, respectively.
(a)(b)(c)(d)(e)(f)Figure 10
Local magnification views of the restoration results for image “Boat” corrupted by type II blur and random-valued impulse noise withp2=50%. From top-left to bottom-right: (a) original image, (b) degraded image and recovered images generated by (c) TV-L1, (d) SA-TV-L1, (e) SA-TV1,2-L1 and (f) NLTV-L1, respectively.
(a)(b)(c)(d)(e)(f)


## 6.4.2. Influence of Spatially Adapted Regularization Parameter Selection

As shown in Figures5 and 6, SA-NLTV-L1 generates better quantitative results than NLTV-L1 in terms of PSNR and SSIM under different image degradation conditions. In this subsection, we will investigate the influence of spatially adapted regularization parameter selection on qualitative visual appearance. In Figure 11, we first assess the image restoration performance for image “Cameraman” corrupted by type I blur and salt-and-pepper noise with p1=50%. In the case of small type I blur, both NLTV-L1 and SA-NLTV-L1 were capable of maintaining a good balance between noise/blur removal and detail preserving. To facilitate the qualitative visual evaluation, we also presented the absolute difference images between original and degraded/recovered images shown in the bottom panel of Figure 11. To achieve high-quality image restoration, the corresponding absolute difference images should contain as little structural details as possible. As observed in Figure 11, the geometrical structures were more noticeable in absolute difference images yielded by NLTV-L1. In contrast, our SA-NLTV-L1 generated smaller absolute differences due to the spatially adapted regularization parameter shown in Figure 12. As can be observed, the values of the adaptive parameter λ are illustrated in a gray scale. Light gray regions refer to larger λ, whereas dark gray belongs to regions where λ is smaller. From the figure, we find that λ Figure 13 becomes large in the regions of the camera and tripod which contain significative structural information. Note that larger regularization parameters can perform well in fine detail preservation, whereas smaller ones in homogeneous regions demonstrate the results of artifact reduction.Figure 11
Example restoration results for image “Cameraman” corrupted by type I blur and salt-and-pepper noise withp1=50%. Top: degraded image and recovered versions yielded by NLTV-L1 and SA-NLTV-L1, respectively. Bottom: absolute difference images of the degraded and recovered versions shown in top figures. As can be observed, SA-NLTV-L1 shows better performance in both noise/blur removal and detail preserving. (a) Degraded image. (b) NLTV-L1. (c) SA-NLTV-L1.Figure 12
The spatially adapted regularization parametersλ for SA-NLTV-L1 in the case of image “Cameraman” corrupted by type I blur and random-valued impulse noise with p1=50%. Larger λ tends to perform well in preservation of significative structures, whereas smaller ones can demonstrate the results of undesirable artifacts suppression.Figure 13
First row: image “Lena” corrupted by type II blur and different levelsp2 of random-valued impulse noise ranged from (a–d) 20% to 50%. From second to third row: restoration results yielded by (e–h) NLTV-L1 and (i–l) SA-NLTV-L1, respectively (the images are best viewed in full-screen mode).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)In the case of large type II blur and random-valued impulse noise withp2 ranged from 20% to 50%, we will assess the restoration performance for image “Lena” under different degradation conditions. The restoration results and their associated magnified views are displayed in Figures 13 and 14, respectively. Similar effects as those obtained in the preceding experiments can be found in this example. Both NLTV-L1 and SA-NLTV-L1 were capable of strong impulse noise removal. However, due to the large type II blur, restoration results generated by NLTV-L1 suffered from the ripple artifacts in homogeneous regions (e.g., face and bare shoulder). In contrast, SA-NLTV-L1 could effectively remove the ripple artifacts and preserved the meaningful geometrical structures. The improvements of the recovered images are due to the spatially adapted regularization parameter λ shown in Figure 15. It can be observed that λ is large in the regions of hair and is small in the homogeneous regions. Thus, our SA-NLTV-L1 could establish a proper balance between detail preservation and artifact reduction and generate better visual appearance. In conclusion, the superior performance of the proposed method benefits from the nonlocal regularizer and spatially adapted regularization parameter selection.Figure 14
Local magnification views of the degraded and recovered images show in. As can be observed, SA-NLTV-L1 yields better recovered geometrical structures compared to NLTV-L1 in the local regions of interest. (a–d)20% to 50%. (e–h) NLTV-L1. (i–l) SA-NLTV-L1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 15
The spatially adapted regularization parameterλ for SA-NLTV-L1 in the case of image “Lena” corrupted by type II blur and different levels p2 of random-valued impulse noise ranged from 20% to 50%. Larger λ tends to perform well in preservation of significative structures; whereas smaller ones can demonstrate the results of artifact reduction. (a)20%. (b) 30%. (c) 40%. (d) 50%.
(a)(b)(c)(d)


## 7. Conclusions and Future Work

In this paper, we have studied the variational method for restoring blurred images with impulse noise. Based on the good features of the nonlocal regularizer and the spatially adapted parameter selection scheme proposed by Hintermüller and Rincon–Camacho in [36], we proposed a spatially adapted NLTV regularized L1 variational image restoration method. This method is capable of restoring degraded images while preserving edges and fine details. The major contributions of this paper were threefold. First, a NLTV regularizer was incorporated into the MAP-based image restoration framework. The main benefit of this nonlocal regularizer is that it can take full advantage of the natural redundancy of patterns within an image. Therefore, edges and fine details can be preserved during image restoration. Second, to guarantee solution stability, an alternating optimization algorithm was proposed to efficiently solve the proposed NLTV-L1 model. Third, to future improve image visual quality, we introduced an automated method for selection of spatially adapted regularization parameter in terms of both salt-and-pepper and random-valued impulse noise. Numerous experiments preformed on several images have demonstrated the superior performance of our proposed method in terms of both quantitative and visual quality evaluations. However, the proposed method still suffers from several potential limitations in current version. To make the proposed method available in practical applications, the work presented in this paper can be extended in the following directions:(i)
The results for computational CPU time are illustrated in Figure16. It can be observed that both NLTV-L1 and SA-NLTV-L1 suffer from high computational cost under different image degradation conditions. The high computational cost is mainly due to the cost of weights calculation for all pixels during image restoration. To significantly accelerate the image restoration, both NLTV-L1 and SA-NLTV-L1 will be implemented using graphics processing unit (GPU). Over the past few decades, GPU has rapidly evolved into a cost-effective hardware platform for parallel computing [59] and has been successfully implemented to accelerate the NLTV-based motion estimation [60] and NLTGV-based optical flow estimation [61]. There is thus a strong incentive to accelerate both NLTV-L1 and SA-NLTV-L1 for real-time applications on GPU-based parallel computing platform.(ii)
In this paper, both blur kernel and noise level are assumed to be known before image restoration. Therefore, the recovered image can be directly induced from the sufficient information provided by the blur kernel and image noise. However, both blur kernel and noise level are commonly unknown in practical applications. In current literature, a significant amount of effort has been dedicated to impulse noise estimation and reduction [9, 62]. By contrast, blind image restoration, which means that the blur kernel needs to be estimated, is a much more challenging ill-posed problem. Many methods have been proposed for blind image restoration [63], but they often tend to suffer from severe ringing artifacts especially in the presence of image noise. In future research, we will propose to first identify image pixels contaminated by impulse noise, as done in [31] and to simultaneously implement noise reduction and blind image restoration based on the estimated location of noise-free pixels. Due to the spatially adapted regularization parameter selection and nonlocal regularizer insensitive to noise, we believe that our SA-NLTV-L1 can be successfully generalized to blind image restoration in the presence of impulse noise.(iii)
Due to the powerful learning capacity, deep learning has significantly improved the image restoration results [64–66]. It is well known that deep learning is an essentially data-driven method, which does not take into consideration the image degradation process. The restored images sometimes suffer from instabilities and unwanted artifacts under complex imaging conditions. In addition, the variational method is essentially related to the image prior which indicates the intrinsic image knowledge. It is able to robustly restore the degraded images. In future study, there is a great potential to incorporate the deep learning strategy into our nonlocal regularized variational model. The learning module is capable of accurately estimating the important variables, which are difficult to handle in traditional variational methods. Therefore, the imaging performance will be further enhanced, learning to more reliable postprocessing steps, such as object detection, recognition, and tracking in practical applications.Figure 16
Computational CPU time (unit: seconds) of different restoration methods for image “Barbara” under different degradation conditions: (a) Type I blur and salt-and-pepper noise. (b) Type II blur and salt-and-pepper noise. (c) Type I blur and random-valued impulse noise and (d) Type II blur and random-valued impulse noise. The noise levelp ranges from 10% to 50% with 10% in step.
(a)(b)(c)(d)Honestly, our current study has some potential limitations, but it is still worthy of consideration since the proposed SA-NLTV-L1 could guarantee superior restoration performance compared to several state-of-the-art methods. We believe there is a great potential to use the proposed method to recover degraded image in practical applications.



---

*Source: 9419410-2022-06-15.xml*