From a79e8078d3adc4481db3b66ce5990b390bddecc7 Mon Sep 17 00:00:00 2001 From: Ian Butterworth Date: Sat, 9 Mar 2024 00:11:06 -0500 Subject: [PATCH] wip --- demo.jl | 24 ++++++++++++++--- src/imfilter.jl | 72 ++++++++++++++++++++++++++++++++++--------------- 2 files changed, 72 insertions(+), 24 deletions(-) diff --git a/demo.jl b/demo.jl index 609dd43..7743112 100644 --- a/demo.jl +++ b/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() \ No newline at end of file +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() \ No newline at end of file diff --git a/src/imfilter.jl b/src/imfilter.jl index f4d69de..7f3e596 100644 --- a/src/imfilter.jl +++ b/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)