renamed to normalize_win, bug fixes and docstring, more tests

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,15 +249,15 @@ 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
     end
     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...)

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