# Constant time Gaussian and Bilateral filter

**URL:** <https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867>\
**Category:** Performance\
**Tags:** images\
**Created:** [November 27, 2020, 3:08pm UTC](https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867 "2020-11-27T15:08:39Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Geoffrey](https://avatars.discourse-cdn.com/v4/letter/g/bc8723/32.png) [@Geoffrey](https://discourse.julialang.org/u/Geoffrey)\
**Post date:** [November 27, 2020, 3:08pm UTC](https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867/1 "2020-11-27T15:08:39Z")

</div>

Hi,

I’m still learning the julia language and to do so I try to implement different published algorithms that I found useful for my work (I’m a biologist and I do a lot of image analysis). So I have an implementation of the constant time Gaussian filter (`O1gaussnd`) from the papers of Sugimoto et al. ([https://doi.org/10.3169/mta.3.12](https://doi.org/10.3169/mta.3.12) ; [Fast Gaussian filter with second-order shift property of DCT-5 | IEEE Conference Publication | IEEE Xplore](https://doi.org/10.1109/ICIP.2013.6738106) ; [Compressive Bilateral Filtering | IEEE Journals & Magazine | IEEE Xplore](https://doi.org/10.1109/TIP.2015.2442916)) that I found to outperform the solution of `imfilter`. The code can be found on here: [GitHub - gschivre/Constant-time-gaussian-filter](https://github.com/gschivre/Constant-time-gaussian-filter)

This code also propose an implementation of the compressive bilateral filter (`CBLF`) (from the same authors) that lead to a constant time bilateral filtering as well as a draft of the implementation of the constant time first and second order Gaussian derivative filter (`Fastdgauss` and `Fastd2gauss`).

I create this topic because this code might be of interest for the images analysis community, and also because I feel like some people more experimented than me with julia could easily improve my code for including them in the `ImageFiltering` package 🙂

![benchmark](https://global.discourse-cdn.com/julialang/original/3X/5/3/53e52525d529128fca4416a0271095d4133d1ae4.png)

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 27, 2020, 10:12pm UTC](https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867/2 "2020-11-27T22:12:18Z")

</div>

Can you show your code? (EDIT: I see, it’s in [https://github.com/gschivre/Constant-time-gaussian-filter/blob/main/benchgauss.jl.](https://github.com/gschivre/Constant-time-gaussian-filter/blob/main/benchgauss.jl.)) And which kernel from ImageFiltering did you benchmark against? You should include `KernelFactors.IIRGaussian` in your comparisons, if you haven’t already. (Since `imfilter` is pretty much constant with `σ`, perhaps you used it already, but I can’t tell.)

---

<div class="post-metadata">

**Author:** ![Geoffrey](https://avatars.discourse-cdn.com/v4/letter/g/bc8723/32.png) [@Geoffrey](https://discourse.julialang.org/u/Geoffrey)\
**Post date:** [November 28, 2020, 12:38am UTC](https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867/3 "2020-11-28T00:38:05Z")

</div>

@tim.holy sorry for the lack of detail for the benchmark. The code for the benchmark are in my github in the file `benchgauss.jl`. However, thank you for pointing out the fact that I didn’t use the IIR Gaussian filter of `imfilter`. I did it just now and effectively there is no more improvement on the computing time for my implementation of the Sugimoto’s code:

```julia
using BenchmarkTools, ImageFiltering, TestImages, Images, Statistics, Plots
include("CstTimeGauss.jl");
img = testimage("mandrill");
# using IIR gaussian
σ = [1; 3; 5; 10; 15; 20; 30; 50; 70; 100]
tcst = []
timf = []
for s in σ
    b = @benchmark O1gaussnd(float64.($img), $s);
    push!(tcst, b.times)
    b = @benchmark imfilter($img, KernelFactors.IIRGaussian(($s, $s)), "reflect");
    push!(timf, b.times)
end
# display mean time and standard deviation for each σ
m = hcat([mean(tcst[i]) for i in 1:length(σ)], [mean(timf[i]) for i in 1:length(σ)]) ./ 1e6
s = hcat([std(tcst[i]) for i in 1:length(σ)], [std(timf[i]) for i in 1:length(σ)]) ./ 1e6
Plots.scatter(σ, m[:, 1]; yerrors = s[:, 1], xlabel = "σ", ylabel = "time (ms)", label = "O1gaussnd", legend = :topleft)
Plots.scatter!(σ, m[:, 2]; yerrors = s[:, 2], label = "imfilter IIR")

```

![benchmarkIIR](https://global.discourse-cdn.com/julialang/original/3X/4/c/4cc20d7894604d4710377d4f16740219b4529b52.png)

But I also wanted to point out another thing: it seems that the precision of the Gaussian convolution with the code of Sugimoto et al. is better than the one from `imfilter`. To evidence this I looked at the difference from `imfilter` and `O1gaussnd` and I see that this 2 give much similar results if I force `imfilter` to use gaussian kernel of size of at least 6 \sigma:

```julia
# default kernel size
maximum(abs.(imfilter(img, KernelFactors.gaussian((10, 10)), "reflect") .- O1gaussnd(float64.(img), 10)))
0.05299627892037459
# kernel size of 6 σ
maximum(abs.(imfilter(img, KernelFactors.gaussian((10, 10), (61, 61)), "reflect") .- O1gaussnd(float64.(img), 10)))
0.0033663324215758017

```

But with the IIR filter I get these difference:

```julia
# IIR filter vs imfilter default
maximum(abs.(imfilter(img, KernelFactors.IIRGaussian((10, 10)), "reflect") .- imfilter(img, KernelFactors.gaussian((10, 10)), "reflect")))
1.025380269215756
# IIR filter vs imfilter 6 σ
maximum(abs.(imfilter(img, KernelFactors.IIRGaussian((10, 10)), "reflect") .- imfilter(img, KernelFactors.gaussian((10, 10), (61, 61)), "reflect")))
1.023126836702493
# IIR filter vs O1gaussnd
maximum(abs.(imfilter(img, KernelFactors.IIRGaussian((10, 10)), "reflect") .- O1gaussnd(float64.(img), 10)))
1.0234280940523726

```

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 28, 2020, 8:17am UTC](https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867/4 "2020-11-28T08:17:28Z")

</div>

Better accuracy would be great. I haven’t looked at your code, but in principle it would be welcomed as a contribution to ImageFiltering. Alternatively, I could add a link to your repo from ImageFiltering’s docs. Let me know.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [November 28, 2020, 11:47am UTC](https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867/5 "2020-11-28T11:47:09Z")

</div>

Usually the IIR methods has lower accuracy than those based on scaling of the kernel (For some reason called Compressive in the last few years).

Bus mostly, the fastest Gaussian Blur is by using the scaling property of the parameter \sigma and use Downscaling and Upscaling. This is usually gives the best performance and the best accuracy.

---

<div class="post-metadata">

**Author:** ![Geoffrey](https://avatars.discourse-cdn.com/v4/letter/g/bc8723/32.png) [@Geoffrey](https://discourse.julialang.org/u/Geoffrey)\
**Post date:** [November 28, 2020, 12:47pm UTC](https://discourse.julialang.org/t/constant-time-gaussian-and-bilateral-filter/50867/6 "2020-11-28T12:47:47Z")

</div>

I think it could be great to having this code in `ImageFiltering` together with the code for the constant time bilateral filter. However, I must also admit that my implementation may need to be improved before fitting with the standard of `ImageFiltering`. For exemple in my implementation I use reflecting boundary conditions by default as I didn’t produce a padded array but rather check during the loop where the entry is and replace it by the mirror extension if it fall before the first or after the last element of the array. I also need to pass the image into `float64` before running the filter which may be also avoided.