Better handling of (N)FFT operators for different NFFT backends (#264)

* Readd FFT params to NFFTOp constructor

* Specialise fft parameter for different AbstractNFFT backend

* Fix backend selection module

* Test NonuniformFFTs as well

* Seperate fft and nfft params

* Fix nfftParams for GPU operators

* Fix fft params for gpu

* Just set fftParams

Author: nHackel

MRIOperators/Project.toml

@@ -1,7 +1,7 @@
 name = "MRIOperators"
 uuid = "fb1137e3-90a6-46ce-a672-6e1e53d120f2"
 author = ["Tobias Knopp <tobias@knoppweb.de>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -34,16 +34,18 @@ Reexport = "1"
 LinearOperators = "2.3"
 LinearOperatorCollection = "2"
 NFFT = "0.14"
+NonuniformFFTs = "0.9"
 Wavelets = "0.9, 0.10"
 julia = "1.6"
 
 [extras]
+NonuniformFFTs = "cd96f58b-6017-4a02-bb9e-f4d81626177f"
 JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 MRISimulation = "8988da37-ea20-4fa6-9af7-8a6f6f9a8970"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "MRISimulation", "JLArrays"]
+test = ["Test", "MRISimulation", "JLArrays", "NonuniformFFTs"]
 
 [extensions]
 MRIOperatorsGPUArraysExt = ["GPUArrays", "KernelAbstractions"]

MRIOperators/ext/MRIOperatorsGPUArraysExt/MRIOperatorsGPUArraysExt.jl

@@ -2,13 +2,15 @@ module MRIOperatorsGPUArraysExt
 
 using MRIOperators, GPUArrays, KernelAbstractions
 using MRIOperators.FFTW, MRIOperators.LinearAlgebra
+using MRIOperators.NFFT
 
 include("ExplicitOp.jl")
 include("Shutter.jl")
 include("SensitivityOp.jl")
 include("FieldmapNFFTOp.jl")
 
+# Flags for GPUs seem to break dispatch
 MRIOperators.fftParams(::Type{<:AbstractGPUArray}) = (;)
-
+MRIOperators.nfftParams(::NFFT.NFFTBackend, ::Type{<:AbstractGPUArray}) = (;)
 
 end # module

MRIOperators/src/EncodingOp.jl

@@ -164,7 +164,7 @@ function fourierEncodingOp(shape::NTuple{D,Int64}, tr::Trajectory{T}, opName::St
         ftOp = SamplingOp(Complex{T}; pattern=idx, shape, S = S) ∘ FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...)
       end
     else
-      ftOp = NFFTOp(Complex{T}; nodes = kspaceNodes(tr), shape, S = S, kargs...)
+      ftOp = NFFTOp(Complex{T}; nodes = kspaceNodes(tr), shape, S = S, nfftParams(S)..., kargs...)
     end
   else
     @error "opName $(opName) is not known"

MRIOperators/src/MRIOperators.jl

@@ -83,7 +83,10 @@ function Base.copy(S::LinearOperator{T}) where T
   deepcopy(S)
 end
 
-fftParams(::Type{<:AbstractArray}) = (;:flags => FFTW.MEASURE)
+fftParams(T::Type{<:AbstractArray}) = (;:flags => FFTW.MEASURE)
+nfftParams(T::Type{<:AbstractArray}) = nfftParams(AbstractNFFTs.active_backend(), T)
+nfftParams(::AbstractNFFTBackend, ::Type{<:AbstractArray}) = (;:fftflags => FFTW.MEASURE)
+
 
 # https://github.com/JuliaLang/julia/issues/35543
 stripParameters(arrT::Type{<:AbstractArray}) = Base.typename(arrT).wrapper

MRIOperators/test/runtests.jl

@@ -1,4 +1,5 @@
 using Test, MRIBase, MRIOperators, MRISimulation, MRIOperators.NFFT, MRIOperators.NFFT.FFTW
+using NonuniformFFTs
 using LinearAlgebra, MRIOperators.LinearOperatorCollection
 using JLArrays
 
@@ -7,5 +8,11 @@ areTypesDefined = @isdefined arrayTypes
 arrayTypes = areTypesDefined ? arrayTypes : [Array, JLArray]
 
 @testset "MRIOperators" begin
-  include("testOperators.jl")
+  for backend in [NFFT.backend(), NonuniformFFTs.backend()]
+    with(nfft_backend => backend) do
+      @testset "Operators with $(string(typeof(backend)))" begin
+        include("testOperators.jl")
+      end
+    end
+  end
 end

src/Reconstruction/RecoParameters.jl

@@ -181,6 +181,7 @@ volumeSize(reconSize::NTuple{3,Int}, numSlice::Int) = reconSize
 
 executor(::Type{<:AbstractArray}) = nothing
 copyOpsFn(::Type{<:AbstractArray}) = copy
+# TODO also set nfftParams with new lower bound to MRIOperators 0.4.2
 normalOpParams(::Type{aT}) where aT <: AbstractArray = (; :copyOpsFn => copyOpsFn(aT), MRIOperators.fftParams(aT)...)
 
 executor(f::Function) = executor(f{Complex{Float32}}(undef, 0))