Output image incorrect when OffsetArray vectors passed as kernel components

[BioTurboNick]

If a 2d kernel is passed as two OffsetArray vectors, the second one transposed or adjoint, the final image is incorrect.

The issue does not occur if the second component is an OffsetArray matrix with one row, or if the vector is transposed /adjointed before being passed to centered.

image = [0.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0]
image_out_typical = imfilter(image, Kernel.gaussian(1))

factor1, factor2 = ImageFiltering.factorkernel(Kernel.gaussian(1))
@info factor1
@info factor2
	
image_out_factors = imfilter(image, (factor1, factor2))

factor1a = factor1 |> vec |> centered
factor2a = (factor2 |> vec |> centered)'

@info factor1a
@info factor2a

image_out_factors_vec = imfilter(image, (factor1a, factor2a))
	
plot(heatmap(image, aspectratio = 1), heatmap(image_out_typical, aspectratio = 1), heatmap(image_out_factors, aspectratio = 1), heatmap(image_out_factors_vec, aspectratio = 1))

[BioTurboNick]

I've traced the first difference down to here:

https://github.com/JuliaImages/ImageFiltering.jl/blob/167611b1ade92d5d4449a6f1cf7b93fb92946a6b/src/imfilter.jl#L326C9-L326C9

With the aberrant case, this line returns Pad(:replicate, (1, 2), (3, 2)); it should be Pad(:replicate, (2, 2), (2, 2))

And further...

image = zeros(3, 3)
kernel1 = factor1, factor2 = ImageFiltering.factorkernel(Kernel.gaussian(1))
factor2a = (factor2 |> vec |> centered)'
kernel2 = (factor1, factor2a)

ImageFiltering.expand(axes(kernel1[1]), axes(kernel1[2]))
# (OffsetArrays.IdOffsetRange(values=-2:2, indices=-2:2), OffsetArrays.IdOffsetRange(values=-2:2, indices=-2:2))

ImageFiltering.expand(axes(kernel2[1]), axes(kernel2[2]))
# (OffsetArrays.IdOffsetRange(values=-1:3, indices=-1:3), OffsetArrays.IdOffsetRange(values=-2:2, indices=-2:2))

And even further:

a = centered(zeros(3))'
b = centered(zeros(3)')

axes(a)
# (Base.OneTo(1), OffsetArrays.IdOffsetRange(values=-1:1, indices=-1:1))
# axes(v::Union{LinearAlgebra.Adjoint{T, var"#s972"}, LinearAlgebra.Transpose{T, var"#s972"}} where {T, var"#s972"<:(AbstractVector)}) in LinearAlgebra at src/adjtrans.jl:298

axes(b)
# (OffsetArrays.IdOffsetRange(values=0:0, indices=0:0), OffsetArrays.IdOffsetRange(values=-1:1, indices=-1:1))
# axes(A::OffsetArrays.OffsetArray) in OffsetArrays at src/OffsetArrays.jl:295

[BioTurboNick]

Since the root issue is that the Julia axes function doesn't account for OffsetArrays, this is a bit awkward. Fixing this would be arguably breaking because old code may have been relying on the behavior, which is how I discovered it.

Perhaps for now, ImageFiltering could manually check for this condition and provide the correct axes. This would be a breaking release. And can add a note that it's a fix for underlying Julia behavior, and that it should be revisited once Julia contains the true fix; at a minimum to make the method conditional on the installed Julia version.

[BioTurboNick]

So looks like there need to be several parts to fix this:

1. Julia needs to not assume that all higher dimensions have a range of 1:1
2. OffsetArrays needs to support specifying the offset for all higher dimensions.
3. OffsetArrays.centered needs to use that feature to make an object that has 0:0 axes in all higher dimensions.

Crossref: JuliaArrays/OffsetArrays.jl#339

[mkitti]

Since the root issue is that the Julia axes function doesn't account for OffsetArrays, this is a bit awkward. Fixing this would be arguably breaking because old code may have been relying on the behavior, which is how I discovered it.

Just because some behavior change does not necessarily make it "breaking". If a documented behavior changed, then it would definitely be breaking. However, if there is a bug, and people have been relying on that bug, but you fixed it then this is a bug fix. That said, it would be nice if you created some facilities to help people adapt to the fix if it disrupts their software.