I Introduction
The bilateral filter was introduced by Tomasi and Manduchi in [1]
as a nonlinear extension of the classical Gaussian filter. The bilateral filter employs a range kernel along with a spatial kernel for performing edgepreserving 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
(1) 
where is the range kernel and is the spatial kernel. The spatial kernel is usually a Gaussian [1],
(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
(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 constanttime 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 pointwiseconvergent 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
(4) 
where . By setting (4) as the range kernel , we can decompose the numerator in (1) as
where
(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 wellknown, the Gaussian convolution in (5) can be efficiently implemented at constanttime 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 complexconjugate 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 raisedcosines 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 onthefly. 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 userdefined 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 runtime 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 halfinterval . We propose to use the shiftable Fourier basis for this purpose. In particular, we fix some order , and consider the shiftable function
(6) 
where . As is wellknown, 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
(7) 
This is the classical linear leastsquares problem, where the unknowns are . Indeed, using matrixnotation, 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 ,
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 continuouslydifferentiable function, e.g., Gaussian and polynomials. Convergence is also guaranteed for functions that are continuous and piecewisedifferentiable [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 userdefined 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 leastsquares 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 backsubstitution, where is the QRdecomposition 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
(9) 
By construction, for all ,
(10) 
Similar to [18], we consider the (worstcase) error
(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)
(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 runtime 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 runtime 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 raisedcosines [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 runtime of the direct implementation of (1) depends only on the image size and . On the other hand, the runtime of the proposed algorithm depends on , tolerance , image size, and . The fact that the runtime is almost independent of (constanttime 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 
The runtime 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 Diraclike distribution [23]. As is wellknown, 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 runtime. 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 
In Table III, we present the variation of runtime with tolerance for a fixed filter setting. It is seen that the order and hence the runtime 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 subpixel 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.
1e5  1e4  1e3  0.01  0.1  

12  11  10  8  7  
Fast ()  910ms  850ms  780ms  715ms  670ms 
1e8  1e5  1e4  1e3  0.01  
15  12  11  10  8  
Actual  2.7e8  1.14  9e4  0.01  0.3 
Predicted  2.4e4  0.2  2.5  29.5  561 
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
and
Similarly, we write (9) as , where
and
Then can be expressed as
(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 .
References
 [1] C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” Proc. IEEE International Conference on Computer Vision, pp. 839846, 1998.
 [2] E. P. Bennett and L. McMillan, “Video enhancement using perpixel virtual exposures,” ACM Transactions on Graphics, vol. 24, no. 3, pp. 845852, Proc. ACM Siggraph, 2005.
 [3] H. Winnemoller, S. C. Olsen, and B. Gooch, “Realtime video abstraction,” Proc. ACM Siggraph, pp. 12211226, 2006.
 [4] J. Xiao, H. Cheng, H. Sawhney, C. Rao, and M. Isnardi, “Bilateral filteringbased optical flow estimation with occlusion detection,” Proc. European Conference on Computer Vision, pp. 211224, 2006.
 [5] E. P. Bennett and L. McMillan, “Video enhancement using perpixel virtual exposures,” ACM Transactions on Graphics, vol. 24, no. 3, pp 845852, 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. 108112, 2015.
 [7] B. M. Oh, M. Chen, J. Dorsey, and F. Durand, “Imagebased modeling and photo editing,” Proc. Annual Conference on Computer Graphics and Interactive Techniques, pp. 433442, 2001.
 [8] J. Xiao, H. Cheng, H. Sawhney, C. Rao, and M. Isnardi, “Bilateral filteringbased optical flow estimation with occlusion detection,” Proc. European Conference on Computer Vision, pp. 211224, 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. 26902696, 2011.

[10]
K. AlIsmaeil, D. Aouada, B. Ottersten, and B. Mirbach, “Bilateral filter evaluation based on exponential kernels,”
International Conference on Pattern Recognition
, pp. 258261, 2012.  [11] Q. Yang, “Hardwareefficient bilateral filtering for stereo matching,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 36, no. 5, pp.10261032, 2014.
 [12] F. Durand and J. Dorsey. “Fast bilateral filtering for the display of highdynamicrange images,” ACM Transactions on Graphics, vol. 21, no. 3, pp. 257266, 2002.
 [13] Q. Yang, K. H. Tan, and N. Ahuja, “Realtime bilateral filtering,” Proc. IEEE Conference on Computer Vision and Pattern Recognition, pp. 557564, 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. 568580, 2006.
 [15] F. Porikli, “Constant time bilateral filtering,” Proc. IEEE Conference on Computer Vision and Pattern Recognition, pp. 18, 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. 33763382, 2011.
 [17] K. N. Chaudhury, “Fast and accurate bilateral filtering using Gausspolynomial decomposition,” Proc. IEEE International Conference on Image Processing, pp. 20052009, 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. 202206, 2015.
 [19] K. N. Chaudhury, “Constanttime 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. 33573369, 2015.
 [21] R. Deriche, “Recursively implementing the Gaussian and its derivatives, Research Report, INRIA00074778, 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. 12911300, 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.
Comments
There are no comments yet.