On Fast Bilateral Filtering using Fourier Kernels

03/26/2016 ∙ by Sanjay Ghosh, et al. ∙ ERNET India 0

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.



There are no comments yet.


page 1

page 2

page 3

page 4

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


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


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


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


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



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.

Data: Image ;
Parameters: Filter , and , , and .
Result: Shiftable approximation of (1).
1 Initialize: Set and for all ;
2 for  do
3       for ;
4       for ;
5       for ;
6       Compute and ;
7       for ;
8       for ;
10 end for
Set for all .
Algorithm 1 Shiftable Bilateral Filtering.

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

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


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


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


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 ,

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


By construction, for all ,


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


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)

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).

Data: Kernel , half-period , and tolerance .
Result: and .
1 ;
2 ;
3 ;
4 Initialize: ;
5 Set

to be the all-ones vector;

6 ;
7 ;
8 ;
9 ;
10 ;
11 while  do
12       Set to be the all-zeros vector;
13       ;
14       ;
15       ;
16       for  do
17             Set to be the -th column of ;
18             ;
19             ;
21       end for
22      ;
23       ;
24       ;
25       ;
26       Add a row of zeros to ;
27       ;
28       Solve using back-substitution;
29       ;
31 end while
Algorithm 2 Progressive Fourier Approximation.

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.

Fig. 1: Left: Target Gaussian and the Fourier approximation () obtained using Algorithm 1. Right: Coefficients .
Fig. 2: Comparison of the normalized binomial coefficients from [16] with that obtained using Algorithm 2 () when .

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.

1 2 5 8 10 12
Fast () 630ms 635ms 638ms 640ms 645ms 650ms
TABLE I: Run-time for the Barbara image at different and fixed . The run-time of the direct implementation is seconds.

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.

10 15 20 30 50 100
Fast () 2.1s 1.5s 840ms 634ms 450ms 200ms
TABLE II: Same as in Table I, except that is varied and . The run-time of the direct implementation is seconds.

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.

1e-5 1e-4 1e-3 0.01 0.1
12 11 10 8 7
Fast () 910ms 850ms 780ms 715ms 670ms
TABLE III: Variation of the run-time with for the Barbara image when and . Also shown is the order .
1e-8 1e-5 1e-4 1e-3 0.01
15 12 11 10 8
Actual 2.7e-8 1.1-4 9e-4 0.01 0.3
Predicted 2.4e-4 0.2 2.5 29.5 561
TABLE IV: Comparison of the predicted bound and the actual error for the Barbara image at and .

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


Similarly, we write (9) as , where


Then can be expressed as


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

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

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


  • [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.