NFCT/NFCT branch (#89)

* adds NFCT interface to AbstractNFFTs and implementation to NFFT3 wrapper #85

* adds NDCT, fixes small bug in wrapper #85

* generalize parts in AbstractNFFTs (less code)

* implement NDST

* implement NDST/NFST (NFFT3 wrapper) + tests

* apodization->deconvolve + adjoint->transpose (where appropriate)

* rename AnyNFFTPlan since the name was ambiguous

Co-authored-by: Michael Schmischke <[email protected]>

Author: tknopp

AbstractNFFTs/src/AbstractNFFTs.jl

@@ -5,14 +5,17 @@ using LinearAlgebra
 using Printf
 
 # interface
-export AnyNFFTPlan, AbstractNFFTPlan, AbstractNNFFTPlan, 
-       plan_nfft, mul!, size_in, size_out, nodes!
+export AbstractFTPlan, AbstractRealFTPlan, AbstractComplexFTPlan,
+       AbstractNFFTPlan, AbstractNFCTPlan, AbstractNFSTPlan, AbstractNNFFTPlan, 
+       plan_nfft, plan_nfct, plan_nfst, mul!, size_in, size_out, nodes!
 
 # optional
-export apodization!, apodization_adjoint!, convolve!, convolve_adjoint!
+export deconvolve!, deconvolve_transpose!, convolve!, convolve_transpose!
 
 # derived
-export nfft, nfft_adjoint, ndft, ndft_adjoint
+export nfft, nfft_adjoint, ndft, ndft_adjoint, 
+       nfct, nfct_transpose, ndct, ndct_transpose,
+       nfst, nfst_transpose
 
 # misc
 export TimingStats, accuracyParams, reltolToParams, paramsToReltol, 

AbstractNFFTs/src/derived.jl

@@ -1,79 +1,102 @@
 
 ##########################
-# plan_nfft constructors
+# plan_* constructors
 ##########################
 
-# The following automatically call the plan_nfft version for type Array
 
-plan_nfft(x::AbstractArray, N::Union{Integer,NTuple{D,Int}}, args...; kargs...) where {D} =
-    plan_nfft(Array, x, N, args...; kargs...)
+for op in [:nfft, :nfct, :nfst]
+planfunc = Symbol("plan_"*"$op")
+@eval begin 
 
-plan_nfft(x::AbstractArray, y::AbstractArray, args...; kargs...) where {D} =
-    plan_nfft(Array, x, y, args...; kargs...)
+# The following automatically call the plan_* version for type Array
 
-# The follow convert 1D parameters into the format required by the NFFT plan
+$(planfunc)(x::AbstractArray, N::Union{Integer,NTuple{D,Int}}, args...; kargs...) where {D} =
+    $(planfunc)(Array, x, N, args...; kargs...)
 
-plan_nfft(Q::Type, x::AbstractVector, N::Integer, rest...; kwargs...) where {D}  =
-    plan_nfft(Q, collect(reshape(x,1,length(x))), (N,), rest...; kwargs...)
+$(planfunc)(x::AbstractArray, y::AbstractArray, args...; kargs...) where {D} =
+    $(planfunc)(Array, x, y, args...; kargs...)
 
-plan_nfft(Q::Type, x::AbstractVector, N::NTuple{D,Int}, rest...; kwargs...) where {D} =
-    plan_nfft(Q, collect(reshape(x,1,length(x))), N, rest...; kwargs...) 
+# The follow convert 1D parameters into the format required by the plan
 
-plan_nfft(Q::Type, x::AbstractMatrix, N::NTuple{D,Int}, rest...; kwargs...) where {D}  =
-    plan_nfft(Q, collect(x), N, rest...; kwargs...)
+$(planfunc)(Q::Type, x::AbstractVector, N::Integer, rest...; kwargs...) where {D}  =
+    $(planfunc)(Q, collect(reshape(x,1,length(x))), (N,), rest...; kwargs...)
 
-plan_nfft(Q::Type, x::AbstractVector, y::AbstractVector, rest...; kwargs...) where {D}  =
-    plan_nfft(Q, collect(reshape(x,1,length(x))), collect(reshape(y,1,length(x))), rest...; kwargs...)
+$(planfunc)(Q::Type, x::AbstractVector, N::NTuple{D,Int}, rest...; kwargs...) where {D} =
+    $(planfunc)(Q, collect(reshape(x,1,length(x))), N, rest...; kwargs...) 
 
+$(planfunc)(Q::Type, x::AbstractMatrix, N::NTuple{D,Int}, rest...; kwargs...) where {D}  =
+    $(planfunc)(Q, collect(x), N, rest...; kwargs...)
+
+end
+end
+
+## NNFFT constructor
+plan_nnfft(Q::Type, x::AbstractVector, y::AbstractVector, rest...; kwargs...) where {D}  =
+    plan_nnfft(Q, collect(reshape(x,1,length(x))), collect(reshape(y,1,length(x))), rest...; kwargs...)
 
-##########################
-# Allocating nfft functions
-##########################
 
+
+###############################################
+# Allocating trafo functions with plan creation
+###############################################
+
+for (op,trans) in zip([:nfft, :nfct, :nfst],
+                      [:adjoint, :transpose, :transpose])
+planfunc = Symbol("plan_$(op)")
+tfunc = Symbol("$(op)_$(trans)")
+@eval begin 
+
+# TODO fix comments (how?)
 """
-nfft(x, f::AbstractArray{T,D}, rest...; kwargs...)
+nfft(x, f, rest...; kwargs...)
 
-calculates the NFFT of the array `f` for the nodes contained in the matrix `x`
+calculates the nfft of the array `f` for the nodes contained in the matrix `x`
 The output is a vector of length M=`size(nodes,2)`
 """
-function nfft(x, f::AbstractArray{T,D}; kargs...) where {T,D}
-  p = plan_nfft(x, size(f); kargs... )
+function $(op)(x, f::AbstractArray; kargs...) 
+  p = $(planfunc)(x, size(f); kargs... )
   return p * f
 end
 
 """
-nfft_adjoint(x, N, fHat::AbstractArray{T,D}, rest...; kwargs...)
+nfft_adjoint(x, N, fHat, rest...; kwargs...)
 
-calculates the adjoint NFFT of the vector `fHat` for the nodes contained in the matrix `x`.
+calculates the adjoint nfft of the vector `fHat` for the nodes contained in the matrix `x`.
 The output is an array of size `N`
 """
-function nfft_adjoint(x, N, fHat::AbstractVector{T};  kargs...) where T
-  p = plan_nfft(x, N;  kargs...)
-  return adjoint(p) * fHat
+function $(tfunc)(x, N, fHat;  kargs...) 
+  p = $(planfunc)(x, N;  kargs...)
+  return $(trans)(p) * fHat
+end
+
+end
 end
 
+############################
+# Allocating trafo functions
+############################
 
 """
         *(p, f) -> fHat
 
-For a **non**-directional `D` dimensional plan `p` this calculates the NFFT of a `D` dimensional array `f` of size `N`.
+For a **non**-directional `D` dimensional plan `p` this calculates the NFFT/NNFFT of a `D` dimensional array `f` of size `N`.
 `fHat` is a vector of length `M`.
 (`M` and `N` are defined in the plan creation)
 
 For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
 dimensional arrays, and the dimension specified in the plan creation is
 affected.
 """
-function Base.:*(p::AnyNFFTPlan{T}, f::AbstractArray{Complex{U},D}; kargs...) where {T,U,D}
+function Base.:*(p::AbstractComplexFTPlan{T}, f::AbstractArray{Complex{U},D}; kargs...) where {T,U,D}
   fHat = similar(f, Complex{T}, size_out(p))
   mul!(fHat, p, f; kargs...)
   return fHat
 end
 
 """
-        *(p::Adjoint{T,<:AnyNFFTPlan{T}}, fHat) -> f
+        *(p::Adjoint{T,<:AbstractFTPlan{T}}, fHat) -> f
 
-For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFFT of a length `M` vector `fHat`
+For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFFT/NNFFT of a length `M` vector `fHat`
 `f` is a `D` dimensional array of size `N`.
 (`M` and `N` are defined in the plan creation)
 
@@ -82,22 +105,69 @@ dimensional arrays, and the dimension specified in the plan creation is
 affected.
 """
 
-function Base.:*(p::Adjoint{Complex{T},<:AnyNFFTPlan{T}}, fHat::AbstractArray{Complex{U},D}; kargs...) where {T,U,D}
+function Base.:*(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T}}, fHat::AbstractArray{Complex{U},D}; kargs...) where {T,U,D}
   f = similar(fHat, Complex{T}, size_out(p))
   mul!(f, p, fHat; kargs...)
   return f
 end
 
 # The following two methods are redundant but need to be defined because of a method ambiguity with Julia Base
-function Base.:*(p::Adjoint{Complex{T},<:AnyNFFTPlan{T}}, fHat::AbstractVector{Complex{U}}; kargs...) where {T,U}
+function Base.:*(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T}}, fHat::AbstractVector{Complex{U}}; kargs...) where {T,U}
   f = similar(fHat, Complex{T}, size_out(p))
   mul!(f, p, fHat; kargs...)
   return f
 end
-function Base.:*(p::Adjoint{Complex{T},<:AnyNFFTPlan{T}}, fHat::AbstractArray{Complex{U},2}; kargs...) where {T,U}
+function Base.:*(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T}}, fHat::AbstractArray{Complex{U},2}; kargs...) where {T,U}
   f = similar(fHat, Complex{T}, size_out(p))
   mul!(f, p, fHat; kargs...)
   return f
 end
 
 
+
+"""
+        *(p, f) -> fHat
+
+For a **non**-directional `D` dimensional plan `p` this calculates the NFCT/NFST of a `D` dimensional array `f` of size `N`.
+`fHat` is a vector of length `M`.
+(`M` and `N` are defined in the plan creation)
+
+For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
+dimensional arrays, and the dimension specified in the plan creation is
+affected.
+"""
+function Base.:*(p::AbstractRealFTPlan{T}, f::AbstractArray{U,D}; kargs...) where {T,U,D}
+  fHat = similar(f, T, size_out(p))
+  mul!(fHat, p, f; kargs...)
+  return fHat
+end
+
+"""
+        *(p::Transpose{T,AbstractRealFTPlan{T}}, fHat) -> f
+
+For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFCT/NFST of a length `M` vector `fHat`
+`f` is a `D` dimensional array of size `N`.
+(`M` and `N` are defined in the plan creation)
+
+For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
+dimensional arrays, and the dimension specified in the plan creation is
+affected.
+"""
+
+function Base.:*(p::Transpose{T,<:AbstractRealFTPlan{T}}, fHat::AbstractArray{U,D}; kargs...) where {T,U,D}
+  f = similar(fHat, T, size_out(p))
+  mul!(f, p, fHat; kargs...)
+  return f
+end
+
+# The following two methods are redundant but need to be defined because of a method ambiguity with Julia Base
+function Base.:*(p::Transpose{T,<:AbstractRealFTPlan{T}}, fHat::AbstractVector{U}; kargs...) where {T,U}
+  f = similar(fHat, T, size_out(p))
+  mul!(f, p, fHat; kargs...)
+  return f
+end
+function Base.:*(p::Transpose{T,<:AbstractRealFTPlan{T}}, fHat::AbstractArray{U,2}; kargs...) where {T,U}
+  f = similar(fHat, T, size_out(p))
+  mul!(f, p, fHat; kargs...)
+  return f
+end