Merge branch 'main' into fractiona_fourier_transform

Author: roflmaostc

Project.toml

@@ -1,7 +1,7 @@
 name = "FourierTools"
 uuid = "b18b359b-aebc-45ac-a139-9c0ccbb2871e"
 authors = ["Felix Wechsler (roflmaostc) <[email protected]>", "rheintzmann <[email protected]>"]
-version = "0.3.2"
+version = "0.3.7"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -11,19 +11,22 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NDTools = "98581153-e998-4eef-8d0d-5ec2c052313d"
 NFFT = "efe261a4-0d2b-5849-be55-fc731d526b0d"
 PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 ShiftedArrays = "1277b4bf-5013-50f5-be3d-901d8477a67a"
 
 [compat]
-ChainRulesCore = "0.9, 0.10, 1.0, 1.1, 1"
+
+ChainRulesCore = "1, 1.0, 1.1"
 FFTW = "1.5"
-ImageTransformations = "0.8.12"
+ImageTransformations = "0.9"
 IndexFunArrays = "0.2"
-NDTools = "0.4.2"
-NFFT = "0.11, 0.12"
-PaddedViews = "0.5.8"
-ShiftedArrays = "1"
+NDTools = "0.5.1"
+NFFT = "0.11, 0.12, 0.13"
+PaddedViews = "0.5"
+Reexport = "1"
+ShiftedArrays = "2"
 Zygote = "0.6"
-julia = "1.6, 1.7"
+julia = "1, 1.6, 1.7, 1.8"
 
 [extras]
 FractionalTransforms = "e50ca838-b4f0-4a10-ad18-4b920bf1ae5c"

README.md

@@ -27,16 +27,25 @@ The main features are:
 * array/image rotation 
 * array/image shifting (including noteworthy subpixel shifts)
 * array/image shearing
+* convenient wrappers of [NFFT.jl](https://github.com/JuliaMath/NFFT.jl)
 * several tools like `ffts`, `ft`, `fftshift_view` etc. allowing simpler use with Fourier transforms
 * Chirp Z-Transform
 * Fractional Fourier Transform
+* reexports [FFTW.jl](https://github.com/JuliaMath/FFTW.jl)
 
 Have a look in the [examples folder](examples/) for interactive examples. The [documentation](https://bionanoimaging.github.io/FourierTools.jl/dev/) offers a quick overview.
 
 ## FFTW Threading
 By default we set 4 Threads. Use `FFTW.set_num_threads(N)` to set `N` threads.
 
 
+
+## Related Packages
+There are numerous packages related to Fourier transforms which offer similar functions or which this package is based on:
+* [FFTW.jl](https://github.com/JuliaMath/FFTW.jl) for FFTs
+* [NFFT.jl](https://github.com/JuliaMath/NFFT.jl) for non-uniform FFTs
+* [FractionalTransforms.jl](https://github.com/SciFracX/FractionalTransforms.jl) offers 1D fractional Fourier transforms
+
 ## Cite
 If you use this package in an academic work, please cite us!
 See on the right side the *Cite this repository*:
@@ -47,6 +56,9 @@ title = {FourierTools.jl - Efficiently Working with Fourier Space},
 url = {https://github.com/bionanoimaging/FourierTools.jl}}
 ```
 
+## Development
+Feel free to file an issue regarding problems, suggestions or improvement ideas for this package!
+
 
 [docs-dev-img]: https://img.shields.io/badge/docs-dev-pink.svg
 [docs-dev-url]: https://bionanoimaging.github.io/FourierTools.jl/dev/

docs/make.jl

@@ -7,11 +7,12 @@ makedocs(modules = [FourierTools],
          pages = Any[
             "FourierTools.jl" => "index.md",
             "FFT Helpers" => "helpers.md",
-            "FFT Based Convolutions" => "convolutions.md",
+            "FFT Based Convolutions and Cross-Correlation" => "convolutions.md",
             "Resampling (sinc Interpolation)" => "resampling.md",
             "Shifting with FFTs" => "shifting.md",
             "Image Shearing with FFTs" => "shear.md",
             "Image Rotation with FFTs" => "rotate.md",
+            "NFFT" => "nfft.md",
             "Utility Functions" => "utils.md",
          ]
         )

docs/src/convolutions.md

@@ -1,6 +1,7 @@
 ## Function References
 ```@docs
 conv
+ccorr
 conv_psf
 plan_conv
 plan_conv_psf

docs/src/helpers.md

@@ -5,8 +5,6 @@ ffts!
 iffts
 rffts
 irffts
-fftshift!
-ifftshift!
 FourierTools.fftshift_view
 FourierTools.ifftshift_view
 FourierTools.rfftshift_view

docs/src/nfft.md

@@ -0,0 +1,7 @@
+# NFFT helpers
+Based on [NFFT.jl](https://github.com/JuliaMath/NFFT.jl) we provide some convenient helper functions:
+
+```@docs
+FourierTools.nfft_nd
+FourierTools.plan_nfft_nd
+```

docs/src/utils.md

@@ -6,8 +6,4 @@ FourierTools.select_region
 FourierTools.center_set!
 FourierTools.get_indices_around_center
 FourierTools.center_extract
-FourierTools.selectsizes
-FourierTools.expanddims
-FourierTools.slice
-FourierTools.slice_indices
 ```

src/FourierTools.jl

@@ -1,26 +1,32 @@
 module FourierTools
 
+
+using Reexport
 using PaddedViews, ShiftedArrays
-using FFTW
-FFTW.set_num_threads(4)
+@reexport using FFTW
 using LinearAlgebra
 using IndexFunArrays
 using ChainRulesCore
 using NDTools
-using NFFT
+@reexport using NFFT
+FFTW.set_num_threads(4)
+
+
 
 include("utils.jl")
 include("nfft_nd.jl")
 include("resampling.jl")
 include("custom_fourier_types.jl")
 include("fourier_resizing.jl")
 include("fourier_shifting.jl")
+include("fourier_filtering.jl")
 include("fourier_resample_1D_based.jl")
 include("fourier_rotate.jl")
 include("fourier_shear.jl")
 include("fftshift_alternatives.jl")
 include("fft_helpers.jl")
 include("convolutions.jl")
+include("correlations.jl")
 include("damping.jl")
 include("czt.jl")
 include("fractional_fourier_transform.jl")

src/convolutions.jl

@@ -136,18 +136,23 @@ function plan_conv(u::AbstractArray{T1, N}, v::AbstractArray{T2, M}, dims=ntuple
         # FFTW.MEASURE flag might overwrite input! Hence copy!
         if (:flags in keys(kwargs) && 
             (getindex(kwargs, :flags) == FFTW.MEASURE || getindex(kwargs, :flags) == FFTW.PATIENT)) 
-            P = plan(copy(u), dims; kwargs...)
+            plan(copy(u), dims; kwargs...)
         else
-            P = plan(u, dims; kwargs...)
+            plan(u, dims; kwargs...)
         end
     end
-    P_inv = inv(P)
 
     v_ft = fft_or_rfft(T1)(v, dims)
     # construct the efficient conv function
     # P and P_inv can be understood like matrices
     # but their computation is fast
-    conv(u, v_ft=v_ft) = p_conv_aux(P, P_inv, u, v_ft)
+    conv = let P = P,
+               P_inv = inv(P),
+               # put a different name here! See https://discourse.julialang.org/t/type-issue-with-captured-variables-let-workaround-failed/85661
+               v_ft = v_ft
+        conv(u, v_ft=v_ft) = p_conv_aux(P, P_inv, u, v_ft)
+    end
+    
     return v_ft, conv
 end
 
@@ -211,9 +216,7 @@ function plan_conv_psf(u::AbstractArray{T, N}, psf::AbstractArray{T, M}, dims=nt
 end
 
 function p_conv_aux(P, P_inv, u, v_ft)
-    tmp = (P_inv.p * ((P * u) .* v_ft))
-    tmp .*= P_inv.scale
-    return tmp
+    return (P_inv.p * ((P * u) .* v_ft .* P_inv.scale))
 end
 
 function ChainRulesCore.rrule(::typeof(p_conv_aux), P, P_inv, u, v)

src/correlations.jl

@@ -0,0 +1,69 @@
+export ccorr
+
+
+"""
+    ccorr(u, v[, dims]; centered=false)
+
+Calculates the cross-correlation between `u` and `v` along `dims`.
+`centered=true` moves the output of the cross-correlation to the Fourier center.
+
+If `u` and `v` are both a real valued array we use `rfft` and hence
+the output is real as well.
+If either `u` or `v` is complex we use `fft` and output is hence complex.
+
+Per default the correlation is performed along `min(ndims(u), ndims(v))`.
+
+
+```jldoctest
+julia> ccorr([1,1,0,0], [1,1,0,0], centered=true)
+4-element Vector{Float64}:
+ 0.0
+ 1.0
+ 2.0
+ 1.0
+
+julia> ccorr([1,1,0,0], [1,1,0,0])
+4-element Vector{Float64}:
+ 2.0
+ 1.0
+ 0.0
+ 1.0
+
+julia> ccorr([1im,0,0,0], [0,1im,0,0])
+4-element Vector{ComplexF64}:
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 1.0 + 0.0im
+
+julia> ccorr([1im,0,0,0], [0,1im,0,0], centered=true)
+4-element Vector{ComplexF64}:
+ 0.0 + 0.0im
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+```
+"""
+function ccorr(u::AbstractArray{T, N}, v::AbstractArray{D, M}, 
+              dims=ntuple(+, min(N, M));
+              centered=false) where {T, D, N, M}
+    out = ifft(fft(u, dims) .* conj.(fft(v, dims)), dims)
+    
+    if centered
+        return fftshift(out)
+    else
+        return out
+    end
+end
+
+function ccorr(u::AbstractArray{<:Real, N}, v::AbstractArray{<:Real, M}, 
+              dims=ntuple(+, min(N, M));
+              centered=false) where {N, M}
+    out = irfft(rfft(u, dims) .* conj.(rfft(v, dims)), size(u, dims[1]), dims)
+    
+    if centered
+        return fftshift(out)
+    else
+        return out
+    end
+end