Fix gaussian integral

src/fourier_filtering.jl

@@ -90,7 +90,7 @@
     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,17 +100,29 @@
         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
     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}}
+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 @@
         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 @@
 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,8 +187,15 @@
             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...)
+
+    # 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)
@@ -310,9 +329,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, keep_integral=true, sigma=sigma, kwargs...)
         return arr
     else
         return fourier_filter!(arr; border_in=border_in, border_out=border_out, kwargs...)

test/fourier_filtering.jl

@@ -14,20 +14,24 @@ 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]