Exchangable NFFT backends

Project.toml

@@ -24,14 +24,15 @@ julia = "1.9"
 GPUArrays = "11"
 KernelAbstractions = "0.9"
 JLArrays = "0.2"
-NFFT = "0.13"
+AbstractNFFTs = "0.9"
 LinearOperators = "2"
 RadonKA = "0.6"
 Wavelets = "0.9, 0.10"
 Reexport = "1.0"
 FFTW = "1.0"
 
 [weakdeps]
+AbstractNFFTs = "7f219486-4aa7-41d6-80a7-e08ef20ceed7"
 GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 NFFT = "efe261a4-0d2b-5849-be55-fc731d526b0d"
@@ -43,7 +44,7 @@ RadonKA = "86de8297-835b-47df-b249-c04e8db91db5"
 test = ["Test", "FFTW", "Wavelets", "NFFT", "JLArrays", "RadonKA"]
 
 [extensions]
-LinearOperatorNFFTExt = ["NFFT", "FFTW"]
+LinearOperatorNFFTExt = ["AbstractNFFTs", "FFTW"]
 LinearOperatorFFTWExt = "FFTW"
 LinearOperatorWaveletExt = "Wavelets"
 LinearOperatorGPUArraysExt = "GPUArrays"

ext/LinearOperatorNFFTExt/LinearOperatorNFFTExt.jl

@@ -1,6 +1,6 @@
 module LinearOperatorNFFTExt
 
-using LinearOperatorCollection, NFFT, NFFT.AbstractNFFTs, FFTW, FFTW.AbstractFFTs
+using LinearOperatorCollection, AbstractNFFTs, FFTW, FFTW.AbstractFFTs
 
 include("NFFTOp.jl")
 

ext/LinearOperatorNFFTExt/NFFTOp.jl

@@ -6,13 +6,16 @@ generates a `NFFTOpImpl` which evaluates the MRI Fourier signal encoding operato
 
 # Arguments:
 * `shape::NTuple{D,Int64}`  - size of image to encode/reconstruct
-* `tr`                      - Either a `Trajectory` object, or a `ND x Nsamples` matrix for an ND-dimenensional (e.g. 2D or 3D) NFFT with `Nsamples` k-space samples
-* (`nodes=nothing`)         - Array containg the trajectory nodes (redundant)
-* (`kargs`)                 - additional keyword arguments
+* `nodes=nothing`         - Array containg the trajectory nodes
+* `toeplitz=false`        - 
+* `oversamplingFactor=1.25`
+* `kernelSize=3`
+* `precompute = AbstractNFFTs.TENSOR` Precompute flag for the NFFT backend
+* (`kargs`)                 - additional keyword arguments for the NFFT plan, 
 """
 function LinearOperatorCollection.NFFTOp(::Type{T};
     shape::Tuple, nodes::AbstractMatrix{U}, toeplitz=false, oversamplingFactor=1.25, 
-   kernelSize=3, kargs...) where {U <: Number, T <: Number}
+   kernelSize=3, precompute = AbstractNFFTs.TENSOR, kargs...) where {U <: Number, T <: Number}
   return NFFTOpImpl(shape, nodes; toeplitz, oversamplingFactor, kernelSize, kargs... )
 end
 
@@ -38,11 +41,11 @@ end
 
 LinearOperators.storage_type(op::NFFTOpImpl) = typeof(op.Mv5)
 
-function NFFTOpImpl(shape::Tuple, tr::AbstractMatrix{T}; toeplitz=false, oversamplingFactor=1.25, kernelSize=3, S = Vector{Complex{T}}, kargs...) where {T}
+function NFFTOpImpl(shape::Tuple, tr::AbstractMatrix{T}; toeplitz, oversamplingFactor, kernelSize, S = Vector{Complex{T}}, kargs...) where {T}
 
   baseArrayType = Base.typename(S).wrapper # https://github.com/JuliaLang/julia/issues/35543
-  plan = plan_nfft(baseArrayType, tr, shape, m=kernelSize, σ=oversamplingFactor, precompute=NFFT.TENSOR,
-		                          fftflags=FFTW.ESTIMATE, blocking=true)
+  plan = plan_nfft(baseArrayType, tr, shape; m=kernelSize, σ=oversamplingFactor,
+		                          fftflags=FFTW.ESTIMATE, blocking=true, kargs...)
 
   return NFFTOpImpl{eltype(S), S, typeof(plan)}(size(tr,2), prod(shape), false, false
             , (res,x) -> produ!(res,plan,x)
@@ -143,7 +146,7 @@ function NFFTToeplitzNormalOp(nfft::NFFTOp{T}, W=nothing; kwargs...) where {T}
   shape_os = 2 .* shape
   baseArrayType = Base.typename(typeof(tmpVec)).wrapper # https://github.com/JuliaLang/julia/issues/35543
   p = plan_nfft(baseArrayType, nfft.plan.k, shape_os; m = nfft.plan.params.m, σ = nfft.plan.params.σ,
-		precompute=NFFT.POLYNOMIAL, fftflags=FFTW.ESTIMATE, blocking=true)
+		precompute=AbstractNFFTs.POLYNOMIAL, fftflags=FFTW.ESTIMATE, blocking=true)
   tmpOnes = similar(tmpVec, size(nfft.plan.k, 2))
   tmpOnes .= one(T)
 

test/testOperators.jl

@@ -218,7 +218,7 @@ function testNFFT2d(N=16;arrayType = Array)
   # Operator
   xop = arrayType(vec(x))
   nodes = [(idx[d] - N÷2 - 1)./N for d=1:2, idx in vec(CartesianIndices((N,N)))]
-  F_nfft = NFFTOp(ComplexF64; shape=(N,N), nodes, symmetrize=false, S = typeof(xop))
+  F_nfft = NFFTOp(ComplexF64; shape=(N,N), nodes, S = typeof(xop))
 
   # test against FourierOperators
   y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N)) )
@@ -249,7 +249,7 @@ function testNFFT2d(N=16;arrayType = Array)
   # test type stability;
   # TODO: Ensure type stability for Trajectory objects and test here
   nodes = Float32.(nodes)
-  F_nfft = NFFTOp(ComplexF32; shape=(N,N), nodes, symmetrize=false, S = typeof(ComplexF32.(xop)))
+  F_nfft = NFFTOp(ComplexF32; shape=(N,N), nodes, S = typeof(ComplexF32.(xop)))
 
   y_nfft = F_nfft * ComplexF32.(xop)
   y_adj_nfft = adjoint(F_nfft) * ComplexF32.(xop)
@@ -273,7 +273,7 @@ function testNFFT3d(N=12;arrayType = Array)
   # Operator
   xop = arrayType(vec(x))
   nodes = [(idx[d] - N÷2 - 1)./N for d=1:3, idx in vec(CartesianIndices((N,N,N)))]
-  F_nfft = NFFTOp(ComplexF64; shape=(N,N,N), nodes=nodes, symmetrize=false, S = typeof(xop))
+  F_nfft = NFFTOp(ComplexF64; shape=(N,N,N), nodes=nodes, S = typeof(xop))
 
   # test agains FourierOperators
   y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N,N)) )