Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix gaussian integral #46

Merged
merged 4 commits into from
Sep 23, 2024
Merged

Fix gaussian integral #46

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
cc71ff5
bug fix with normalization in real-space Gaussian filtering
RainerHeintzmann Sep 16, 2024
9e9f30d
test update
RainerHeintzmann Sep 16, 2024
ec61ecb
bug fix
RainerHeintzmann Sep 16, 2024
6f60872
renamed to `normalize_win`, bug fixes and docstring, more tests
RainerHeintzmann Sep 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
109 changes: 94 additions & 15 deletions src/fourier_filtering.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,19 @@
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}}
"""
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}}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe transform_win and normalize_win.
I think being more verbose would be better here since win is confusing

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry. Don't understand your comment. It is called "normalize_win" and you suggest the same?
We should definitely not change transform_win to stay backwards compatible. So I agree that "win" should be more verbose but consistency is also important.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok didn't see this was already in the old version


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
Expand All @@ -100,17 +112,42 @@
fft!(arr, d)
arr .*= let
if transform_win
fft(wins[d], d)
tmp = fft(wins[d], 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]

Check warning on line 125 in src/fourier_filtering.jl

View check run for this annotation

Codecov / codecov/patch

src/fourier_filtering.jl#L124-L125 

Added lines #L124 - L125 were not covered by tests
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this line is not covered by tests

end
else
wins[d]
end
end
end
end
ifft!(arr, dims)
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}}
"""
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
Expand All @@ -124,10 +161,22 @@
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, 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}}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

windows instead of wins

if isempty(dims)
return arr
end
Expand All @@ -138,23 +187,48 @@
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]

Check warning on line 200 in src/fourier_filtering.jl

View check run for this annotation

Codecov / codecov/patch

src/fourier_filtering.jl#L198-L200 

Added lines #L198 - L200 were not covered by tests
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, 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
Expand All @@ -175,8 +249,15 @@
win .= fct(real(eltype(arr)), select_sizes(arr_ft,d), offset=CtrRFFT, scale=2 ./size(arr,d); kwargs...)
end
end
if (normalize_win)
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...)

# 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, 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)
Expand Down Expand Up @@ -310,9 +391,7 @@
"""
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, normalize_win=true, sigma=sigma, kwargs...)
return arr
else
return fourier_filter!(arr; border_in=border_in, border_out=border_out, kwargs...)
Expand Down
14 changes: 11 additions & 3 deletions test/fourier_filtering.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,23 +14,31 @@ 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, (1.0,), 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
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
@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
Loading