{'fraction_non_alphanumeric': 0.05126424006668519, 'fraction_numerical': 0.016879688802445125, 'mean_word_length': 4.306671581275341, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimized Fourier Bilateral Filtering

We consider the problem of approximating a truncated Gaussian kernel using Fourier (trigonometric) functions. The computation-intensive bilateral filter can be expressed using fast convolutions by applying such an approximation to its range kernel, where the truncation in question is the dynamic range of the input image. The error from such an approximation depends on the period, the number of sinusoids, and the coefficient of each sinusoid. For a fixed period, we recently proposed a model for optimizing the coefficients using least-squares fitting. Following the Compressive Bilateral Filter (CBF), we demonstrate that the approximation can be improved by taking the period into account during the optimization. The accuracy of the resulting filtering is found to be at least as good as CBF, but significantly better for certain cases. The proposed approximation can also be used for non-Gaussian kernels, and it comes with guarantees on the filtering accuracy.


I. INTRODUCTION
The bilateral filter is popularly used in image processing and computer vision for edge-preserving smoothing [1], [2]. Unlike classical convolutional filters, it uses an additional kernel for measuring proximity in range (intensity) space. In particular, the bilateral filtering of an image f is given by where ω and ϕ are the spatial and range kernels. Both kernels are typically Gaussian [1]: where θ and σ are the respective standard deviations. Notice that the range kernel acts on the intensity difference between the pixel of interest i and its neighbor i−j. If |f (i)− f (i − j)| 3σ (e.g. the pixels are on opposite sides of an edge), then the weight assigned to pixel i − j is small, whereby it is excluded from the aggregation. This helps in preserving sharp edges [1]. However, the range kernel also makes the filter non-linear and computation intensive. In particular, the bruteforce computation of (1) requires 2(6θ + 1) 2 operations per pixel. This makes the real-time implementation challenging for practical values of θ. Several fast approximations [3]- [18] have been proposed to accelerate the brute-force computation of (1). The present focus is on approximations that use trigonometric functions [8], [10], [14], [18]. We briefly describe their relation to the present work and summarize our contributions.
Previous work. The idea of fast bilateral filtering using trigonometric functions was first proposed in [8]. It was observed that, by approximating ϕ using trigonometric polynomials, we can compute (1) using fast convolutions. This was later refined and extended in [10]- [18]. In particular, state-ofthe-art results were obtained using Fourier series approximation of a truncated ϕ in [10], where the truncation is simply the intensity range of the input image. An important observation from [10], which was overlooked in prior work [14], is that the period of the sinusoids plays a critical role in the approximation. In [14], the period was set as [−R, R], where R is the intensity range. While this works well when σ is small, it results in poor approximations when 3σ R. The kernel does not flatten out sufficiently over [−R, R] in such cases. This induces (higherorder) discontinuities at the boundaries, causing the coefficients to decay slowly (e.g. see Figure 4 in [10]).
Contributions. In this work, we combine some of the ideas from [10] and [14] to develop an approximation with improved filtering quality. By filtering quality, we mean the error between (1) and the output obtained using the kernel approximation in question. The highlights of our approximation model and its key differences with [10] are as follows: • Unlike [10], we do not use a continuous approximation of ϕ over the interval [−R, R]. Instead, we just approximate the discrete samples ϕ(t) corresponding to the intensity levels t ∈ {−R, . . . , R}. Notice that it is precisely these samples that appear in (1) and (2). As a result, we can control the filtering quality (see Proposition 1) by adjusting the kernel error at will.
• In [10], the Gaussian kernel is approximated using Fourier series. However, since there is no known analytical formula for the Fourier coefficients in this case, the coefficients are further approximated using the Fourier transform of a Gaussian (cf. [10,Lemma 8]). In contrast, the coefficients in our model are computed exactly using least-squares optimization.
• The approximation error for our method provably decreases with the increase in the number of sinusoids (K), and eventually vanishes (see Proposition 2). Such a simple guarantee is not offered by [10]. This is because the period changes with K in [10]. Therefore, we cannot use standard convergence guarantees from Fourier analysis (where the period is assumed to be fixed).
• The approximation in [10] is specialized for Gaussians, whereas our approximation can be used for any range kernel. As is well-known, it is rather difficult to analytically compute the Fourier coefficients (or the Fourier transform) of an arbitrary kernel. As a result, it is difficult to use the integration-based method in [10] for computing the coefficients for non-Gaussian kernels [3], [11], [12]. In other words, the issue with [10] for non-Gaussian kernels is the fast and accurate computation of Fourier coefficients.
Organization. The rest of the paper is organized as follows.  The proposed approximation, along with some numerical results, is explained in Section II. In Section III, we compare our filtering results with [10], [14], both visually and in terms of PSNR. We summarize the results in Section IV.

II. PARAMETER OPTIMIZATION
Following [10], consider the K-term Fourier series approximation of ϕ over the period [−T, T ]: Since the original kernel ϕ is symmetric, we consider only the cosine terms in (3). Notice that, unlike the standard definition, we divide by 2T + 1 and not by 2T in (3). The reason for this deviation will be made precise later. The important observation is that, by usingφ in place of ϕ, we can express the numerator and denominator of (1) using spatial convolutions; see Section II in [14] for a detailed account. This observation is at the heart of the fast algorithms in [10], [14]. For fixed σ, the design problem is that of fixing the intrinsic parameters K, T and (c k ). More precisely, since the number of convolutions required in the fast algorithm is proportional to K, the goal is to find the smallest K (and the corresponding T and (c k )) such that the error between ϕ andφ is within a specified tolerance [10], [14].
A key difference between the proposals in [10], [14] is the definition of approximation error. Notice that the argument t in ϕ(t) assumes the values The error (up to a normalization) was defined in [10] However, notice that the domain of t is not the full interval [−R, R], but rather just the integers Λ = {−R, . . . , R}. Based on this observation, the following error was considered in [14]: We choose to work with (4) for couple of reasons. First, it was shown in [14] that a bound on (4) automatically translates into a bound on the filtering accuracy.
The other point is that (4) is a quadratic function of the coefficients. Consider the vector b of length 2R + 1 consisting of the samples {ϕ(t) : t ∈ Λ}, and matrix A of size (2R+1)× K whose columns are the discretized sinusoids in (3). That is, b(i) = ϕ(i−R−1) and A(i, j) = cos(ν(i−R−1)(j −1)) for i = 1, . . . , 2R + 1 and j = 1, . . . , K, where ν = 2π/(2T + 1). We can then simply write (4) as Ac − b 2 , where c = (c k ). Importantly, for fixed K and T , we can exactly minimize (4) with respect to c using linear algebra. In particular, let which is the smallest possible error for fixed K and T . Notice that the size of A and its components depend on K and T . The next question is how does E(K, T ) behave with T for some fixed K? Intuitively, it is clear that the error is large if T is too small or too large with respect to σ. If T is too large, then we see from (3) that the sinusoids effectively degenerate to constant functions over [−3σ, 3σ]. On the other hand, as mentioned previously, if T 3σ, then the Gaussian does not flatten out sufficiently within the period [−T, T ]. This induces a discontinuity at the boundary after the periodization, which causes the Fourier coefficients to decay slowly. This adversely affects the approximation for a fixed K. We noticed that this problem persists even after optimizing the coefficients (see Figure 1). Thus, following [10], we propose to optimize (5) with T . In particular, let Given some user-defined tolerance ε > 0, the goal is to find the smallest K such that e(K) ≤ ε. The existence of such a K is guaranteed by the following observation (see supplement for the proof; this is exactly where we require the division by 2T + 1 in (3)). Proposition 2: The error given by (6) is non-increasing in K, i.e., e(K + 1) ≤ e(K), and it vanishes when K = 2R + 1.
The overall optimization procedure is summarized in Algorithm 1, where the optimal order and period are denoted by K and T . Proposition 2 guarantees that the "condition e ≤ ε" in line 2 is satisfied for some K. As with the example in Figure  Unlike [10], [14], we propose to perform the optimization offline for practical values of σ and ε, and store the corresponding K and T in a lookup table. At run time, we simply need to read K and T from the table, and perform the leastsquares optimization in (5) to get the optimal coefficients. We empirically found that for fixed σ, K and T scale almost linearly with log(1/ε). Whereas, for fixed ε, K (resp. T ) decreases (resp. increases) almost linearly with σ. Additional plots are provided in the supplement. As shown in Figure 2, K and T varies smoothly with σ and log(1/ε), and hence the off-grid values can accurately be estimated using bilinear interpolation. The loss in filtering accuracy owing to the suboptimality of the interpolated K and T is at most 1-2 dB. For reference, the values of K and T for different values of σ and ε are provided in the file LUT.mat in the supplementary material. Of course, we can use a large table with more entries depending on the application at hand. On the other hand, for hardware implementations of the filter, the values of σ and ε would typically be hardcoded [19], and a lookup table would not be required.

III. RESULTS AND COMPARISONS
The simulations in this section were performed using Matlab on a 3.40 GHz quad-core machine with 32 GB memory. We used the Matlab code of [10], [14] provided by the authors. The Gaussian convolutions involved in the fast algorithm were implemented using the Matlab routine "imfilter". We have used 8-bit grayscale images [20], for which R = 255.    (b) Error between the target kernel and its approximation. The order is K = 4 in both cases; for our method, we tuned ε in Algorithm 1 to fix K.
We first compare the proposed approximation with Fourier Bilateral Filtering (FBF) [14]. The order for the former was fixed by adjusting ε in Algorithm 1; the same order was then used in FBF. A particular result is shown in Figure 3, where we compare different approximations (of identical order) with the target kernel. Notice that our approximation is much better than that of FBF. In particular, by optimizing the period, we are able to suppress the oscillation on the tail appearing in the FBF kernel. The approximation around the origin is also better for our method. In Figure 4, we show the decay of kernel error with the approximation order for both methods. As expected, the error is consistently lower with our method for different values of σ.
Following existing works [3], [7], [13], we measure the filtering accuracy using the peak-signal-to-noise ratio: PSNR = 10 log 10 (255 2 /MSE), where MSE is the mean-squared error between the brute-force and the fast approximation of (1). The timing and PSNR of the proposed approximation for different values of σ and ε are shown in Table I. Notice that the timing scales linearly with log(1/ε). For fixed ε, the timing (order) is more when σ is small. This is generally the case with Fourier approximations [8], [10], [14]-it it is difficult to approximate a narrow Gaussian pulse using sinusoids. On the other hand,    [14], CBF [10], and the proposed method for fixed order K = 5. The error between the brute-force implementation of (1) and the respective approximations are also shown. The parameters used are σ = 55 and θ = 5. The corresponding PSNR value is mentioned in the caption.   our PSNR is consistently larger than the acceptable threshold of 40 dB [6], [13], even when ε is as large as 0.1. We next compare the filtering performance with CBF [10]. For a fair comparison, we have used the same order (number of convolutions) for all three methods. That is, we kept the timings same and compared the PSNRs. First, we set the order for our method using ε in Algorithm 1. We next tuned the tolerance parameter τ in CBF to obtain the same order. The order in FBF was set directly, since T is not optimized in this case. A couple of comparisons are shown in Figures 5 and 6 on the grayscale images Montage and Cameraman [20]. For the example in Figure 5, the PSNR improvement over CBF is about 10 dB. While it is somewhat difficult to access this improvement by directly comparing the filtered images, it is evident from the respective error images that the approximation is better near edges in our method. Considering the state-ofthe-art performance of CBF, the 3 dB improvement in Figure  6 is significant. An additional visual comparison at σ = 25 is provided in the supplement, where the increment is by 4.2 dB. Finally, in Table II, we have compared the average PSNRs over a set of images [20] for practical settings of θ and σ. Notice that our method is at least as good as CBF in terms of filtering accuracy. The PSNR improvement is in the range 0.7 to 1.3 dB.

IV. CONCLUSION
We showed that by jointly optimizing the order, period, and coefficients in [14], we can outperform the state-of-the-art CBF in certain cases. In general, we showed that the accuracy of our method is at least as good as CBF. As mentioned in the introduction, unlike CBF, the proposed approximation can also be used for non-Gaussian range kernels [9], [11], [12]. Also, we were able to establish convergence and provide a bound on the filtering error.