adds NFCT interface to AbstractNFFTs and implementation to NFFT3 wrap…

…per #85

AbstractNFFTs/src/AbstractNFFTs.jl

@@ -5,14 +5,14 @@ using LinearAlgebra
 using Printf
 
 # interface
-export AnyNFFTPlan, AbstractNFFTPlan, AbstractNNFFTPlan, 
-       plan_nfft, mul!, size_in, size_out, nodes!
+export AnyNFFTPlan, AbstractNFFTPlan, AbstractNFCTPlan, AbstractNNFFTPlan, 
+       plan_nfft, plan_nfct, mul!, size_in, size_out, nodes!
 
 # optional
 export apodization!, apodization_adjoint!, convolve!, convolve_adjoint!
 
 # derived
-export nfft, nfft_adjoint, ndft, ndft_adjoint
+export nfft, nfft_adjoint, ndft, ndft_adjoint, nfct, nfct_transposed
 
 # misc
 export TimingStats, accuracyParams, reltolToParams, paramsToReltol, 

AbstractNFFTs/src/derived.jl

@@ -24,6 +24,32 @@ plan_nfft(Q::Type, x::AbstractMatrix, N::NTuple{D,Int}, rest...; kwargs...) wher
 
 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...)
+
+##########################
+# plan_nfct constructors
+##########################
+
+# The following automatically call the plan_nfft version for type Array
+
+plan_nfct(x::AbstractArray, N::Union{Integer,NTuple{D,Int}}, args...; kargs...) where {D} =
+    plan_nfct(Array, x, N, args...; kargs...)
+
+plan_nfct(x::AbstractArray, y::AbstractArray, args...; kargs...) where {D} =
+    plan_nfct(Array, x, y, args...; kargs...)
+
+# The follow convert 1D parameters into the format required by the NFFT plan
+
+plan_nfct(Q::Type, x::AbstractVector, N::Integer, rest...; kwargs...) where {D}  =
+    plan_nfct(Q, collect(reshape(x,1,length(x))), (N,), rest...; kwargs...)
+
+plan_nfct(Q::Type, x::AbstractVector, N::NTuple{D,Int}, rest...; kwargs...) where {D} =
+    plan_nfct(Q, collect(reshape(x,1,length(x))), N, rest...; kwargs...) 
+
+plan_nfct(Q::Type, x::AbstractMatrix, N::NTuple{D,Int}, rest...; kwargs...) where {D}  =
+    plan_nfct(Q, collect(x), N, rest...; kwargs...)
+
+plan_nfct(Q::Type, x::AbstractVector, y::AbstractVector, rest...; kwargs...) where {D}  =
+    plan_nfct(Q, collect(reshape(x,1,length(x))), collect(reshape(y,1,length(x))), rest...; kwargs...)
 
 
 ##########################
@@ -64,7 +90,7 @@ 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::AbstractNFFTPlan{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
@@ -82,22 +108,94 @@ 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},<:AbstractNFFTPlan{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},<:AbstractNFFTPlan{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},<:AbstractNFFTPlan{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
 
+##########################
+# Allocating nfct functions
+##########################
+
+"""
+nfct(x, f::AbstractArray{T,D}, rest...; kwargs...)
+
+calculates the NFCT of the array `f` for the nodes contained in the matrix `x`
+The output is a vector of length M=`size(nodes,2)`
+"""
+function nfct(x, f::AbstractArray{T,D}; kargs...) where {T,D}
+  p = plan_nfct(x, size(f); kargs... )
+  return p * f
+end
+
+"""
+nfct_transposed(x, N, fHat::AbstractArray{T,D}, rest...; kwargs...)
+
+calculates the transposed NFCT of the vector `fHat` for the nodes contained in the matrix `x`.
+The output is an array of size `N`
+"""
+function nfct_transposed(x, N, fHat::AbstractVector{T};  kargs...) where T
+  p = plan_nfct(x, N;  kargs...)
+  return transpose(p) * fHat
+end
+
 
+"""
+        *(p, f) -> fHat
+
+For a **non**-directional `D` dimensional plan `p` this calculates the NFCT 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::AbstractNFCTPlan{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,AbstractNFCTPlan{T}}, fHat) -> f
+
+For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFCT 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,<:AbstractNFCTPlan{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,<:AbstractNFCTPlan{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,<:AbstractNFCTPlan{T}}, fHat::AbstractArray{U,2}; kargs...) where {T,U}
+  f = similar(fHat, T, size_out(p))
+  mul!(f, p, fHat; kargs...)
+  return f
+end

AbstractNFFTs/src/interface.jl

@@ -18,6 +18,16 @@ Abstract type for an NFFT plan.
 """
 abstract type AbstractNFFTPlan{T,D,R} <: AnyNFFTPlan{T,D,R} end
 
+"""
+  AbstractNFCTPlan{T,D,R}
+
+Abstract type for an NFCT plan.
+* T is the element type (Float32/Float64)
+* D is the number of dimensions of the input array.
+* R is the number of dimensions of the output array.
+"""
+abstract type AbstractNFCTPlan{T,D,R} <: AnyNFFTPlan{T,D,R} end
+
 """
   AbstractNNFFTPlan{T,D,R}
 
@@ -32,16 +42,38 @@ abstract type AbstractNNFFTPlan{T,D,R} <: AnyNFFTPlan{T,D,R} end
 # Function needed to make AnyNFFTPlan an operator 
 #####################
 
-Base.eltype(p::AnyNFFTPlan{T}) where T = Complex{T}
+Base.eltype(p::AbstractNFFTPlan{T}) where T = Complex{T}
+Base.eltype(p::AbstractNNFFTPlan{T}) where T = Complex{T}
+Base.eltype(p::AbstractNFCTPlan{T}) where T = T
+
 Base.size(p::AnyNFFTPlan) = (prod(size_out(p)), prod(size_in(p)))
-Base.size(p::Adjoint{Complex{T}, U}) where {T, U<:AnyNFFTPlan{T}} = 
+
+Base.size(p::Adjoint{Complex{T}, U}) where {T, U<:AbstractNFFTPlan{T}} = 
+  (prod(size_in(p.parent)), prod(size_out(p.parent)))
+
+Base.size(p::Adjoint{Complex{T}, U}) where {T, U<:AbstractNNFFTPlan{T}} = 
   (prod(size_in(p.parent)), prod(size_out(p.parent)))
 
-LinearAlgebra.adjoint(p::AnyNFFTPlan{T,D,R}) where {T,D,R} = 
+Base.size(p::Transpose{T, U}) where {T, U<:AbstractNFCTPlan{T}} = 
+  (prod(size_in(p.parent)), prod(size_out(p.parent)))
+
+LinearAlgebra.adjoint(p::AbstractNFFTPlan{T,D,R}) where {T,D,R} = 
 Adjoint{Complex{T}, typeof(p)}(p)
 
-size_in(p::Adjoint{Complex{T},<:AnyNFFTPlan{T}}) where {T} = size_out(p.parent)
-size_out(p::Adjoint{Complex{T},<:AnyNFFTPlan{T}}) where {T} = size_in(p.parent)
+LinearAlgebra.adjoint(p::AbstractNNFFTPlan{T,D,R}) where {T,D,R} = 
+Adjoint{Complex{T}, typeof(p)}(p)
+
+LinearAlgebra.transpose(p::AbstractNFCTPlan{T,D,R}) where {T,D,R} = 
+Transpose{T, typeof(p)}(p)
+
+size_in(p::Adjoint{Complex{T},<:AbstractNFFTPlan{T,D,R}}) where {T,D,R} = size_out(p.parent)
+size_out(p::Adjoint{Complex{T},<:AbstractNFFTPlan{T,D,R}}) where {T,D,R} = size_in(p.parent)
+
+size_in(p::Adjoint{Complex{T},<:AbstractNNFFTPlan{T,D,R}}) where {T,D,R} = size_out(p.parent)
+size_out(p::Adjoint{Complex{T},<:AbstractNNFFTPlan{T,D,R}}) where {T,D,R} = size_in(p.parent)
+
+size_in(p::Transpose{T,<:AbstractNFCTPlan{T,D,R}}) where {T,D,R} = size_out(p.parent)
+size_out(p::Transpose{T,<:AbstractNFCTPlan{T,D,R}}) where {T,D,R} = size_in(p.parent)
 
 #####################
 # define interface
@@ -58,14 +90,32 @@ Both `f` and `fHat` must be complex arrays of element type `Complex{T}`.
 """
 @mustimplement LinearAlgebra.mul!( fHat::AbstractArray, p::AbstractNFFTPlan{T}, f::AbstractArray) where {T}
 
+"""
+    mul!(fHat, p, f) -> fHat
+
+Inplace NFCT transforming the `D` dimensional array `f` to the `R` dimensional array `fHat`.
+The transformation is applied along `D-R+1` dimensions specified in the plan `p`.
+Both `f` and `fHat` must be complex arrays of element type `Complex{T}`.
+"""
+@mustimplement LinearAlgebra.mul!( fHat::AbstractArray, p::AbstractNFCTPlan{T}, f::AbstractArray) where {T}
+
 """
     mul!(f, p, fHat) -> f
 
 Inplace adjoint NFFT transforming the `R` dimensional array `fHat` to the `D` dimensional array `f`.
 The transformation is applied along `D-R+1` dimensions specified in the plan `p`.
 Both `f` and `fHat` must be complex arrays of element type `Complex{T}`.
 """
-@mustimplement LinearAlgebra.mul!(f::AbstractArray, p::Adjoint{Complex{T},<:AnyNFFTPlan{T}}, fHat::AbstractArray) where {T}
+@mustimplement LinearAlgebra.mul!(f::AbstractArray, p::Adjoint{Complex{T},<:AbstractNFFTPlan{T}}, fHat::AbstractArray) where {T}
+
+"""
+    mul!(f, p, fHat) -> f
+
+Inplace transposed NFCT transforming the `R` dimensional array `fHat` to the `D` dimensional array `f`.
+The transformation is applied along `D-R+1` dimensions specified in the plan `p`.
+Both `f` and `fHat` must be arrays of element type `T`.
+"""
+@mustimplement LinearAlgebra.mul!(f::AbstractArray, p::Transpose{T,<:AbstractNFCTPlan{T}}, fHat::AbstractArray) where {T}
 
 """
     size_in(p)
@@ -75,6 +125,14 @@ Note that this will be the output array for an adjoint NFFT.
 """
 @mustimplement size_in(p::AbstractNFFTPlan{T,D,R}) where {T,D,R}
 
+"""
+    size_in(p)
+
+Size of the input array for an NFCT operation. The returned tuple has `D` entries. 
+Note that this will be the output array for a transposed NFCT.
+"""
+@mustimplement size_in(p::AbstractNFCTPlan{T,D,R}) where {T,D,R}
+
 """
     size_out(p)
 
@@ -83,12 +141,21 @@ Note that this will be the input array for an adjoint NFFT.
 """
 @mustimplement size_out(p::AbstractNFFTPlan{T,D,R}) where {T,D,R}
 
+"""
+    size_out(p)
+
+Size of the output array for an NFCT operation. The returned tuple has `R` entries. 
+Note that this will be the input array for a transposed NFCT.
+"""
+@mustimplement size_out(p::AbstractNFCTPlan{T,D,R}) where {T,D,R}
+
 """
     nodes!(p, x) -> p
 
 Change nodes `x` in the plan `p` operation and return the plan.
 """
 @mustimplement nodes!(p::AbstractNFFTPlan{T}, x::Matrix{T}) where {T}
+@mustimplement nodes!(p::AbstractNFCTPlan{T}, x::Matrix{T}) where {T}
 
 
 ## Optional Interface ##