Replicating behavior of scipy.ndimate.gaussian_filter

General Usage
,
1

I am trying to convert an old Python script to Julia. At one point, the script makes a call to scipy.ndimage.gaussian_filter, for which I presume the Julia equivalent should be ImageFiltering.imfilter with Kernel.gaussian. However, I cannot get the latter to replicate the output of the former.

Naively:

julia> using ImageFiltering, PyCall

julia> a = reshape(1:25, (5,5)) .|> Float64
5×5 Array{Float64,2}:
 1.0   6.0  11.0  16.0  21.0
 2.0   7.0  12.0  17.0  22.0
 3.0   8.0  13.0  18.0  23.0
 4.0   9.0  14.0  19.0  24.0
 5.0  10.0  15.0  20.0  25.0

julia> py_filt = pyimport("scipy.ndimage").gaussian_filter(a, (1,0))
5×5 Array{Float64,2}:
 1.42704  6.42704  11.427   16.427   21.427 
 2.06782  7.06782  12.0678  17.0678  22.0678
 3.0      8.0      13.0     18.0     23.0   
 3.93218  8.93218  13.9322  18.9322  23.9322
 4.57296  9.57296  14.573   19.573   24.573 

julia> jl_naive = imfilter(a, Kernel.gaussian((1,0)))
5×5 Array{Float64,2}:
 1.35318  6.35318  11.3532  16.3532  21.3532
 2.05449  7.05449  12.0545  17.0545  22.0545
 3.0      8.0      13.0     18.0     23.0   
 3.94551  8.94551  13.9455  18.9455  23.9455
 4.64682  9.64682  14.6468  19.6468  24.6468

In this example, the first element is different by ~5%! Looking a bit deeper into the documentation for both functions, I notice that the default padding mode for the scipy function is “reflect” whereas Julia’s is “replicate”, however, changing this still does not yield the same array:

julia> jl_reflect = imfilter(a, Kernel.gaussian((1,0)), "reflect")
5×5 Array{Float64,2}:
 1.70636  6.70636  11.7064  16.7064  21.7064
 2.10898  7.10898  12.109   17.109   22.109 
 3.0      8.0      13.0     18.0     23.0   
 3.89102  8.89102  13.891   18.891   23.891 
 4.29364  9.29364  14.2936  19.2936  24.2936

Also Julia’s ImageFiltering reflect(Kernel.gaussian((1,0))) has no effect.

These differences seem way too big to be due to floating point errors, so what else am I doing wrong?

2

From what I understand, there isn’t a single standard implementation of a fast discrete Gaussian filter, so neither version is “wrong” – they’re just different approximations:

3

To see if there’s a difference in kernel, filter

a = zeros(5, 5); a[3, 3] = 1

If those are the same, see if the difference can be explained by padding. Pad manually, filter the pre-padded image, and crop back.

3 Likes
4

@GunnarFarneback Thank you very much! Playing around with padding and cropping, I found that what Scipy calls “reflect” padding is what ImageFiltering calls “symmetric”. Thus, setting this and the filter length to 9, I get the same result.

Cheers!

julia> py_filt = pyimport("scipy.ndimage").gaussian_filter(a, (1,0))
5×5 Array{Float64,2}:
 1.42704  6.42704  11.427   16.427   21.427 
 2.06782  7.06782  12.0678  17.0678  22.0678
 3.0      8.0      13.0     18.0     23.0   
 3.93218  8.93218  13.9322  18.9322  23.9322
 4.57296  9.57296  14.573   19.573   24.573 

julia> jl_filt = imfilter(a, Kernel.gaussian((1,0), (9,1)), "symmetric")
5×5 Array{Float64,2}:
 1.42704  6.42704  11.427   16.427   21.427 
 2.06782  7.06782  12.0678  17.0678  22.0678
 3.0      8.0      13.0     18.0     23.0   
 3.93218  8.93218  13.9322  18.9322  23.9322
 4.57296  9.57296  14.573   19.573   24.573
1 Like