Acceleration of the shiftable O(1) algorithm for bilateral filtering and non-local means

Introduction

The bilateral filter is an edge-preserving diffusion filter, which was introduced by Tomasi et al. in . The edge-preserving property comes from the use of a range kernel (along with the spatial kernel) that is used to control the diffusion in the vicinity of edges. In this work, we will focus on the Gaussian bilateral filter where both the spatial and range kernels are Gaussian . This is given by

Here, $g_{\sigma_{s}}(\boldsymbol{x})$ is the centered Gaussian distribution on the plane with variance $\sigma_{s}^{2}$, and $g_{\sigma_{r}}(s)$ is the one-dimensional Gaussian distribution with variance $\sigma_{r}^{2}$; $\Omega$ is the support of $g_{\sigma_{s}}(\boldsymbol{x})$ over which the averaging takes place. We call $g_{\sigma_{s}}(\boldsymbol{x})$ and $g_{\sigma_{r}}(s)$ the spatial and the range kernel.

The range kernel is controlled by the local distribution of intensity. Sharp discontinuities (jumps) in intensity typically occur in the vicinity of edges. This is picked up by the range kernel, which is then used to inhibit the spatial diffusion. On the other hand, the range kernel becomes inoperative in regions with smooth variations in intensity. The spatial kernel then takes over, and the bilateral filter behaves as a standard diffusion filter. Together, the spatial and range kernels perform smoothing in homogeneous regions, while preserving edges at the same time .

The bilateral filter has found widespread use in several image processing, computer graphics, and computer vision applications ; see for further applications. More recently, the bilateral filter was extended by Baudes et al. in the form of the non-local means filter, where the similarity between pixels is measured using patches centered around the pixel.

The direct implementation of (1) is computationally intensive, especially when $\sigma_{s}$ is large ($\sigma_{r}$ has no effect on the run time in this case). In particular, the direct implementation requires $O(\sigma_{s}^{2})$ operations per pixel. This makes the filter slow for real-time applications. Several efficient algorithms have been proposed in the past for implementing the filter in real time, e.g., see . In , Porikli demonstrated for the first time that the bilateral filter could be implemented using $O(1)$ operations per pixel (with respect to $\sigma_{s}$). This was done for two different settings: $(a)$ Spatial box filter and arbitrary range filter, and $(b)$ Arbitrary spatial filter and polynomial range filter. The author extended $(b)$ to the Gaussian bilateral filter in (1) by approximating $g_{\sigma_{r}}(s)$ with its Taylor polynomial. The run time of this approximation was linear in the order of the polynomial. The problem with Taylor polynomials, however, is that they provide good approximations of $g_{\sigma_{r}}(s)$ only locally around the origin. In particular, they have the following drawbacks:

Taylor polynomials are not guaranteed to be positive and monotonic away from the origin, where the approximation is poor. Moreover, they tend to blow up at the tails.

It is difficult to approximate $g_{\sigma_{r}}(s)$ using the Taylor expansion when $\sigma_{r}$ is small. In particular, a large order polynomial is required to get a good approximation of a narrow Gaussian, and this considerably increases the run time of the algorithm.

The first of these problems was addressed in . In this paper, the authors observed that it is important that the kernel used to approximate $g_{\sigma_{r}}(s)$ be positive, monotonic, and symmetric. While it is easy to ensure symmetry, the other two properties are hard to enforce using Taylor approximations. It was noticed that, in the absence of these properties, the bilateral filter in created strange artifacts in the processed image (cf. Figure 3 in ). The authors proposed to fix this problem using the family of raised cosines, namely, functions of the form

Here $N$ is the order of the kernel, which controls the width of $\phi(s)$. The kernel can be made narrow by increasing $N$.

The key parameter in (2) is the quantity $T$. The idea here is that $[\cos(s)]^{N}$ is guaranteed to be positive and monotonic provided that $s$ is restricted to the interval $[-\pi/2,\pi/2]$. Note that the argument $s$ in (2) takes on the values $|f(\boldsymbol{x}-\boldsymbol{y})-f(\boldsymbol{x})|$ as $\boldsymbol{x}$ and $\boldsymbol{y}$ varies over the image. Therefore, by letting

one could guarantee $\phi(s)$ to be positive and monotonic over $[-T,T]$. In , $T$ was simply set to the maximum dynamic range, for example, $255$ for grayscale images. We refer the readers to Figure 2 in for a comparison of (2) and the polynomial kernels used in . It was later observed in that polynomials could also be used for the same purpose. The polynomials suggested were of the form

2 Fast O​(1)𝑂1O(1) implementation using shiftable kernels

For completeness, we now explain how the above kernels can be used to compute (1) using $O(1)$ operations. As observed in , (2) and (3) are essentially the simplest kernels that have the so-called property of shiftability. This means that, for a given $N$, we can find a fixed set of basis functions $\phi_{1}(s),\ldots,\phi_{N}(s)$ and coefficients $c_{1},\ldots,c_{N}$, so that for any translation $\tau$, we can write

The coefficients depend continuously on $\tau$, but the basis functions have no dependence on $\tau$. For (2), both the basis functions and coefficients are cosines, while they are polynomials for (3). This shiftability property is at the heart of the $O(1)$ algorithm. Let $\overline{f}(\boldsymbol{x})$ denote the output of the Gaussian filter $g_{\sigma_{s}}(\boldsymbol{x})$ with neighborhood $\Omega$,

Note that, by replacing $g_{\sigma_{r}}(s)$ with $\phi(s)$, we can write (1) as

where we have set $F_{i}(\boldsymbol{x})=f(\boldsymbol{x})\phi_{i}(f(\boldsymbol{x}))$. Similarly, by setting $G_{i}(\boldsymbol{x})=\phi_{i}(f(\boldsymbol{x}))$, we can write

Now, it is well-known that certain approximation of (5) can be computed using just $O(1)$ operations per pixel. These recursive algorithms are based on specialized kernels, such as the box and the hat function , and the more general class of box splines . Putting all these together, we arrive at the following $O(1)$ algorithm for approximating (1):

Fix $N$, and approximate $g_{\sigma_{r}}(s)$ using (2) or (3).

For $i=1,2,\ldots,N$, set up the images $F_{i}(\boldsymbol{x})=f(\boldsymbol{x})\phi_{i}(f(\boldsymbol{x}))$ and $G_{i}(\boldsymbol{x})=\phi_{i}(f(\boldsymbol{x}))$, and the coefficients $c_{i}(f(\boldsymbol{x}))$ .

Use a recursive $O(1)$ algorithm to compute each $\overline{F_{i}}(\boldsymbol{x})$ and $\overline{G_{i}}(\boldsymbol{x})$.

Plug these into (6) and (7) to get the filtered image.

It is clear that better approximations are obtained when $N$ is large. On the other hand, the run time scales linearly with $N$. One key advantage of the above algorithm, however, is that the $\overline{F_{i}}(\boldsymbol{x})$ and $\overline{G_{i}}(\boldsymbol{x})$ can be computed in parallel. For small orders ($N<10$), the serial implementation is found to be comparable, and often better, than the state-of-the-art algorithms. The parallel implementation, however, turns out to be much faster than the competing algorithms, at least for $N<50$. Henceforth, we will refer to the above algorithm as SHIFTABLE-BF, the shiftable bilateral filter.

This brings us to the question as to whether we can always work with, say, $N<50$ basis functions, in SHIFTABLE-BF? To answer this question, we must explain in some detail step (1) of the algorithm, where we approximate the Gaussian range kernel,

on the interval $[-T,T]$. This could be done either using (2),

These approximations were proposed in . Note that, we have to rescale (2) by $\sqrt{N}$, and (3) by $N$, to get to the right limit. On the other hand, to ensure positivity and monotonicity, we need to guarantee that the arguments of (8) and (9) are in the intervals $[-\pi/2,\pi/2]$ and $$. A simple calculation shows that this is the case provided that$N$is larger than$N_{0}=4T^{2}/\pi^{2}\sigma_{r}^{2}=0.405(T/\sigma_{r})^{2}$for the former, and$N_{0}=0.5(T/\sigma_{r})^{2}$for the latter. In other words, it is not sufficient to set$N$large – it must be at least be as large as$N_{0}$. In Table 1, we give the values of$N_{0}$for different values of$\sigma_{r}$when$T=255$. It is seen that$N_{0}$gets impracticable large for$\sigma_{r}<30$. This does not come as a surprise since it is well-known that one requires a large number of trigonometric functions (or polynomials) to closely approximate a narrow Gaussian on a large interval. As pointed out earlier, this was also one of the problems in .

4 Present Contributions

In this paper, we address the above problem, namely that $N_{0}$ grows as $O(T^{2}/\sigma_{r}^{2})$ with $\sigma_{r}$. In Section 2, we propose a fast algorithm for determining $T$ exactly. Besides cutting down $N_{0}$, this is essential for determining the (local) dynamic range of a grayscale image that has been deformed, e.g., by additive noise. Setting $T=255$ in this case can lead to artifacts in the processed image. Next, in Section 3, we provide a simple and practical means of reducing the order, which leads to quite dramatic reductions in the run time of SHIFTABLE-BF. These modifications are also applicable to the shiftable algorithms proposed in . Finally, in Section 4, we provide some experimental results to demonstrate the acceleration that is achieved using these modifications. We also compare our algorithm with the Porikli’s algorithms , both in terms of speed and accuracy.

Fast algorithm for finding T𝑇T

For the rest of the discussion, we work with finite-sized images (bounded $\Omega$) on the Cartesian grid. We continue to use $\boldsymbol{x}$ and $\boldsymbol{y}$ to denote points on the grid. The integral in (1) is simply replaced by a finite sum over $\Omega$. We will use the norm $\lVert\boldsymbol{x}\rVert=|x_{1}|+|x_{2}|$, where $\boldsymbol{x}=(x_{1},x_{2})$. Without loss of generality, we assume that $\Omega$ is a square neighborhood, that is, $\Omega=\{\boldsymbol{x}:\lVert\boldsymbol{x}\rVert\leq R\}$ where, say, $R=3\sigma_{s}$ (if $\Omega$ is not rectangular, we take the smallest rectangle containing $\Omega$).

Note that, for a given $\sigma_{r}$, we can cut down $N_{0}$ by using a tight estimate for

The smaller the estimate, the lower is the threshold $N_{0}$. The point is that the worst-case estimate $T=255$ is often rather loose for grayscale images. For example, we give the exact values of $T$ for a test image in Table 2, computed at different values of $\sigma_{s}$. We also give the time required to compute $T$.

It is seen that the exact values of $T$ are indeed much less than the worst-case estimate, particularly for small $\sigma_{s}$. For $\sigma_{s}=3$, $T$ is only about $205$. Even for $\sigma_{s}$ as large as $30$, $T$ is about $215$. Consider the bilateral filter with $\sigma_{s}=10$ and $\sigma_{r}=10$. From Table 1, $N_{0}=263$ using $T=255$. However, using the exact value $T=210$, we can bring this down to $263\cdot(210/255)^{2}\approx 178$, a reduction by almost $100$. For smaller values of $\sigma_{r}$, this gain is even more drastic. However, notice the time required to compute $T$ in Table 1. This increases quickly with the increase in $\sigma_{s}$ (in fact, scales as $O(R^{2})$). Experiments show us that, for large $\sigma_{s}$, this is comparable to the time required to compute the bilateral filter. It would thus help to have an $O(1)$ algorithm for computing $T$. Motivated by our previous work on filtering using running sums , we recently devised an algorithm that does exactly this. We later found that the algorithm had already been discovered two decades back in a different context .

Our algorithm is based on the following observations. First, note that we can take out the modulus from (10) using symmetry.

This follows from the observations that $|t|=\max(t,-t)$, and that $\lVert\boldsymbol{x}-\boldsymbol{y}\rVert\leq R$ is symmetric in $\boldsymbol{x}$ and $\boldsymbol{y}$. Moreover, note that the operation that takes two numbers $a$ and $b$ and returns $\max(a,b)$ is associative. Using associativity, we can write (10) as $T=\max(T_{+},T_{-})$, where

We claim that $T_{+}=T_{-}$, so that we need not compute them separately. Indeed, suppose that the first maximum is attained at $\boldsymbol{x}_{0}$ and $\boldsymbol{y}_{0}$, that is, $T_{+}=f(\boldsymbol{x}_{0})-f(\boldsymbol{x}_{0}-\boldsymbol{y}_{0})$. Taking $\boldsymbol{x}=\boldsymbol{x}_{0}-\boldsymbol{y}_{0}$ and $\boldsymbol{y}=-\boldsymbol{y}_{0}$, and noting that $\lVert\boldsymbol{x}-\boldsymbol{y}\rVert\leq R$, we must have

By an identical argument, $T_{+}\geq T_{-}$, and the proposition follows. ∎

The problem is now reduced to that of computing the windowed maximums in (11). A direct computation would still require $O(R^{2})$ comparisons. It turns out that we can do this very fast (no matter how large is $R$) by exploiting the overlap between adjacent windows. This is done by adapting the so-called MAX-FILTER algorithm.

There is an $O(1)$ algorithm for computing

This algorithm was first proposed by van Herk, and Gil and Werman . It is clear that since the search domain $\Omega$ is separable, it suffices to solve the problem in one dimension. The problem in two-dimensions can be solved simply by iterating the one-dimensional MAX-FILTER along each dimension. From (11), we arrive at Algorithm 1 for computing $T$ with $O(1)$ operations. We note that Algorithm 1 does not actually compute $\max\ \{|f(\boldsymbol{x}-\boldsymbol{y})-f(\boldsymbol{x})|:\lVert\boldsymbol{y}\rVert\leq R\}$ at every $\boldsymbol{x}$. It only has access to the distribution of the maximums. For completeness, we have explained the MAX-FILTER algorithm in the Appendix. For further details, we refer the readers to .

Acceleration using truncations

We have seen that, by using the exact value of $T$, we can bring down the run time by $10-20\%$. Unfortunately, Table 1 tells us that this alone is not sufficient in the regime $\sigma_{r}<15$. For example, $N_{0}$ is of the order $10^{3}$ in the regime $\sigma_{r}<5$. So why do we require so many terms in (8) and (9) to approximate a narrow Gaussian? This is exactly because we are forcing $\phi(s)$ to positive and monotonic on its broad tails, where $g_{\sigma_{r}}(s)$ is close to zero. For example, consider the approximation in (8):

By requiring $N>N_{0}$, we can guarantee that (1) $\phi(s)$ is close to $g_{\sigma_{r}}(s)$, and (2) $\phi(s)$ is positive and monotonic over $[-T,T]$. Note, however, that $g_{\sigma_{r}}(s)$ falls off very fast, and almost vanishes outside $\pm 3\sigma_{r}$. It turns out that only a few significant terms in (14) contribute to the approximation in the $\pm 3\sigma_{r}$ region. The rest of the terms have a negligible contribution, and are required only to force positivity at the tails.

Thanks to expression (14), it is now straightforward to determine which are the significant terms. Note that the coefficients $2^{-N}\binom{N}{n}$ in (14) are positive and sum up to one. In fact, they are unimodal and closely follow the shape of the target Gaussian. The smallest coefficients are at the tails, and the largest coefficients are at the center. In particular, the smallest coefficient is $1/2^{N}$, while the largest one is $2^{-N}[(N/2)!]^{-2}N!$ (assuming $N$ to be even). For large $N$, the latter is approximately $\sqrt{2/\pi N}$ using Stirling’s formula. Thus, as $N$ gets large, the coefficients get smaller. What is perhaps significant is that the ratio of the smallest to the largest coefficient is $(\pi N2^{2N-1})^{-1/2}$, and this keeps shrinking at an exponential rate with $N$. On the other hand, the cosine functions (which act as the interpolating function) are always bounded between $$.

The above observation suggests dropping the small terms on the tail. In particular, for a given tolerance $\varepsilon>0$, let $M=M(\varepsilon)$ be the smallest term for which

It follows that the error $|\phi(s)-\phi_{\varepsilon}(s)|$ is within $\varepsilon$ for all $-T\leq s\leq T$. Note that $\phi_{\varepsilon}(s)$ is symmetric, but is no longer guaranteed to be positive on the tails, where oscillations begin to set in. However, what we can guarantee is that the negative overshoots are within $-\varepsilon$. In fact, the quality of the final approximation turns out to be quite satisfactory. This is illustrated with an example in Figure 1. The main point is that, in the regime $\sigma_{r}<15$, we are now able to bring down the order to well within $100$. We list some of them in Table 3. Notice that we can drop more terms for a given accuracy as the kernel gets narrow. For $\sigma_{r}<5$, we can drop almost $95\%$ of the terms, while keeping the error within $0.5\%$ of the peak value.

Note that (15) actually requires us to compute a large number of tails coefficients, which are eventually not used in (16). It is thus better to estimate $M$ when $N$ is large. A good estimate of (15) is provided by the Chernoff bound for the binomial distribution , namely,

It can be verified the estimate is quite tight for $N>100$. By setting the bound to $\varepsilon/2$, we get

The final algorithm obtained by combining the proposed modifications is given in Algorithm 2. Henceforth, we will continue to refer to this as the SHIFTABLE-BF. The Matlab implementation of SHIFTABLE-BF can be found here .

Experiments

We now provide some results on synthetic and natural images to understand the improvements obtained used our proposal. While all the experiments were done on Matlab, we took the opportunity to report the run time of a multithreaded Java implementation of Algorithm 2. All experiments were run on an Intel quad core $2.83$ GHz processor.

First, we tested the speedup obtained using Algorithm 1 for computing $T$. We used a Matlab implementation of this algorithm . It is expected that the run time remain roughly the same for different $\sigma_{s}$. As seen in table 4, this is indeed the case. The run time of the direct method, for the same image and the same settings of $\sigma_{s}$, was already provided in table 2. Note that we have been able to cut down the time by a few orders using our fast algorithm.

We then compared the run times of multithreaded Java implementations of SHIFTABLE-BF proposed in and its present refinement. For this, we used the test image Checker shown in Fig. 2. We have used small values of $\sigma_{r}$, and a fixed $\sigma_{s}=15$. The average run times are shown in Table 5. The tolerance $\varepsilon$ used for the truncation are also given. We use a smaller $\varepsilon$ (larger truncation) as $\sigma_{r}$ gets small. Notice how we have been able to cut down the run time by more than $70\%$. This is not surprising, since we have discarded more than $85\%$ of terms. Notice that the run times are now well within $1$ second. The run time of the direct implementation (which does not depend on $\sigma_{r}$) was around $10$ seconds. The main point is that we can now implement the filter in a reasonable amount of time for small $\sigma_{r}$, which could not be done previously in .

2 Accuracy

We next studied the effect of truncation. To get an idea of the noise that is injected into the filter due to the truncation, we used the Checker image. This particular image allowed us to test both the diffusive and the edge-preserving properties of the filter at the same time. We used the setting $\sigma_{s}=30$ and $\sigma_{r}=10$. First, we tried the direct implementation of (1), using a very fine discretization. Then we tried Algorithm 2. The difference between the two outputs is shown in Figure 3. The artifacts shown in the image are actually quite insignificant, within $10^{-5}$ times the peak value. Notice that most of the artifacts are around the edges. This comes from the oscillations induced at the tails of kernel by the truncation. To compare the filter outputs (with and without truncation), we extracted two horizontal scan profiles from the respective outputs. These are shown in Figure 4. Notice that it is rather hard to distinguish the two.

We then applied the filters on the standard grayscale image of Lena. We first applied the direct implementation followed by Algorithm 2. In this case, $T$ was computed to be $215$. The difference image is shown in Figure 5. It is again seen that the small artifacts are cluttered near the edges. We have also tried measuring the mean-squared-error (MSE) for different $\sigma_{r}$. The results are given in Table 6. Note that relatively larger MSEs are obtained at small $\sigma_{r}$. This is because we are forced to use a large truncation to speed up the filter at small $\sigma_{r}$. The above results show that we can drastically cut down the run time of filter using the proposed modifications, without incurring significant errors.

3 Comparison with a benchmark algorithm

We next compared the performance of the improved SHIFTABLE-BF algorithm with those proposed in . The latter algorithms are considered as benchmark in the literature on fast bilateral filtering. Porikli proposed a couple of algorithms in – one using a variable spatial filter and a polynomial range filter (we call this BF1), and the other using a constant spatial filter and a variable range filter (we call this BF2). The difficulty with BF1 is that it is rather difficult to control the width of the polynomial range filter. In particular, as was already mentioned in the introduction, it is difficult to approximate narrow Gaussian range kernels using BF1. We refer the interested readers to the experimental results in , where a comparison was already made between BF1 and SHIFTABLE-BF. For completeness, we perform a single experiment to compare these filters when $\sigma_{r}$ is small (a rather large value of $\sigma_{r}$ was used in the experiments in ). For this, and the remaining experiments, we will consider the standard test image of Barbara of size $512\times 512$. This image has several texture patterns, and is well-suited for comparing the performance of bilateral filters with narrow range kernelswe thank one of the reviewers for suggesting this example.. In particular, we consider the Gaussian bilateral filter with $\sigma_{s}=20$ and $\sigma_{r}=20$. The results obtained using SHIFTABLE-BF and BF1 are shown in Figure 6. The error between the direct implementation of the bilateral filter and SHIFTABLE-BF was within $10^{-3}$. On the other hand, note how BF1 completely breaksdown. The reason for this was already mentioned in the introduction. A similar breakdown, with a larger $\sigma_{r}$, was also observed in Figure 3 in .

We next considered BF2, which does not suffer from the above problem. However, we note that BF2 cannot be used to perform Gaussian bilater filtering – it only works with constant spatial filters (box filters). This is because BF2 uses fast integral histograms, and this only works for box filters. To make the comparison even, we considered bilateral filters with constant spatial filters and Gaussian range kernels. We note SHIFTABLE-BF can be trivially modified to work with arbitrary spatial filters.

First, we compared the run times of the Matlab implementations of SHIFTABLE-BF and BF2. The results obtained at particular settings of $\sigma_{s}$ (radius of box filter) and $\sigma_{r}$ are shown in Figure 7. We see that the run time of SHIFTABLE-BF is consistently better than that of BF2. The difference is particularly large when $\sigma_{r}>10$, and it closes down as $\sigma_{r}$ gets small. All these can be perfectly explained. Note that, as per the design, the computational complexity of BF2 is $O(1)$ both with respect to $\sigma_{s}$ and $\sigma_{r}$. This is indeed seen to be the case from the run times. On the other hand, the computational complexity of SHIFTABLE-BF is $O(1)$ with respect to $\sigma_{s}$ (this is again clear from the plots in Figure 7). However, for a given $\sigma_{s}$, the complexity of the the original algorithm scales as $O(1/\sigma_{r}^{2})$. The complexity, in fact, remains roughly the same even after the proposed truncation. This explains the step rise in the run time for small values of $\sigma_{r}$, as shown in Figure 7. However, the actual run time goes down substantially as a result of the truncations (cf. Table 5). In particular, we have noticed that the worst case run time of SHIFTABLE-BF is less than the average run time of BF2 for $\sigma_{r}$ as low as $3$.

We note that the run time of BF2 depends on the number of bins used for the integral histogram. In the above experiments, we used as many bins as the grayscale levels of the image. It is thus possible to reduce the run time by cutting down the resolution of the histogram. However, this comes at the cost of the quality of the filtered image. This lead us to compare the outputs of SHIFTABLE-BF and BF2. In Figure 8, we compared the MSEs of the two algorithms for different $\sigma_{s}$ and $\sigma_{r}$ for the image Barbara. We note that the MSE for SHIFTABLE-BF is significantly lower than BF2. The gap is around $30$ dB for $\sigma_{r}>10$, and around $90$ dB when $\sigma_{r}\leq 10$. We noticed that this difference becomes even more pronounced if we use a smaller number of bins. As expected, note that the individual MSEs do not vary much with $\sigma_{s}$. The reader will notice that the MSE of SHIFTABLE-BF suddenly drops by $60$ dB when $\sigma_{r}\leq 10$. To explain this, we plot the effective order $N-2M$ for different $\sigma_{r}$ in Figure 9. We see that the order (of approximation) suddenly jumps up when $\sigma_{r}$ goes below $10$, which explains the jump in the MSE in Figure 8. This is simply due to the rule (17) used in Algorithm 2

In Figure 10, the filtered outputs of the two algorithms are compared with the direct implementation, for a small value of $\sigma_{r}$. Note that the pointwise error between the direct implementation and SHIFTABLE-BF is of the order $10^{-3}$. On the other hand, the corresponding error between the direct implementation and BF2 is substantial, a few orders larger than that for SHIFTABLE-BF. This explains the large gap between the MSEs in Figure 8.

We close this section by commenting on the memory usage of SHIFTABLE-BF and BF2. The former requires us to compute and store a total of $2(2N-M)$ images, while the latter requires us to store a histogram with $B$ bins per pixel (equivalent of $B$ images). It is clear from Figure 9 that even for a half-resolution histogram ($B=128$) and for $\sigma_{r}>10$, the memory requirement of BF2 is comparable to that of SHIFTABLE-BF. For smaller values of $\sigma_{r}$, SHIFTABLE-BF clearly requires more memory than BF2.

Discussion

In this paper, we proposed some simple ways of accelerating the bilateral filtering algorithm proposed in . This, in particular, opened up the possibility of implementing the algorithm in real-time for small $\sigma_{r}$. We note that the problem of determining the optimal $\sigma_{s}$ and $\sigma_{r}$ for a given application is extrinsic to our algorithm. We are only required to determine the parameter $T$, which is intrinsic to our algorithm. A fast algorithm was proposed in the paper for this purpose. However, we note that having a fast algorithm does make it easier to determine the optimal parameters. In this regard, we note that Kishan et al. have recently shown how our fast algorithm can be used to tune the parameters for image denoising, under different noise models . One crucial observation used in these papers is that the a certain unbiased estimator of the MSE can be efficiently computed for our fast bilater filter, using the linear expansions in (6) and (7). The “best” parameters are choosen by optimizing this MSE estimator. While this can also be done for the polynomial-based bilateral filter in , this trick cannot be used for other fast implementations of the bilateral filter, at least to the best of our knowledge.

Finally, we note that the ideas proposed here can also be extended to the $O(1)$ algorithm for non-local means given in . In non-local means , the range kernel operates on patches centered around the pixel of interest. A coarse non-local means was considered in , where a small patch neighborhood consisting of the pixels $\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{p}$ (where, say, $\boldsymbol{u}_{1}=0$) was used. In this case, the main observation was that formula for the non-local means can be written in terms of following sums:

where $g(s_{1},\ldots,s_{p})$ is an anisotropic Gaussian in $p$ variables, and has a diagonal covariance. This looks very similar to (1), except that we now have a multivariate range kernel. By using the separability of $g(s_{1},\ldots,s_{p})$, and by approximating each Gaussian component by either (8) or (9), a $O(1)$ algorithm for computing (18) and (19) was developed. We refer the readers to for further details. The key consideration with this algorithm is that the overall order scales as $N^{p}$, where $N$ is the order of the Gaussian approximation for each component. In this case, it is thus important to keep $N$ as low as possible for a given covariance. After examining (18) and (19), it is clear that the interval over which the one dimensional Gaussians need to be approximated is $[-T,T]$, where $T$ is as defined in (10). We can compute this using Algorithm 1. Moreover, we can further reduce the order by truncation, especially when the covariance is small. However, $N^{p}$ can still be large (even for $p=3$ or $4$), and hence a parallel implementation must be used for real-time implementation.

Appendix : MAX-FILTER algorithm

We explain how the MAX-FILTER algorithm works in one dimension. Let $f_{1},f_{2},\ldots,f_{N}$ be given, and we have to compute $\max(f_{i-R},\ldots,f_{i+R})$ at every interior point $i$. Assume $R$ is an integer, and $N$ is a multiple of the window size $W=2R+1$, say, $N=pW$ (padding is used if this not the case). The idea is to compute the local maximums using running maximums, similar to running sums used for local averaging . The difference here is that, unlike averaging, the max operation is not linear. This can be fixed using “local” running maximums.

We begin by dividing $f_{1},f_{2},\ldots,f_{N}$ into $p$ equal partitions. The $k$th partition ($k=0,1,\ldots,p-1$) is composed of $f_{1+kW},\ldots,f_{W+kW}$. For a given partition $k$, we recursively compute the two running maximums (of length $W$), one from from the left and one from the right. Let $l^{(k)}$ and $r^{(k)}$ be the left and right running maximums for the $k$th partition. The sequence $l^{(k)}$ start at the left of the partition with $l^{(k)}_{1}=f_{1+kW}$, and is recursively given by $l^{(k)}_{i}=\max(\ l^{(k)}_{i-1},f_{i+kW}\ )$ for $i=2,3,\ldots,W$. It ends on the right end of the partition. On the other hand, $r^{(k)}$ start at the right and ends on the left: $r^{(k)}_{W}=f_{W+kW}$, and $r^{(k)}_{W-i}=\max(\ r^{(k)}_{W-i+1},f_{W-i+kW}\ )$ for $i=1,2,\ldots,W-1$. This is done for every partition to get $l^{(0)},\ldots,l^{(p-1)}$ and $r^{(0)},\ldots,r^{(p-1)}$.

We now concatenate the left maximums into a single function $l_{1},\ldots,l_{N}$, that is, we set $l=(l^{(0)},\ldots,l^{(p-1)})$. Similarly, we concatenate the right maximums in order, $r=(r^{(0)},\ldots,r^{(p-1)})$. In practice, we just need to recursively compute $l_{1},\ldots,l_{N}$ and $r_{1},\ldots,r_{N}$, resetting the recursion at the boundary of every partition. We now split $f_{i-R},\ldots,f_{i+R}$ into two segments, which either belong to the same partition or two adjacent partitions. Then, from the associativity of the max operation, it is seen that $\max(f_{i-R},\ldots,f_{i+R})=\max(r_{i-R},l_{i+R})$. Note that, we need just $3$ max operations per point to get the result, independent of the window size $R$.

Acknowledgments

This work was partly supported by the Swiss National Science Foundation under grant PBELP2-$135867$. The author thanks M. Unser and D. Sage for interesting discussions, and the anonymous referees for their helpful comments and suggestions. The author also thanks A. Singer and the Program in Applied and Computational Mathematics at Princeton University for hosting him during this work.

References