Optimal O(1) Bilateral Filter with Arbitrary Spatial and Range Kernels Using Sparse Approximation

Abstract

A number of acceleration schemes for speeding up the time-consuming bilateral filter have been proposed in the literature. Among these techniques, the histogram-based bilateral filter trades the flexibility for achieving O(1) computational complexity using box spatial kernel. A recent study shows that this technique can be leveraged for O(1) bilateral filter with arbitrary spatial and range kernels by linearly combining the results of multiple-box bilateral filters. However, this method requires many box bilateral filters to obtain sufficient accuracy when approximating the bilateral filter with a large spatial kernel. In this paper, we propose approximating arbitrary spatial kernel using a fixed number of boxes. It turns out that the multiple-box spatial kernel can be applied in many O(1) acceleration schemes in addition to the histogram-based one. Experiments on the application to the histogram-based acceleration are presented in this paper. Results show that the proposed method has better accuracy in approximating the bilateral filter with Gaussian spatial kernel, compared with the previous histogram-based methods. Furthermore, the performance of the proposed histogram-based bilateral filter is robust with respect to the parameters of the filter kernel.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Academic Editor:Chung-Hao Chen

College of Mechatronics Engineering and Automation, National University of Defence Technology (NUDT), Changsha, Hunan 410073, China

Received 15 October 2013; Accepted 9 January 2014; 18 February 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The bilateral filter is proposed by Tomasi and Manduchi in [1]. Before this name, it is called SUSAN filter [2]. The general idea of the bilateral filter stems from the earlier work of the neighborhood filter [3] and the sigma filter [4]. The bilateral filter is a noniterative, nonlinear spatial filter which imposes both the geometric affinity and the intensity similarity when aggregating the contribution of each neighboring pixel within the spatial support. Hence the shape of the filter kernel is data adaptive, which enables the filter to preserve abrupt edges when smoothing out small variations. Such an edge-avoiding smoothing property is later demonstrated to be of much potential use in many applications in computer vision and graphics. It has become a general tool in image processing literature, where it is employed in noise cancelation [1, 5-8], high-dynamic-range tone mapping [9-11], image enhancement [12-14], and many other applications [15-17].

On the flip side of its convenient properties for image- and video-based applications, the direct implementation of the standard bilateral filter requires O ( ^{σ s 2} ) operations per pixel, where _{σ s} is the radius of the effective support of the spatial kernel. The computational complexity is too intensive for time-critical applications. Consequently, a plenty of studies on its simplification and acceleration can be found in the literature. Durand and Dorsey proposed a piecewise linear bilateral filter in [9]. They quantize the intensity into several segments and perform FFT-based linear filtering for each segment. The final results are pooled using a linear interpolation on these linearly filtered images. The complexity of Durand's piecewise linear bilateral filter is O ( log ... _{σ s} ) . Pham and van Vliet [18] decomposed the 2D bilateral convolution using two 1D bilateral convolutions and reduced the complexity to O ( _{σ s} ) . Paris et al. [19] and later Chen et al. [11] generalize the idea of the piecewise linear bilateral filter by Durand et al. to form a 3-D bilateral grid. Then the bilateral filtering can be interpreted as a linear convolution of a vector-valued image. Based on an equipollent subsampling in the augmented data space, the complexity of this method is O ( 1 + | ... | / _{σ s σ r} ) , where | ... | denotes the number of grids of the intensity and _{σ r} is the bandwidth of the range kernel. Thus the bilateral grid method runs faster with larger spatial kernel. Weiss [20] develops an O ( log ... _{σ s} ) algorithm for local histogram calculation, which is later used to derive an O ( log ... _{σ s} ) bilateral filter with box spatial kernel. Porikli [21] proposes an O ( 1 ) box bilateral filter by virtue of the O ( 1 ) integral histogram [22]. Porikli [21] and Yang et al. [23] and later Chaudhury et al. [24, 25] suggest the use of some series of the intensity to approximate arbitrary range kernel and employ O ( 1 ) algorithms to efficiently compute the spatial filtering for each term of the series.

Based on Porikli's O ( 1 ) single-box bilateral filter proposed in [21], Gunturk proposes an O ( 1 ) bilateral filter with arbitrary spatial and range kernels by linearly combining multiple single-box bilateral filters [26]. The accuracy can be improved compared with the single-box bilateral filter in approximating the bilateral filter with arbitrary spatial kernel. However, when the radius of the effective spatial support is large, a large number of single-box bilateral filters are required to guarantee the accuracy. Furthermore, no instruction for the selection over all the possible box spatial kernels is presented when the number of single-box bilateral filter is limited. Consequently, Gunturk's method does not produce sufficiently good approximation with only a few additional single-box bilateral filters when the spatial kernel is large. More details about the limitations of Gunturk's method are shown in the following text.

In this paper, we focus our attention on approximating arbitrary spatial kernel with a fixed number of boxes. With this constraint, the MSE-based optimization becomes a well-known sparse representation problem. We then employ the Batch-OMP method published recently in [27] to solve the problem for the radiuses and the coefficients of the boxes. In addition to the application to the histogram-based bilateral filter, the proposed method can be employed in other acceleration schemes by substituting the spatial filtering with a fixed times of box filtering. Hence a number of O ( 1 ) bilateral filters are produced. Afterwards, the application to the histogram-based bilateral filter is discussed with some representative experiments. Results show that the proposed histogram-based bilateral filter is robust and sufficiently accurate over a large range of the size of the spatial kernel and outperforms the other congeneric methods in both the accuracy and the robustness.

The rest of this paper is arranged as follows. Brief descriptions of the original bilateral filter and its acceleration schemes are provided in Section 2, where the previous accelerations are categorized into several types and analyzed, respectively. Then in Section 3, the proposed method for approximating arbitrary spatial kernel with a fixed number of boxes is deduced in detail. The application to some popular acceleration schemes is also provided in this section. Afterwards, the application to the histogram-based bilateral filter is tested in Section 4, together with some comparisons with the previous histogram-based methods. Finally, the conclusions are drawn in Section 5, followed by some further discussions.

2. Related Works

The bilateral filter is proposed by Tomasi and Manduchi in [1]. It is a normalized convolution in which the contribution of a neighboring pixel depends on both the geometric distance and the intensity difference with regard to the center pixel. Generally, arbitrary kernels can be applied to both the spatial and range filtering to measure the affinity between one pixel and its neighbors. The discrete form of the bilateral filter with arbitrary spatial and range kernels can be formulated as follows: [figure omitted; refer to PDF] where v and u are, respectively, the original image and the filtered image and N ( x ) denotes the effective spatial support of the filtering kernel, which is centered at x . _{K s} and _{K r} are, respectively, the spatial and range kernels. The original bilateral filter employs Gaussian kernels for both the spatial and range filtering, which provides intuitive control over the similarity measure with the variances, respectively, for these two kernels. However, the brute-force implementation of the Gaussian bilateral filter is rather time consuming and is thus too slow for real-time applications. Consequently, many research papers have been devoted to accelerating the computation of the bilateral filter, and many publications can be found during the past decade.

2.1. Kernel Separation

Pham and Vliet proposed a method to approximate the 2-D data-adaptive convolution using two 1-D data-adaptive convolutions [18]. The computational complexity is reduced from O ( ^{σ r 2} ) operations per pixel to O ( _{σ r} ) operations per pixel, where _{σ r} denotes the radius of the effective spatial support. This method performs well in regions with simple structure. However, it does not produce satisfactory results in textured regions.

2.2. Piecewise Linear Approximation

Durand and Dorsey proposed a piecewise linear approximation of the bilateral filter in [9]. The authors first quantize the original intensity into several segments. Then for each segment _{i k} , a range-weighted image pair ( _{J k} , _{W k} ) is calculated. Afterwards, a linearly filtered result _{u k} for each _{i k} is generated by linear filtering with the corresponding range-weighted image pair, which can be formulated as follows: [figure omitted; refer to PDF] where _{W k} ( y ) = _{K r} ( _{i k} - v ( y ) ) and _{J k} ( y ) = ^{W k} ( y ) v ( y ) . Then the final result is obtained by linearly interpolating between two closest filtered images; that is, for v ( x ) ∈ [ _{i k} , _{i k + 1} ] , the filtered output at x is given as follows: [figure omitted; refer to PDF]

The computational complexity is dramatically reduced since the implementation requires only O ( log ... _{σ r} ) operations per pixel using fast Fourier transform (FFT) to compute each _{u k} . It is further accelerated by performing a subsampling in the spatial domain.

Recently, based on Durand's formulation, Yang et al. [23] employed Deriche's recursive approximation of Gaussian kernel [28] to calculate the spatial filtering, which can be implemented in constant time and can be realized in real time in some graphic processing units (GPU).

2.3. 3-D Bilateral Grid Manipulations

The piecewise linear bilateral filter is later generalized and further accelerated by Paris et al. in [19]. The authors express the filtering in a higher-dimensional space where another dimension for the intensity is added to the original spatial dimensions. Then the bilateral filtering can be expressed as simple linear convolutions with a vector-valued image in this augmented space, followed by a point-by-point division. With the new representation of the data, simple criteria are derived for subsampling the data to achieve significant acceleration of the bilateral filter. With the same computational time, this method is more accurate than the piecewise linear bilateral filter by Durand and Dorsey in [9].

Built upon the technique by Paris et al. in [19], Chen et al. [11] extend the higher-dimensional approach and introduce a new compact data structure--the bilateral grid. This enables a variety of fast edge-preserving image processing applications. The authors parallelize the bilateral-grid algorithm in modern graphics hardware with graphic processing units (GPU). The GPU-based implementation is two orders of magnitude faster than the equivalent CPU-based implementations. This enables the real-time edge-preserving manipulations on high-definition images. The authors demonstrate the use of the bilateral-grid method on a variety of applications including image editing, high-dynamic-range tone mapping, and image enhancement.

The computational complexity of the bilateral-grid methods is O ( 1 + | ... | / _{σ s σ r} ) operations per pixel, where ... is the dynamic range of the intensity and | ... | is the number of grids occupied by the range dimension. This results in a paradoxical property that the algorithm proceeds faster when the size of the spatial kernel is larger, due to larger subsampling rate. However, the exact output is dependent on the phase of subsampling. Furthermore, the operations with small _{σ s} and _{σ r} require fine subsampling. This requires larger memory and more computation time. Furthermore, the aggregation of the final result relies on a trilinear interpolation among the results on the grids. Better accuracy can be obtained using higher-order interpolations. Consequently, there is a tradeoff between the quality and the computational cost using this type of bilateral filters.

2.4. Histogram-Based Approximation Using Box Spatial Kernel

Weiss developed an iterative method to compute local histogram with O ( log ... _{σ r} ) complexity in [20]. He later demonstrated that such fast local histogram computation can be applied for the O ( log ... _{σ r} ) bilateral filter when the spatial kernel is a uniform box. Similar to Weiss' work, Porikli proposed an O ( 1 ) bilateral filter [21] based on his earlier work on O ( 1 ) integral histogram [22]. The bilateral filter with single-box spatial kernel is formulated as follows: [figure omitted; refer to PDF] where B ( x ) denotes the box with uniform weight centered at x . Then the above formulation can be further expressed in terms of local histograms as follows: [figure omitted; refer to PDF] where P is the number of histogram bin, _{I p} denotes the intensity level of the p th histogram bin, and _{H x} ( _{I p} ) indicates the number of pixels in the p th histogram bin within the spatial neighborhood centered at x . Then the computational time is constant with O ( 1 ) integral histogram technique [22]. And the calculation can be further accelerated by the quantizing the intensity into a small number of histogram bins.

Based on Porikli's single-box bilateral filter, Gunturk proposed a method to approximate the bilateral filter with arbitrary spatial and range kernels using a weighted sum of multiple single-box bilateral filters [26], which can be formulated as follows: [figure omitted; refer to PDF] The coefficients { _{k ^ m} } are then obtained by solving a set of linear equations which minimize the squared error between the discrete forms of [figure omitted; refer to PDF]

This improvement provides better approximation of the bilateral filter with arbitrary spatial and range kernels by costing a little computation time compared with the single-box bilateral filter. However, further study shows that the approximation given by (6) only performs well when the size of spatial kernel is relatively small. As the radius of the spatial support increases, more single-box bilateral filters should be combined to guarantee the approximation accuracy according to Gunturk's scheme. Thus the computational complexity will be unbearable for efficient applications if the size of the spatial kernel is too large.

2.5. Accelerations Based on Series Expansion of the Range Kernel

In [21], Porikli proposes another O ( 1 ) bilateral filter with arbitrary spatial kernel and Gaussian range kernel, which employs Taylor series expansion of the range kernel. Then the bilateral filtering is decomposed into a set of spatial filtering steps with the image series computed prior to the convolution. To guarantee a constant-time processing, the author proposes to subsample the separable 1D linear spatial filters to a constant number of taps asymmetrically. However, the accuracy is poor using a low order Taylor expansion when the variance of Gaussian range kernel is small. Moreover, the range weight blows up when the intensity difference is too large with respect to the variance. Hence, the decaying monotonicity for an admissible kernel is broken.

As an extension of Porikli's polynomial bilateral filter in [21], Chaudhury et al. propose a constant-time bilateral filter with trigonometric range kernels in [24]. Based on the fact that the raised cosine series converges to a Gaussian, the trigonometric range kernel is applied to approximate Gaussian range kernel. The authors show that the bilateral filter using trigonometric series is more accurate than the one using Taylor series with the same number of terms. The computation time of the trigonometric bilateral filter is constant with respect to the variance of the spatial kernel. However, the complexity is O ( ^{| ... | 2} / ^{σ r 2} ) with respect to the variance of the range kernel, where | ... | denotes the number of quantization levels of the intensity and _{σ r} is the variance of the range kernel. Consequently, the bilateral filtering with small _{σ r} will be time consuming using trigonometric range kernel.

More recently in [25], the trigonometric bilateral filter is accelerated using truncations when _{σ r} is small and further extended to a larger class of data-dependent filters including the well-known nonlocal means filter.

3. Fast Implementation of the Bilateral Filter Using Sparse Approximation with Fixed Number of Boxes

In this paper, we propose a method to approximate arbitrary spatial kernel with multiple boxes, which can then be leveraged for the constant-time implementation of the bilateral filter with arbitrary spatial and range kernels. Let L be the radius of the spatial support of the given spatial kernel _{K s} applied in the bilateral filter. Then all the candidate boxes together form a series { _{B l} } , where l is the radius of the box _{B l} and l = 0,1 , 2 , ... , L . For arbitrary _{K s} , it can be approximated using the weighted sum of all the candidate boxes, which is formulated as follows: [figure omitted; refer to PDF]

Because of the symmetry and monotonicity of any admissible spatial kernel, it is possible to find a real positive series { _{k l} } that minimizes the following squared error: [figure omitted; refer to PDF] where the spatial dependency is omitted for simplicity.

3.1. Approximating Arbitrary Spatial Kernel with Fixed Number of Boxes

For that the number of boxes is great when L is large, the computational cost will become unbearable for real-time applications. In order to handle this problem, we add a constraint that the number of boxes used in the approximation should not be larger than a preset number N , such that the constrained minimization of the squared error can be further formulated as follows: [figure omitted; refer to PDF] The left part of the above formulation is again the least squared error norm. And the right part is employed to limit the number of boxes used in the approximation, which denotes the number of boxes that have nonzero coefficients. For any l ∈ [ 0 , L ] , we align the center of the corresponding box with that of _{K s} and pad it with zeros up to the same size as _{K s} . Then the columns of each padded box are concatenated to form a column vector _{b l} . We then put these column vectors together to form a matrix B of size S × ( L + 1 ) , where S = ( 2 L + 1 ^{) 2} is the number of elements in _{K s} . Correspondingly, we concatenate the columns of _{K s} to form a column vector q . By defining k as a column vector containing all the coefficients _{k l} , the optimization problem given in (10) can be reformulated as follows: [figure omitted; refer to PDF] where _{|| k || 0} is the _{L 0} norm of the vector k , which denotes the number of non-zero elements in k .

One should be very familiar with the formulation given by (11), since it is a well-known sparse representation problem. This formulation searches the sparse representation of the signal q using the known dictionary B . This sparse approximation problem, which is known to be NP-hard, can be efficiently solved using several available techniques, including Orthogonal Matching Pursuit (OMP) [27, 29, 30], Basis Pursuit (BP) [31, 32], and FOCUSS [33].

In this paper, we employ the efficient OMP algorithm by Rubinstein et al. in [27], namely, Batch-OMP, to solve (11). After we find out the radiuses and coefficients of the boxes corresponding to the given spatial kernel and the given number of boxes, the further computation of the bilateral filter can be given as follows: [figure omitted; refer to PDF] With this respect, we can benefit from the integral map technique to efficiently calculate the spatial convolution in constant time. It is worth mentioning that N in (10) and (11) is an input parameter defined by the users. For a given _{σ s} , the more the number of boxes used is, the better the approximation accuracy can be achieved. But more computational time is needed for more boxes. Thus there is a tradeoff between the approximation accuracy and the computational efficiency with our method. Generally, we can find a minimum N that guarantees the acceptable quality of the approximation, which will be discussed in detail in the next section.

3.2. Application to the Histogram-Based Bilateral Filters

Now having that an arbitrary spatial kernel is approximated with a preset number of boxes, we first exploit its use in the histogram-based bilateral filters, which is originally developed by Weiss in [20], later accelerated by Porikli in [21], and recently extended by Gunturk in [26].

Firstly, let us denote the solution of (11) using the notation { _{k n} } , which contains N non-zero coefficients for N boxes with different radiuses selected by the Batch-OMP algorithm. The subscripts n 's now are not the radiuses but the indices of the selected boxes; that is, n = 1,2 , ... , N . Then the proposed bilateral filter with multiple-box spatial kernel given in (12) can be reformulated as follows: [figure omitted; refer to PDF]

Following the single-box histogram-based bilateral filter by Weiss [20] and Porikli [21], (13) can be further formulated in terms of local histogram as follows: [figure omitted; refer to PDF]

Prior to the filtering process, the intensity is quantized into P bins, resulting in an intensity series { _{I p} ^{} 1 P} . Then the corresponding P frames of range weight { _{K r} ( v ( x ) - _{I p} ) } are calculated. Based on the work in [22], an integral histogram is established. Then N local histogram maps are calculated in constant time using the same integral histogram. Compared with the single-box bilateral filter, the proposed method requires N - 1 additional local histogram calculations, which can be efficiently obtained using the integral histogram technique. The calculations of N times of box filtering for both the numerator and dominator in (14) can be parallelized for more efficient implementation.

3.3. Application to the Interpolation-Based Bilateral Filters

Since Durand's piecewise linear bilateral filter [9], Yang's extension [23], and also Paris' method [19] employ an interpolation in pooling the filtering results, we call them together the interpolation-based bilateral filters for convenience. As the origin of this type of bilateral filters, Durand's method is discussed in this paper for demonstrating the application of the proposed approximation to this type of bilateral filters. Similar adaptation can be easily made for the application to Yang's method.

Let us replace the spatial kernel _{K s} in (2) by Durand and Dorsey using the weighted sum of multiple boxes, and the calculation of _{u k} is formulated as follows: [figure omitted; refer to PDF]

Given an image, it is first divided into a preset number of intensity segments, forming an image series { _{i k} } . Correspondingly, a range-weighted image series { _{J k} ( y ) } and a weight series { _{W k} ( y ) } are computed. Afterwards, the corresponding integral maps of these two image series are established. The filtered result _{u k} for each image segment _{i k} is computed according to (15) with O ( 1 ) complexity. The final result is obtained using an interpolation, as given in (3) for an example. The running time of this bilateral filter is constant with respect to the size of the spatial kernel.

3.4. Applications to the Series-Based Bilateral Filters

Since a convolution with spatial kernels is also involved in the series-based bilateral filters, the proposed sparse approximation is also possible to be applied in this type of accelerations. The detailed adaptation for the series-based bilateral filters is omitted in this paper, which can be easily derived by the readers.

4. Experiments and Comparisons

The proposed algorithm is implemented in MATLAB 2010b under Windows XP SP3 (32 bits), in a PC computer with Pentium Dual-Core T4200 @ 2.00 GHz each and 3 GB RAM. The proposed histogram-based bilateral filter with multiple-box spatial kernel is tested and discussed in this section. The experiments on approximating the bilateral filter with Gaussian spatial and range kernels provide the readers with a glance at the validity and the effectiveness of the proposed method. The application to other types of the acceleration of the bilateral filter can be easily implemented and evaluated under the guidance provided in the previous section. But it is not given here for saving the length of this paper. In the following text of this section, we evaluate the performance of the proposed histogram-based bilateral filter in both qualitative and quantitative aspects.

Before the evaluation, we briefly introduce the basic settings of the experiments discussed in the following text. Four standard test images are employed, which are 8 bit grayscale images of size 512 × 512 , as shown in Figure 1. The numerical results given in the following texts are all produced using the standard test images, unless otherwise noticed.

Four standard test images used in the experiments in this paper. They are 8-bit grayscale images of size 512 × 512 .

(a) Barbara

[figure omitted; refer to PDF]

(b) Boat

[figure omitted; refer to PDF]

(c) Lena

[figure omitted; refer to PDF]

(d) Pirate

[figure omitted; refer to PDF]

4.1. Qualitative Demonstration

For the qualitative demonstration, the proposed histogram-based method, together with Porikli's single-box bilateral filter and Gunturk's bilateral filter, is employed to approximate the standard bilateral filter with Gaussian spatial and range kernels. The results are illustrated in Figure 2. As shown in Figure 2, the images in the first, third, and fifth columns are, respectively, the results of the proposed histogram-based bilateral filter, the results of Porikli's single-box bilateral filter [21], and the results of Gunturk's bilateral filter [26]. The color-coded maps in the second, fourth, and sixth columns are the corresponding absolute error between the images, respectively, in the first, third, and fifth columns and the results of the standard Gaussian bilateral filter, respectively, with the same input parameters. The variances of the spatial kernels are _{σ s} = 1.2,1.8,3.0 , and 6.0 , respectively, for the four rows. The variances of the range kernels are the same; that is, _{σ r} = 50 . The number of boxes for the proposed method is N = 5 in producing the images in the first column. Correspondingly, the number of box bilateral filters used in Gunturk's method is M = 5 in producing the results in the fifth column. For Porikli's method, the results are obtained with optimal spatial parameter derived by Gunturk in [26]; that is, r = [ 1.4 _{σ s} ] . For producing all the results shown in Figure 2, P = 16 . The ranges of the absolute error maps shown here are uniformly limited to [ 0,50 ] for better illustrating the difference of the performance of the three methods.

Figure 2: Barbara: a demonstration of the visual quality of the approximation results. First column: results obtained using the proposed method, where N = 5 . Second column: the color-coded absolute error maps corresponding to the images in the first column. Third column: results obtained using Porikli's single-box bilateral filter with optimal spatial parameter. Fourth column: the color-coded absolute error maps corresponding to the images in the third column. Fifth column: results obtained using Gunturk's bilateral filter, where M = 5 . Last column: the color-coded absolute error maps corresponding to the images in the fifth column. First row: images are obtained with _{σ s} = 1.2 . Second row: images are obtained with _{σ s} = 1.8 . Third row: images are obtained with _{σ s} = 3.0 . Fourth row: images are obtained with _{σ s} = 6.0 . For all these results, _{σ r} = 50 and P = 16 . The ranges of the color-coded maps are uniformly limited to [ 0,50 ] to facilitate visual comparison.

[figure omitted; refer to PDF]

It is shown in Figure 2 that the proposed histogram-based bilateral filter produces less error than the other two methods when varying _{σ s} . Hence, better accuracy is achieved in approximating the standard bilateral filter. Note that the error maps in the fourth column show that Porikli's method with optimal box performs better than Gunturk's method when _{σ s} is sufficiently large. Especially noted in the fifth column, Gunturk's method performs well with small _{σ s} , but the result images gradually become darker as _{σ s} increases with fixed M . The performance of Gunturk's bilateral filter can be improved by increasing M . However, the computational cost will increase as well. The numerical results on this experiment are also included in Tables 1 and 2.

Table 1: Comparison of the PSNR values.

σ s	Barbara	Boat	Lena	Pirate
Ours	[21]	[26]	Ours	[21]	[26]	Ours	[21]	[26]	Ours	[21]	[26]
0.90	41.90	40.34	41.62	43.15	41.89	42.92	41.45	40.89	41.43	42.68	41.70	42.60
1.20	43.28	39.21	43.10	44.70	40.10	44.37	42.58	41.06	42.53	44.15	40.83	44.04
1.50	44.27	41.68	43.50	45.70	43.26	43.92	43.35	42.61	42.85	45.17	43.52	44.36
1.80	45.05	41.98	37.89	46.44	40.95	37.08	43.90	42.01	37.61	45.90	41.98	38.48
2.40	46.14	43.89	24.81	47.35	44.05	24.11	44.75	43.37	24.55	46.91	44.27	25.35
3.00	46.84	44.79	18.25	47.85	44.29	17.60	45.51	44.10	17.97	47.69	44.94	18.78
3.60	47.51	45.28	14.56	48.41	44.39	13.93	46.09	44.57	14.27	48.36	45.25	15.09
4.50	48.07	45.00	11.52	48.85	44.43	10.90	46.75	44.53	11.23	48.93	44.93	12.05
6.00	48.88	44.74	9.16	49.33	44.40	8.53	47.58	44.52	8.86	49.66	44.56	9.69
7.50	49.34	44.25	8.06	49.62	44.11	7.42	48.19	44.47	7.75	50.21	44.25	8.59
9.00	49.67	44.89	7.46	49.92	44.27	6.83	48.61	45.22	7.16	50.51	45.23	7.99

Table 2: Comparison of the FSIM values.

σ s	Barbara	Boat	Lena	Pirate
Ours	[21]	[26]	Ours	[21]	[26]	Ours	[21]	[26]	Ours	[21]	[26]
0.90	0.9874	0.9859	0.9873	0.9928	0.9921	0.9927	0.9830	0.9810	0.9829	0.9911	0.9902	0.9911
1.20	0.9902	0.9893	0.9902	0.9943	0.9911	0.9941	0.9862	0.9864	0.9862	0.9931	0.9908	0.9930
1.50	0.9924	0.9914	0.9924	0.9954	0.9949	0.9952	0.9888	0.9879	0.9889	0.9946	0.9940	0.9946
1.80	0.9939	0.9922	0.9939	0.9962	0.9901	0.9959	0.9907	0.9896	0.9909	0.9956	0.9913	0.9956
2.40	0.9960	0.9942	0.9929	0.9972	0.9949	0.9941	0.9934	0.9911	0.9910	0.9969	0.9948	0.9944
3.00	0.9973	0.9953	0.9817	0.9978	0.9941	0.9834	0.9955	0.9931	0.9820	0.9977	0.9950	0.9849
3.60	0.9979	0.9956	0.9589	0.9981	0.9931	0.9629	0.9964	0.9939	0.9629	0.9982	0.9947	0.9658
4.50	0.9984	0.9954	0.9137	0.9984	0.9927	0.9235	0.9973	0.9939	0.9248	0.9985	0.9944	0.9280
6.00	0.9988	0.9956	0.8427	0.9985	0.9931	0.8631	0.9979	0.9940	0.8660	0.9988	0.9942	0.8697
7.50	0.9991	0.9956	0.7928	0.9987	0.9932	0.8198	0.9984	0.9943	0.8246	0.9990	0.9944	0.8287
9.00	0.9992	0.9962	0.7602	0.9989	0.9935	0.7916	0.9987	0.9958	0.7979	0.9990	0.9952	0.8021

4.2. Quantitative Evaluation

We quantitatively evaluate the proposed method in three aspects, including the running time on our system with nonoptimized MATLAB codes, the quality in approximating the standard bilateral filter according to some previously established criteria, and the minimum number of boxes to guarantee sufficient accuracy compared with the minimum number of box bilateral filters in Gunturk's method.

4.2.1. Computational Efficiency

We test our method by varying values of N and P and record the running time. Experiments show that the proposed method costs about 2.36 s for a 512 × 512 grayscale image when N = 5 and P = 16 . And every additional histogram bin costs about 0.04 s when N = 5 , and every additional box costs about 0.36 s when P = 16 . The local histogram calculation costs 0.34 s each time on our system, which is the most time-consuming part of the nonoptimized codes for every additional box. Further acceleration of the local histogram calculation should be done for more efficient implementation of the proposed method.

Experiments show that the running time of the proposed method is about 0.1 s more than that of Gunturk's method with N = M and P = 16 . It is because the sparse approximation algorithm in the proposed method costs about 0.1 s more than the linear problem solver in Gunturk's algorithm. However, the computation of any additional box costs as much time as an additional single-box bilateral filter. As a result, the running time of the rest codes of the two algorithms is almost the same.

4.2.2. Approximation Accuracy

The two previously established criteria are employed in quantitatively evaluating the approximation quality in this paper. Following the previous work, the PSNR is employed. Recently, Zhang et al. propose a feature similarity index (FSIM) in [34] for image quality assessment with reference. It is one of the state-of-the-art methods for predicting the similarity between the test and reference images. In order to provide a more comprehensive evaluation, the FSIM is also applied for evaluating the quality of the results.

Given the result images of the standard bilateral filter with different parameter settings, the PSNR and FSIM values of the results by the methods to be evaluated are calculated. The PSNR and FSIM values on four standard test images of the proposed histogram-based method together with that of Porikli's and Gunturk's methods are, respectively, listed in Tables 1 and 2, where _{σ r} = 50 , P = 16 , and N = M = 5 . The best scores among the three methods for each individual _{σ s} and each test image are shown in boldface. We can see from these two tables that our method outperforms the other two according to both criteria. And it is more robust with respect to the change of _{σ s} with a fixed N . Besides, Gunturk's bilateral filter performs well when _{σ s} is small. However, both PSNR and FSIM values of Gunturk's method decrease quickly as _{σ s} increases when M is fixed. It is because the result images become darker as _{σ s} increases for a fixed M , as can be illustrated in Figure 2. Larger M is required for better quality in approximating Gaussian bilateral filter with larger _{σ s} using Gunturk's method.

In order to provide a more intuitive illustration, we plot the PSNR and FSIM values on "Lena," respectively, in Figures 3 and 4 with respect to _{σ s} . As can be seen in these two figures, the proposed histogram-based bilateral filter has higher PSNR and FSIM values compared with the other two methods over a large range of _{σ s} . Besides, the value of N does not affect the performance too much. As a result, only a small number of boxes are sufficient for a vast range of _{σ s} using the proposed method. On the contrary, for a fixed M , the scores given by both principles of Gunturk's method decay quickly as _{σ s} increases. Using larger M is helpful in producing satisfactory results; however, as discussed before, the computational efficiency will be degraded as a consequence. Although small _{σ s} is sufficient for most denoising applications, large _{σ s} is usually required in many other applications including high-dynamic-range tone mapping [9], image enhancement [14], stereo vision [17], and dehazing [35]. Consequently, the computational cost will be unbearable for efficient implementation combining a large number of single-box bilateral filters.

The PSNR values of the results with respect to _{σ s} . These results are calculated on "Lena" in approximating the Gaussian bilateral filter, where P = 16 .

(a) _{σ r} = 25

[figure omitted; refer to PDF]

(b) _{σ r} = 50

[figure omitted; refer to PDF]

(c) _{σ r} = 100

[figure omitted; refer to PDF]

The FSIM values of the results with respect to _{σ s} . These results are calculated on "Lena" in approximating the Gaussian bilateral filter, where P = 16 .

(a) _{σ r} = 25

[figure omitted; refer to PDF]

(b) _{σ r} = 50

[figure omitted; refer to PDF]

(c) _{σ r} = 100

[figure omitted; refer to PDF]

Furthermore, as can be seen in Figures 3 and 4, our histogram-based bilateral filter is apparently more robust than the other two methods with respect to _{σ r} . On the contrary, as the value of _{σ r} increases, the performance of Porikli's single-box bilateral filter becomes less stable with respect to the change of _{σ s} .

4.2.3. Minimum Number of Box

For Gunturk's method, larger M leads to better approximation accuracy. Thus sufficient single-box bilateral filters should be combined to guarantee an acceptable accuracy. Correspondingly, as mentioned in the previous section, sufficient boxes should also be used in the proposed histogram-based bilateral filter. For the acceptable quality threshold, it is often considered as visually undistinguishable when PSNR ...5; 40 dB [21]. We experiment on a standard test image, namely, "Lena," shown as Figure 1(c). In order to find out the minimum N for the proposed histogram-based bilateral filter and the corresponding minimum M for Gunturk's method with respect to the scale of the spatial kernel, we vary the value of _{σ s} and keep other parameters unchanged. For each given _{σ s} , we gradually increase both numbers and calculate the PSNR values of the results, with respect to the results of the standard Gaussian bilateral filter. Then the minimum numbers which enable the quality to exceed the acceptable threshold for both algorithms can be obtained. The minimum N and the minimum M with respect to _{σ s} are shown in Figure 5. It is clearly seen that the minimum M for Gunturk's method increases as _{σ s} increases. On the contrary, the minimum N for the proposed method is much smaller than the minimum M for Gunturk's method. And it is nearly constant with respect to _{σ s} . P = 16 and _{σ r} = 50 for both methods in this experiment. It is shown that the proposed histogram-based bilateral filter is more efficient to obtain sufficient accuracy, and the performance is more robust to the change of the spatial scale of the filter.

Figure 5: Minimum number of boxes required to guarantee PSNR > 40 dB, where P = 16 and _{σ r} = 50 .

[figure omitted; refer to PDF]

Finally, it is worth noting that, when N = 1 , the proposed histogram-based bilateral filter is exactly the same as Porikli's single-box bilateral filter with optimal spatial parameter.

5. Conclusions and Discussions

We have presented a method for approximating arbitrary spatial kernel using a fixed number of boxes based on sparse approximation techniques. With applications to the acceleration of the bilateral filter, the proposed method can be leveraged for a broad class of constant-time bilateral filters with arbitrary spatial and range kernels. Once the parameter of the spatial kernel and the number of boxes are given, the radiuses and the coefficients of the boxes are determined by the Batch-OMP algorithm. Hence they can be computed ahead of the filtering process. The application to the histogram-based O ( 1 ) bilateral filter is demonstrated in this paper, followed by a number of convictive experiments. Results tell that the proposed histogram-based bilateral filter has better accuracy compared with other histogram-based bilateral filters in approximating the standard bilateral filter. Meanwhile, the performance of the proposed histogram-based bilateral filter is robust with respect to the parameters of the filter kernels. And a small number of boxes are sufficient to achieve satisfactory accuracy for a large range of the spatial parameters. With a very little more computational cost, the accuracy of the proposed histogram-based bilateral filter is better than the previous histogram-based ones.

The proposed O ( 1 ) bilateral filter with arbitrary spatial and range kernels can be parallelized for real-time implementation in GPUs. Moreover, it is possible for the proposed method to be extended to the bilateral grid. 3-D cubic boxes can be used to approximate 3-D arbitrary spherical filter kernels in the augmented data space. With the proposed method, higher-dimensional manipulations are also possible in this direction.

Acknowledgments

The authors would like to thank Dr. Jian Li for his insightful discussion over the mathematical deductions and the implementation of the algorithms. This work is supported by the National Natural Science Foundation of China, under Grant 90820302.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

References

[1] C. Tomasi, R. Manduchi, "Bilateral filtering for gray and color images," in Proceedings of the 6th IEEE International Conference on Computer Vision, pp. 839-846, January 1998.

[2] S. M. Smith, J. M. Brady, "SUSAN--a new approach to low level image processing," International Journal of Computer Vision , vol. 23, no. 1, pp. 45-78, 1997.

[3] L. P. Yaroslavsky Digital Picture Processing--An Introduction , Springer, Berlin, Germany, 1985.

[4] J.-S. Lee, "Digital image smoothing and the sigma filter," Computer Vision, Graphics and Image Processing , vol. 24, no. 2, pp. 255-269, 1983.

[5] M. Elad, "On the origin of the bilateral filter and ways to improve it," IEEE Transactions on Image Processing , vol. 11, no. 10, pp. 1141-1151, 2002.

[6] E. P. Bennett, J. L. Mason, L. McMillan, "Multispectral bilateral video fusion," IEEE Transactions on Image Processing , vol. 16, no. 5, pp. 1185-1194, 2007.

[7] M. Zhang, B. K. Gunturk, "Multiresolution bilateral filtering for image denoising," IEEE Transactions on Image Processing , vol. 17, no. 12, pp. 2324-2333, 2008.

[8] A. Wong, "Adaptive bilateral filtering of image signals using local phase characteristics," Signal Processing , vol. 88, no. 6, pp. 1615-1619, 2008.

[9] F. Durand, J. Dorsey, "Fast bilateral filtering for the display of high-dynamic-range images," in Proceedings of the 29th Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH '02), pp. 257-266, San Antonio, Tex, USA, July 2002.

[10] S. Bae, S. Paris, F. Durand, "Two-scale tone management for photographic look," in Proceedings of the ACM SIGGRAPH (SIGGRAPH '06), pp. 637-645, Boston, Mass, USA, August 2006.

[11] J. Chen, S. Paris, F. Durand, "Real-time edge-aware image processing with the bilateral grid," ACM Transactions on Graphics , vol. 26, no. 3, pp. 1-9, 2007.

[12] E. Eisemann, F. Durand, "Flash photography enhancement via intrinsic relighting," in Proceedings of the ACM SIGGRAPH, pp. 673-678, Los Angeles, Calif, USA, August 2004.

[13] M. Elad, "Retinex by two bilateral filters," in Proceedings of Scale-Space and PDE Methods in Computer Vision, vol. 3459, pp. 217-229, Hofgeismar, Germany, April 2005.

[14] B. Zhang, J. P. Allebach, "Adaptive bilateral filter for sharpness enhancement and noise removal," IEEE Transactions on Image Processing , vol. 17, no. 5, pp. 664-678, 2008.

[15] R. Ramanath, W. E. Snyder, "Adaptive demosaicking," Journal of Electronic Imaging , vol. 12, no. 4, pp. 633-642, 2003.

[16] J. Xiao, H. Cheng, H. Sawhney, C. Rao, M. Isnardi, "Bilateral filtering-based optical flow estimation with occlusion detection," in Proceedings of the European Conference on Computer Vision (ECCV '06), pp. 211-224, May 2006.

[17] Q. Yang, L. Wang, R. Yang, H. Stewénius, D. Nistér, "Stereo matching with color-weighted correlation, hierarchical belief propagation, and occlusion handling," IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 31, no. 3, pp. 492-504, 2009.

[18] T. Q. Pham, L. J. van Vliet, "Separable bilateral filtering for fast video preprocessing," in Proceedings of the IEEE International Conference on Multimedia and Expo (ICME '05), pp. 1-4, Amsterdam, Netherlands, July 2005.

[19] S. Paris, F. Dur, A. fast, "A fast approximation of the bilateral filter using a signal processing approach," in Proceedings of the European Conference on Computer Vision (ECCV '06), pp. 829-836, Graz, Austria, May 2006.

[20] B. Weiss, "Fast median and bilateral filtering," in Proceedings of the ACM SIGGRAPH, pp. 519-526, Boston, Mass, USA, August 2006.

[21] F. Porikli, "Constant time O (1) bilateral filtering," in Proceedings of the IEEE International Conference on Computer Vision and Pattern Recognition (CVPR '08), pp. 1-8, Anchorage, Alaska, USA, June 2008.

[22] F. Porikli, "Integral histogram: a fast way to extract histograms in Cartesian spaces," in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR '05), pp. 829-836, San Diego, Calif, USA, June 2005.

[23] Q. Yang, K.-H. Tan, N. Ahuja, "Real-time O (1) bilateral filtering," in Proceedings of IEEE International Conference on Computer Vision and Pattern Recognition (CVPR '09), pp. 557-564, June 2009.

[24] K. N. Chaudhury, D. Sage, M. Unser, "Fast O (1) bilateral filtering using trigonometric range kernels," IEEE Transactions on Image Processing , vol. 20, no. 12, pp. 3376-3382, 2011.

[25] K. N. Chaudhury, "Acceleration of the shiftable O (1) algorithm for bilateral filtering and nonlocal means," IEEE Transactions on Image Processing , vol. 22, no. 4, pp. 1291-1300, 2013.

[26] 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.

[27] R. Rubinstein, M. Zibulevsky, M. Elad, "Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit,", no. CS-2008-08, Technion-Computer Science Department, 2008.

[28] R. Deriche, "Recursively implementing the gaussian and its derivatives," in Proceedings of IEEE International Conference on Image Processing (ICIP '92), pp. 263-267, 1992.

[29] Y. C. Pati, R. Rezaiifar, P. S. Krishnaprasad, "Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition," in Proceedings of the 27th Asilomar Conference on Signals, Systems & Computers, pp. 40-44, November 1993.

[30] G. Davis, S. Mallat, M. Avellaneda, "Adaptive greedy approximations," Constructive Approximation , vol. 13, no. 1, pp. 57-98, 1997.

[31] S. S. Chen, D. L. Donoho, M. A. Saunders, "Atomic decomposition by basis pursuit," SIAM Review , vol. 43, no. 1, pp. 129-159, 2001.

[32] D. L. Donoho, M. Elad, "Optimally sparse representation in general (nonorthogonal) dictionaries via ^{L 1} minimization," Proceedings of the National Academy of Sciences of the United States of America , vol. 100, no. 5, pp. 2197-2202, 2003.

[33] I. F. Gorodnitsky, B. D. Rao, "Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm," IEEE Transactions on Signal Processing , vol. 45, no. 3, pp. 600-616, 1997.

[34] L. Zhang, L. Zhang, X. Mou, D. Zhang, "FSIM: a feature similarity index for image quality assessment," IEEE Transactions on Image Processing , vol. 20, no. 8, pp. 2378-2386, 2011.

[35] P. Carr, R. Hartley, "Improved single image dehazing using geometry," in Proceedings of the Digital Image Computing: Techniques and Applications (DICTA '09), pp. 103-110, December 2009.