adds NDCT, fixes small bug in wrapper #85

Author: mischmi96

AbstractNFFTs/src/AbstractNFFTs.jl

@@ -12,7 +12,7 @@ export AnyNFFTPlan, AbstractNFFTPlan, AbstractNFCTPlan, AbstractNNFFTPlan,
 export apodization!, apodization_adjoint!, convolve!, convolve_adjoint!
 
 # derived
-export nfft, nfft_adjoint, ndft, ndft_adjoint, nfct, nfct_transposed
+export nfft, nfft_adjoint, ndft, ndft_adjoint, nfct, nfct_transposed, ndct, ndct_transposed
 
 # misc
 export TimingStats, accuracyParams, reltolToParams, paramsToReltol, 

Wrappers/NFFT3.jl

@@ -51,9 +51,9 @@ function NFFT3Plan(x::Matrix{T}, N::NTuple{D,Int};
   f2 = UInt32(fftflags | NFFT3.FFTW_DESTROY_INPUT)
 
   n = ntuple(d -> (ceil(Int,σ*N[d])÷2)*2, D) # ensure that n is an even integer 
-  σ = n[1] / N[1]
+  #σ = n[1] / N[1]
 
-  parent = NFFT3.NFFT(Int32.(reverse(N)), size(x,2), Int32.(n), Int32(m), f1, f2)
+  parent = NFFT3.NFFT(Int32.(reverse(N)), size(x,2), Int32.(reverse(n)), Int32(m), f1, f2)
 
   xx = Float64.(reverse(x, dims=1))
   parent.x = D == 1 ? vec(xx) : xx
@@ -99,9 +99,9 @@ function NFCT3Plan(x::Matrix{T}, N::NTuple{D,Int};
   f2 = UInt32(fftflags | NFFT3.FFTW_DESTROY_INPUT)
 
   n = ntuple(d -> (ceil(Int,σ*N[d])÷2)*2, D) # ensure that n is an even integer 
-  σ = n[1] / N[1]
+  #σ = n[1] / N[1]
 
-  parent = NFFT3.NFCT(Int32.(reverse(N)), size(x,2), Int32.(n), Int32(m), f1, f2)
+  parent = NFFT3.NFCT(Int32.(reverse(N)), size(x,2), Int32.(reverse(n)), Int32(m), f1, f2)
 
   xx = Float64.(reverse(x, dims=1))
   parent.x = D == 1 ? vec(xx) : xx
@@ -135,6 +135,19 @@ function LinearAlgebra.mul!(fHat::StridedArray, p::NFFT3Plan{D}, f::AbstractArra
   return fHat
 end
 
+function LinearAlgebra.mul!(f::StridedArray, pl::Adjoint{Complex{T},<:NFFT3Plan{T}}, fHat::AbstractArray;
+                     verbose=false, timing::Union{Nothing,TimingStats} = nothing) where T
+  #consistencyCheck(p, f, fHat)
+  p = pl.parent
+
+  p.parent.f = vec(fHat)
+  p.parent.fhat = vec(f)
+  tadjoint = fApprox = NFFT3.nfft_adjoint(p.parent)
+  f[:] .= p.parent.fhat
+
+  return f
+end
+
 function LinearAlgebra.mul!(fHat::StridedArray, p::NFCT3Plan{D}, f::AbstractArray;
              verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {D}
   #consistencyCheck(p, f, fHat)

src/direct.jl

@@ -8,6 +8,15 @@ end
 AbstractNFFTs.size_in(p::NDFTPlan) = p.N
 AbstractNFFTs.size_out(p::NDFTPlan) = (p.M,)
 
+mutable struct NDCTPlan{T,D} <: AbstractNFCTPlan{T,D,1}
+  N::NTuple{D,Int64}
+  M::Int64
+  x::Matrix{T}
+end
+
+AbstractNFFTs.size_in(p::NDCTPlan) = p.N
+AbstractNFFTs.size_out(p::NDCTPlan) = (p.M,)
+
 mutable struct NNDFTPlan{T} <: AbstractNNFFTPlan{T,1,1}
   N::Int64
   M::Int64
@@ -31,6 +40,17 @@ function NDFTPlan(x::Matrix{T}, N::NTuple{D,Int}; kwargs...) where {T,D}
   return NDFTPlan{T,D}(N, M, x)
 end
 
+function NDCTPlan(x::Matrix{T}, N::NTuple{D,Int}; kwargs...) where {T,D}
+
+  if D != size(x,1)
+    throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $D"))
+  end
+
+  M = size(x, 2)
+
+  return NDCTPlan{T,D}(N, M, x)
+end
+
 function NNDFTPlan(x::Matrix{T}, y::Matrix{T}; kwargs...) where {T}
 
   M = size(x, 2)
@@ -47,6 +67,13 @@ ndft(x::AbstractArray, f::AbstractArray; kargs...) =
 ndft_adjoint(x, N, fHat::AbstractVector; kargs...) =
     adjoint(NDFTPlan(x, N; kargs...)) * fHat
 
+ndct(x::AbstractArray, f::AbstractArray; kargs...) =
+    NDCTPlan(x, size(f); kargs...) * f
+
+ndct_transposed(x, N, fHat::AbstractVector; kargs...) =
+    transpose(NDCTPlan(x, N; kargs...)) * fHat
+
+
 
 
 function LinearAlgebra.mul!(g::AbstractArray{Tg}, plan::NDFTPlan{Tp,D}, f::AbstractArray{T,D}) where {D,Tp,T,Tg}
@@ -99,7 +126,55 @@ function LinearAlgebra.mul!(g::AbstractArray{Tg,D}, pl::Adjoint{Complex{Tp},<:ND
   return g
 end
 
+function LinearAlgebra.mul!(g::AbstractArray{Tg}, plan::NDCTPlan{Tp,D}, f::AbstractArray{T,D}) where {D,Tp,T,Tg}
+
+    plan.N == size(f) ||
+        throw(DimensionMismatch("Data f is not consistent with NDFTPlan"))
+    plan.M == length(g) ||
+        throw(DimensionMismatch("Output g is inconsistent with NDFTPlan"))
+
+    g .= zero(Tg)
+
+    for l=1:prod(plan.N)
+        idx = CartesianIndices(plan.N)[l]
+
+        for k=1:plan.M
+            arg = one(T)
+            for d=1:D
+                arg *= cos( 2 * pi * plan.x[d,k] * ( idx[d] - 1 ) )
+            end
+            g[k] += f[l] * arg
+        end
+    end
+
+    return g
+end
+
+
+function LinearAlgebra.mul!(g::AbstractArray{Tg,D}, pl::Transpose{Tp,<:NDCTPlan{Tp,D}}, fHat::AbstractVector{T}) where {D,Tp,T,Tg}
+  p = pl.parent
 
+  p.M == length(fHat) ||
+      throw(DimensionMismatch("Data f inconsistent with NDFTPlan"))
+  p.N == size(g) ||
+      throw(DimensionMismatch("Output g inconsistent with NDFTPlan"))
+
+  g .= zero(Tg)
+
+  for l=1:prod(p.N)
+      idx = CartesianIndices(p.N)[l]
+
+      for k=1:p.M
+          arg = one(T)
+          for d=1:D
+              arg *= cos( 2 * pi * plan.x[d,k] * ( idx[d] - 1 ) )
+          end
+          g[l] += fHat[k] * arg
+      end
+  end
+
+  return g
+end
 
 
 function LinearAlgebra.mul!(g::AbstractArray{Tg}, plan::NNDFTPlan{Tp}, f::AbstractArray{T}) where {Tp,T,Tg}