Bug in 3d cascade

[timholy]

From https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/12?u=tim.holy (with some typos fixed):

using ImageFiltering, Test

a = rand(10,10,10);
b = similar(a);
kf2 = kernelfactors((ones(3)/3, ones(3)/3))
kf3 =  kernelfactors((ones(3)/3, ones(3)/3, [1.0]))
for i in 1:10
       b[:,:,i] = imfilter(a[:,:,i], kf2)
end
c = imfilter(a, kf3);
b == c

gives false. And it's not just floating-point roundoff: the only slice they agree on is c[:,:,10].

[timholy]

Oh, duh, this isn't a bug. Watch what happens with a smaller and easier-to-understand example:

julia> a = reshape(1:24, 3, 4, 2)
3×4×2 reshape(::UnitRange{Int64}, 3, 4, 2) with eltype Int64:
[:, :, 1] =
 1  4  7  10
 2  5  8  11
 3  6  9  12

[:, :, 2] =
 13  16  19  22
 14  17  20  23
 15  18  21  24

julia> b = similar(a);

julia> kf2 = kernelfactors(([1.0], [1.0]))
(ImageFiltering.KernelFactors.ReshapedOneD{Float64,2,0,Array{Float64,1}}([1.0]), ImageFiltering.KernelFactors.ReshapedOneD{Float64,2,1,Array{Float64,1}}([1.0]))

julia> kf3 =  kernelfactors(([1.0], [1.0], [1.0]))
(ImageFiltering.KernelFactors.ReshapedOneD{Float64,3,0,Array{Float64,1}}([1.0]), ImageFiltering.KernelFactors.ReshapedOneD{Float64,3,1,Array{Float64,1}}([1.0]), ImageFiltering.KernelFactors.ReshapedOneD{Float64,3,2,Array{Float64,1}}([1.0]))

julia> for i in 1:size(b, 3)
              b[:,:,i] = imfilter(a[:,:,i], kf2, Fill(0))
       end

julia> imfilter(a, kf3, Fill(0))
3×4×2 Array{Float64,3}:
[:, :, 1] =
 17.0  20.0  23.0  0.0
 18.0  21.0  24.0  0.0
  0.0   0.0   0.0  0.0

[:, :, 2] =
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> b
3×4×2 Array{Int64,3}:
[:, :, 1] =
 5  8  11  0
 6  9  12  0
 0  0   0  0

[:, :, 2] =
 17  20  23  0
 18  21  24  0
  0   0   0  0

Both of these are just following the formula for a convolution. A key point is that ImageFiltering pays attention to the axes of the kernels (https://juliaimages.org/latest/imagefiltering.html): in particular, [1.0] is indexed starting at 1 so it causes a translation by 1 pixel; in contrast, centered([1.0]) is indexed starting at 0 so there is no translation.

julia> axes([1.0])
(Base.OneTo(1),)

julia> axes(centered([1.0]))
(Base.Slice(0:0),)

You want this:

kf2 = kernelfactors(centered.((ones(3)/3, ones(3)/3)))
kf3 =  kernelfactors(centered.((ones(3)/3, ones(3)/3, [1.0])))

and then they will agree.

[aramirezreyes]

Interesting. Thanks for taking a look.