From cc71ff5c5dc71d0ada674c43764b223c2bcbc64a Mon Sep 17 00:00:00 2001 From: RainerHeintzmann Date: Mon, 16 Sep 2024 19:51:54 +0200 Subject: [PATCH 1/4] bug fix with normalization in real-space Gaussian filtering --- src/fourier_filtering.jl | 33 ++++++++++++++++++++++++--------- test/fourier_filtering.jl | 6 +++++- 2 files changed, 29 insertions(+), 10 deletions(-) diff --git a/src/fourier_filtering.jl b/src/fourier_filtering.jl index 98454b1..291a3a4 100644 --- a/src/fourier_filtering.jl +++ b/src/fourier_filtering.jl @@ -90,7 +90,7 @@ function fourier_filter(arr, fct=window_gaussian; kwargs...) return fourier_filter!(copy(arr), fct; kwargs...) end -function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} +function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=false, keep_integral=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} if isempty(dims) return arr end @@ -100,9 +100,21 @@ function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=f fft!(arr, d) arr .*= let if transform_win - fft(wins[d], d) + tmp = fft(wins[d], d) + if (keep_integral) + if (tmp[1] != 0 && tmp[1] != 1) + tmp ./= tmp[1] + end + end + tmp else - wins[d] + if (keep_integral) + if (tmp[1] != 0 && tmp[1] != 1) + wins[d] ./ wins[d][1] + end + else + wins[d] + end end end end @@ -110,7 +122,7 @@ function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=f return arr end -function fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, kwargs...) where {N, TA <: AbstractArray{<:Complex, N}} +function fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, keep_integral=false, kwargs...) where {N, TA <: AbstractArray{<:Complex, N}} if isempty(dims) return arr end @@ -124,7 +136,7 @@ function fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; dims=(1:ndims(ar win .= fct(real(eltype(arr)), select_sizes(arr, d); kwargs...) wins[d] = ifftshift(win) end - return fourier_filter_by_1D_FT!(arr, wins; transform_win=transform_win, dims=dims) + return fourier_filter_by_1D_FT!(arr, wins; transform_win=transform_win, keep_integral=keep_integral, dims=dims) end function fourier_filter_by_1D_RFT!(arr::TA, wins::AbstractVector; dims=(1:ndims(arr)), transform_win=false, kwargs...) where {T<:Real, N, TA<:AbstractArray{T, N}} @@ -154,7 +166,7 @@ function fourier_filter_by_1D_RFT!(arr::TA, wins::AbstractVector; dims=(1:ndims( end # transforms the first dim as rft and then hands over to the fft-based routines. -function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, kwargs...) where {T<:Real, N, TA<:AbstractArray{T, N}} +function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, keep_integral=false, kwargs...) where {T<:Real, N, TA<:AbstractArray{T, N}} if isempty(dims) return arr end @@ -175,6 +187,11 @@ function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(a win .= fct(real(eltype(arr)), select_sizes(arr_ft,d), offset=CtrRFFT, scale=2 ./size(arr,d); kwargs...) end end + if (keep_integral) + if (win[1] != 0 && win[1] != 1) + win ./= win[1] + end + end arr_ft .*= win fourier_filter_by_1D_FT!(arr_ft, fct; dims=dims[2:end], transform_win=transform_win, kwargs...) # go back to real space now and return because shift_by_1D_FT processed @@ -310,9 +327,7 @@ See also `filter_gaussian()` and `fourier_filter!()`. """ function filter_gaussian!(arr, sigma=eltype(arr)(1); real_space_kernel=true, border_in=(real(eltype(arr)))(0), border_out=(real(eltype(arr))).(2 ./ (pi .* sigma)), kwargs...) if real_space_kernel - mysum = sum(arr) - fourier_filter!(arr, gaussian; transform_win=true, sigma=sigma, kwargs...) - arr .*= (mysum/sum(arr)) + fourier_filter!(arr, gaussian; transform_win=true, keep_integral=true, sigma=sigma, kwargs...) return arr else return fourier_filter!(arr; border_in=border_in, border_out=border_out, kwargs...) diff --git a/test/fourier_filtering.jl b/test/fourier_filtering.jl index a351946..ea95c17 100644 --- a/test/fourier_filtering.jl +++ b/test/fourier_filtering.jl @@ -14,7 +14,9 @@ Random.seed!(42) gfc = conv_psf(x, k) @test ≈(gf,gfc, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places gfr = filter_gaussian(x, sigma, real_space_kernel=true) - @test ≈(gfr,gfc) # it can be debated how to best normalize a Gaussian filter + @test ≈(gfr, gfc) # it can be debated how to best normalize a Gaussian filter + gfr = filter_gaussian(zeros(5).+1im, sigma, real_space_kernel=true) + @test ≈(gfr, zeros(5).+1im) # it can be debated how to best normalize a Gaussian filter end @testset "Gaussian filter real" begin @@ -28,6 +30,8 @@ Random.seed!(42) k = k./sum(k) # different than "normal". gf2 = conv_psf(x, k) @test ≈(gf,gf2, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places + gf2 = filter_gaussian(zeros(5), sigma, real_space_kernel=true) + @test ≈(gf2, zeros(5)) # it can be debated how to best normalize a Gaussian filter end @testset "Other filters" begin @test filter_hamming(FourierTools.delta(Float32, (3,)), border_in=0.0, border_out=1.0) ≈ [0.23,0.54, 0.23] From 9e9f30d3654a6bd5833ad5a45cc8fd1b5744cbdc Mon Sep 17 00:00:00 2001 From: RainerHeintzmann Date: Mon, 16 Sep 2024 19:59:03 +0200 Subject: [PATCH 2/4] test update --- test/fourier_filtering.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/fourier_filtering.jl b/test/fourier_filtering.jl index ea95c17..19fc583 100644 --- a/test/fourier_filtering.jl +++ b/test/fourier_filtering.jl @@ -15,7 +15,7 @@ Random.seed!(42) @test ≈(gf,gfc, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places gfr = filter_gaussian(x, sigma, real_space_kernel=true) @test ≈(gfr, gfc) # it can be debated how to best normalize a Gaussian filter - gfr = filter_gaussian(zeros(5).+1im, sigma, real_space_kernel=true) + gfr = filter_gaussian(zeros(5).+1im, (1.0,), real_space_kernel=true) @test ≈(gfr, zeros(5).+1im) # it can be debated how to best normalize a Gaussian filter end @@ -30,8 +30,8 @@ Random.seed!(42) k = k./sum(k) # different than "normal". gf2 = conv_psf(x, k) @test ≈(gf,gf2, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places - gf2 = filter_gaussian(zeros(5), sigma, real_space_kernel=true) - @test ≈(gf2, zeros(5)) # it can be debated how to best normalize a Gaussian filter + gf2 = filter_gaussian(zeros(sz), sigma, real_space_kernel=true) + @test ≈(gf2, zeros(sz)) # it can be debated how to best normalize a Gaussian filter end @testset "Other filters" begin @test filter_hamming(FourierTools.delta(Float32, (3,)), border_in=0.0, border_out=1.0) ≈ [0.23,0.54, 0.23] From ec61ecbb3ee901893542817f47cc27f6f232c789 Mon Sep 17 00:00:00 2001 From: RainerHeintzmann Date: Mon, 16 Sep 2024 20:19:55 +0200 Subject: [PATCH 3/4] bug fix --- src/fourier_filtering.jl | 4 +++- test/fourier_filtering.jl | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/fourier_filtering.jl b/src/fourier_filtering.jl index 291a3a4..f41dc16 100644 --- a/src/fourier_filtering.jl +++ b/src/fourier_filtering.jl @@ -193,7 +193,9 @@ function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(a end end arr_ft .*= win - fourier_filter_by_1D_FT!(arr_ft, fct; dims=dims[2:end], transform_win=transform_win, kwargs...) + + # hand over to the fft-based routines for all other dimensions + fourier_filter_by_1D_FT!(arr_ft, fct; dims=dims[2:end], transform_win=transform_win, keep_integral=keep_integral, kwargs...) # go back to real space now and return because shift_by_1D_FT processed # the other dimensions already mul!(arr, inv(p), arr_ft) diff --git a/test/fourier_filtering.jl b/test/fourier_filtering.jl index 19fc583..db7c0ee 100644 --- a/test/fourier_filtering.jl +++ b/test/fourier_filtering.jl @@ -22,14 +22,14 @@ Random.seed!(42) @testset "Gaussian filter real" begin sz = (21, 22) x = randn(Float32, sz) - sigma = (1.1,2.2) + sigma = (1.1, 2.2) gf = filter_gaussian(x, sigma, real_space_kernel=true) # Note that this is not the same, since one kernel is generated in real space and one in Fourier space! # with sizes around 10, the difference is huge! k = gaussian(sz, sigma=sigma) k = k./sum(k) # different than "normal". gf2 = conv_psf(x, k) - @test ≈(gf,gf2, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places + @test ≈(gf, gf2, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places gf2 = filter_gaussian(zeros(sz), sigma, real_space_kernel=true) @test ≈(gf2, zeros(sz)) # it can be debated how to best normalize a Gaussian filter end From 6f60872124ee9d756ce20d862892f6e5cc2b9b73 Mon Sep 17 00:00:00 2001 From: RainerHeintzmann Date: Tue, 17 Sep 2024 07:56:34 +0200 Subject: [PATCH 4/4] renamed to `normalize_win`, bug fixes and docstring, more tests --- src/fourier_filtering.jl | 92 ++++++++++++++++++++++++++++++++------- test/fourier_filtering.jl | 4 ++ 2 files changed, 81 insertions(+), 15 deletions(-) diff --git a/src/fourier_filtering.jl b/src/fourier_filtering.jl index f41dc16..28b9375 100644 --- a/src/fourier_filtering.jl +++ b/src/fourier_filtering.jl @@ -90,7 +90,19 @@ function fourier_filter(arr, fct=window_gaussian; kwargs...) return fourier_filter!(copy(arr), fct; kwargs...) end -function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=false, keep_integral=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} +""" + fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=false, normalize_win=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} + +filters an array by sequential multiplication in Fourierspace using one-dimensional FFTs with one-dimensional kernels as specified in `wins`. + +#Arguments ++ `arr`: the array to filter by separable multiplication with windows. ++ `wins`: the window as a vector of vectors. Each of the original array dimensions corresponds to an entry in the outer vector. ++ `transform_win`: specifies whether the directional windows each need to be FFTed before applying. ++ `normalize_win`: specifies whether the directional windows (after potential FFT) each need to be normalized to a value of 1 at the zero frequency coordinate. ++ `dims`: the dimensions to apply this separable multiplication to. If empty, the array will not be filtered. +""" +function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=false, normalize_win=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} if isempty(dims) return arr end @@ -101,15 +113,15 @@ function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=f arr .*= let if transform_win tmp = fft(wins[d], d) - if (keep_integral) + if (normalize_win) if (tmp[1] != 0 && tmp[1] != 1) tmp ./= tmp[1] end end tmp else - if (keep_integral) - if (tmp[1] != 0 && tmp[1] != 1) + if (normalize_win) + if (wins[d][1] != 0 && wins[d][1] != 1) wins[d] ./ wins[d][1] end else @@ -122,7 +134,20 @@ function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=f return arr end -function fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, keep_integral=false, kwargs...) where {N, TA <: AbstractArray{<:Complex, N}} +""" + fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; transform_win=false, normalize_win=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} + +filters an array by sequential multiplication in Fourierspace using one-dimensional FFTs with one-dimensional kernels as specified in `wins`. + +#Arguments ++`arr`: the array to filter by separable multiplication with windows. ++`fct`: the window as a fuction. It needs to accept a datatype as the first argument and a size as a second argument. + If specified in real space (`transform_win=true`), it should interpret the zero position near the center of the array. ++`transform_win`: specifies whether the directional windows each need to be FFTed before applying. ++`normalize_win`: specifies whether the directional windows (after potential FFT) each need to be normalized to a value of 1 at the zero frequency coordinate. ++`dims`: the dimensions to apply this separable multiplication to. If empty, the array will not be filtered. +""" +function fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, normalize_win=false, kwargs...) where {N, TA <: AbstractArray{<:Complex, N}} if isempty(dims) return arr end @@ -136,10 +161,22 @@ function fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; dims=(1:ndims(ar win .= fct(real(eltype(arr)), select_sizes(arr, d); kwargs...) wins[d] = ifftshift(win) end - return fourier_filter_by_1D_FT!(arr, wins; transform_win=transform_win, keep_integral=keep_integral, dims=dims) + return fourier_filter_by_1D_FT!(arr, wins; transform_win=transform_win, normalize_win=normalize_win, dims=dims) end -function fourier_filter_by_1D_RFT!(arr::TA, wins::AbstractVector; dims=(1:ndims(arr)), transform_win=false, kwargs...) where {T<:Real, N, TA<:AbstractArray{T, N}} +""" + fourier_filter_by_1D_RFT!(arr::TA, wins::AbstractVector; transform_win=false, normalize_win=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} + +filters an array by sequential multiplication in Fourierspace using one-dimensional RFFTs (first transform dimension) and FFTS (other transform dimensions) with one-dimensional kernels as specified in `wins`. + +#Arguments ++`arr`: the array to filter by separable multiplication with windows. ++`wins`: the window as a vector of vectors. Each of the original array dimensions corresponds to an entry in the outer vector. ++`transform_win`: specifies whether the directional windows each need to be FFTed before applying. ++`normalize_win`: specifies whether the directional windows (after potential FFT) each need to be normalized to a value of 1 at the zero frequency coordinate. ++`dims`: the dimensions to apply this separable multiplication to. If empty, the array will not be filtered. +""" +function fourier_filter_by_1D_RFT!(arr::TA, wins::AbstractVector; dims=(1:ndims(arr)), transform_win=false, normalize_win=false, kwargs...) where {T<:Real, N, TA<:AbstractArray{T, N}} if isempty(dims) return arr end @@ -150,23 +187,48 @@ function fourier_filter_by_1D_RFT!(arr::TA, wins::AbstractVector; dims=(1:ndims( arr_ft .*= let if transform_win pw = plan_rfft(wins[d], d) - pw * wins[d] + tmp = pw * wins[d] + if (normalize_win) + if (tmp[1] != 0 && tmp[1] != 1) + tmp ./= tmp[1] + end + end + tmp else - wins[d] + if (normalize_win) + if (wins[d][1] != 0 && wins[d][1] != 1) + wins[d] ./ wins[d][1] + end + else + wins[d] + end end end # since we now did a single rfft dim, we can switch to the complex routine # new_limits = ntuple(i -> i ≤ d ? 0 : limits[i], N) # workaround to mimic in-place rfft - fourier_filter_by_1D_FT!(arr_ft, wins; dims=dims[2:end], transform_win=transform_win, kwargs...) + fourier_filter_by_1D_FT!(arr_ft, wins; dims=dims[2:end], transform_win=transform_win, normalize_win=normalize_win, kwargs...) # go back to real space now and return because shift_by_1D_FT processed # the other dimensions already mul!(arr, inv(p), arr_ft) return arr end -# transforms the first dim as rft and then hands over to the fft-based routines. -function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, keep_integral=false, kwargs...) where {T<:Real, N, TA<:AbstractArray{T, N}} +""" + fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; transform_win=false, normalize_win=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}} + +filters an array by sequential multiplication in Fourierspace using one-dimensional RFFTs (first transform dimension) and FFTS (other transform dimensions) with one-dimensional kernels as specified in `wins`. + +#Arguments ++`arr`: the array to filter by separable multiplication with windows. ++`fct`: the window as a fuction. It needs to accept a datatype as the first argument and a size as a second argument. + If specified in real space (`transform_win=true`), it should interpret the zero position near the center of the array. ++`transform_win`: specifies whether the directional windows each need to be FFTed before applying. ++`normalize_win`: specifies whether the directional windows (after potential FFT) each need to be normalized to a value of 1 at the zero frequency coordinate. ++`dims`: the dimensions to apply this separable multiplication to. If empty, the array will not be filtered. +""" +function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(arr)), transform_win=false, normalize_win=false, kwargs...) where {T<:Real, N, TA<:AbstractArray{T, N}} + # transforms the first dim as rft and then hands over to the fft-based routines. if isempty(dims) return arr end @@ -187,7 +249,7 @@ function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(a win .= fct(real(eltype(arr)), select_sizes(arr_ft,d), offset=CtrRFFT, scale=2 ./size(arr,d); kwargs...) end end - if (keep_integral) + if (normalize_win) if (win[1] != 0 && win[1] != 1) win ./= win[1] end @@ -195,7 +257,7 @@ function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; dims=(1:ndims(a arr_ft .*= win # hand over to the fft-based routines for all other dimensions - fourier_filter_by_1D_FT!(arr_ft, fct; dims=dims[2:end], transform_win=transform_win, keep_integral=keep_integral, kwargs...) + fourier_filter_by_1D_FT!(arr_ft, fct; dims=dims[2:end], transform_win=transform_win, normalize_win=normalize_win, kwargs...) # go back to real space now and return because shift_by_1D_FT processed # the other dimensions already mul!(arr, inv(p), arr_ft) @@ -329,7 +391,7 @@ See also `filter_gaussian()` and `fourier_filter!()`. """ function filter_gaussian!(arr, sigma=eltype(arr)(1); real_space_kernel=true, border_in=(real(eltype(arr)))(0), border_out=(real(eltype(arr))).(2 ./ (pi .* sigma)), kwargs...) if real_space_kernel - fourier_filter!(arr, gaussian; transform_win=true, keep_integral=true, sigma=sigma, kwargs...) + fourier_filter!(arr, gaussian; transform_win=true, normalize_win=true, sigma=sigma, kwargs...) return arr else return fourier_filter!(arr; border_in=border_in, border_out=border_out, kwargs...) diff --git a/test/fourier_filtering.jl b/test/fourier_filtering.jl index db7c0ee..71fbd2b 100644 --- a/test/fourier_filtering.jl +++ b/test/fourier_filtering.jl @@ -36,5 +36,9 @@ Random.seed!(42) @testset "Other filters" begin @test filter_hamming(FourierTools.delta(Float32, (3,)), border_in=0.0, border_out=1.0) ≈ [0.23,0.54, 0.23] @test filter_hann(FourierTools.delta(Float32, (3,)), border_in=0.0, border_out=1.0) ≈ [0.25,0.5, 0.25] + @test FourierTools.fourier_filter_by_1D_FT!(ones(ComplexF64, 6), [ones(ComplexF64, 6)]; transform_win=true, normalize_win=false) ≈ 6 .* ones(ComplexF64, 6) + @test FourierTools.fourier_filter_by_1D_FT!(ones(ComplexF64, 6), [ones(ComplexF64, 6)]; transform_win=true, normalize_win=true) ≈ ones(ComplexF64, 6) + @test FourierTools.fourier_filter_by_1D_RFT!(ones(6), [ones(6)]; transform_win=true, normalize_win=false) ≈ 6 .* ones(6) + @test FourierTools.fourier_filter_by_1D_RFT!(ones(6), [ones(6)]; transform_win=true, normalize_win=true) ≈ ones(6) end end