wip

Author: IanButterworth

demo.jl

@@ -27,10 +27,12 @@ function test(mats)
         r_noncached = CPU1(Algorithm.FFT())
         for i in axes(mat, 3)
             frame = @view mat[:, :, i]
+            @info "imfilter! noncached"
             imfilter!(r_noncached, f2, frame, kernel)
+            @info "imfilter! cached"
             imfilter!(r_cached, f1, frame, kernel)
-            @show f1[1:2, 1:2] f2[1:2, 1:2]
-            all(f1 .≈ f2) || error("f1 !≈ f2")
+            @show f1[1:4] f2[1:4]
+            f1 ≈ f2 || error("f1 !≈ f2")
         end
         return
     end
@@ -58,4 +60,20 @@ function profile()
     # GC.gc(true)
 end
 
-profile()
+profile()
+
+using ImageFiltering
+using ImageFiltering.RFFT
+
+function mwe()
+    a = rand(Float64, 10, 10)
+    out1 = rfft(a)
+
+    buf = RFFT.RCpair{Float64}(undef, size(a))
+    rfft_plan = RFFT.plan_rfft!(buf)
+    copy!(buf, a)
+    out2 = complex(rfft_plan(buf))
+
+    return out1 ≈ out2
+end
+mwe()

src/imfilter.jl

@@ -829,65 +829,95 @@ function _imfilter_fft!(r::AbstractCPU{FFT},
     Af = filtfft(A, krn, r.settings.plan1, r.settings.plan2, r.settings.plan3)
     if map(first, axes(out)) == map(first, axes(Af))
         R = CartesianIndices(axes(out))
+        # @info "here" axes(out) axes(Af)
         copyto!(out, R, Af, R)
     else
         # Exploit the periodic boundary conditions of FFTView
         dest = FFTView(out)
         # src = OffsetArray(view(FFTView(Af), axes(dest)...), axes(dest))
         src = view(FFTView(Af), axes(dest)...)
+        @info "there"
         copyto!(dest, src)
     end
-    out
+    return out
 end
 
 function buffered_planned_rfft(a::AbstractArray{T}) where {T<:AbstractFloat}
     buf = RFFT.RCpair{T}(undef, size(a))
-    plan = RFFT.plan_rfft!(buf)
-    return function (a::AbstractArray)
-        copy!(buf, OffsetArrays.no_offset_view(a))
+    plan = RFFT.plan_rfft!(buf; flags=FFTW.MEASURE)
+    return function (arr::AbstractArray{T}) where {T<:AbstractFloat}
+        copy!(buf, FFTW.AbstractFFTs.to1(arr))
         return plan(buf)
     end
 end
-function buffered_planned_irfft(a::AbstractArray{T}) where {T<:AbstractFloat}
+function buffered_planned_irfft(a::AbstractArray{T}) where {T}
+    # @show size(a)
     buf = RFFT.RCpair{T}(undef, size(a))
-    plan = RFFT.plan_irfft!(buf)
-    return function (a::AbstractArray)
-        copy!(buf, a)
+    plan = RFFT.plan_irfft!(buf; flags=FFTW.MEASURE)
+    return function (arr::AbstractArray{T}) where {T<:Complex}
+        # @info "inner" size(buf.R) size(buf.C) size(a)
+        copy!(buf, FFTW.AbstractFFTs.to1(arr))
         return plan(buf)
     end
 end
 
 function planned_fft(A::AbstractArray{T,N},
-                    kernel::Tuple{AbstractArray,Vararg{AbstractArray}},
-                    border::BorderSpecAny=Pad(:replicate)) where {T,N}
+    kernel::Tuple{AbstractArray,Vararg{AbstractArray}},
+    border::BorderSpecAny=Pad(:replicate)) where {T,N}
     bord = border(kernel, A, Algorithm.FFT())
     _A = padarray(T, A, bord)
+    # @info "bfp1" size(_A)
     bfp1 = buffered_planned_rfft(_A)
-
+    # B = complex(bfp1(_A))
     kern = samedims(_A, kernelconv(kernel...))
     krn = FFTView(zeros(eltype(kern), map(length, axes(_A))))
+    for I in CartesianIndices(axes(kern))
+        krn[I] = kern[I]
+    end
+    # @info "bfp2" size(B)
+    @show krn[1:3] sum(krn)
+    @show size(_A) size(krn)
     bfp2 = buffered_planned_rfft(krn)
+    # B .*= conj!(complex(bfp2(real(krn))))
+    # @info "bfp3" size(B) size(krn) size(A) size(_A) size(conj!(_A))
     bfp3 = buffered_planned_irfft(_A)
+    # @info "return"
     return Algorithm.FFT(bfp1, bfp2, bfp3)
 end
 
 function filtfft(A, krn, planned_rfft1::Function, planned_rfft2::Function, planned_irfft::Function)
-    B = complex(planned_rfft1(A)) * FFTW.AbstractFFTs.to1(A)
-    _B = rfft(A)
-    @show B[1:4] _B[1:4]
+    B = complex(planned_rfft1(A))
+    # @info "filtfft" size(A) size(FFTW.AbstractFFTs.to1(A)) size(B)
+    # @info "filtfft planned" FFTW.AbstractFFTs.to1(A)[1:2] B[1:2] typeof(krn) axes(krn)
+    @show krn[1:3] sum(krn)
+    @info "bad one 1" complex(planned_rfft2(krn))[1:3]
     B .*= conj!(complex(planned_rfft2(krn)))
-    out = irfft(B, length(axes(A, 1)))
-    @show complex(B)[1:2,1:2] axes(complex(B)) out[1:2,1:2]
-    return out
-    # return real(planned_irfft(complex(B)))
+    Af = real(planned_irfft(complex(B)))
+
+    #
+    buf = RFFT.RCpair{Float64}(undef, size(krn))
+    _rfft_plan = RFFT.plan_rfft!(buf)
+    copy!(buf, FFTW.AbstractFFTs.to1(krn))
+    _k = complex(_rfft_plan(buf))
+    @info "bad one 2" _k[1:3]
+    #
+
+    @info "filtfft planned" B[1:2] Af[1:2]
+    # TODO: convert Af back to a padded view so that the center of the image is copied in _imfilter_fft!
+    # @info "filtfft planned return" typeof(A) size(A) size(B) typeof(Af) size(Af)
+    return Af
 end
 filtfft(A, krn, ::Nothing, ::Nothing, ::Nothing) = filtfft(A, krn)
 function filtfft(A, krn)
     B = rfft(A)
+    # @info "filtfft unplanned" FFTW.AbstractFFTs.to1(A)[1:2] B[1:2] typeof(krn) axes(krn)
+    # @show krn[1:3] sum(krn)
+    @info "good one" rfft(krn)[1:3]
     B .*= conj!(rfft(krn))
-    out = irfft(B, length(axes(A, 1)))
-    @show B[1:2,1:2] axes(B) out[1:2,1:2]
-    return out
+    Af = irfft(B, length(axes(A, 1)))
+    @info "filtfft unplanned" B[1:2] Af[1:2]
+    # @info "filtfft unplanned return" typeof(A) size(A) size(B) typeof(Af) size(Af)
+    return Af
 end
 function filtfft(A::AbstractArray{C}, krn) where {C<:Colorant}
     Av, dims = channelview_dims(A)