Constant time Gaussian and Bilateral filter

General Usage Performance
1

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 ; Fast Gaussian filter with second-order shift property of DCT-5 | IEEE Conference Publication | IEEE Xplore ; Compressive Bilateral Filtering | IEEE Journals & Magazine | IEEE Xplore) that I found to outperform the solution of imfilter. The code can be found on here: GitHub - 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 :slight_smile:

benchmark

2 Likes
2

Can you show your code? (EDIT: I see, it’s in 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.)

3

@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:

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

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:

# 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:

# 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
1 Like
4

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.

5

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.

6

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.