fix freqkernel (#207)

Author: JeffFessler

Project.toml

@@ -1,7 +1,7 @@
 name = "ImageFiltering"
 uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 author = ["Tim Holy <[email protected]>", "Jan Weidner <[email protected]>"]
-version = "0.6.19"
+version = "0.6.20"
 
 [deps]
 CatIndices = "aafaddc9-749c-510e-ac4f-586e18779b91"
@@ -11,7 +11,6 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 FFTViews = "4f61f5a4-77b1-5117-aa51-3ab5ef4ef0cd"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
-ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"

src/ImageFiltering.jl

@@ -3,6 +3,7 @@ module ImageFiltering
 using FFTW
 using ImageCore, FFTViews, OffsetArrays, StaticArrays, ComputationalResources, TiledIteration
 # Where possible we avoid a direct dependency to reduce the number of [compat] bounds
+# using FixedPointNumbers: Normed, N0f8 # reexported by ImageCore
 using ImageCore.MappedArrays
 using Statistics, LinearAlgebra
 using ColorVectorSpace  # for filtering RGB arrays

src/mapwindow.jl

@@ -133,6 +133,9 @@ resolve_window(window::AbstractArray) = resolve_window((window...,))
 resolve_window(window::AbstractUnitRange) = (window,)
 resolve_window(window::Indices) = window
 
+# avoid method ambiguity between ::Dims and ::Indices
+resolve_window(window::Tuple{}) = throw(ArgumentError("empty window"))
+
 resolve_border(border::AbstractString) = borderinstance(border)
 resolve_border(border::BorderSpecAny) = border
 

src/utils.jl

@@ -1,9 +1,8 @@
 """
     shiftedkernel = centered(kernel)
 
-Shift the origin-of-coordinates to the center of `kernel`. The
-center-element of `kernel` will be accessed by `shiftedkernel[0, 0,
-...]`.
+Shift the origin-of-coordinates to the center of `kernel`.
+The center-element of `kernel` will be accessed by `shiftedkernel[0, 0, ...]`.
 
 This function makes it easy to supply kernels using regular Arrays,
 and provides compatibility with other languages that do not support
@@ -21,20 +20,35 @@ end
     kernfft = freqkernel([T::Type], kern, sz=size(kern); rfft=false)
 
 Return a frequency-space representation of `kern`.
-This embeds `kern` in an array of size `sz`, in a manner that implicitly imposes periodic boundary
-conditions, and then returns the fourier transform. This is sometimes called the optical transfer
-function, and known in some frameworks as `psf2otf`. If `rfft` is `true`, the fft for real-valued
-arrays (`rfft`) is returned instead and the first dimension size will be approximately half of `sz[1]`.
-
-`kern` should be zero-centered, i.e., `kern[0, 0]` should reference the center of your kernel.
-See [`centered`](@ref). Optionally specify the numeric type `T` (which must be one of the
-types supported by FFTW, either `Float32` or `Float64`).
+This embeds `kern` in an array of size `sz`,
+in a manner that implicitly imposes periodic boundary conditions,
+and then returns the Fourier transform (frequency response).
+This is sometimes called the optical transfer function,
+and is known in some frameworks as `psf2otf`.
+If `rfft` is `true`, the FFT for real-valued arrays (`rfft`) is returned instead
+and the first dimension size will be approximately half of `sz[1]`.
+
+`kern` should be zero-centered, i.e.,
+`kern[0, 0]` should reference the center of your kernel,
+and `sz` must be large enough to support `kern`.
+See [`centered`](@ref).
+Optionally specify the numeric type `T`
+(which must be one of the types supported by FFTW,
+either `Float32` or `Float64`).
 
 The inverse of `freqkernel` is [`spacekernel`](@ref).
 """
 function freqkernel(::Type{T}, kern::AbstractArray, sz::Dims=size(kern); rfft=false) where T<:Union{Float32,Float64}
-    wrapindex(i, s) = i<1 ? i+s : i
-    all(size(kern) .<= sz) || throw(DimensionMismatch("kernel size $(size(kern)) is bigger than supplied size $sz"))
+    wrapindex(i, s) = 1 + (i<0 ? i+s : i)
+    all(size(kern) .<= sz) ||
+        throw(DimensionMismatch("kernel size $(size(kern)) exceeds supplied size $sz"))
+    rhs = Tuple(last(CartesianIndices(kern)))
+    lhs = Tuple(first(CartesianIndices(kern)))
+    limit = collect((sz .+ 1) .÷ 2) # handle odd and even sizes
+    all(rhs .< limit) ||
+        throw(DimensionMismatch("kernel last index $rhs >= limit $limit"))
+    all(-limit .<= lhs) ||
+        throw(DimensionMismatch("kernel first index $lhs < limit $(-limit)"))
     kernw = zeros(T, sz...)
     for I in CartesianIndices(kern)
         J = CartesianIndex(map(wrapindex, Tuple(I), sz))
@@ -48,17 +62,22 @@ freqkernel(kern::AbstractArray{T}, args...; rfft=false) where T =
 """
     kern = spacekernel(kernfft, axs; rfftsz=0)
 
-Return a real-space representation of `kernfft`, a frequency-space representation of a kernel.
-This performs an inverse fourier transform, implicitly imposes periodic boundary conditions,
-and then trims & truncates axes of the output to `axs`. By default `kernfft` is assumed to have
-been generated by `fft`; if it was instead generated by `rfft`, specify the original size
-of the first dimension. (If `kernfft` was generated by [`freqkernel`](@ref), this is just `sz[1]`.)
+Return a real-space representation of `kernfft`,
+the frequency-space representation of a kernel.
+This performs an inverse Fourier transform,
+implicitly imposes periodic boundary conditions,
+and then trims & truncates axes of the output to `axs`.
+By default `kernfft` is assumed to have been generated by `fft`;
+if it was instead generated by `rfft`,
+the specify the original size of the first dimension.
+(If `kernfft` was generated by [`freqkernel`](@ref), this is just `sz[1]`.)
 
 The inverse of `spacekernel` is [`freqkernel`](@ref).
 """
 function spacekernel(kernfft::AbstractArray, axs::Indices; rfftsz=0)
-    wrapindex(i, s) = i<1 ? i+s : i
+    wrapindex(i, s) = 1 + (i<0 ? i+s : i)
     kernw = rfftsz > 0 ? irfft(kernfft, rfftsz) : ifft(kernfft)
+    # there could be some checking of axs vs size(kernfft) here
     kern = zeros(eltype(kernw), axs...)
     sz = size(kernw)
     for I in CartesianIndices(kern)

test/basic.jl

@@ -109,6 +109,15 @@ end
     k2 = spacekernel(kfft, axes(k); rfftsz=31)
     @test k2 ≈ k
 
+    k = centered([0, 1, 0] * [0, 1, 0]') # Kronecker impulse
+    @test freqkernel(k, (8, 8)) ≈ ones(8,8) # flat frequency response
+
+    @test_throws DimensionMismatch freqkernel(ones(5), (3,)) # kernel too big
+    k = OffsetArray(ones(3), (-5,)) # index -4:-2 too far left
+    @test_throws DimensionMismatch freqkernel(k, (5,))
+    k = OffsetArray(ones(3), (+1,)) # index 2:4 too far right
+    @test_throws DimensionMismatch freqkernel(k, (5,))
+
     for T in (Float64, Float32, Float16, N0f16)
         k = T.(Kernel.gaussian(3))
         kfft = freqkernel(k)

test/mapwindow.jl

@@ -71,6 +71,9 @@ using ImageFiltering: IdentityUnitRange
         Amin = groundtruth(min, A, w)
         @test minval == Amin
     end
+
+    @test_throws ArgumentError mapwindow(sum, ones(5,5), ()) # resolve_window
+
     # offsets
     @testset "offsets" begin
         n = 5