A robust iterative algorithm for image restoration

We present a new image restoration method by combining iterative VanCittert algorithm with noise reduction modeling. Our approach enables decoupling between deblurring and denoising during the restoration process, so allows any well-established noise reduction operator to be implemented in our model, independent of the VanCittert deblurring operation. Such an approach has led to an analytic expression for error estimation of the restored images in our method as well as simple parameter setting for real applications, both of which are hard to attain in many regularization-based methods. Numerical experiments show that our method can achieve good balance between structure recovery and noise reduction, and perform close to the level of the state of the art method and favorably compared to many other methods.


Introduction
Image restoration aims to compensate for or undo the defects that degrade an image. Degradation can come in many forms such as motion blur, noise, and camera defocus. In optical microscopes, there are predominately two sources for degradation in the imaging systems, blurring and noise, which can be described by the general imaging model where I, J are the ground truth and the corresponding observation, respectively, P is a point spread function (PSF), and N is a noise which is assumed to be independent to the ground truth. The simplest way to estimate the ground truth from the observation is by minimizing the residual, which can lead to the least square solution. Unfortunately, a unbounded noise will be introduced into the solution because the PSF matrix always has small eigenvalues even *Correspondence: W.Lu@hw.ac.uk 2 School of Engineering and Physical Sciences, Heriot Watt University, Edinburgh, UK Full list of author information is available at the end of the article it is invertible. This is not surprising as (2) is well known to be ill-posed [1].
There are now vast literatures to tackle the problem of image restoration. A recent trend is concentrated on a sparse block matching 3-D (BM3D)-based restoration technique. BM3D algorithms are initially developed for collaborative filtering through a non-local modeling of images by collecting similar image patches in 3D arrays [2]. They have recently been incorporated into image restoration for solving regularized inverse problems for image denoising as well as deblurring [3]. Another development based on BM3D is sparse representation for image restoration, where the image is considered to be a combination of a few atomic functions taken from a certain dictionary and can be parameterized and approximated locally or non-locally by these functions [4]. The dictionary is usually considered as an over-complete system in order to better describe all variety of images. There are now many published works on the sparsity-based models and methods [5]. For example, the formulation of IDD-BM3D image modeling in terms of the overcomplete sparse frame representation for image reconstruction has led to impressive restoration performance [6]. This approach allows decoupling between deblurring and denoising by considering the optimization problem as a generalized Nash equilibrium balance of two objective functions. A distinct advantage of this approach is that various denoising algorithms can be selected independently with respect to deblurring algorithms, which have demonstrated better performance than those where deblurring and denoising are jointly performed in many cases. However, for the decoupled algorithms such as [6], [7], and [8], the parameters setting for optimal performance of regularization is usually complicated and reasons for the best setting are often not explained. Another shortage of these methods is the lack of error analysis for the solutions because of the complexity of the regularization factors.
In this paper, we present a new image restoration method based on the inverse operator theory. As we know, the inverse operator theory [9] gives the solution of P −1 J = I +P −1 N for the general imaging model (1), where P −1 is the inverse or pseudo-inverse of the PSF matrix. Due to small eigenvalues of P, P −1 N leads to significant noise amplification to the ground truth. To overcome this problem, we propose a new approach that combines iterative VanCittert algorithm with noise reduction modeling, the latter enables to reliably estimate the gradient in the presence of noise so that the VanCittert iteration can converge to the ground truth even when the observation is noise contaminated. This work has several contributions to the research area of image restoration. Firstly, it extends the inverse operator theory to image restoration in the presence of noise, which offers a different approach to that of the present popular regularization methods. Secondly, our method enables decoupling between deblurring and denoising, so any well-established noise reduction operator can be selected in our model, independent of the VanCittert deblurring operation. Thirdly, our approach allows error analysis of the solutions because the structure recovery and noise amplification in the VanCittert iterations can be separated analytically, which is an advantage over many regularization methods for which errors are difficult to be estimated due to complicated regularization factors. Finally, parameter setting in our method is simple and robust to the performance. There are only two parameters in our method: σ , as the noise reduction strength, and s, the interval between two neighboring denoising operations. We have further developed an automated parameter setting procedure for our method, which has no need to set the parameters manually. The above points have been verified by numerical experiments, which also show that our method performs close to the level of the state of the art method and favorably compared to many other methods.

Methods
Our method is motivated by the iterative VanCittert algorithm, which has a long history as a simple and efficient approach for image restoration. The algorithm is formulated for spatially invariant or variant restoration problems with neglect of noise contribution in (1). Originally, it is a steepest descent method but the solution does not converge if the step parameter is assumed to be real values. To overcome this shortage, an iterative procedure was proposed [10], which converges to the ground truth only if noise in an observation is negligible, where P T is the transpose of P. When an observation comprises noise, VanCittert iteration (3) can be expressed as [9,11,12] where {ζ uv : u = 1, 2 . . . , R, v = 1, 2, . . . , C} and Z uv are the eigenvalues and eigenvectors of P, and R, C are image size, and u, v the indices of image pixels and β the step parameter. The first term involving I describes structure recovery while the second term involving N shows noise amplification, so structures and noise are separated in (4). For a noisy observation, however, small eigenvalues ζ uv can lead to significant noise amplification in the second term of (4) so the inverse problem becomes ill-posed. Therefore, we have to suppress the noise in the second term if the iteration gives any hope to converge to the ground truth. To tackle the above problem, we first introduce a noise reduction operator, , which minimizes the estimation error of a cost function. Letting I be the ground truth, N be a white noise and V = I + N, we define the cost function for the noise reduction operator, 2 2 , where E{·} is the expectation taken over the noise distribution. The error is measured by L 2 norm and averaged over the noise distribution. For the general imaging model (1), we propose our method as is the first term in (4). As seen from (5), the VanCittert iteration (5a) tends to recover the structures by searching a solution along the gradient of (2). However, the gradient is contaminated by noise, leading to noise amplification in the iterative solution. A noise reduction operation (5b) and (5c) is then applied to remove noise and to optimize the gradient for the next iteration.
As for noise reduction of (5b), our method do not expect an ideal operator removing all noise [13]. Instead, it can be any denoising algorithm as long as the operator satisfies the following condition, where σ 2 = Var(N) is the variances of noise N and 1 is the noise reduction factor. The condition (6) implies that remaining noise o(N) has a variance far less than the initial noise N after applying the operator R N on a noisy image. We will show in error analysis below that when the condition (6) is satisfied, the iterative solution of (5) converges to the ground truth with a higher order small noise term, i.e., I k→∞ = I + o(N). A necessary condition for the iterative process (5) to converge is that |1 − β|ζ uv | 2 | must fall within [ 0, 1) for all the eigenvalues of P, which leads to Since most of PSFs act as a low-pass filter, the maximum absolute value of their eigenvalues is about 1. This means that it is easy to set a value for β that satisfies condition (7) and the following condition: For example, we set β = 1 in our experiments.
To implement the method (5), we can apply any wellestablished noise reduction algorithms to combine with the VanCittert iteration, for example, the wavelet domain shrinking filter T SW (V , δ) =Î = Ww, whereÎ is the estimated image, W is a group of wavelet bases and w is a vector of shrinking coefficients depending on the smooth parameter δ [13]. The smooth parameter δ can be determined in a similar form to (5b) by which has a noise shrinkage strength of = (2logRC + 1)(logRC + 1)/RC, where R and C are the image size. For images of modest size, 1 so the wavelet algorithm satisfies (6). Another popular denoising method is the state of the art BM3D method. BM3D improves from wavelet domain shrinking by incorporating the concept of image patches and non-local mean (NLM) [14] into a transformed domain and has shown the highest peak signal-to-noise ratio in its performance compared to the wavelet domain and other algorithms. Moreover, BM3D has simple parameter setting and is easy to use. Mathematically, BM3D can be expressed as where Z is the stacked noisy blocks, T 3D is the transformation from spatial domain to frequency domain with discrete cosine bases, and W wie is the Wiener shrinkage operator and A is the aggregation operator, all defined in [2]. In view of the advantage and excellent performance of BM3D [2], we choose the operator R N = O BM3D in our method (5) for the numerical experiments below. We note that while we follow the same decoupled approach for deblurring and denoising as IDD-BM3D [6], our method has two advantages. Firstly, structure restoration by the VanCittert algorithm in (5) has a simple step parameter of β = 1, while the regularization factors for optimal deblurring in [6] and [7] are much more complex to set. This leads to overall simple parameter setting of our method compared to many regularization methods. Secondly, structure and noise can be separated analytically in (5), which allows us to perform error analysis for a restored image, while error analysis for regularization methods is generally hard to attain due to the complexity of the regularization factors.

Error analysis
Given (4), the noise amplification is separated from structures so VanCittert algorithm allows error analysis for our method. We begin with the following lemma and then give theorem (1).
Theorem 1 For any operator R N satisfying condition (6) and β satisfying (8), the iterative solution of model (5) leads to where I is the noise free solution(ground truth) and o(N) denotes the remaining noise with the variance far less than the variance of the noise N.
be the noise amplification factors for k = 1, 2 . . . . Known from lemma (1), F k is infinite so there must be a minimum number k satisfying when noise N is bounded. Here, we suppose k = 1 for convenience though this number depends the eigenvalues of P. Thus, we start with the iteration solution(4) for k = 1, where I 0 is the initial image and D 1 is the first sum on right side. The first term in (12) is noise-free iterative solution while the second term is noise contribution with the factor F 1 = βζ uv . Since most of the eigenvalues have absolute values small than 1, Var(N 1 ) ∝ Var(N). Then we apply the filter (6) to (12) and have The noise is now reduced by a factor 1 according to (6), i.e., Var(o(N 1 )) ∝ Var(N 1 ) ∝ σ 2 , where σ 2 is the variance of noise N.
The intensity of image I 1 can be rewritten where o( uv )Z uv is used and F 1 is defined at the beginning of the proof.
From (14) and (4), we give the iteration solution for k = 2 where N 2 is the sum of the two noise terms on the right side. By (11) and lemma (1), 1 ≤ F 2 /F 1 < 2, Var(N 2 ) ∝ Var(N) = σ 2 is obtained. Repeating in the operation (13) we have where the variance of o(N 2 ) is far less than that of N 2 , also far less than that of N.
In general, we have By applying R N to I k , we get the iterative image where Var(o(N k )) ∝ Var(N) = σ 2 . Therefore, the iterative solutions converge to the real scene and the noise is controlled to the order of σ 2 in the iterative process.
It can be seen from the above derivation that the separation of structure recovery from noise amplification in the VanCittert expression is the key that enables us to express noise amplification factor in the kth iteration to be F k /F k−1 , which is always between 1 and 2 for all k. This makes it possible to control the noise amplification over a finite number of iterations by a noise reduction operator satisfying (6). Consequently, the iterative solution converge to the ground truth with higher-order infinitesimal noise.

Results and discussion
such as lines, buttons, and textures. These images are the subjects of a recent extensive investigation by an iterative decoupled deblurring BM3D algorithm (IDD-BM3D) [6], which is formulated based on the Nash equilibrium balance of two objective functions undertaking separate denoising and deblurring operations. IDD-BM3D has showed state of the art restoration performance compared to seven other existing methods, which include Fourier-Wavelet regularized deconvolution (ForWaRD) [15], space-variant Gaussian scale mixtures (SV-GCM) [16], shape-adaptive discrete cosine transform(SA-DCT) [17], BM3D deblurring (BM3DDEB) [3], analysis-based sparsity (L0-Abs) [18], adaptive total variation image deblurring by a majorization minimization approach (TVMM) [19], and finally a method based on spatially weighted total variation (CGMK) [20]. We test on the same six scenarios in [6], which have different PSF shapes and blurring strengths as well as noise levels listed in Table 1. Comparisons with all the eight methods are made quantitatively through the measurement of peak signalto-noise ratio (PSNR).
As discussed earlier, we choose BM3D (http://www. cs.tut.fi/~foi/GCF-BM3D) as our noise reduction filter because it combines the transform-domain filter [13] with non-local mean filter [14] and has shown improved performance over the both methods individually on their own. Our method is easy to operate, requiring only two parameters: noise standard deviation, σ , as an input for BM3D denoising and the step interval, s, between two neighboring denoising operations in the iteration (denoising is not necessary for each iteration step for efficient computing). In general, the two parameters depend on the levels of blur and noise in an observation (input  image). We have found that our method can produce good performance in a large area in the two parameter space, showing the robustness of the method against the setting of the two parameters. The solutions converge around 1200 iterations for all scenarios except scenario 3 which requires 10,000 iterations because of severe blur in this scenario. Due to high noise levels in scenarios 4 and 6, BM3D is applied to the observations before our iterative algorithm is implemented. To investigate the effects of different denoising algorithms on the performance of our method, we have also implemented the wavelet domain shrinking filter T SW as an alternative denoising operator in our method. Table 2 shows the results of PSNR for our algorithm, both with BM3D and T SW , and the eight existing methods, the latter are from [6]. From the table, we can conclude that our method with BM3D outperforms the existing seven methods for both images under different scenarios and is not far behind the state-of-the-art IDD-BM3D. As expected, the algorithm with T SW performs not as good as that with BM3D, because the latter is better than the former as a denoising method. We have further investigated the above results by using the structural similarity (SSIM) index matric, which is a method for measuring the structural similarity between two images. The results are shown in Table 3. As seen from the table, L0-Abs, BM3DDEB, and our methods both with BM3D and T SW have all performed better than IDD-BM3D in terms of SSIM for the first five scenarios, although IDD-BM3D gives the highest PSNR values as discussed above. By comparing the results in Tables 2  and 3, it is pleasing to see that our method with BM3D gives very good and balanced performance in terms of both noise reduction and structure preservation. Figures 1  and 2 show the restored images by four methods used. As seen by visual inspection, BM3DDEB produces some obvious artifacts around the edge of cameraman and L0-Abs cannot restore some details of the eye in the image of Lena because of noise. In comparison, IDD-BM3D and our method with BM3D denoise very well and our method shows better recovery of elaborative features in these images.
We have further developed an automated parameter setting procedure for our method. Here, we first estimate the noise standard deviation in the observation and set this value as the denoising threshold, σ thr . During the iteration, the noise level of the image is estimated at each iteration step and when it exceeds σ thr , then BM3D denoising is applied. The procedure is simple to operate, which can be important for real scene applications. We Fig. 1 Deblurring of Cameraman image, scenario 3. From left to right and from top to bottom are presented zoomed fragments of the following images: original, blurred, and noisy, reconstructed by BM3DDEB, L0-Abs, IDD-BM3D, and our method. In our method, the two input parameters used, (σ , s), are (7.5, 550) for this scenario, (7.5, 85) for scenario 1, (7.5, 25) for scenario 2, (7.5, 10) for scenario 4, (7.5, 50) for scenario 5, and (7.5, 5) for scenario 6

Conclusions
In summary, we have developed a new robust iterative method for image restoration in which an iterative cost function is utilized to optimize the gradient in the steepest descent by adaptively adjusting to the current state in the iterative process. We show that the iterative solution converges to the real scene despite noise contamination in an observation, and the restoration error can be controlled by an order of magnitude smaller than the noise level in the observation. Different from the well-established regularization methods, which introduce a penalty to solve the ill-posed problem, we directly apply VanCittert algorithm to minimize residual along gradient for structure restoration and to suppress noise through a noise reduction operator. It turns out to be a denoising problem in an iterative manner, and the noise can be removed judiciously by applying existing noise reduction methods. We have undertaken two numerical experiments to investigate the performance of this method and compare to existing regularization methods. We show that our method performs to the level close to the best of the methods currently available in terms of recovering elaborate structures and reducing noise, and favorably compared to the many other existing methods. Moreover, our method requires simple parameter setting, particularly in the second experiment where a single parameter is estimated from the observation. This could be a great advantage for real-world applications. We note finally that we have only considered additive noise in this paper, for images contaminated by multiplicative noise, some newly developed noise reduction filters, such as [22] may be applied, which can be our future work in this area.

Availability of data and materials
The web links to the sources of the data ( namely, images) used for our experiments and comparisons in this work have been provided in this article.

Authors' contributions
The method development and test are mainly carried out by Liu, with extensive discussion and supervision from Lu. All authors read and approved the final manuscript.