Skip to content

Commit 3b5f9d0

Browse files
bug fixes and added the possibility to specify dims
1 parent 8d2ba49 commit 3b5f9d0

File tree

3 files changed

+215
-39
lines changed

3 files changed

+215
-39
lines changed

src/fourier_filtering.jl

Lines changed: 182 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -2,71 +2,214 @@ export fourier_filter!, fourier_filter
22
export filter_gaussian # , filter_gaussian2
33

44
"""
5+
fourier_filter!(arr::AbstractArray, fct=window_gaussian; kwargs...)
6+
7+
filters an array by multiplication in Fourierspace. This version uses in-place Fourier-transforms and multiplication whereever possible.
8+
The filter window function is assumed to be separable. Depending on the array type either full-complex or real-to-complex FFTs will be used.
9+
Note that some array types cannot directly be processed due to the type inference mechanism, in which case you should transform such types to
10+
a standard array type (usually via calling "collect").
11+
High-pass filters can be used by exchanging the border_in and border_out argument values.
12+
13+
#Arguments
14+
+`arr`: the array to filter
15+
+`fct`: the separable window function to use. The function has an ND-array size as the first parameter but wil be called with only one non-singleton dimension at a time.
16+
+`kwargs...`: keyword arguments which will be passed to fct. These typically are `border_in=0.8` and `border_out=1.0` for the window-functions as defined in the `IndexFunArrays.jl` toolbox.
17+
18+
#Example
19+
```jdoctest
20+
julia> arr = zeros(10,10); arr[3,4]=1.0;
21+
22+
julia> fourier_filter!(arr)
23+
10×10 Matrix{Float64}:
24+
-0.00747645 0.00747645 -0.00747645 -0.07899 -0.00747645 … -0.00747645 0.00747645 -0.00747645 0.00747645
25+
0.00747645 -0.00747645 0.00747645 0.07899 0.00747645 0.00747645 -0.00747645 0.00747645 -0.00747645
26+
0.07899 -0.07899 0.07899 0.834544 0.07899 0.07899 -0.07899 0.07899 -0.07899
27+
0.00747645 -0.00747645 0.00747645 0.07899 0.00747645 0.00747645 -0.00747645 0.00747645 -0.00747645
28+
-0.00747645 0.00747645 -0.00747645 -0.07899 -0.00747645 -0.00747645 0.00747645 -0.00747645 0.00747645
29+
0.00747645 -0.00747645 0.00747645 0.07899 0.00747645 … 0.00747645 -0.00747645 0.00747645 -0.00747645
30+
-0.00747645 0.00747645 -0.00747645 -0.07899 -0.00747645 -0.00747645 0.00747645 -0.00747645 0.00747645
31+
0.00747645 -0.00747645 0.00747645 0.07899 0.00747645 0.00747645 -0.00747645 0.00747645 -0.00747645
32+
-0.00747645 0.00747645 -0.00747645 -0.07899 -0.00747645 -0.00747645 0.00747645 -0.00747645 0.00747645
33+
0.00747645 -0.00747645 0.00747645 0.07899 0.00747645 0.00747645 -0.00747645 0.00747645 -0.00747645
34+
```
535
"""
6-
function fourier_filter!(arr::AbstractArray{<:Complex, N}, fct=window_gaussian; border_in=0.0, border_out=1.0) where {N}
7-
return fourier_filter_by_1D_FT!(arr, fct; border_in=border_in, border_out=border_out)
36+
function fourier_filter!(arr::AbstractArray{<:Complex, N}, fct=window_gaussian; kwargs...) where {N}
37+
return fourier_filter_by_1D_FT!(arr, fct; kwargs...)
838
end
939

1040
function fourier_filter!(arr::AbstractArray{<:Int, N}, fct=window_gaussian; kwargs...) where {N}
1141
throw(ArgumentError("FFTW.jl does not accept AbstractArrays{<:Int, N}. Convert array to a Real/Complex."))
1242
end
1343

14-
function fourier_filter!(arr::AbstractArray{<:Real, N}, fct=window_gaussian; border_in=0.0, border_out=1.0) where {N}
15-
return fourier_filter_by_1D_RFT!(arr, fct; border_in=border_in, border_out=border_out)
44+
function fourier_filter!(arr::AbstractArray{<:Real, N}, fct=window_gaussian; kwargs...) where {N}
45+
return fourier_filter_by_1D_RFT!(arr, fct; kwargs...)
46+
end
47+
48+
function fourier_filter!(arr::AbstractArray{<:Complex, N}, wins::AbstractVector; transform_win=true, kwargs...) where {N}
49+
return fourier_filter_by_1D_FT!(arr, wins; transform_win=transform_win, kwargs...)
50+
end
51+
52+
function fourier_filter!(arr::AbstractArray{<:Real, N}, wins::AbstractVector; transform_win=true, kwargs...) where {N}
53+
return fourier_filter_by_1D_RFT!(arr, fct; transform_win=transform_win, kwargs...)
1654
end
1755

1856
"""
19-
fourier_filter(arr, limits)
57+
fourier_filter(arr::AbstractArray, fct=window_gaussian; kwargs...)
58+
59+
filters an array by multiplication in Fourierspace. This version uses in-place Fourier-transforms and multiplication whereever possible.
60+
The filter window function is assumed to be separable. Depending on the array type either full-complex or real-to-complex FFTs will be used.
61+
Note that some array types cannot directly be processed due to the type inference mechanism, in which case you should transform such types to
62+
a standard array type (usually via calling "collect").
63+
High-pass filters can be used by exchanging the border_in and border_out argument values.
64+
65+
#Arguments
66+
+`arr`: the array to filter
67+
+`fct`: the separable window function to use. The function has an ND-array size as the first parameter but wil be called with only one non-singleton dimension at a time.
68+
+`kwargs...`: keyword arguments which will be passed to fct. These typically are `border_in=0.8` and `border_out=1.0` for the window-functions as defined in the `IndexFunArrays.jl` toolbox.
69+
#Example
70+
```jdoctest
71+
julia> using IndexFunArrays
72+
73+
julia> arr = zeros(10,10); arr[3,4]=1.0;
74+
75+
julia> fourier_filter(arr, window_hanning; border_in=0.5, border_out=1.0)
76+
10×10 Matrix{Float64}:
77+
-0.00528825 0.0132184 -0.0217829 -0.0994384 … -0.00528825 0.00102161 0.000100793 0.00102161
78+
0.00896055 -0.0223976 0.0369096 0.168491 0.00896055 -0.00173105 -0.000170787 -0.00173105
79+
0.0317295 -0.0793102 0.130698 0.59663 0.0317295 -0.00612969 -0.00060476 -0.00612969
80+
0.00896055 -0.0223976 0.0369096 0.168491 0.00896055 -0.00173105 -0.000170787 -0.00173105
81+
-0.00528825 0.0132184 -0.0217829 -0.0994384 -0.00528825 0.00102161 0.000100793 0.00102161
82+
0.00186562 -0.00466326 0.00768472 0.0350805 … 0.00186562 -0.000360412 -3.55585e-5 -0.000360412
83+
-1.38778e-18 1.11022e-17 5.55112e-18 -2.22045e-17 0.0 1.73472e-19 -4.29344e-18 -6.41848e-18
84+
-0.000499354 0.00124817 -0.0020569 -0.00938969 -0.000499354 9.64682e-5 9.51763e-6 9.64682e-5
85+
-5.55112e-18 0.0 1.11022e-17 8.88178e-17 -5.55112e-18 -1.21431e-17 -2.47198e-18 -2.08167e-18
86+
0.00186562 -0.00466326 0.00768472 0.0350805 0.00186562 -0.000360412 -3.55585e-5 -0.000360412```
2087
21-
Returning a fourier_filtered array.
22-
See [`fourier_filter!`](@ref fourier_filter!) for more details
2388
"""
24-
function fourier_filter(arr, fct=window_gaussian; border_in=0.0, border_out=1.0)
25-
return fourier_filter!(copy(arr), fct; border_in=border_in, border_out=border_out)
89+
function fourier_filter(arr, fct=window_gaussian; kwargs...)
90+
return fourier_filter!(copy(arr), fct; kwargs...)
2691
end
2792

28-
function fourier_filter_by_1D_FT!(arr::TA, fct=window_gaussian; border_in=0.0, border_out=1.0, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}}
93+
function fourier_filter_by_1D_FT!(arr::TA, wins::AbstractVector; transform_win=false, dims=(1:ndims(arr))) where {N, TA <: AbstractArray{<:Complex, N}}
94+
if isempty(dims)
95+
return arr
96+
end
2997
# iterates of the dimension d using the corresponding shift
30-
sz = size(arr)
3198
for d in dims #(d, limit) in pairs(limits)
32-
w = TA(ifftshift(fct(real(eltype(arr)), select_sizes(arr,d), border_in=border_in, border_out=border_out)))
33-
#w = TA(collect(ifftshift(fct(real(eltype(arr)), select_sizes(arr,d), border_in=border_in, border_out=border_out))))
34-
# go to fourier space and apply w
99+
# go to fourier space and apply window
35100
fft!(arr, d)
36-
arr .*= w
101+
arr .*= let
102+
if transform_win
103+
fft(wins[d], d)
104+
else
105+
wins[d]
106+
end
107+
end
37108
end
38-
39109
ifft!(arr, dims)
40110
return arr
41111
end
42112

43-
# the idea is the following:
44-
# rfft(x, 1) -> exp shift -> fft(x, 2) -> exp shift -> fft(x, 3) -> exp shift -> ifft(x, [2,3]) -> irfft(x, 1)
45-
# So once we did a rft to shift something we can call the routine for complex arrays to shift
46-
function fourier_filter_by_1D_RFT!(arr::TA, fct=window_gaussian; border_in=0.0, border_out=1.0) where {T<:Real, N, TA<:AbstractArray{T, N}}
47-
sz = size(arr)
48-
for d = 1:length(sz) #(d, limit) in pairs(limits)
49-
p = plan_rfft(arr, d)
50-
51-
arr_ft = p * arr
52-
w = TA(fct(real(eltype(arr)), select_sizes(arr_ft,d), offset=CtrRFFT, scale=2 ./size(arr,d), border_in=border_in, border_out=border_out))
53-
# w = TA(collect(fct(real(eltype(arr)), select_sizes(arr_ft,d), offset=CtrRFFT, scale=2 ./size(arr,d), border_in=border_in, border_out=border_out)))
54-
arr_ft .*= w
55-
# since we now did a single rfft dim, we can switch to the complex routine
56-
# new_limits = ntuple(i -> i ≤ d ? 0 : limits[i], N)
57-
# workaround to mimic in-place rfft
58-
fourier_filter_by_1D_FT!(arr_ft; border_in=border_in, border_out=border_out, dims=(2:ndims(arr_ft)))
59-
# go back to real space now and return because shift_by_1D_FT processed
60-
# the other dimensions already
61-
mul!(arr, inv(p), arr_ft)
62-
return arr # this breaks the for loop and finishes the algorithm
113+
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}}
114+
if isempty(dims)
115+
return arr
116+
end
117+
# iterates of the dimension d using the corresponding shift
118+
TR = real_arr_type(TA)
119+
wins = Vector{TR}(undef, N)
120+
# only calculate the necessary windows
121+
for d in dims
122+
# these will possibly be transformed later
123+
wins[d] = TR(ifftshift(fct(real(eltype(arr)), select_sizes(arr, d); kwargs...)))
124+
end
125+
return fourier_filter_by_1D_FT!(arr::TA, wins; transform_win=transform_win, dims=dims)
126+
end
127+
128+
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}}
129+
if isempty(dims)
130+
return arr
63131
end
132+
d = dims[1]
133+
p = plan_rfft(arr, d)
134+
135+
arr_ft = p * arr
136+
arr_ft .*= let
137+
if transform_win
138+
pw = plan_rfft(wins[d], d)
139+
pw * wins[d]
140+
else
141+
wins[d]
142+
end
143+
end
144+
# since we now did a single rfft dim, we can switch to the complex routine
145+
# new_limits = ntuple(i -> i ≤ d ? 0 : limits[i], N)
146+
# workaround to mimic in-place rfft
147+
fourier_filter_by_1D_FT!(arr_ft, wins; dims=dims[2:end], transform_win=transform_win, kwargs...)
148+
# go back to real space now and return because shift_by_1D_FT processed
149+
# the other dimensions already
150+
mul!(arr, inv(p), arr_ft)
64151
return arr
65152
end
66153

67-
function filter_gaussian(arr, sigma=eltype(arr)(1))
68-
sigma_k = 2 ./ (pi .* sigma)
69-
fourier_filter(arr, border_out=sigma_k)
154+
# transforms the first dim as rft and then hands over to the fft-based routines.
155+
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}}
156+
if isempty(dims)
157+
return arr
158+
end
159+
TR = real_arr_type(TA)
160+
d = dims[1]
161+
p = plan_rfft(arr, d)
162+
arr_ft = p * arr
163+
win = TR(fct(real(eltype(arr)), select_sizes(arr_ft,d), offset=CtrRFFT, scale=2 ./size(arr,d); kwargs...))
164+
if transform_win
165+
pw = plan_rfft(win, d)
166+
win = pw*win
167+
end
168+
arr_ft .*= win
169+
fourier_filter_by_1D_FT!(arr_ft, fct; dims=dims[2:end], transform_win=transform_win, kwargs...)
170+
# go back to real space now and return because shift_by_1D_FT processed
171+
# the other dimensions already
172+
mul!(arr, inv(p), arr_ft)
173+
return arr # this breaks the for loop and finishes the algorithm
174+
end
175+
176+
"""
177+
filter_gaussian(arr, sigma=eltype(arr)(1); border_in=(real(eltype(arr)))(0), border_out=(real(eltype(arr))).(2 ./ (pi .* sigma)), kwargs...)
178+
179+
performs Gaussian filtering by multiplying a Gaussian function in Fourier space.
180+
Note that this is not identical to the Fourier-transform of a real-space Gaussian, especially for small kernel sizes and small arrays.
181+
See also `filter_gaussian!()` and `fourier_filter()`.
182+
183+
#Arguments
184+
+`arr`: the array to filter
185+
+`sigma`: the real-space standard deviation to filter with. From this the Fourier-space standard deviation will be calculated.
186+
+`kwargs...`: additional arguments to be passed to `window_gaussian`, which is the underlying function from `IndexFunArray.jl`. This can be useful to create Fourier-shifted (Gabor-) filtering.
187+
"""
188+
function filter_gaussian(arr, sigma=eltype(arr)(1); real_space_kernel=false, border_in=(real(eltype(arr)))(0), border_out=(real(eltype(arr))).(2 ./ (pi .* sigma)), kwargs...)
189+
filter_gaussian!(copy(arr), sigma; real_space_kernel=real_space_kernel, border_in=border_in, border_out=border_out, kwargs...)
190+
end
191+
192+
"""
193+
filter_gaussian!(arr, sigma=eltype(arr)(1); border_in=(real(eltype(arr)))(0), border_out=(real(eltype(arr))).(2 ./ (pi .* sigma)), kwargs...)
194+
195+
performs in-place Gaussian filtering by multiplying a Gaussian function in Fourier space.
196+
Note that this is not identical to the Fourier-transform of a real-space Gaussian, especially for small kernel sizes and small arrays.
197+
See also `filter_gaussian()` adn `fourier_filter!()`.
198+
199+
#Arguments
200+
+`arr`: the array to filter
201+
+`sigma`: the real-space standard deviation to filter with. From this the Fourier-space standard deviation will be calculated.
202+
+`kwargs...`: additional arguments to be passed to `window_gaussian`, which is the underlying function from `IndexFunArray.jl`. This can be useful to create Fourier-shifted (Gabor-) filtering.
203+
"""
204+
function filter_gaussian!(arr, sigma=eltype(arr)(1); real_space_kernel=false, border_in=(real(eltype(arr)))(0), border_out=(real(eltype(arr))).(2 ./ (pi .* sigma)), kwargs...)
205+
if real_space_kernel
206+
mysum = sum(arr)
207+
fourier_filter!(arr, gaussian; transform_win=true, sigma=sigma, kwargs...)
208+
arr .*= (mysum/sum(arr))
209+
return arr
210+
else
211+
return fourier_filter!(arr; border_in=border_in, border_out=border_out, kwargs...)
212+
end
70213
end
71214

72215
# Suggested code by Felix. But it is not as fast as the code above:

test/fourier_filtering.jl

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
Random.seed!(42)
2+
3+
@testset "Fourier filtering" begin
4+
5+
@testset "Gaussian filter complex" begin
6+
sz = (21, 22)
7+
x = randn(ComplexF32, sz)
8+
sigma = (1.1,2.2)
9+
gf = filter_gaussian(x, sigma)
10+
# Note that this is not the same, since one kernel is generated in real space and one in Fourier space!
11+
# with sizes around 10, the difference is huge!
12+
k = gaussian(Float32, sz, sigma=sigma)
13+
k = k./sum(k) # different than "normal".
14+
gfc = conv_psf(x, k)
15+
@test (gf,gfc, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places
16+
gfr = filter_gaussian(x, sigma, real_space_kernel=true)
17+
@test (gfr,gfc) # it can be debated how to best normalize a Gaussian filter
18+
end
19+
20+
@testset "Gaussian filter real" begin
21+
sz = (21, 22)
22+
x = randn(Float32, sz)
23+
sigma = (1.1,2.2)
24+
gf = filter_gaussian(x, sigma)
25+
# Note that this is not the same, since one kernel is generated in real space and one in Fourier space!
26+
# with sizes around 10, the difference is huge!
27+
k = gaussian(sz, sigma=sigma)
28+
k = k./sum(k) # different than "normal".
29+
gf2 = conv_psf(x, k)
30+
@test (gf,gf2, rtol=1e-2) # it is realatively inaccurate due to the kernel being generated in different places
31+
end
32+
end

test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,5 +21,6 @@ include("custom_fourier_types.jl")
2121
include("damping.jl")
2222
include("czt.jl")
2323
include("nfft_tests.jl")
24+
include("fourier_filtering.jl")
2425

2526
return

0 commit comments

Comments
 (0)