diff --git a/Project.toml b/Project.toml index caf4722..41e899a 100644 --- a/Project.toml +++ b/Project.toml @@ -24,7 +24,7 @@ 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" @@ -32,6 +32,7 @@ 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" diff --git a/ext/LinearOperatorNFFTExt/LinearOperatorNFFTExt.jl b/ext/LinearOperatorNFFTExt/LinearOperatorNFFTExt.jl index 8140e24..6818041 100644 --- a/ext/LinearOperatorNFFTExt/LinearOperatorNFFTExt.jl +++ b/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") diff --git a/ext/LinearOperatorNFFTExt/NFFTOp.jl b/ext/LinearOperatorNFFTExt/NFFTOp.jl index 9662688..0e0d47d 100644 --- a/ext/LinearOperatorNFFTExt/NFFTOp.jl +++ b/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) diff --git a/test/testOperators.jl b/test/testOperators.jl index 31f41ef..0cb150d 100644 --- a/test/testOperators.jl +++ b/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)) )