Two iterative methods for reverse image filtering

We consider the problem of recovering an original image 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 its filtered version 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}, assuming that the internal structure of 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} is unknown to us (i.e., we can only query the filter as a black-box and, for example, cannot invert it). We present two new iterative methods to attack the problem, analyze, and evaluate them on various smoothing and edge-preserving image filters.


Introduction and contribution
In this paper, we deal with a reverse imaging problem introduced recently in [19]. Given an image filter y = f (x) with input x ∈ R d and output y ∈ R d , the problem consists of restoration of x from y under the condition that we can apply the filter f (·) as many times as needed but without any knowledge of its internal structure. In particular, it means that we cannot directly invert the filter. Following [19], we call such a restoration process reverse filtering or defiltering.
The reverse image filtering schemes introduced in [19] and in the subsequent works [4,12] are aimed at removing slight filtering effects. In contrast, our work is focused on removing, at least partially, severe filtering effects. Three examples are shown in Fig. 1, where reverse filtering of images filtered by a Gaussian, Laplace-of-Gaussian (LoG), and the adaptive manifold filter (AMF) [6] is demonstrated.
Mathematically, the defiltering problem consists of numerically solving the equation B Alexander G. Belyaev a.belyaev@hw.ac.uk Pierre-Alain Fayolle fayolle@u-aizu.ac.jp 1 School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, UK 2 Computer Graphics Laboratory, University of Aizu, Aizu-Wakamatsu, Japan for x, given y. Here, x is an (unknown) input image, f (·) is a (black-box) filter, and y is the (given) output image. Numerical methods for solving (1) include the Newton's method and its relatives, various gradient descent schemes, and fixed-point iterations. Straightforward uses of the Newton's and gradient descent methods require an estimation of the Jacobian of f (·) and are not computationally feasible, as a typical image consists of tens of thousands to millions of pixels.
A simple and efficient fixed-point iteration scheme for image defiltering was recently proposed by Tao et al. [19]. The scheme consists of the following iterations and can be considered as an iterative unsharp masking with filter f (·). Mathematically, (2) corresponds to the Picard's iterative method x n+1 = f (x n ), a simple and popular method for numerically solving nonlinear equation (1) (see, for example, [1] for a comprehensive study of fixed-point iteration methods). A rigorous convergence analysis of (2) for image filters satisfying certain conditions is given in [19]. As shown in [19], (2) demonstrates excellent results for removing light filtering effects but, according to our experiments, often fails when dealing with severe image deformations. So in this paper, we introduce two iterative methods and (4) which are inspired by a gradient descent scheme studied by Polyak [15] and Steffensen's method [17], respectively (hence the notations P(·) and S(·)). We propose a brief analysis of (2), (3) and (4) in the next section. Given x n the reconstructed image at the n-th step of the iterative process, we consider two error functions in order to measure the reconstruction accuracy Although we use both error functions to evaluate the T, P, and S-iterative processes, only E y can be used in practice, as x is not known a priori. Figure 2 presents the E x and E y error graphs corresponding to the reverse filtering examples shown in Fig. 1.

Analysis of the previous and proposed methods
Numerical methods for solving (1) include Newton's method and its relatives, as well as various gradient descent schemes, and fixed-point iterations [13]. Unfortunately, the vast majority of these methods require knowledge of the gradient of f (·) or impose too many restrictions on the behavior of f (·). Tao et al. [19] use a fixed-point iteration method to analyze convergence properties of their T-iterations (2). Unfortunately, their analysis can be applied only to linear mappings.
A similar scheme (6) was proposed by Milanfar [12] who used different considerations. Here, γ and μ are positive parameters chosen such that both γ < 1 and γ μ is small. Although the iterative scheme (6) is supported by a solid analysis of its convergence properties [12], our experiments have not revealed practical advantages of (6) over (2). In Fig. 3 we present the E x -error graphs for defiltering six linear and nonlinear filters using (2) and (6) with γ μ = 0.001 and μ = 0.5 k , k = 0, 1, 2, 3. One can see that in the majority of cases, (2) slightly outperforms (6). In [4], an interesting iterative scheme for solving (1) was proposed Here, F stands for the Fourier transform and element-wise multiplication and division are used. Interestingly, (7) can be considered as a variant of (6) employed in the frequency domain and used with variable parameter γ n instead of γ and  (2) and (6). Top-left: Gaussian; top-right: AMF [6]; middle-left: bilateral filtering; middle-right: RTV [24]; bottom-left: WMF [26]; bottom-right: RGF [25]. In the majority of these tests (2), corresponding to the curves in black, slightly outperforms (6) μ = 0. Indeed, the iterative process yields (7). This iterative scheme shows an excellent performance for linear filters f (·), provided that periodic boundary conditions are imposed, as demonstrated by the top-left and top-middle images of Fig. 4. Unfortunately, if other types of boundary conditions are chosen, (7) does not lead to good defiltering results, as shown by the bottom-left and bottommiddle images. In addition, this method does not show a good performance for some nonlinear filters, as seen in the right images of Fig. 4. It is interesting to note that the problem of reverse image filtering is closely related to stochastic approximation [2,9], an active research area initiated by Robbins and Monro (RM) almost seventy years ago [16]. Remarkably, the RM algorithm looks very similar to (2). Given noisy measurements where {η n } are measurement errors, the algorithm implements a simple iterative procedure for solving the equation defiltering Gaussian defiltering LoG defiltering RTV [24] defiltering Gaussian defiltering LoG defiltering AMF [6] Fig. 4 Defiltering imgaussfilt(x,5,'Padding','circular') (top-left) and imgaussfilt(x,5) (bottom-left) by (7). Similar defiltering results for imfilter(x,fspecial('log',7,0.4)) with (top-middle) and without (bottom-middle) periodic boundary conditions by (7). Top-right: defiltering RTV [24] with default settings by (7). Bottom-right: defiltering AMF [6] adaptive_manifold_filter(x,20,0.4) by (7) y = f (x) where {a n } is a sequence of nonnegative real numbers satisfying a n = ∞ and a 2 n < ∞. Our P-iterations (3) can be considered as an approximation of the gradient descent with a variable step-size for minimizing least-square energy whose anti-gradient is approximated by a symmetric finite difference where, as before, h n = y−f (x n ) is assumed to be sufficiently small. It turns out that in the one-dimensional case (single-pixel images), our iterative scheme (3) approximates the standard Newton iterations for numerically solving y = f (x). Indeed, the Newton iterations for solving F(x) ≡ y − f (x) = 0 are given by WMF [26] RTV [24] RGF [25] ILS [11] Now we can observe that in the one-dimensional case (3) is reduced to It is also worth mentioning a link between (3) and the Kiefer-Wolfowitz (KW) stochastic approximation method [8]. The KW algorithm for a univariate function y = f (x) consists of an iterative process where c n → 0 in addition to some other conditions imposed on the sequences {a n } and {c n }. A multivariate version of the KW algorithm can be found, for example, in [9] (initially both the RM and KW algorithms were proposed by their authors for univariate functions).
In the one-dimensional case, our S-iterations (4) become Steffensen's method [17] x which delivers a numerical solution to y = f (x) and has quadratic convergence like Newton's method.

Numerical experiments and discussion
In this section, we test our methods (3) and (4), and compare them against (2) using a wide selection of edge-preserving image filters, in addition to the two linear filters (Gaussian and LoG) and the adapted manifold filter (AMF) [6] considered in Sect. 1. Results obtained are shown in Fig. 5 for the image filters: WMF [26], RTV [24], RGF [25], ILS [11], WLS [5] and L 0 [23], and in Fig. 6 for the bilateral filter (BF), LLF [14] and LE [18]. For each image filter, we show an example of a filtered test image (top row), followed by the results obtained from T-iterations of Tao et al. [19] (2), P-iterations (3) and S-iterations (4) in the second, third and fourth row (from the top). For each method, the minimal E x -error is indicated under the image. The bottom row provides error graphs for E x (5) for the three methods. The best method is decided based on the error function E x (5) and indicated by a thick color line.  The decision is made based on the error function E x . P stands for P-iterations (3), S for S-iterations (4), and T for T-iterations (2) Surprisingly, reverse filtering of linear filters (in this paper, we consider Gaussian and LoG filters) turns out to be a difficult task. In the case of periodic boundary conditions, perfect defiltering is delivered by the method of Dong et al. (7) proposed in [4]. However, (7) is very sensitive to small perturbations and fails to recover images resulting from linear filters with non-periodic boundary conditions (see, e.g., Fig. 4).
The method of Tao et al. (2) [19] works well when it is used for defiltering Gaussians with small variances. In our example of Gaussian defiltering (see the left images in Figs. 1 and 2), the variance is σ 2 = 5 and (2) fails to produce a reasonably good result. Our experiments show that (2) also fails to recover the LoG filtering results (see, for example, Figs. 1 and 2, middle column).
Our P-iteration scheme (3) is capable of defiltering the Gaussian and LoG filters, as demonstrated in Figs. 1 and 2. However, the convergence can be slow. For example, restoring the left and middle images of Fig. 1 requires 100K iterations.
In addition to the Gaussian and LoG filters, we test ten state-of-the-art edge-preserving filters, see Figs. 5 and 6 and the right images of Figs. 1 and 2. For some filters, e.g., AMF, excellent defiltering results are obtained. For some others, e.g., LE, only very modest results are achieved. In general, we see that our P-iteration and S-iteration schemes together outperform the T-iteration method (2) [19]. However, for zero-order reverse filtering, (2) is the fastest method, as it requires only one call of function f (·) per iteration and, if it approaches an error minimum value, does it quickly.
We find it surprising that for the majority of the tested filters, small-scale image details can be accurately recovered if a proper defiltering method is applied. This means that typically edge-preserving filters only suppress those smallscale details but do not remove them completely.
In our experiments, we observe only a weak correlation between the two error functions E x and E y (5) for an iteratively defiltered image x n . Namely, if E y does not grow, then it is very likely that the defiltered image has its main features preserved, as demonstrated by Fig. 7. In real-life scenarios when only E y is available, the number of reverse filtering iterations has to be selected manually. This is the only userspecified parameter in our approach.
While, in this paper, we deal with deterministic reverse filtering methods, we feel that randomized algorithms could be very successful at achieving high-quality defiltering results. Fig. 7 Overrestoration examples. Left: WMF defiltering with 100 iterations of (4) amplifies certain middle-range frequencies. Right: BF defiltering with 1K iterations of (4) amplifies some high frequencies.
Compare these results with those for WMF and BF from Figures 5 and 6, respectively LE [18] WLS [5] Fig. 8 R-iterations (10) achieve the smallest E x -errors for the LE filter [18] (left) and WLS filter [5] (right). Top: images defiltered by (10). Bottom: the corresponding E x -error graphs As an example, let us consider the following simple stochastic iterative process where quadratic energy E is defined by (8), r is a random image whose components are random numbers uniformly distributed on the interval [0, 1] and c n = 1/ √ n. One can see that (10) combines (9) with stochastic gradient descent ideas. Pillbox filtering and defiltering with (3). c Motion blur filtering and defiltering with (3). d

Image-guided filtering
and defiltering with (4). e Image-guided filtering using a smoothed guidance and defiltering with (4). See the main text for details. Zoom in the images to see how well small-scale details are restored While (10) shows a modest performance for the majority of the considered filters, it allows us to achieve some progress with defiltering the LE filter [18] for which (10) shows the lowest E x error. In Fig. 8, we show the result obtained by (10) for defiltering the LE and WLS filters. The top row shows the best obtained results, while the bottom row compares the result of the R-iterations (10) against the results of the T, S and P-iterations for the E x error function. Visually (10) does a better work for retrieving fine details and high frequency texture in the image.
When dealing with color images, the simplest approach consists of defiltering each RGB channel. Obviously, more sophisticated strategies should lead to better results. Figure 9 shows several examples of color images, defiltered with our approach, from different domains (nature, architecture, humans, animals). In these examples, we test defiltering capabilities of our schemes (3) and (4) using popular image filters which have not yet been considered so far in this paper: (a) Wiener filter, (b) circular averaging filter (pillbox), (c) motion blur, (d) image guided filter [7], and (e) image guided filter, where Gaussian smoothing is applied to the guidance image: (a) wiener2(x, [5,5] For the linear filters (first three filters), reverse filtering was done by applying (3). The remaining nonlinear filters were defiltered by (4). In these five examples, for each color channel, we stop iterative processes (3) and (4) immediately after the relative change in successive iterations becomes less than a user-specified threshold t > 0 Better results could be obtained if an optimal number of iterations is chosen for each filter. Achieving reasonably good threshold t in (11) (for all the examples in Fig. 9 defilering is done with t = 0.0005) demonstrates the robustness of our schemes (3) and (4). One seemingly promising application of defiltering schemes lies in image enhancing. It seems natural to use the proposed reverse image filtering schemes as follows: The image x in (1) can be considered as an enhanced version of the input image y if f (·) is an appropriate filter. In Fig. 10, we combine AMF [6] with P-iterations (3) for image sharpening and with S-iterations (4) for image dehazing.

Conclusion and directions for future work
In this work, we have considered the problem of reverse filtering, or defiltering, a severely filtered image y = f (x). We assume no knowledge of the internal structure of the filter f (·) (and thus cannot compute its inverse). Two proposed methods (3) and (4) are compared against the existing method (2) introduced in [19]. The comparison is done on severely filtered images.
Possible directions of future work on defiltering include the use of stochastic approximation methods for image defiltering and adapting image reverse filtering schemes for geometric modeling applications (see, for example, the recent work [22] devoted to reverse geometry filtering).
Using reverse filtering schemes for image enhancing and sharpening constitutes another direction for future work. See a very recent work [3] where single-step polynomial defiltering is used for fast and natural-looking sharpening large-size images.
Iterative schemes (3) and (4) can be considered as zeroorder optimization methods. Problems requiring zero-order optimization appear frequently in signal processing and machine learning [10].
As (3) and (4) can be used to estimate how much of highfrequency image content is partially suppressed but remains preserved by an image filtering method, it would be interesting to establish a link with the ability of the method to protect deep learning schemes against adversarial attacks which often target high-frequency image components invisible to the human eye [20,21].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.