Backward Filtering of Images by Iterating Forward Filters

Given a filtering operation for images, we consider a problem for inverting it. If the algorithm of the original, i.e. forward filter is explicitly expressed, its inversion is tractable. However, if its algorithm is unknown exactly, it is hard to be reversed. We present, in this paper, a simple solution for this problem only by iterating the forward filter. This method can be used for a wide class of image filters including ones built in photo-retouch softwares. Examples are shown for smoothing, sharpening, halftoning and resizing filters.


Introduction
Major classes of image filtering operation are smoothing and sharpening. The smoothing filter includes Gaussian, bilateral [1], non-local means [2], median filters and so on. Typical example of sharpening filter is the unsharp masking. In this paper, we consider the inverse of these filters. Except for linear filters such as the Gaussian filter, their inversion is not so simple. Furthermore, many filters implemented in general-purpose photo-retouch softwares such as the Photoshop and the Paintshop are executed by built-in programs of which algorithm is unknown for users. This closedness inhibits their inversion.
A method [3] has been presented for inverting some filters such as the bilateral filter. We, in this paper, present a similar simple method for general image filters and analyzed the behavior of the method. In this paper, a scalar is denoted by a lower-case letter, vectors by thick lower-case letters and matrices by upper-case letters.

Fundamental Theory
Let x be a vector of pixel values in an image and y be a vector outputted by a filter as is illustrated in Fig. 1 where we call the ordinary filter "forward filter" in contrast to its reverse called "backward filter". The vectors x and y satisfy an equation a(x, y) = 0 where 0 is the zero vector. In this paper, we classify a(x, y) to the following three classes: type 1: a(x, y) is nonlinear with respect to both of x and y. type 2: a(x, y) = y − b(x), i.e. linear with respect to y. The function b(x) takes nonnegative values because y is pixel values. type 3: a(x, y) = x − c(y) where c(y) ≥ 0.

Iterative Methods
Iterative methods are commonly used for solving nonlinear equations. Iterations are generally expressed by x (ξ+1) = f(x (ξ) ) where ξ is the iteration counter ξ = 0, 1, 2, ... and x (ξ) is called ξth iterant. General theorems for their convergence are the Liapunov's theorem and the contractive mapping theorem [4], in which the latter one is mainly concerned in this paper except for the mode filter [5] for which the Liapunov's theorem is resorted. We use simple iterative methods for solving the equation a(x, y) = 0 as follows. type 1: we convert a(x, y) = 0 to y = y − ka(x, y) where k > 0. The iteration is written by backward filter : We convert a(x, y) = 0 to x = x + ka(x, y) for which iteration is type 2: forward filter : The equation a(x, y) = y − b(x) = 0 leads to y = b(x) which needs no iteration.
backward filter : and solve it with the iteration type 3: In this type, x and y are exchanged from the above type 2. The forward filter needs the iteration , while the backward filter is explicitly expressed as x = c(y).
Some examples of these three types and their convergence are presented in the following sections where the role of the coefficient k is explained.

Smoothing Filters
Let a value of pixel (i, j) in an input image be x i j , and an output of a filter be y i j . Arrays of these pixel values are x = [x i j ] and y = [y i j ]. Some filters are derived from a statistical estimation as

Mode Filter
In the mode filter [5], the objective function in eq.(4) is where s lm = e −α(l 2 +m 2 ) . There exists a unique solution for eq.(4) because this d(x, y i j ) is convex downward with respect to y i j and is given by This a(x, y i j ) is nonlinear with respect to both of x and y i j , hence this mode filter is of tyle 1 mentioned above.
forward filter : If we set k = 1, eq.(1) leads to which converges for arbitral initial values y (0) i j because d(x, y (ξ) i j ) decreases monotonically [4], i.e. d(x, y (ξ) i j ) is the Liapunov's function for this iteration.
backward filter : Similarly, by setting k = 1, eq.(2) becomes which is written in a vector form as is the coefficient matrix which is positive definite and its eigenvalues satisfy 0 < λ ≤ 1. The eigenvalues of I − H(x (ξ) , y) satisfy 0 ≤ 1 − λ < 1, hence the right hand side of eq.(8) is a contractive function of x (ξ) , which ensures the convergence of eq.(8) for arbitral initial values x (0) [4].

Bilateral Filter
In the bilateral filter (BF) [1], the objective function in eq.(4) is forward filter : This function a(x, y i j ) is linear with respect to y i j , hence this BF is of type 2. From eq.(10), we obtain backward filter : If we set k = 1, eq.(3) leads to which converges because its right hand side is a contractive function of x (ξ) . This convergence is ensured for an arbitral initial value, while we set x (0) i j = y i j for avoiding wasteful computation. Equation (12) is written in a vector form as We analyze the convergence process of this iteration. For simplifying the analysis, we approximate H(x) to constant H and simplifying eq.(13) to with substitution of y = Hx. From the graph spectral analysis [6], the frequency of signals is given by ω = 1 − λ where λ is the eigenvalue of H. From a straightforward manipulation of eq. (14), the frequency response of x (ξ) is derived as 1 − ω ξ+1 which is plotted for first three iterants x (0) , x (1) and x (2) in Fig. 2. These graphs reveal that each x (ξ) is a low-pass filter which converges to an all-pass filter except for the highest frequency ω = 1 where the response is always fixed to zero. Thus it is confirmed that the iteration of eq. (14) converges to the original x inputted to the forward filter in Fig. 1, i.e. x (∞) = x. Note the difference between the iteration of eq.(13) and the repeated unsharp masking: which diverges because the eigenvalue of 2I − H(x (ξ) ) is 2 − λ which is larger than 1. Fixation of y i j in eq.(12) is prerequisite for the convergence of x (ξ) . We will consider the unsharp masking again in section 4.
Equation (11) is the standard BF. If we set the denominator in eq.(11) to a constant value, we obtain a new type of BF called the unnormalized BF [7].

Regularized Nonlinear Diffusion Filter
The regularized nonlinear diffusion filter (RNDF) [8] is defined by which is linear with respect to x i j , hence this RNDF is of type 3.
forward filter : From eq.(17), a(x i j , y) = 0 leads to the iteration where ∑ l,m 0 means the sum for l and m skipping the center pixel where l = 0 and m = 0. An output of this RNDF for image in Fig. 3(a) is shown in Fig. 3(b) where we set p = 5, α = 0.01, β = 0.01, ϕ = 0.02.
by which we can get the image in Fig. 3(a) from Fig. 3 a weighted unsharp masking induced from the above mentioned unnormalized BF. Standard sharpening filters are analyzed below in section 4.

Built-in Filters
In addition to the above filters, many filters such as the median filter and the non-local means filter are of type 2. We can also use similar filters implemented in photo-retouch software or online at a webpage. Among them, we experimented a nonlinear smoothing filter in a photo-retouch software "Paintshop". Its output y is shown in Fig. 4(b) for the input image x in Fig. 4(a) which is an ear of Koala. Iterants in the iteration of eq.(3) are shown in Fig. 4(c) and (d). The iterant x (ξ) sufficiently converges at the fifth iteration. Another example of smoothing filter is shown in Fig. 5 where Fig. 5(b) is an output of the anisotropic filter [9] for Fig. 5(a). The algorithm in this filter is too complex to be analytically inverted. Iterants in the iteration of eq.(3) are shown in Fig. 5(c) and (d). This filter smoothes images strongly, hence its inversion needs more iteration than Fig. 4. Figure 7: Inversion of image enhancer.

Finally Converged Value
The iteration converges for all of these filters. The iterant x (ξ) almost converges to the original x, thus the filters can be inverted. In some highly nonlinear filters, however, the final iterant x (∞) deviates from the original x. This is because such filters remove some components in x completely. If we decompose x into a robust signal component x r and a frag- The iteration for inverting the forward filter y = b(x) converges to x (∞) = x r . An image with addition to 30% random noise to Fig. 3(a) is shown in Fig. 6(a) for which output of the median filter of 5 × 5 window is shown in Fig. 6(b). Noise in Fig. 6(a) is almost removed while some textures also disappear. Iterants are shown in Fig. 6(c),(d),(e) where the lost textures are gradually restored while noise still remains removed. Thus this iteration method is useful for removal of noise from images while maintaining textures.

Sharpening Filters
As is seen above, the coefficient k in the iterations in section 2.1 can be safely set to 1 for ensuring their convergence for smoothing filters. This is intrinsically owing to the contractiveness of these filters. Sharpening filters, however, does not satisfy this condition as is explained in this section. Typical sharpening filter is the unsharp masking (UM), e.g. if we write eq.(11) in a vector form as y = H(x)x, the unsharp masking derived from it is expressed by with µ > 0. The eigenvalue of (1 + µ)I − µH(x) is 1 + µ − µλ where λ is the eigenvalue of H(x). This eigenvalue 1+µ−µλ exceeds 1 because 0 < λ ≤ 1. this eigenvalue larger than 1 means that this filter magnifies high frequency components in x, i.e. UM is a sharpening filter. This UM is of type 2 We cannot stop such oscillating iteration at its midway due to errors. Hence 0 ≤ 1 − k(1 + µ − µλ) ≤ 1 is required for monotonic convergence which can be safely stopped at any midpoint. Therefore, 0 < k ≤ 1/(1 + µ) is desired for smooth convergence. Various sharpening filters are also included in photo- Figure 10: Iterative regularization.
retouch software. The magnification factor is unknown for those built-in filters, for such filters k should be safely set to a small value, typically k = 0.5. An example of foggy image as x is shown in Fig. 7(a) and its enhanced output y by an online software "One-Shot-Enhancer" in the webpage http://www.picturetopeople.org is shown in Fig. 7(b). Iterants x (ξ) are shown in Fig. 7(c) and (d) where we set k = 0.5.

Relaxed Iterative Methods
The above smoothing and sharpening are fundamental processing of images, hence its inversion is not so difficult. However, if filters are more advanced, the iteration formula should be extended.

Quantizing Filters
By quantizing x to discrete levels, y takes discrete values, extremely only two levels, black and white. For such quantization, halftoning techniques such as the error diffusion are used. Inversion of such quantizing filters becomes an integer programming (IP) problem which is hard to solve. We approximate it by a nonlinear programming (NP) problem by smoothing the discrete image y.
Let a quantizing operation be denoted by y = q(x) which is rewritten as q(x) − y = 0. We relax this equation to s(q(x) − y) = 0 where s(·) is a smoothing filter. Next we convert s(q( The iteration needs no complex computation demanded for many inverse halftoning methods [10] [11] even without explicit programming by utilizing a photo retouch software. We experimented this iteration for the images in Fig. 8 where (a) is an input grayscale image x and (b) is y with the use of the Floyd-Steinberg error diffusion as q(·). Some iterants are shown in Fig. 8 (c),(d) where we used the Gaussian filter as s(·).

Resizing Filters
When r(·) is a resizing operation, y = r(x) cannot be directly inverted because image size is different for x and y. Similarly to the above quantization, we relax r( . This iteration resembles the back-projection method for image superresolution [12] [13]. An example image shown in Fig. 9(a) is shrunk to Fig. 9(b) of which width and height are 1/4 of those of Fig. 9(a). Jaggies obstruct the perceptual quality of Fig. 9(b). Enlargement of Fig. 9(b) to Fig. 9(a) is an example called super-resolution. Iterants x (ξ) are shown in Fig. 9(c) and (d) where we set k = 0.5. We used the flow filter in Fig. 5 for interpolating enlarged images. Jaggies are weakened in Fig. 9(d).

Iterative Regularization
Equation (22) is effective even when q(x) = x. Image denoising is an operation for recovering x from y = x + n where n is noise. A conventional way of denoising is smoothing of y and outputting s(y) where s(·) is a smoothing filter. This smoothing removes noise but also smoothes out signal components slightly.
A method called the iterative regularization [14] restores lost signal components by iterating where its initial value is set to x (0) = s(y). Note that this iteration should be stopped at its midway because x approaches gradually to y and eventually x (∞) = y. Addition of Gaussian noise of standard deviation 10 to the image in Fig. 3(a) yields Fig. 10(a). Blurring this Fig. 10(a) by a Gaussian filter, we get Fig. 10(b) where noise is reduced while textures disappear slightly. The iterant x (1) and x (5) of eq.(23) are shown in Fig. 10(c),(d) where lost textures are partially recovered. We set s(·) as the Gaussian filter of the standard deviation 3.

Conclusion
We have presented a method for inverting image filters by simply iterating them. An outstanding benefit of the proposed method lies in its applicability to wide class of filters including built-in softwares of which algorithm is unknown exactly. We have experimented this method for smoothing and enhancing filters, image halftoning and resizing operations. Automatic setting rule of the coefficient in the iteration is under study for image sharpening and enhancing filters.