# On Fast Bilateral Filtering using Fourier Kernels

It was demonstrated in earlier work that, by approximating its range kernel using shiftable functions, the non-linear bilateral filter can be computed using a series of fast convolutions. Previous approaches based on shiftable approximation have, however, been restricted to Gaussian range kernels. In this work, we propose a novel approximation that can be applied to any range kernel, provided it has a pointwise-convergent Fourier series. More specifically, we propose to approximate the Gaussian range kernel of the bilateral filter using a Fourier basis, where the coefficients of the basis are obtained by solving a series of least-squares problems. The coefficients can be efficiently computed using a recursive form of the QR decomposition. By controlling the cardinality of the Fourier basis, we can obtain a good tradeoff between the run-time and the filtering accuracy. In particular, we are able to guarantee sub-pixel accuracy for the overall filtering, which is not provided by most existing methods for fast bilateral filtering. We present simulation results to demonstrate the speed and accuracy of the proposed algorithm.

Comments

There are no comments yet.

## Authors

• 4 publications
• 20 publications
• ### Fast and Accurate Bilateral Filtering using Gauss-Polynomial Decomposition

The bilateral filter is a versatile non-linear filter that has found div...
05/01/2015 ∙ by Kunal N. Chaudhury, et al. ∙ 0

read it

• ### Fast and Provably Accurate Bilateral Filtering

The bilateral filter is a non-linear filter that uses a range filter alo...
03/26/2016 ∙ by Kunal N. Chaudhury, et al. ∙ 0

read it

• ### Constant-time filtering using shiftable kernels

It was recently demonstrated in [5] that the non-linear bilateral filter...
07/22/2011 ∙ by Kunal Narayan Chaudhury, et al. ∙ 0

read it

• ### Index Set Fourier Series Features for Approximating Multi-dimensional Periodic Kernels

Periodicity is often studied in timeseries modelling with autoregressive...
05/14/2018 ∙ by Anthony Tompkins, et al. ∙ 0

read it

• ### Malliavin-Mancino estimators implemented with non-uniform fast Fourier transforms

We implement and test kernel averaging Non-Uniform Fast-Fourier Transfor...
03/05/2020 ∙ by Patrick Chang, et al. ∙ 0

read it

• ### Recursive Filter for Space-Variant Variance Reduction

We propose a method to reduce non-uniform sample variance to a predeterm...
11/12/2019 ∙ by Alexander Zamyatin, et al. ∙ 15

read it

• ### Optimal Envelope Approximation in Fourier Basis with Applications in TV White Space

Lowpass envelope approximation of smooth continuous-variable signals are...
06/03/2017 ∙ by Animesh Kumar, et al. ∙ 0

read it

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

The bilateral filter was introduced by Tomasi and Manduchi in [1]

as a non-linear extension of the classical Gaussian filter. The bilateral filter employs a range kernel along with a spatial kernel for performing edge-preserving smoothing of images. Since its introduction, the bilateral filter has found widespread applications in image processing, computer graphics, computer vision, and computational photography

[2] - [8].

In this paper, we will consider a general form of the bilateral filter where an arbitrary kernel is used for the range filtering, and a box or Gaussian kernel is used for the spatial filtering [1]. In particular, consider an image , where is a finite rectangular lattice. The output of the bilateral filter is given by

 fBF(i)=∑j∈Ωw(j) ϕ(f(i−j)−f(i)) f(i−j)∑j∈Ωw(j) ϕ(f(i−j)−f(i)), (1)

where is the range kernel and is the spatial kernel. The spatial kernel is usually a Gaussian [1],

 w(i)=exp(−∥i∥2/2σ2s). (2)

The window of the spatial kernel is a local neighbourhood of the origin. For example, for the Gaussian kernel, where . The original proposal in [1] was to use a Gaussian range kernel given by

 ϕ(t)=exp(−t2/2σ2r). (3)

In more recent work, exponential range kernels have been used [9, 10, 11].

The direct computation of (1) requires operations per pixel. In fact, the direct computation is slow for practical settings of [12]. To address this issue, researchers have come up with various fast algorithms [12] - [17]. While some of these algorithms can reduce the complexity to operations per pixel for any arbitrary , there is, however, no available guarantee on the approximation quality that can be achieved using these algorithms. In fact, as reported in [16], a poor approximation can lead to visible distortions in the filtered image. Only recently, a quantitative analysis of Yang’s fast algorithm was presented in [18].

In Section II, we recall the idea of constant-time bilateral filtering using Fourier (complex exponential) kernels [19]. In this work, we build on this idea to propose a new algorithm for approximating (1) using the shiftable Fourier basis. The contribution of this work is not the fast algorithm itself, but rather the approximation scheme in Section III, and the subsequent approximation guarantee in Section IV. The approximation scheme can be applied to any arbitrary range kernel that has a pointwise-convergent Fourier series. In this respect, we note that all previous approaches based on shiftable approximation were restricted to Gaussian range kernels [16, 19, 20]. We provide some representative results concerning the speed and accuracy of the resulting algorithm in Section V, where we also compare the empirical accuracy of the filtering with the bounds predicted by our analysis.

## Ii Shiftable Bilateral Filtering

It was demonstrated in [16, 19] that the bilateral filter can be decomposed into a series of Gaussian convolutions using shiftable functions. In particular, since our present interest is in the shiftable complex exponential, consider the function

 φN(t)=N∑n=−Ncnexp(ιnωt), (4)

where . By setting (4) as the range kernel , we can decompose the numerator in (1) as

 N∑n=−Ncnexp(−ιnwf(i))Fn(i),

where

 Fn(i)=∑j∈Ωw(j)f(i−j)exp(ιnωf(i−j)). (5)

It is clear that a similar decomposition can be obtained for the denominator of (1). We readily recognize (5) to be a Gaussian convolution. As is well-known, the Gaussian convolution in (5) can be efficiently implemented at constant-time complexity (with respect to ) using separability and recursion [21]. In summary, we can decompose the bilateral filtering into a series of Gaussian filtering. The fast shiftable algorithm resulting from this decomposition is summarized in Algorithm 1. We use in line 1 to denote the complex-conjugate of . In line 1, we use and to denote the Gaussian filtering of the images and . To avoid confusion, we note that the formal structure of Algorithm 1 is somewhat different from that of the shiftable algorithms in [16, 19]. While the cosine and sine components of the complex exponential were used in [16, 19], we work directly with the complex exponential in Algorithm 1. Note that we have abused notation in using to denote the shiftable approximation of (1) in Algorithm 1.

If the range kernel is not shiftable, one can approximate it using a shiftable function. For example, the shiftable raised-cosines were used in [16] to approximate the Gaussian kernel. Shiftable approximation using polynomials was later presented in [19]. More recently, the classical Fourier basis was used for this purpose in [20]. The above approximations, however, come with the following shortcomings:
They are customized to work with the Gaussian kernel, and cannot be extended to general range kernels, such as the exponential kernel [10, 11]. Even for the Gaussian kernel, the proposal in [20] requires one to compute the coefficients of the Fourier series. This is computationally intensive (e.g., requires numerical integration, or some analytical properties particular to the kernel), and cannot be done on-the-fly. Indeed, the authors in [20] work with an approximation of the Fourier coefficients, which is only valid for small .
Notice that, in most applications of the bilateral filter, the argument in (3) assumes discrete values. This should be taken into consideration while designing the shiftable approximation. The approximations in [16, 20], however, do not necessarily guarantee that the approximation error at these discrete points are within some user-defined tolerance. This makes it difficult to quantify the overall filtering accuracy. In this paper, we propose a rather simple optimization principle, which has an efficient implementation. This provides us with the desired control on the numerical accuracy of the overall filtering.

## Iii Progressive Fourier Approximation

We now explain how the above shortcomings can be fixed. As noted above, the argument in (3) takes on the values as and varies over the image. In particular, takes values in , where

 T=maxi∈I maxj∈Ω |f(i−j)−f(j)|.

Thus, is the dynamic range of the image measured over the window , which is typically smaller than the full dynamic range. We can compute using the fast algorithm in [22]; the run-time of the algorithm does not depend on the size of . Without loss of generality, we assume that the range kernel is symmetric. The problem is that of approximating using a shiftable function over the half-interval . We propose to use the shiftable Fourier basis for this purpose. In particular, we fix some order , and consider the shiftable function

 φN(t)=d0+N∑n=1dncos(nωt), (6)

where . As is well-known, using the identity , we can write (6) as in (4), where , and for . The key difference with [20] is with respect to the rule used to set the coefficients in (6). These are set to be the standard Fourier coefficients of in [20]. In keeping with the arguments presented in earlier, we take a different approach and instead try to minimize the error at the discrete points . In particular, we consider the problem of finding that minimizes the gross error

 ∑t∈ΛT (ϕ(t)−φN(t))2. (7)

This is the classical linear least-squares problem, where the unknowns are . Indeed, using matrix-notation, we can write (7) as , where , is the discretization of at the points , and the columns of are the corresponding discretization of the basis functions in (6). In particular, let us denote

 (8)

The following fact is the basis of our approximation algorithm to be discussed next.

###### Proposition III.1 (Decay of Error)

Assume that the Fourier series of the range kernel converges pointwise on the interval . That is, for ,

 limN→∞φN(t)=ϕ(t),

where in (6) are the Fourier coefficients of . Then decays to zero as .

Indeed, let be the error in (7) when is taken to be the -th order Fourier approximation of . Then, by optimality, we have . Since, by assumption, as , the proposition follows. We note that the Fourier series converges pointwise for any continuously-differentiable function, e.g., Gaussian and polynomials. Convergence is also guaranteed for functions that are continuous and piecewise-differentiable [23], such as the exponential. Thus, the assumption in Proposition III.1 covers the commonly used kernels [1, 10, 11].

Proposition III.1 suggests the following numerical scheme: We fix some user-defined tolerance . We begin with , and solve (8) to get . If , we stop. Else, we increase by one and proceed, until . In other words, we solve a series of least-squares problems, where the basis matrix at each step is obtained by augmenting the in the previous step. The whole process can be efficiently implemented using a recursive version of the modified QR algorithm [24]. The main idea is that (8) can be computed by solving using back-substitution, where is the QR-decomposition of . In the recursive computation, , and at each iteration is computed from the corresponding quantities in the previous iteration using cheap operations. An adaptation of this recursive algorithm to our problem is provided in Algorithm 2. In steps 2 and 2, we discretize the kernel and the incoming Fourier basis. In step 2, denotes the -th component of .

## Iv Filtering Accuracy

Suppose we are given a range kernel and tolerance . We compute the approximation order and the corresponding coefficients using Algorithm 2. This gives us the corresponding kernel in (4), which is used to approximate (1) using Algorithm 1. In particular, the approximation provided by Algorithm 1 is given by

 ^fBF(i)=∑j∈Ωw(j)φN(f(i−j)−f(i))f(i−j)∑j∈Ωw(j)φN(f(i−j)−f(i)). (9)

By construction, for all ,

 |ϕ(t)−φN(t)|≤ε. (10)

Similar to [18], we consider the (worst-case) error

 ∥fBF−^fBF∥∞=max{|fBF(i)−^fBF(i)|: i∈I}. (11)

Our goal is to bound (11

), which provides us with an estimate of the pixelwise difference between the outputs of the exact and the approximate bilateral filter. In fact, a simple analysis (cf. Appendix) give us the following result.

###### Proposition IV.1 (Filtering Accuracy)
 ∥fBF−^fBF∥∞≤2Tεw(0)−ε. (12)

In other words, the filtering error is essentially within a certain factor of the kernel approximation error . To arrive at (12), we have assumed that the weights of the spatial filter add up to unity. Indeed, this assumption can be made since the spatial filter appears in both the numerator and denominator of (1) and (9).

## V Simulation and Conclusion

All simulations reported here were performed using Matlab 8.4 on a MacBook Air with 1.3 GHz Intel Core i5 processor and 4 GB memory. The typical run-time of Algorithm 2 was between - milliseconds (depending on the order ) for the simulations reported in this section. This is a small fraction of the overall run-time of Algorithm 1. Indeed, the time required to filter a single image with a Gaussian kernel is already about milliseconds. In Figure 1, we give an example of the approximation result obtained using Algorithm 2 with . In Figure 2, we compare the coefficients obtained using Algorithm 2 with that obtained by expanding the raised-cosines [16] into (4). Notice that the former decays much more rapidly and hence requires fewer terms.

We present some results on the Barbara image for which was computed to be . We note that the run-time of the direct implementation of (1) depends only on the image size and . On the other hand, the run-time of the proposed algorithm depends on , tolerance , image size, and . The fact that the run-time is almost independent of (constant-time algorithm) is evident from the results in Table I

. The small fluctuations are essentially due to the variable padding required to handle the boundary conditions for the spatial filtering.

The run-time of the proposed algorithm scales inversely with , which was also observed for the shiftable filtering in [16, 22, 20]. In particular, as gets small, the Gaussian range kernel tends to a Dirac-like distribution [23]. As is well-known, the Dirac distribution is formally composed of all frequencies. The implication of this fact is that a large is required to approximate the kernel for small , and hence the increase in run-time. This is demonstrated with an example in Table II. However, notice that even for small , the proposed algorithm is much faster than the exact implementation. For , the speedup is by a couple of orders.

In Table III, we present the variation of run-time with tolerance for a fixed filter setting. It is seen that the order and hence the run-time changes rather slowly with (almost logarithmically in ). We have however not been able to establish this empirical fact, which is deeply tied to the working of Algorithm 2. We next compare the bound in (12) with the actual error for the Barbara image in Table IV. We note that the error is within the predicted bound. In fact, we are able to predict sub-pixel accuracy when . The bounds are, however, far from being tight. One of the reasons for this is that we have not incorporated any information about the local intensity distribution into our analysis. Derivation of a tighter bound will require a more sophisticated analysis. The present work is a first step in that direction. To best of our knowledge, with the exception of [13], this is the only approximation algorithm that comes with a provable guarantee on the filtering accuracy.

## Vi Acknowledgement

This work was supported by the Startup Grant awarded by the Indian Institute of Science. The authors would like to thank the anonymous reviewers for their comments and suggestions.

## Vii Appendix

In this section, we outline the main steps in the derivation of (12). We write (1) as , where

 P1(i)=∑j∈Ωw(j)ϕ(f(i−j)−f(i))f(i−j),

and

 Q1(i)=∑j∈Ωw(j)ϕ(f(i−j)−f(i)).

Similarly, we write (9) as , where

 P2(i)=∑j∈Ωw(j)φN(f(i−j)−f(i))f(i−j),

and

 Q2(i)=∑j∈Ωw(j)φN(f(i−j)−f(i)).

Then can be expressed as

 1Q2(i)[fBF(i)(Q2(i)−Q1(i))+P1(i)−P2(i)]. (13)

From (10), we have . On the other hand, note that . This is because is given by the convex combination of . Therefore, from (10), we get .

To obtain a lower-bound for in (13), we note that

 Q1(i)=w(0)φ(0)+positive terms≥w(0),

where we have used the non-negativity of the range and spatial kernels. Therefore, using the inverse triangle inequality, we get

 |Q2(i)|≥Q1(i)−|Q2(i)−Q1(i)|≥w(0)−ε.

By incorporating the above bounds into (13), we get (12).

## References

• [1] C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” Proc. IEEE International Conference on Computer Vision, pp. 839-846, 1998.
• [2] E. P. Bennett and L. McMillan, “Video enhancement using per-pixel virtual exposures,” ACM Transactions on Graphics, vol. 24, no. 3, pp. 845-852, Proc. ACM Siggraph, 2005.
• [3] H. Winnemoller, S. C. Olsen, and B. Gooch, “Real-time video abstraction,” Proc. ACM Siggraph, pp. 1221-1226, 2006.
• [4] J. Xiao, H. Cheng, H. Sawhney, C. Rao, and M. Isnardi, “Bilateral filtering-based optical flow estimation with occlusion detection,” Proc. European Conference on Computer Vision, pp. 211-224, 2006.
• [5] E. P. Bennett and L. McMillan, “Video enhancement using per-pixel virtual exposures,” ACM Transactions on Graphics, vol. 24, no. 3, pp 845-852, 2005.
• [6] K. N. Chaudhury and K. Rithwik, “Image denoising using optimally weighted bilateral filters: A SURE and fast approach,” Proc. IEEE International Conference on Image Processing, pp. 108-112, 2015.
• [7] B. M. Oh, M. Chen, J. Dorsey, and F. Durand, “Image-based modeling and photo editing,” Proc. Annual Conference on Computer Graphics and Interactive Techniques, pp. 433-442, 2001.
• [8] J. Xiao, H. Cheng, H. Sawhney, C. Rao, and M. Isnardi, “Bilateral filtering-based optical flow estimation with occlusion detection,” Proc. European Conference on Computer Vision, pp. 211-224, Springer, 2006.
• [9] B. K. Gunturk, “Fast bilateral filter with arbitrary range and domain kernels,” IEEE Transactions on Image Processing, vol. 20, no. 9, pp. 2690-2696, 2011.
• [10] K. Al-Ismaeil, D. Aouada, B. Ottersten, and B. Mirbach, “Bilateral filter evaluation based on exponential kernels,”

International Conference on Pattern Recognition

, pp. 258-261, 2012.
• [11] Q. Yang, “Hardware-efficient bilateral filtering for stereo matching,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 36, no. 5, pp.1026-1032, 2014.
• [12] F. Durand and J. Dorsey. “Fast bilateral filtering for the display of high-dynamic-range images,” ACM Transactions on Graphics, vol. 21, no. 3, pp. 257-266, 2002.
• [13] Q. Yang, K. H. Tan, and N. Ahuja, “Real-time bilateral filtering,” Proc. IEEE Conference on Computer Vision and Pattern Recognition, pp. 557-564, 2009.
• [14] S. Paris and F. Durand, “A fast approximation of the bilateral filter using a signal processing approach,” Proc. European Conference on Computer Vision, pp. 568-580, 2006.
• [15] F. Porikli, “Constant time bilateral filtering,” Proc. IEEE Conference on Computer Vision and Pattern Recognition, pp. 1-8, 2008.
• [16] K. N. Chaudhury, D. Sage, and M. Unser, “Fast bilateral filtering using trigonometric range kernels,” IEEE Transactions on Image Processing, vol. 20, no. 12, pp. 3376-3382, 2011.
• [17] K. N. Chaudhury, “Fast and accurate bilateral filtering using Gauss-polynomial decomposition,” Proc. IEEE International Conference on Image Processing, pp. 2005-2009, 2015.
• [18] S. An, F. Boussaid, M. Bennamoun, and F. Sohel, “Quantitative error analysis of bilateral filtering,” IEEE Signal Processing Letters, vol. 22, no. 2, pp. 202-206, 2015.
• [19] K. N. Chaudhury, “Constant-time filtering using shiftable kernels,” IEEE Signal Processing Letters, vol. 18, no. 11, pp. 651 - 654, 2011.
• [20] K. Sugimoto and S. I. Kamata, “Compressive bilateral filtering,” IEEE Transactions on Image Processing, vol. 24, no. 11, pp. 3357-3369, 2015.
• [21] R. Deriche, “Recursively implementing the Gaussian and its derivatives, Research Report, INRIA-00074778, 1993.
• [22] K. N. Chaudhury, “Acceleration of the shiftable algorithm for bilateral filtering and nonlocal means,” IEEE Transactions on Image Processing, vol. 22, no. 4, pp. 1291-1300, 2013.
• [23] L. Grafakos, Classical Fourier Analysis, vol. 2, New York: Springer, 2008.
• [24] J. Demmel, Applied Numerical Linear Algebra, Society for Industrial and Applied Mathematics, 1997.