Reverse image filtering with clean and noisy filters

Given an image filter y=f(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{y}}}={{\varvec{f}}}\,({{\varvec{x}}})$$\end{document}, where x\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{x}}}$$\end{document} and y\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{y}}}$$\end{document} are input and output images, respectively, reverse image filtering consists of rendering an approximation to x\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{x}}}$$\end{document} from y\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{y}}}$$\end{document} using the filter f(·)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{f}}}\,(\cdot )$$\end{document} itself as a black box, without knowing the internal structure of the filter. In this paper, we propose to use modified Landweber iterations for reverse image filtering, evaluate the performance of our approach, and present applications to image deblurring and super-resolution. An important advantage of our approach over the existing reverse image filtering methods is high robustness to noise.


Introduction and contributions
Reverse image filtering is a recent and promising research direction in image processing [2,9,20,30]. Given an image filter, the goal of reverse image filtering consists of removing or suppressing of the filter effects using only the filter itself. Furthermore, the filter internal structure is not supposed to be known, the filter is allowed to be used as many times as needed but only as a black box.
The idea of multiple use of the same filter for designing a new filter or improving existing filter properties is not new, of course. It goes back to Tukey's twicing scheme [12,32] and was further developed by Kaiser and Hamming [14]. See also an expository paper [8]. Similar approaches have been recently proposed in [23,25]. Reverse image filtering is a related but different problem, as our aim is to invert the filter, not to enhance filter properties.
The problem setting for reverse image filtering is simple and consists of restoring an input image x from the observed output image y = f (x) where f (·) is an image filter. It can be formulated as an inverse problem x = f −1 (y). We assume that filter f (·) is only available to us as a black-box, i.e., it can be queried, but it cannot be directly inverted, as we don't know its internal structure. Thus, the goal of reverse filtering consists of estimating the filter inverse f −1 (·) using iterations of the filter f (·) itself.
The previous studies on reverse image filtering [2,9,20,30] mostly work with symmetric and clean filters. However, in real-life scenarios, we often deal with non-symmetric filters (e.g., motion blur filters) which can also be contaminated by noise. In Fig. 1, we present a simple experiment with the convolutional filter defined by where h is a motion blur kernel and n is a zero-mean Gaussian noise. The reader can see that both the reverse filtering methods of [30] and [9] fail to achieve a satisfactory reconstruction of original image x: the method of [30] does not cope with non-symmetric filters (the same can be said about the schemes proposed in [20] and [2]) and the method of [9] demonstrates a noise-amplification behavior. In contrast, a Landweber-like iteration scheme proposed in this paper yields a very satisfactory reconstruction of the original image.  (1) is shown. c Blurred and noisy image generated by (1) where n is additive Gaussian noise with zero mean and variance equal to 10 −4 . d The method of Tao et al. [30] is unable to cope with non-symmetric filters.
e The method of Dong et al. [9] amplifies noise. f A phase-corrected Landweber-like iteration scheme (to be discussed later in the paper) yields a satisfactory reconstruction of the original image Despite their seemingly simple nature, the Landweber iteration method and its relatives constitute a family of powerful and popular image restoration tools [3,27].
Our approach to reverse image filtering is based on combining Landweber-like iterations with an iterative convolutionbased approximation of a given filter. We demonstrate advantages of the approach over the existing reverse image filtering methods and consider its applications to image deblurring and super-resolution. We demonstrate that total variation regularization, Nesterov acceleration and some image processing techniques can be incorporated in our reverse image filtering framework.

Previous work
Mathematically we deal with the equation where y is an observed filtered image, f (·) is a given blackbox filter, and the goal is to achieve an accurate estimation of the original image x. At the first glance, solving (2) numerically seems to be a simple task. For example, estimating the gradient (the Jacobian matrix) of f (x) allows us to use standard equation solving methods including Newton's iterations.
The problem here is that even for a very small 100 × 100 image, the Jacobian matrix has size of 10 4 × 10 4 and that makes gradient-based methods computationally infeasible.
This calls for developing and using gradient-free methods for numerically solving (2). Gradient-free methods, also known as zeroth-order optimization methods, are currently a subject of intensive research in connection to their applications in Machine Learning [18]. Below we review gradient-free methods introduced for solving (2) in the context of reverse image filtering.
Zero-order Reverse Filtering. Tao et al. [30] propose the following simple iterative reverse filtering scheme where x 0 = y. The term y − f (y) contains details of y removed by filter f (·). The idea is to iteratively use this residual to approximate the lost details of x and add them back to estimate the original version. The whole procedure can be considered as an iterative version of unsharp masking, a simple and popular image sharpening tool. The authors of [30] consider (3) as a fixed-point map and establish its convergence for linear filters satisfying certain conditions. In practice, (3) shows a very good performance for a number of popular nonlinear image filters. However, as demonstrated in Fig.1, (3) fails to invert non-symmetric filters including linear motion blur filters. Rendition. A generalization of (3) proposed by Milanfar [20] consists of running the following iterative process where γ and τ are user-specified parameters. In theory, (4) has better convergence properties than (3). In practice, however, both methods demonstrate similar performance. First-order Reverse Filtering. Dong et al. [9] propose a gradient-free defiltering scheme which they call an iterative first-order reverse image filtering method, as Newton's iterations are used [9] to derive the method. The method performs its iterations in the frequency domain where point-wise division is used and where, here and everywhere below, the capital letters indicate the Fourier transforms of the images denoted by the corresponding lower case letters. For example, F(x k ) = F{f (x k )} and X k = F{x k }, where F{·} denotes the Fourier transform.
We can interpret (5) as follows. If f (·) is a linear convolution filter then its corresponding frequency response function is given by Thus, in each iteration, (5) approximates f (·) by a linear convolution.
The method of Dong et al. demonstrates excellent performance when dealing with linear filters but is very sensitive to boundary conditions (see examples in [2]) and, in addition, cannot cope with image noise, as shown in Fig.1.
Other reverse image filtering schemes. Finally two reverse filtering methods proposed in [2] are designed to suppress severe filtering effects produced by clean and symmetric image filters and are not appropriate to deal with noisy and/or non-symmetric filters.

Proposed approach
The reverse image filtering approach we pursue in this paper is based on using Landweber-like iterations applied to a quadratic error energy function. Namely let us consider the following energy minimization problem If f (·) is a linear filter, f (x) = Ax, then it seems natural to restore x by applying the standard Landweber iterations [3,15] where A * is the adjoint (Hermitian conjugate) of operator A and τ > 0 is a sufficiently small step-size parameter.
Although, for linear f (·), this iterative process represents a simple gradient descent for quadratic energy (8), (9) and its generalizations and modifications remain popular tools for various linear and nonlinear inverse problems including image deconvolution [27]. If filter f (·) is a linear convolution (6), then, due to Parseval's identity, (8) can be transferred to the frequency domain where the transfer function H is defined by (7). Assume that we deal with a linear convolution filter and H is fixed. Then, in the frequency space, Landweber iterations (9), a special case of gradient descent for (10), are given by where H * denotes the complex conjugate of H and τ is the step size. Now combining (11) with estimating the transfer function H at each iteration leads us to Note that (13) works in the frequency space and, at each step, approximates nonlinear filter f (·) by a linear convolution filter whose transfer function is determined by (12). This simple Landweber-like iterative procedure is summarized in Algorithm 1. Similar to the other reverse filtering methods, the algorithm only requires access to the filter as a black-box (without knowledge of its internal structure).

Algorithm 1 Landweber-like iteration method
One advantage of using the variational approach (8) for reverse image filtering is that (8) can easily be extended by adding additional terms incorporating various image priors. Here we are interested in using the total variation prior [26] which leads to a better preservation of salient image edges. Namely, instead of (8) we consider where γ is a positive parameter controlling the strength of the prior. Our Landweber-like iterations for (14) are given by The results of our numerical experiments presented in Sect. 4 demonstrate that adding the total variation term improves reverse filtering results for deblurring and superresolution tasks.
123 Table 1 Evaluation results in PSNR (dB) for our method compared with the methods of Tao   Italic, bold and bold-italic values indicate the best approach, the second best and the third best respectively

Defiltering
To evaluate the methods developed in this work, we selected 16 popular image filters. The first four, Gaussian filter (GS), motion blur (MT), Laplacian of Gaussian (LOG), and disk filter (DK), are simple linear filters, whose operations on the image can be described by the convolution with a kernel. Except the unsharp masking filter (UMF), the remaining filters are nonlinear. Adaptive Wiener filter (AWF) [19] and weighted median filter (WMF) [34] are used for noise removal. Guided filter (GF) [13], L0 Smoothing (L0) [17], bilateral filter (BF) [31], weighted least squares filter (WLS) [10], domain transform filter (RF) [11], adaptive manifolds filter (AMF) [19], rolling guidance filter (RGF) [24] are used in edge-preserving image smoothing. Local extrema filter (LE) [29] and relative total variation (RTV) [33] are often used for image decomposition purposes. Table 1 shows reverse filtering PSNR results for the cameraman image (see Fig. 2) for the above-mentioned filters. We compare the proposed method with the methods of Tao et al. [30] and Dong et al. [9]. The parameters for our Landweberlike iteration method (15) are set to N = 10K , τ = 0.1, and γ = 0 (i.e., the total variation regularization is not used).
We start from considering the noise-free linear filters from the above list. For them, the Dong et al. method [9]   For nonlinear filters, Tao's method usually performs best, while Dong's method often has difficulty giving reasonable results, or even fails. In contrast, the proposed method demonstrates always good results, though it does not always deliver the best performance. See Fig. 3

Deblurring
In the previous section, we have demonstrated the effectiveness of the proposed method for reversing Gaussian and motion blur. This motivates us to explore applying the proposed method to the deblurring task. We model a blurred image by where y is the observed blurred image contaminated by noise, h is a blurring kernel (we assume that h is not known, but f (·) is available), is the convolution operation, x is the latent image to be recovered, and n is additive noise. In addition to (16), we also consider the corresponding clean filter obtained from (16) by setting n = 0. As illustrated in Fig. 1, the methods of Tao et al. and Dong et al. [9,30] do not work well when dealing with noisy filters. In contrast, as seen in Fig. 4 c, d, our approach is robust to noise. In this example, imnoise(imfilter(x,fspecial('motion', 21,11),'circular'), 'gaussian',0,0.001) is used as a noisy filter (with MATLAB notation). Using PSNR as a quality measure for defiltering, the best result is obtained at around 850 iterations of our Landweber-like iterations, and subsequent iterations seem to degrade the quality of the image (see Li [16] and Pan [1] are the state-of-the-art blind deblurring methods. Chan [4] is the state-of-the-art non-blind deblurring method. The proposed approach corresponds to iterations of (15) 123 Fig. 4). Continuing to iterate the method after the PSNR peak, the resulting images seem to improve very little visually. This can be explained by the fact that the image has already recovered boundaries and textures, but the noise of the image gets also amplified. To better understand this phenomenon, the result of 2K iterations is shown in Fig. 4 e.
In Fig. 4 f, we use several iterations of the proposed method with total variation regularization (15) to reverse the same noisy filter. Visual inspection and numerical evaluation indicate that the total variation acts as a powerful regularization to the proposed method, which can significantly remove the noise and lead to a much higher PSNR value.
To verify the generality of our approach, we try it with different kinds of blur and noise. We used Gaussian and Poisson noise, respectively, as they are the most common noise one deals with. Experimental results show that the proposed method can significantly remove image blurring (Gaussian blur, Motion blur, or even blur described by an unknown kernel) under these two types of noise. We compare our results with the state-of-the-art deblurring methods [1,4,16] in Fig. 5. The performance of the proposed method is between nonblind and blind methods. This observation is not surprising, as compared with blind methods, we can query the filter as a black-box. However, compared to non-blind methods, we do not have exact information about the kernel. From this point of view, the method developed in this work can be treated as a semi-blind deblurring method.

Image super-resolution
Tao et al. [30] demonstrated that reverse filtering could be used to improve results obtained by image super-resolution methods. In this section, we demonstrate that (14) is also capable to improve image super-resolution results and performs better than [30].
We consider a low-resolution image as a down-sampled (filtered) version of a high-resolution image. Tao et al. proposed in [30] that the down-sampling operation and image super-resolution process could be regarded as a black-box filter, so as to ensure the image size consistency. Let us take bicubic interpolation as an example, one can consider for the black-box filter imresize(imresize(x,1/2,' bicubic'),2,' bicubic') using MATLAB notation. In Fig. 6 a-c, we demonstrate this idea by combining iterations of (13), i.e., Algorithm 1, with bicubic interpolation. The PSNR value shows that the proposed method further improves the quality of the restored image compared to the approach of Tao et al. [30]. However, we observe that the interpolation-based method may produce ringing artifacts in the position where the gradient changes significantly. We note that edge artifacts increase the image gradient. Therefore, it is expected that iterations of (15), with the total variation prior, could help removing this artifact. We verified this fact experimentally in Fig. 6d.   (a) (b) (c) (d) Fig. 6 Image super-resolution enhanced by iterations of (13) and (15) With the regularization from the total variation term, the artifacts are removed, and the quality of the restored image has been significantly improved qualitatively and quantitatively. For a more rigorous evaluation, we combine the approach developed in this work with several commonly used image interpolation methods. The best results (selected within 3K iterations for interpolation-based approaches and 300 for SRCNN [5]) are compared to the method of Tao et al. [30] in Table 2. PSNR is used as a quantitative measure for the evaluation. The restored images are shown in Fig. 7 for visual comparison. These results show that the proposed method can enhance various interpolation-based methods, and that it performs better than that of Tao et al. For interpolationbased methods, Tao's method produces edge artifacts at the edges of the image, such as the junction of the beak with the background. In contrast, this artifact is suppressed in the results obtained via the method proposed in this work, see Fig. 7e and f. Moreover, the proposed method also shows some improvements when the advanced SRCNN approach is used, but the improvements are more nuanced than when interpolation-based methods are used.

Accelerated Landweber-like iterations
Although the proposed method (13) demonstrates good reverse filtering results, it requires many iterations. Below we consider two acceleration schemes for (13). Acceleration by phase correction. The van Cittert and Landweber methods are the simplest and most fundamental iterative methods for linear inverse problems [3]. One can easily observe that if image filter f (·) is linear, then the van Cittert method corresponds to (4) and (3), while (11) constitutes Landweber iterations in the frequency domain. As noted in [3,Chapter 6], the method of Landweber can be applied to any imaging system, while the van Cittert method requires an imaging system to have some particular properties. In our experiments, we have already seen that (4) and (3) failed to deal with non-symmetric blurring filters. On the other hand, if the van Cittert method converges, it does it faster than the Landweber one, as each iteration of the latter involves reblurring.
As H * k /|H k | does only a phase correction, it seems natural to call (17) by phase-corrected Landweber iterations. Figure 8 demonstrates that (17) accelerates (13) and is capable to deliver better results for reversing linear and nonlinear filter (the MT and AMF filters are used). Setting τ = 1 in (17) leads to a simplified version of (17) which can be considered as a phase-corrected version of the method of Tao et al. [30], as it transfers (3) into the frequency domain and corrects it by a pure phase filter H * k /|H k |. This phase-corrected version of (17) is used to produce Fig. 1 (f) result. Nesterov-accelerated Landweber-like iterations. Another way to accelerate our basic method (13) is to combine it with the Nesterov gradient descent acceleration scheme [21].
Figure 8 a and g shows that the Nesterov-accelerated iterations converge faster than (13) for linear filters. However, some instability may be developed if reverse filtering with nonlinear filters is considered.

6 Conclusion
In this work, we have considered the reverse image filtering problem and proposed to use Landweber-like iterations for image defiltering purposes. The proposed method is robust to noise and produces state-of-the-art defiltering results for a number of popular linear and nonlinear image filters. In addition to image defiltering, we have verified the effectiveness of our method for image deblurring and super-resolution. For the image super-resolution problem, the proposed method can be applied as a post-processing step and, according to our experiments, is capable of achieving a significant improvement of state-of-the-art image super-resolution methods.
The main limitation of the proposed method is that it demonstrates a slow convergence and too many iterations are required to achieve good defiltering results. We have addressed this limitation by considering and evaluating two acceleration schemes. One of the schemes combines Landweber-like iterations with phase correction steps and, in addition to producing good reverse image filtering results, seems useful for other image processing tasks.