diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml new file mode 100644 index 00000000..8163db81 --- /dev/null +++ b/.buildkite/pipeline.yml @@ -0,0 +1,99 @@ +steps: + - label: "Nvidia GPUs -- MRIReco.jl" + plugins: + - JuliaCI/julia#v1: + version: "1.10" + agents: + queue: "juliagpu" + cuda: "*" + command: | + julia --color=yes --project -e ' + using Pkg + Pkg.develop([PackageSpec(path=pwd(), subdir="MRIBase") + , PackageSpec(path=pwd(), subdir="MRIFiles") + , PackageSpec(path=pwd(), subdir="MRISampling") + , PackageSpec(path=pwd(), subdir="MRISimulation") + , PackageSpec(path=pwd(), subdir="MRIOperators") + , PackageSpec(path=pwd(), subdir="MRICoilSensitivities")]) + Pkg.add("TestEnv") + using TestEnv + TestEnv.activate(); + Pkg.add("CUDA") + Pkg.instantiate() + include("test/gpu/cuda.jl")' + timeout_in_minutes: 30 + + #- label: "AMD GPUs -- MRIReco.jl" + # plugins: + # - JuliaCI/julia#v1: + # version: "1.10" + # agents: + # queue: "juliagpu" + # rocm: "*" + # rocmgpu: "*" + # command: | + # julia --color=yes --project -e ' + # using Pkg + # Pkg.develop([PackageSpec(path=pwd(), subdir="MRIBase") + # , PackageSpec(path=pwd(), subdir="MRIFiles") + # , PackageSpec(path=pwd(), subdir="MRISampling") + # , PackageSpec(path=pwd(), subdir="MRISimulation") + # , PackageSpec(path=pwd(), subdir="MRIOperators") + # , PackageSpec(path=pwd(), subdir="MRICoilSensitivities")]) + # Pkg.add("TestEnv") + # using TestEnv + # TestEnv.activate(); + # Pkg.add("AMDGPU") + # Pkg.instantiate() + # include("test/gpu/rocm.jl")' + # timeout_in_minutes: 30 + + - label: "Nvidia GPUs -- MRIOperators.jl" + plugins: + - JuliaCI/julia#v1: + version: "1.10" + agents: + queue: "juliagpu" + cuda: "*" + command: | + julia --color=yes --project -e ' + using Pkg + Pkg.develop([PackageSpec(path=pwd(), subdir="MRIBase") + , PackageSpec(path=pwd(), subdir="MRIFiles") + , PackageSpec(path=pwd(), subdir="MRISampling") + , PackageSpec(path=pwd(), subdir="MRISimulation") + , PackageSpec(path=pwd(), subdir="MRIOperators") + , PackageSpec(path=pwd(), subdir="MRICoilSensitivities")]) + Pkg.add("TestEnv") + using TestEnv + TestEnv.activate("MRIOperators"); + Pkg.add("CUDA") + Pkg.instantiate() + include("MRIOperators/test/gpu/cuda.jl")' + timeout_in_minutes: 30 + + + - label: "AMD GPUs -- MRIOperators.jl" + plugins: + - JuliaCI/julia#v1: + version: "1.10" + agents: + queue: "juliagpu" + rocm: "*" + rocmgpu: "*" + command: | + julia --color=yes --project -e ' + using Pkg + Pkg.develop([PackageSpec(path=pwd(), subdir="MRIBase") + , PackageSpec(path=pwd(), subdir="MRIFiles") + , PackageSpec(path=pwd(), subdir="MRISampling") + , PackageSpec(path=pwd(), subdir="MRISimulation") + , PackageSpec(path=pwd(), subdir="MRIOperators") + , PackageSpec(path=pwd(), subdir="MRICoilSensitivities")]) + Pkg.add("TestEnv") + using TestEnv + TestEnv.activate("MRIOperators"); + Pkg.add("AMDGPU") + Pkg.instantiate() + include("MRIOperators/test/gpu/rocm.jl")' + timeout_in_minutes: 30 \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 113af779..9a7b759c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -44,12 +44,12 @@ jobs: shell: julia --color=yes --project=. {0} run: | using Pkg - Pkg.develop(PackageSpec(path=pwd(), subdir="MRIBase")) - Pkg.develop(PackageSpec(path=pwd(), subdir="MRIFiles")) - Pkg.develop(PackageSpec(path=pwd(), subdir="MRISampling")) - Pkg.develop(PackageSpec(path=pwd(), subdir="MRISimulation")) - Pkg.develop(PackageSpec(path=pwd(), subdir="MRIOperators")) - Pkg.develop(PackageSpec(path=pwd(), subdir="MRICoilSensitivities")) + Pkg.develop([PackageSpec(path=pwd(), subdir="MRIBase") + , PackageSpec(path=pwd(), subdir="MRIFiles") + , PackageSpec(path=pwd(), subdir="MRISampling") + , PackageSpec(path=pwd(), subdir="MRISimulation") + , PackageSpec(path=pwd(), subdir="MRIOperators") + , PackageSpec(path=pwd(), subdir="MRICoilSensitivities")]) - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - name: Run subpackage tests @@ -78,12 +78,12 @@ jobs: shell: julia --color=yes --project=. {0} run: | using Pkg - Pkg.add(PackageSpec(path=pwd(), subdir="MRIBase")) - Pkg.add(PackageSpec(path=pwd(), subdir="MRIFiles")) - Pkg.add(PackageSpec(path=pwd(), subdir="MRISampling")) - Pkg.add(PackageSpec(path=pwd(), subdir="MRISimulation")) - Pkg.add(PackageSpec(path=pwd(), subdir="MRIOperators")) - Pkg.add(PackageSpec(path=pwd(), subdir="MRICoilSensitivities")) + Pkg.develop([PackageSpec(path=pwd(), subdir="MRIBase") + , PackageSpec(path=pwd(), subdir="MRIFiles") + , PackageSpec(path=pwd(), subdir="MRISampling") + , PackageSpec(path=pwd(), subdir="MRISimulation") + , PackageSpec(path=pwd(), subdir="MRIOperators") + , PackageSpec(path=pwd(), subdir="MRICoilSensitivities")]) - run: | julia --project=docs -e ' using Pkg diff --git a/MRIBase/Project.toml b/MRIBase/Project.toml index 55f8b753..71d595f0 100644 --- a/MRIBase/Project.toml +++ b/MRIBase/Project.toml @@ -1,7 +1,7 @@ name = "MRIBase" uuid = "f7771a9a-6e57-4e71-863b-6e4b6a2f17df" author = ["Tobias Knopp "] -version = "0.4.3" +version = "0.4.4" [deps] AbstractNFFTs = "7f219486-4aa7-41d6-80a7-e08ef20ceed7" diff --git a/MRIBase/src/Trajectories/3D/Cartesian3D.jl b/MRIBase/src/Trajectories/3D/Cartesian3D.jl index bcdf8947..b37e5f74 100644 --- a/MRIBase/src/Trajectories/3D/Cartesian3D.jl +++ b/MRIBase/src/Trajectories/3D/Cartesian3D.jl @@ -23,7 +23,7 @@ function CartesianTrajectory3D(::Type{T}, numProfiles, numSamplingPerProfile , kargs...) where T nodes = cartesian3dNodes(T, numProfiles, numSamplingPerProfile; numSlices=numSlices) times = readoutTimes(T, numProfiles, numSamplingPerProfile, numSlices; TE=TE, AQ=AQ) - return Trajectory("Cartesian3D", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, numSlices, true, false) + return Trajectory("Cartesian3D", nodes, times, T(TE), T(AQ), numProfiles, numSamplingPerProfile, numSlices, true, false) end function cartesian3dNodes(::Type{T}, numProfiles, numSamplingPerProfile diff --git a/MRIOperators/Project.toml b/MRIOperators/Project.toml index 9fc8f981..d5af1311 100644 --- a/MRIOperators/Project.toml +++ b/MRIOperators/Project.toml @@ -1,9 +1,10 @@ name = "MRIOperators" uuid = "fb1137e3-90a6-46ce-a672-6e1e53d120f2" author = ["Tobias Knopp "] -version = "0.2.1" +version = "0.3.0" [deps] +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" MRIBase = "f7771a9a-6e57-4e71-863b-6e4b6a2f17df" @@ -17,15 +18,20 @@ LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" +[weakdeps] +GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" [compat] +Adapt = "3, 4" StatsBase = "0.33, 0.34" Distributions = "0.25" +GPUArrays = "8, 9, 10" +JLArrays = "0.1" julia = "1.6" MRIBase = "0.4" Reexport = "1" LinearOperators = "2.3.3" -LinearOperatorCollection = "1.1" +LinearOperatorCollection = "2" NFFT = "0.13" FLoops = "0.2" LowRankApprox = "0.5" @@ -34,11 +40,10 @@ Wavelets = "0.9, 0.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" MRISimulation = "8988da37-ea20-4fa6-9af7-8a6f6f9a8970" +JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb" [targets] -test = ["Test", "MRISimulation"] - - - - +test = ["Test", "MRISimulation", "JLArrays"] +[extensions] +MRIOperatorsGPUArraysExt = "GPUArrays" \ No newline at end of file diff --git a/MRIOperators/ext/MRIOperatorsGPUArraysExt/ExplicitOp.jl b/MRIOperators/ext/MRIOperatorsGPUArraysExt/ExplicitOp.jl new file mode 100644 index 00000000..85820bd9 --- /dev/null +++ b/MRIOperators/ext/MRIOperatorsGPUArraysExt/ExplicitOp.jl @@ -0,0 +1,110 @@ +function MRIOperators.produ!(out::vecTc, x::vecTc, shape::NTuple{2,Int64}, + nodes::matT, times::vecT, echoOffset::T, + disturbanceTerm::matTc) where {T, Tc <: Union{T, Complex{T}}, vecTc <: AbstractGPUArray{Tc}, matT <: AbstractGPUArray{T}, vecT <: AbstractGPUArray{T}, matTc <: AbstractGPUArray{Tc}} + + factor = Tc(-2 * 1im * pi) + limit = prod(shape) + fill!(out, zero(Tc)) + gpu_call(out, reshape(x, shape), shape, nodes, times, echoOffset, disturbanceTerm; elements = size(nodes, 2)) do ctx, out_, x_, shape_, nodes_, times_, echoOffset_, disturbanceTerm_ + k = linear_index(ctx) + if !(1 <= k <= limit) + return nothing + end + + for nx = 1:shape_[1] + for ny = 1:shape_[2] + phi = (nodes_[1, k] * (nx - shape_[1] / 2 - 1) + nodes_[2, k] * (ny - shape_[2] / 2 - 1)) + e = exp(factor * phi - (times_[k] - echoOffset_) * disturbanceTerm_[nx, ny]) + out_[k] += x_[nx, ny] * e + end + end + return nothing + end + + return out +end + +function MRIOperators.produ!(out::vecTc, x::vecTc, shape::NTuple{3,Int64}, + nodes::matT, times::vecT, echoOffset::T, + disturbanceTerm::arrTc) where {T, Tc <: Union{T, Complex{T}}, vecTc <: AbstractGPUArray{Tc}, matT <: AbstractGPUArray{T}, vecT <: AbstractGPUArray{T}, arrTc <: AbstractGPUArray{Tc, 3}} + + factor = Tc(-2 * 1im * pi) + limit = prod(shape) + fill!(out, zero(Tc)) + gpu_call(out, reshape(x, shape), shape, nodes, times, echoOffset, disturbanceTerm; elements = size(nodes, 2)) do ctx, out_, x_, shape_, nodes_, times_, echoOffset_, disturbanceTerm_ + k = linear_index(ctx) + if !(1 <= k <= limit) + return nothing + end + + for nx = 1:shape_[1] + for ny = 1:shape_[2] + for nz = 1:shape_[3] + phi = (nodes_[1, k] * (nx - shape_[1] / 2 - 1) + nodes_[2, k] * (ny - shape_[2] / 2 - 1) + nodes_[3, k] * (nz - shape_[3] / 2 - 1)) + e = exp(factor * phi - (times_[k] - echoOffset_) * disturbanceTerm_[nx, ny, nz]) + out_[k] += x_[nx, ny, nz] * e + end + end + end + return nothing + end + + return out +end + +function MRIOperators.ctprodu!(out::vecTc, x::vecTc, shape::NTuple{2,Int64}, + nodes::matT, times::vecT, echoOffset::T, + disturbanceTerm::matTc) where {T, Tc <: Union{T, Complex{T}}, vecTc <: AbstractGPUArray{Tc}, matT <: AbstractGPUArray{T}, vecT <: AbstractGPUArray{T}, matTc <: AbstractGPUArray{Tc}} + + factor = Tc(-2 * 1im * pi) + limit = prod(shape) + fill!(out, zero(Tc)) + gpu_call(reshape(out, shape), x, shape, nodes, times, echoOffset, disturbanceTerm; elements = limit) do ctx, out_, x_, shape_, nodes_, times_, echoOffset_, disturbanceTerm_ + linearIndex = linear_index(ctx) + if !(1 <= linearIndex <= limit) + return nothing + end + + cartIndex = CartesianIndices(shape)[linearIndex] + + nx = cartIndex[1] + ny = cartIndex[2] + for k = 1:size(nodes_, 2) + phi = (nodes_[1,k]*(nx-shape_[1]/2-1)+nodes_[2,k]*(ny-shape_[2]/2-1)) + e = exp(factor * phi - (times_[k]-echoOffset_)*disturbanceTerm_[cartIndex]) + out_[cartIndex] += x_[k] * conj(e) + end + return nothing + end + + return out +end + +function MRIOperators.ctprodu!(out::vecTc, x::vecTc, shape::NTuple{3,Int64}, + nodes::matT, times::vecT, echoOffset::T, + disturbanceTerm::arrTc) where {T, Tc <: Union{T, Complex{T}}, vecTc <: AbstractGPUArray{Tc}, matT <: AbstractGPUArray{T}, vecT <: AbstractGPUArray{T}, arrTc <: AbstractGPUArray{Tc, 3}} + + factor = Tc(-2 * 1im * pi) + limit = prod(shape) + fill!(out, zero(Tc)) + gpu_call(reshape(out, shape), x, shape, nodes, times, echoOffset, disturbanceTerm; elements = limit) do ctx, out_, x_, shape_, nodes_, times_, echoOffset_, disturbanceTerm_ + linearIndex = linear_index(ctx) + if !(1 <= linearIndex <= limit) + return nothing + end + + cartIndex = CartesianIndices(shape)[linearIndex] + + nx = cartIndex[1] + ny = cartIndex[2] + nz = cartIndex[3] + for k = 1:size(nodes_, 2) + phi = (nodes_[1,k]*(nx-shape_[1]/2-1)+nodes_[2,k]*(ny-shape_[2]/2-1)+nodes_[3,k]*(nz-shape_[3]/2-1)) + e = exp(factor * phi - (times_[k]-echoOffset_)*disturbanceTerm_[cartIndex]) + out_[cartIndex] += x_[k] * conj(e) + end + return nothing + end + + return out +end \ No newline at end of file diff --git a/MRIOperators/ext/MRIOperatorsGPUArraysExt/FieldmapNFFTOp.jl b/MRIOperators/ext/MRIOperatorsGPUArraysExt/FieldmapNFFTOp.jl new file mode 100644 index 00000000..36a47372 --- /dev/null +++ b/MRIOperators/ext/MRIOperatorsGPUArraysExt/FieldmapNFFTOp.jl @@ -0,0 +1,41 @@ +function MRIOperators.produ_inner!(K, C::matT, A::matT, shape, d::Vector{vecT}, s::vecT, sp, plan, idx, x_::vecT, p::Vector{arrT}) where {T, vecT <: AbstractGPUVector{T}, matT <: AbstractGPUMatrix{T}, arrT <: AbstractGPUArray{T}} + for κ=1:K + if !isempty(idx[κ]) + p[κ][:] .= C[κ,:] .* x_ + mul!(d[κ], plan[κ], p[κ]) + # Assumption: l is unique per idx[κ] + gpu_call(idx[κ], d[κ], view(A, :, κ), s) do ctx, indices, d_, A_, s_ + k = @linearidx(indices) + l = indices[k] + s_[l] += d_[k]*A_[l] + return nothing + end + end + end + + return +end + + +function MRIOperators.ctprodu_inner!(K, C::matT, A::matT, shape, d::Vector{vecT}, y::vecT, sp, plan, idx, x::vecT, p::Vector{arrT}) where {T, vecT <: AbstractGPUVector{T}, matT <: AbstractGPUMatrix{T}, arrT <: AbstractGPUArray{T}} + + for κ=1:K + if !isempty(idx[κ]) + gpu_call(idx[κ], d[κ], view(A, :, κ), x) do ctx, indices, d_, A_, x_ + k = @linearidx(indices) + l = indices[k] + d_[k] = conj(A_[l]) * x_[l] + return nothing + end + mul!(p[κ], adjoint(plan[κ]), d[κ]) + + gpu_call(p[κ], view(C, κ, :), y) do ctx, p_, C_, y_ + k = @linearidx(p_) + y_[k] += conj(C_[k]) * p_[k] + return nothing + end + end + end + + return +end \ No newline at end of file diff --git a/MRIOperators/ext/MRIOperatorsGPUArraysExt/MRIOperatorsGPUArraysExt.jl b/MRIOperators/ext/MRIOperatorsGPUArraysExt/MRIOperatorsGPUArraysExt.jl new file mode 100644 index 00000000..d20bbdfa --- /dev/null +++ b/MRIOperators/ext/MRIOperatorsGPUArraysExt/MRIOperatorsGPUArraysExt.jl @@ -0,0 +1,13 @@ +module MRIOperatorsGPUArraysExt + +using MRIOperators, GPUArrays, MRIOperators.FFTW, MRIOperators.LinearAlgebra + +include("ExplicitOp.jl") +include("Shutter.jl") +include("SensitivityOp.jl") +include("FieldmapNFFTOp.jl") + +MRIOperators.fftParams(::Type{<:AbstractGPUArray}) = (;) + + +end # module diff --git a/MRIOperators/ext/MRIOperatorsGPUArraysExt/SensitivityOp.jl b/MRIOperators/ext/MRIOperatorsGPUArraysExt/SensitivityOp.jl new file mode 100644 index 00000000..3497c09d --- /dev/null +++ b/MRIOperators/ext/MRIOperatorsGPUArraysExt/SensitivityOp.jl @@ -0,0 +1,31 @@ +function MRIOperators.prod_smap!(y::vecT, smaps::matT, x::vecT, numVox, numChan, numContr=1) where {T, vecT <: AbstractGPUVector{T}, matT <: AbstractGPUMatrix{T}} + x_ = reshape(x,numVox,numContr) + y_ = reshape(y,numVox,numChan,numContr) + + @assert size(smaps) == (size(y_,1), size(y_,2)) + + gpu_call(y_, x_, smaps) do ctx, ygpu, xgpu, smapsgpu + cart = @cartesianidx(ygpu) + ygpu[cart] = xgpu[cart[1],cart[3]] * smapsgpu[cart[1],cart[2]] + return nothing + end + return y +end + +function MRIOperators.ctprod_smap!(y::vecT, smapsC::matT, x::vecT, numVox, numChan, numContr=1) where {T, vecT <: AbstractGPUVector{T}, matT <: AbstractGPUMatrix{T}} + x_ = reshape(x,numVox,numChan,numContr) + y_ = reshape(y,numVox,numContr) + + @assert size(smapsC) == (size(x_,1), size(x_,2)) + + y_ .= 0 + # Inner loop to avoid race condition + gpu_call(y_, x_, smapsC, numChan) do ctx, ygpu, xgpu, smapsCgpu, limit + cart = @cartesianidx(ygpu) + for i = 1:limit + ygpu[cart] += xgpu[cart[1], i, cart[2]] * smapsCgpu[cart[1],i] + end + return nothing + end + return y +end \ No newline at end of file diff --git a/MRIOperators/ext/MRIOperatorsGPUArraysExt/Shutter.jl b/MRIOperators/ext/MRIOperatorsGPUArraysExt/Shutter.jl new file mode 100644 index 00000000..148558a4 --- /dev/null +++ b/MRIOperators/ext/MRIOperatorsGPUArraysExt/Shutter.jl @@ -0,0 +1,31 @@ +function MRIOperators.circularShutter!(I::AbstractGPUArray, radiusFactor::Number=1.0) + imgSize = size(I) + center = imgSize./2.0 + radius = maximum(center) * radiusFactor + # Applying filtering + gpu_call(I, center, radius) do ctx, I_, center_, radius_ + cart = @cartesianidx(I_) + if sqrt(sum((Tuple(cart) .-center_).^2)) > radius_ + I[cart] = 0 + end + return nothing + end + return I +end + + +function MRIOperators.circularShutterFreq!(I::AbstractGPUArray, radiusFactor::T=1.0) where T<:Number + imgSize = size(I) + fftI = fftshift(fft(I)) + center = imgSize./2.0 + radius = maximum(center) * radiusFactor + # Applying filtering + gpu_call(fftI, center, radius) do ctx, fftI_, center_, radius_ + cart = @cartesianidx(fftI_) + if sqrt(sum((Tuple(cart) .-center_).^2)) > radius_ + fftI_[cart] = 0 + end + return nothing + end + return MRIOperators.preserveType(I, ifft(ifftshift(fftI))) +end diff --git a/MRIOperators/src/Composition.jl b/MRIOperators/src/Composition.jl deleted file mode 100644 index 0e209466..00000000 --- a/MRIOperators/src/Composition.jl +++ /dev/null @@ -1,146 +0,0 @@ -""" - `mutable struct CompositeOp{T}` - - struct describing the result of a composition/product of operators. - Describing the composition using a dedicated type has the advantage - that the latter can be made copyable. This is particularly relevant for - multi-threaded code -""" -mutable struct CompositeOp{T,U,V} <: AbstractLinearOperator{T} - nrow :: Int - ncol :: Int - symmetric :: Bool - hermitian :: Bool - prod! :: Function - tprod! :: Function - ctprod! :: Function - nprod :: Int - ntprod :: Int - nctprod :: Int - args5 :: Bool - use_prod5! :: Bool - allocated5 :: Bool - Mv5 :: Vector{T} - Mtu5 :: Vector{T} - isWeighting :: Bool - A::U - B::V - tmp::Vector{T} -end - -""" - CompositeOp(ops :: AbstractLinearOperator...) - -composition/product of two Operators. Differs with * since it can handle normal operator -""" -function CompositeOp(A,B;isWeighting=false) - nrow = A.nrow - ncol = B.ncol - S = promote_type(eltype(A), eltype(B)) - tmp_ = Vector{S}(undef, B.nrow) - - function produ!(res, x::AbstractVector{T}, tmp) where T<:Union{Real,Complex} - mul!(tmp, B, x) - return mul!(res, A, tmp) - end - - function tprodu!(res, y::AbstractVector{T}, tmp) where T<:Union{Real,Complex} - mul!(tmp, transpose(A), y) - return mul!(res, transpose(B), tmp) - end - - function ctprodu!(res, y::AbstractVector{T}, tmp) where T<:Union{Real,Complex} - mul!(tmp, adjoint(A), y) - return mul!(res, adjoint(B), tmp) - end - - Op = CompositeOp( nrow, ncol, false, false, - (res,x) -> produ!(res,x,tmp_), - (res,y) -> tprodu!(res,y,tmp_), - (res,y) -> ctprodu!(res,y,tmp_), - 0, 0, 0, false, false, false, S[], S[], - isWeighting, A, B, tmp_) - - return Op -end - -LinearOperators.storage_type(op::CompositeOp) = typeof(op.Mv5) - -""" -∘(A::T1, B::T2; isWeighting::Bool=false) where {T1<:AbstractLinearOperator, T2<:AbstractLinearOperator} - - composition of two operators. -""" -function Base.:∘(A::T1, B::T2; isWeighting::Bool=false) where {T1<:AbstractLinearOperator, T2<:AbstractLinearOperator} - return CompositeOp(A,B;isWeighting=isWeighting) -end - -function Base.copy(S::CompositeOp{T}) where T - A = copy(S.A) - B = copy(S.B) - return CompositeOp(A,B; isWeighting=S.isWeighting) -end - - -### Normal Matrix Code ### -# Left matrix can be build into a normal operator - -mutable struct CompositeNormalOp{T,S,U,V} <: AbstractLinearOperator{T} - nrow :: Int - ncol :: Int - symmetric :: Bool - hermitian :: Bool - prod! :: Function - tprod! :: Nothing - ctprod! :: Nothing - nprod :: Int - ntprod :: Int - nctprod :: Int - args5 :: Bool - use_prod5! :: Bool - allocated5 :: Bool - Mv5 :: Vector{T} - Mtu5 :: Vector{T} - opOuter::S - normalOpInner::U - tmp::V -end - -LinearOperators.storage_type(op::CompositeNormalOp) = typeof(op.Mv5) - - -function CompositeNormalOp(opOuter, normalOpInner, tmp::Vector{T}) where T - - function produ!(y, opOuter, normalOpInner, tmp, x) - mul!(tmp, opOuter, x) - mul!(tmp, normalOpInner, tmp) # This can be dangerous. We might need to create two tmp vectors - return mul!(y, adjoint(opOuter), tmp) - end - - return CompositeNormalOp(size(opOuter,2), size(opOuter,2), false, false - , (res,x) -> produ!(res, opOuter, normalOpInner, tmp, x) - , nothing - , nothing - , 0, 0, 0, false, false, false, T[], T[] - , opOuter, normalOpInner, tmp) -end - - -function LinearOperatorCollection.normalOperator(S::CompositeOp, W=opEye(eltype(S),size(S,1))) - T = promote_type(eltype(S.A), eltype(S.B), eltype(W)) - if S.isWeighting #&& typeof(W) <: opEye - # In this case we are converting the left argument into a - # weighting matrix, that is passed to normalOperator - normalOperator(S.B, S.A) - else - tmp = Vector{T}(undef, size(S.A, 2)) - return CompositeNormalOp(S.B, normalOperator(S.A, W), tmp) - end -end - -function Base.copy(S::CompositeNormalOp) - opOuter = copy(S.opOuter) - opInner = copy(S.normalOpInner) - tmp = copy(S.tmp) - return CompositeNormalOp(opOuter, opInner, tmp) -end diff --git a/MRIOperators/src/EncodingOp.jl b/MRIOperators/src/EncodingOp.jl index 9a4e3b1a..09858480 100644 --- a/MRIOperators/src/EncodingOp.jl +++ b/MRIOperators/src/EncodingOp.jl @@ -35,16 +35,16 @@ in terms of their sensitivities * `senseMaps::Array{Complex{T}}` - coil sensitivities """ function encodingOps_parallel(acqData::AcquisitionData{T,D}, shape::NTuple{D,Int64} - , senseMaps::Array{Complex{T},4} - ; slice=1, kargs...) where {T,D} + , senseMaps::AbstractArray{Complex{T},4} + ; slice=1, S = Vector{Complex{T}}, copyOpsFn = copy, kargs...) where {T,D} smaps = ( D==2 ? senseMaps[:,:,slice,:] : senseMaps ) numContr, numChan = numContrasts(acqData), numChannels(acqData) # fourier operators - ft = encodingOps_simple(acqData, shape; slice=slice, kargs...) - S = SensitivityOp(reshape(smaps,:,numChan),1) - Op = [ DiagOp(ft[i], numChan) ∘ S for i=1:numContr] + ft = encodingOps_simple(acqData, shape; slice=slice, S = S, kargs...) + SOp = SensitivityOp(reshape(smaps,:,numChan),1) + Op = [ DiagOp(ft[i], numChan; copyOpsFn = copyOpsFn) ∘ SOp for i=1:numContr] return Op end @@ -81,17 +81,17 @@ in terms of their sensitivities * `senseMaps::Array{Complex{T}}` - coil sensitivities """ function encodingOp_multiEcho_parallel(acqData::AcquisitionData{T,D}, shape::NTuple{D,Int64} - , senseMaps::Array{Complex{T}} - ; slice::Int64=1, kargs...) where {T,D} + , senseMaps::AbstractArray{Complex{T}} + ; slice::Int64=1, copyOpsFn = copy, kargs...) where {T,D} smaps = ( D==2 ? senseMaps[:,:,slice,:] : senseMaps ) numChan = numChannels(acqData) # fourier operators ft = encodingOps_simple(acqData, shape; kargs...) - S = SensitivityOp(reshape(smaps,:,numChan),numContrasts(acqData)) - ops2 = [copy(ft[n]) for j=1:numChan,n=eachindex(ft)] - return DiagOp(ops2...) ∘ S + SOp = SensitivityOp(reshape(smaps,:,numChan),numContrasts(acqData)) + ops2 = [copyOpsFn(ft[n]) for j=1:numChan,n=eachindex(ft)] + return DiagOp(ops2...) ∘ SOp end ################################### @@ -139,7 +139,7 @@ return Fourier encoding operator (either Explicit or NFFT) """ function fourierEncodingOp(shape::NTuple{D,Int64}, tr::Trajectory{T}, opName::String; subsampleIdx::Vector{Int64}=Int64[], slice::Int64=1, correctionMap::Array{Complex{T}}=Complex{T}[], - echoImage::Bool=true, kargs...) where {T,D} + echoImage::Bool=true, S = Vector{Complex{T}}, kargs...) where {T,D} # extract proper portion of correctionMap if !isempty(correctionMap) @@ -148,21 +148,21 @@ function fourierEncodingOp(shape::NTuple{D,Int64}, tr::Trajectory{T}, opName::St # Fourier transformations if opName=="explicit" @debug "ExplicitOp" - ftOp = ExplicitOp(shape, tr, cmap, echoImage=echoImage) + ftOp = ExplicitOp(shape, tr, cmap, echoImage=echoImage, S = S) elseif opName=="fast" @debug "NFFT-based Op" if !isempty(correctionMap) && correctionMap!=zeros(Complex{T},size(correctionMap)) - ftOp = FieldmapNFFTOp(shape, tr, cmap, echoImage=echoImage; kargs...) + ftOp = FieldmapNFFTOp(shape, tr, cmap, echoImage=echoImage; S = S, fftParams(S)..., kargs...) elseif isCartesian(tr) @debug "FFTOp" if !MRIBase.isUndersampledCartTrajectory(shape,tr) - ftOp = FFTOp(Complex{T}; shape, unitary=false) + ftOp = FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...) else idx = MRIBase.cartesianSubsamplingIdx(shape,tr) - ftOp = SamplingOp(Complex{T}; pattern=idx, shape) ∘ FFTOp(Complex{T}; shape, unitary=false) + ftOp = SamplingOp(Complex{T}; pattern=idx, shape, S = S) ∘ FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...) end else - ftOp = NFFTOp(tr; shape, kargs...) + ftOp = NFFTOp(Complex{T}; nodes = kspaceNodes(tr), shape, S = S, fftParams(S)..., kargs...) end else @error "opName $(opName) is not known" @@ -171,7 +171,7 @@ function fourierEncodingOp(shape::NTuple{D,Int64}, tr::Trajectory{T}, opName::St # subsampling if !isempty(subsampleIdx) && (subsampleIdx != collect(1:size(tr,2))) && isCartesian(tr) β = (D==2) ? (tr.numSamplingPerProfile, tr.numProfiles) : (tr.numSamplingPerProfile, tr.numProfiles, tr.numSlices) - S = SamplingOp(Complex{T}; pattern = subsampleIdx, shape=β) + S = SamplingOp(Complex{T}; pattern = subsampleIdx, shape=β, S = S) return S ∘ ftOp else return ftOp diff --git a/MRIOperators/src/ExplicitOp.jl b/MRIOperators/src/ExplicitOp.jl index c58f9755..807a5be2 100644 --- a/MRIOperators/src/ExplicitOp.jl +++ b/MRIOperators/src/ExplicitOp.jl @@ -1,6 +1,6 @@ export ExplicitOp -mutable struct ExplicitOp{T,F1,F2} <: AbstractLinearOperator{T} +mutable struct ExplicitOp{T, vecT <: AbstractVector{T}, F1, F2} <: AbstractLinearOperator{T} nrow :: Int ncol :: Int symmetric :: Bool @@ -14,8 +14,8 @@ mutable struct ExplicitOp{T,F1,F2} <: AbstractLinearOperator{T} args5 :: Bool use_prod5! :: Bool allocated5 :: Bool - Mv5 :: Vector{T} - Mtu5 :: Vector{T} + Mv5 :: vecT + Mtu5 :: vecT end LinearOperators.storage_type(op::ExplicitOp) = typeof(op.Mv5) @@ -34,39 +34,46 @@ generates a `ExplicitOp` which explicitely evaluates the MRI Fourier signal enco this results in complex valued image even for real-valued input. * `kargs` - additional keyword arguments """ -function ExplicitOp(shape::NTuple{D,Int64}, tr::Trajectory, correctionmap::Array{ComplexF64,D} - ; echoImage::Bool=false, kargs...) where D +function ExplicitOp(shape::NTuple{D,Int64}, tr::Trajectory{T}, correctionmap::AbstractArray{Tc,D} + ; echoImage::Bool=false, S = storage_type(correctionmap), kargs...) where {T, Tc <: Union{Complex{T}, T}, D} nodes,times = kspaceNodes(tr), readoutTimes(tr) nrow = size(nodes,2) ncol = prod(shape) + if isempty(correctionmap) + disturbanceTerm = zeros(Complex{T}, shape...) + else + disturbanceTerm = correctionmap + end + + tmp = S(undef, 0) + baseArrayType = stripParameters(S) + nodes = adapt(baseArrayType, nodes) + times = adapt(baseArrayType, times) + disturbanceTerm = adapt(baseArrayType, disturbanceTerm) + # if echo image is desired echo time is needed as an offset if echoImage echoOffset = echoTime(tr) else echoOffset = 0.0 end + echoOffset = T(echoOffset) - return ExplicitOp{ComplexF64,Nothing,Function}(nrow, ncol, false, false - , (res,x)->(res .= produ(x, shape, nodes, times, echoOffset, correctionmap)) + return ExplicitOp{Complex{T}, S, Nothing, Function}(nrow, ncol, false, false + , (res,x)->(produ!(res, x, shape, nodes, times, echoOffset, disturbanceTerm)) , nothing - , (res,y)->(res .= ctprodu(y, shape, nodes, times, echoOffset, correctionmap)) - , 0,0,0, false, false, false, ComplexF64[], ComplexF64[]) + , (res,y)->(ctprodu!(res, y, shape, nodes, times, echoOffset, disturbanceTerm)) + , 0,0,0, false, false, false, tmp, tmp) end -function produ(x::Vector{T}, shape::NTuple{2,Int64}, - nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64, - correctionmap::Matrix{ComplexF64}) where T<:Union{Real,Complex} - - if isempty(correctionmap) - disturbanceTerm = zeros(ComplexF64,shape...) - else - disturbanceTerm = correctionmap - end +function produ!(out::Vector{Tc}, x::Vector{Tc}, shape::NTuple{2,Int64}, + nodes::Matrix{T}, times::Vector{T}, echoOffset::T, + disturbanceTerm::Matrix{Tc}) where {T, Tc <: Union{T, Complex{T}}} x= reshape(x,shape) - out = zeros(ComplexF64,size(nodes,2)) + fill!(out, zero(Tc)) for nx=1:shape[1] for ny=1:shape[2] for k=1:size(nodes,2) @@ -77,21 +84,16 @@ function produ(x::Vector{T}, shape::NTuple{2,Int64}, end end end - return vec(out) + return out end -function produ(x::Vector{T}, shape::NTuple{3,Int64}, - nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64, - correctionmap::Array{ComplexF64,3}) where T<:Union{Real,Complex} +function produ!(out::Vector{Tc}, x::Vector{Tc}, shape::NTuple{3,Int64}, + nodes::Matrix{T}, times::Vector{T}, echoOffset::T, + disturbanceTerm::Array{Tc,3}) where {T, Tc <: Union{T, Complex{T}}} - if isempty(correctionmap) - disturbanceTerm = zeros(ComplexF64,shape...) - else - disturbanceTerm = correctionmap - end x = reshape(x,shape) - out = zeros(ComplexF64,size(nodes,2)) + fill!(out, zero(Tc)) for nx=1:shape[1] for ny=1:shape[2] for nz=1:shape[3] @@ -103,20 +105,16 @@ function produ(x::Vector{T}, shape::NTuple{3,Int64}, end end end - return vec(out) + return out end -function ctprodu(x::Vector{T}, shape::NTuple{2,Int64}, - nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64, - correctionmap::Matrix{ComplexF64}) where T<:Union{Real,Complex} +function ctprodu!(out::Vector{Tc}, x::Vector{Tc}, shape::NTuple{2,Int64}, + nodes::Matrix{T}, times::Vector{T}, echoOffset::T, + disturbanceTerm::Matrix{Tc}) where {T, Tc <: Union{T, Complex{T}}} - if isempty(correctionmap) - disturbanceTerm = zeros(ComplexF64,shape...) - else - disturbanceTerm = correctionmap - end - out = zeros(ComplexF64,shape) + out = reshape(out, shape) + fill!(out, zero(Tc)) for nx=1:shape[1] for ny=1:shape[2] for k=1:size(nodes,2) @@ -126,20 +124,15 @@ function ctprodu(x::Vector{T}, shape::NTuple{2,Int64}, end end end - return vec(out) + return out end -function ctprodu(x::Vector{T}, shape::NTuple{3,Int64}, - nodes::Matrix{Float64}, times::Vector{Float64}, echoOffset::Float64, - correctionmap::Array{ComplexF64,3}) where T<:Union{Real,Complex} - - if isempty(correctionmap) - disturbanceTerm = zeros(ComplexF64,shape...) - else - disturbanceTerm = correctionmap - end +function ctprodu!(out::Vector{Tc}, x::Vector{Tc}, shape::NTuple{3,Int64}, + nodes::Matrix{T}, times::Vector{T}, echoOffset::T, + disturbanceTerm::Array{Tc,3}) where {T, Tc <: Union{T, Complex{T}}} - out = zeros(ComplexF64,shape) + out = reshape(out, shape) + fill!(out, zero(Tc)) for nx=1:shape[1] for ny=1:shape[2] for nz=1:shape[3] @@ -151,10 +144,10 @@ function ctprodu(x::Vector{T}, shape::NTuple{3,Int64}, end end end - return vec(out) + return out end function Base.adjoint(op::ExplicitOp{T}) where T return LinearOperator{T}(op.ncol, op.nrow, op.symmetric, op.hermitian, - op.ctprod!, nothing, op.prod!) + op.ctprod!, nothing, op.prod!; S = LinearOperators.storage_type(op)) end diff --git a/MRIOperators/src/FieldmapNFFTOp.jl b/MRIOperators/src/FieldmapNFFTOp.jl index cf097653..64e2afac 100644 --- a/MRIOperators/src/FieldmapNFFTOp.jl +++ b/MRIOperators/src/FieldmapNFFTOp.jl @@ -2,17 +2,19 @@ export FieldmapNFFTOp, InhomogeneityData, createInhomogeneityData_ include("ExpApproximation.jl") -mutable struct InhomogeneityData{T} - A_k::Matrix{Complex{T}} - C_k::Matrix{Complex{T}} - times::Vector{T} - Cmap::Matrix{Complex{T}} +mutable struct InhomogeneityData{T, matT <: AbstractArray{Complex{T}, 2}, vecT <: AbstractArray{T, 1}} + A_k::matT + C_k::matT + times::vecT + Cmap::matT t_hat::T z_hat::Complex{T} method::String end -mutable struct FieldmapNFFTOp{T,F1,F2,D} <:AbstractLinearOperator{Complex{T}} +Adapt.adapt_structure(::Type{arrT}, data::InhomogeneityData) where arrT = InhomogeneityData(adapt(arrT, data.A_k), adapt(arrT, data.C_k), adapt(arrT, data.times), adapt(arrT, data.Cmap), data.t_hat, data.z_hat, data.method) + +mutable struct FieldmapNFFTOp{T, vecT <: AbstractVector{Complex{T}},F1,F2,D, vecI, matT, vecTR} <:AbstractLinearOperator{Complex{T}} nrow :: Int ncol :: Int symmetric :: Bool @@ -26,14 +28,13 @@ mutable struct FieldmapNFFTOp{T,F1,F2,D} <:AbstractLinearOperator{Complex{T}} args5 :: Bool use_prod5! :: Bool allocated5 :: Bool - Mv5 :: Vector{Complex{T}} - Mtu5 :: Vector{Complex{T}} + Mv5 :: vecT + Mtu5 :: vecT plans - idx::Vector{Vector{Int64}} + idx::Vector{vecI} circTraj::Bool shape::NTuple{D,Int64} - cparam::InhomogeneityData{T} - + cparam::InhomogeneityData{T, matT, vecTR} end LinearOperators.storage_type(op::FieldmapNFFTOp) = typeof(op.Mv5) @@ -73,8 +74,8 @@ function FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory, K=20, K_tol::Float64=1.e-3, step::Int64=10, + S = Vector{Complex{T}}, kargs...) where {T,D} - nodes,times = kspaceNodes(tr), readoutTimes(tr) if echoImage times = times .- echoTime(tr) @@ -82,31 +83,46 @@ function FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory, nrow = size(nodes,2) ncol = prod(shape) - x_tmp = zeros(Complex{T}, ncol) - y_tmp = zeros(Complex{T}, nrow) - # create and truncate low-rank expansion cparam = createInhomogeneityData_(vec(times), correctionmap; K=K, alpha=alpha, m=m, method=method, K_tol=K_tol, numSamp=numSamplingPerProfile(tr),step=step) K = size(cparam.A_k,2) @debug "K = $K" - plans = Vector{NFFT.NFFTPlan{T,D,1}}(undef,K) - idx = Vector{Vector{Int64}}(undef,K) + plans = [] # Dont fully specify type yet + idx = [] + baseArrayType = stripParameters(S) for κ=1:K - idx[κ] = findall(x->x!=0.0, cparam.A_k[:,κ]) - plans[κ] = plan_nfft(nodes[:,idx[κ]], shape, m=3, σ=1.25, precompute = NFFT.POLYNOMIAL) + posIndices = findall(x->x!=0.0, cparam.A_k[:,κ]) + if !isempty(posIndices) + push!(idx, posIndices) + push!(plans, plan_nfft(baseArrayType, nodes[:,idx[κ]], shape, m=3, σ=1.25, precompute = NFFT.POLYNOMIAL)) + end end + plans = identity.(plans) # This gives the vectors (more) concrete types + idx = identity.(idx) + K=length(plans) - d = [zeros(Complex{T}, length(idx[κ])) for κ=1:K ] - p = [zeros(Complex{T}, shape) for κ=1:K] + cparam = adapt(baseArrayType, cparam) + idx = adapt.(baseArrayType, idx) + + if !isa(first(idx), Array) + idx = map(i -> Int32.(i), idx) + end + + tmp = S(undef, 0) + x_tmp = fill!(similar(tmp, Complex{T}, ncol), zero(Complex{T})) + y_tmp = fill!(similar(tmp, Complex{T}, nrow), zero(Complex{T})) + + d = [fill!(similar(tmp, Complex{T}, length(idx[κ])), zero(Complex{T})) for κ=1:K ] + p = [fill!(similar(tmp, Complex{T}, shape), zero(Complex{T})) for κ=1:K] circTraj = isCircular(tr) - return FieldmapNFFTOp{T,Nothing,Function,D}(nrow, ncol, false, false + return FieldmapNFFTOp{T, typeof(tmp), Nothing, Function, D, eltype(idx), typeof(cparam.A_k), typeof(cparam.times)}(nrow, ncol, false, false , (res,x) -> produ!(res,x,x_tmp,shape,plans,idx,cparam,circTraj,d,p) , nothing - , (res,y) -> ctprodu!(res,y,y_tmp,shape,plans,idx,cparam,circTraj,d,p), 0, 0, 0, false, false, false, Complex{T}[], Complex{T}[] + , (res,y) -> ctprodu!(res,y,y_tmp,shape,plans,idx,cparam,circTraj,d,p), 0, 0, 0, false, false, false, tmp, tmp , plans, idx, circTraj, shape, cparam) end @@ -116,7 +132,7 @@ function Base.copy(cparam::InhomogeneityData{T}) where T end -function Base.copy(S::FieldmapNFFTOp{T,Nothing,Function,D}) where {T,D} +function Base.copy(S::FieldmapNFFTOp{T}) where {T} shape = S.shape @@ -124,12 +140,12 @@ function Base.copy(S::FieldmapNFFTOp{T,Nothing,Function,D}) where {T,D} plans = [copy(S.plans[i]) for i=1:K] idx = copy(S.idx) - x_tmp = zeros(Complex{T}, S.ncol) - y_tmp = zeros(Complex{T}, S.nrow) + x_tmp = fill!(similar(S.Mv5, Complex{T}, S.ncol), zero(Complex{T})) + y_tmp = fill!(similar(S.Mv5, Complex{T}, S.nrow), zero(Complex{T})) cparam = copy(S.cparam) - d = [zeros(Complex{T}, length(idx[κ])) for κ=1:K] - p = [zeros(Complex{T}, shape) for κ=1:K] + d = [fill!(similar(S.Mv5, Complex{T}, length(idx[κ])), zero(Complex{T})) for κ=1:K ] + p = [fill!(similar(S.Mv5, Complex{T}, shape), zero(Complex{T})) for κ=1:K] D_ = length(shape) circTraj = S.circTraj @@ -137,19 +153,19 @@ function Base.copy(S::FieldmapNFFTOp{T,Nothing,Function,D}) where {T,D} mul! = (res,x) -> produ!(res,x,x_tmp,shape,plans,idx,cparam,circTraj,d,p) ctmul! = (res,y) -> ctprodu!(res,y,y_tmp,shape,plans,idx,cparam,circTraj,d,p) - return FieldmapNFFTOp{T,Nothing,Function,D_}(S.nrow, S.ncol, false, false + return FieldmapNFFTOp{T, typeof(x_tmp), Nothing, Function, D_, eltype(idx), typeof(cparam.A_k), typeof(cparam.times)}(S.nrow, S.ncol, false, false , mul! , nothing , ctmul!, 0, 0, 0, false, false, false, Complex{T}[], Complex{T}[] , plans, idx, circTraj, shape, cparam) end -function produ!(s::AbstractVector{T}, x::AbstractVector{T}, x_tmp::Vector{T},shape::Tuple, plan, - idx::Vector{Vector{Int64}}, cparam::InhomogeneityData, - shutter::Bool, d::Vector{Vector{T}}, p::Vector{Array{T,D}}) where {T,D} +function produ!(s::AbstractVector{T}, x::AbstractVector{T}, x_tmp::AbstractVector{T},shape::Tuple, plan, + idx, cparam::InhomogeneityData, + shutter::Bool, d, p) where {T} s .= zero(T) - K = size(cparam.A_k,2) + K=length(plan) if shutter circularShutter!(reshape(x, shape), 1.0) @@ -187,11 +203,11 @@ function produ_inner!(K, C, A, shape, d, s, sp, plan, idx, x_, p) end -function ctprodu!(y::AbstractVector{Complex{T}}, x::AbstractVector{Complex{T}}, x_tmp::Vector{Complex{T}}, shape::Tuple, plan, idx::Vector{Vector{Int64}}, - cparam::InhomogeneityData{T}, shutter::Bool, d::Vector{Vector{Complex{T}}}, p::Vector{Array{Complex{T},D}}) where {T,D} +function ctprodu!(y::AbstractVector{Complex{T}}, x::AbstractVector{Complex{T}}, x_tmp::AbstractVector{Complex{T}}, shape::Tuple, plan, idx, + cparam::InhomogeneityData{T}, shutter::Bool, d, p) where {T} y .= zero(Complex{T}) - K = size(cparam.A_k,2) + K=length(plan) x_tmp .= x diff --git a/MRIOperators/src/MRIOperators.jl b/MRIOperators/src/MRIOperators.jl index 8503dab0..7d1f23fa 100644 --- a/MRIOperators/src/MRIOperators.jl +++ b/MRIOperators/src/MRIOperators.jl @@ -2,6 +2,7 @@ module MRIOperators using Base: hcat, vcat, \ +using Adapt using Reexport using MRIBase @reexport using LinearOperatorCollection @@ -15,14 +16,10 @@ using LowRankApprox: psvd using Distributions using StatsBase -export DiagOp - include("Shutter.jl") -include("Composition.jl") include("ExplicitOp.jl") include("SensitivityOp.jl") include("MapSliceOp.jl") -include("NFFTOp.jl") include("FieldmapNFFTOp.jl") include("EncodingOp.jl") include("SparseOp.jl") @@ -54,169 +51,41 @@ function LinearOperators.vcat(A::AbstractLinearOperator, n::Int) return op end -function diagOpProd(y::AbstractVector{T}, x::AbstractVector{T}, nrow::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T +function LinearOperatorCollection.diagOpProd(y::Vector{T}, x::Vector{T}, nrow::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T @floop for i=1:length(ops) mul!(view(y,yIdx[i]:yIdx[i+1]-1), ops[i], view(x,xIdx[i]:xIdx[i+1]-1)) end return y end -function diagOpTProd(y::AbstractVector{T}, x::AbstractVector{T}, ncol::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T +function LinearOperatorCollection.diagOpTProd(y::Vector{T}, x::Vector{T}, ncol::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T @floop for i=1:length(ops) mul!(view(y,yIdx[i]:yIdx[i+1]-1), transpose(ops[i]), view(x,xIdx[i]:xIdx[i+1]-1)) end return y end -function diagOpCTProd(y::AbstractVector{T}, x::AbstractVector{T}, ncol::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T +function LinearOperatorCollection.diagOpCTProd(y::Vector{T}, x::Vector{T}, ncol::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T @floop for i=1:length(ops) mul!(view(y,yIdx[i]:yIdx[i+1]-1), adjoint(ops[i]), view(x,xIdx[i]:xIdx[i+1]-1)) end return y end -function Base.copy(S::LinearOperator{T}) where T - deepcopy(S) -end - -mutable struct DiagOp{T} <: AbstractLinearOperator{T} - nrow :: Int - ncol :: Int - symmetric :: Bool - hermitian :: Bool - prod! :: Function - tprod! :: Function - ctprod! :: Function - nprod :: Int - ntprod :: Int - nctprod :: Int - args5 :: Bool - use_prod5! :: Bool - allocated5 :: Bool - Mv5 :: Vector{T} - Mtu5 :: Vector{T} - ops - equalOps :: Bool - xIdx :: Vector{Int} - yIdx :: Vector{Int} +function LinearOperatorCollection.diagNormOpProd!(y::Vector{T}, normalOps, idx, x::Vector{T}) where T + @floop for i=1:length(normalOps) + mul!(view(y,idx[i]:idx[i+1]-1), normalOps[i], view(x,idx[i]:idx[i+1]-1)) + end + return y end - -LinearOperators.storage_type(op::DiagOp) = typeof(op.Mv5) - - -""" - DiagOp(ops::AbstractLinearOperator...) - DiagOp(ops::Vector{AbstractLinearOperator}) - DiagOp(ops::NTuple{N,AbstractLinearOperator}) - -create a bloc-diagonal operator out of the `LinearOperator`s contained in ops -""" -DiagOp(ops::AbstractLinearOperator...) = DiagOp(ops) - -function DiagOp(ops) - nrow = 0 - ncol = 0 - S = eltype(ops[1]) - for i = 1:length(ops) - nrow += ops[i].nrow - ncol += ops[i].ncol - S = promote_type(S, eltype(ops[i])) - end - - xIdx = cumsum(vcat(1,[ops[i].ncol for i=1:length(ops)])) - yIdx = cumsum(vcat(1,[ops[i].nrow for i=1:length(ops)])) - - Op = DiagOp{S}( nrow, ncol, false, false, - (res,x) -> (diagOpProd(res,x,nrow,xIdx,yIdx,ops...)), - (res,y) -> (diagOpTProd(res,y,ncol,yIdx,xIdx,ops...)), - (res,y) -> (diagOpCTProd(res,y,ncol,yIdx,xIdx,ops...)), - 0, 0, 0, false, false, false, S[], S[], - [ops...], false, xIdx, yIdx) - - return Op -end - -function DiagOp(op::AbstractLinearOperator, N=1) - nrow = N*op.nrow - ncol = N*op.ncol - S = eltype(op) - ops = [copy(op) for n=1:N] - - xIdx = cumsum(vcat(1,[ops[i].ncol for i=1:length(ops)])) - yIdx = cumsum(vcat(1,[ops[i].nrow for i=1:length(ops)])) - - Op = DiagOp{S}( nrow, ncol, false, false, - (res,x) -> (diagOpProd(res,x,nrow,xIdx,yIdx,ops...)), - (res,y) -> (diagOpTProd(res,y,ncol,yIdx,xIdx,ops...)), - (res,y) -> (diagOpCTProd(res,y,ncol,yIdx,xIdx,ops...)), - 0, 0, 0, false, false, false, S[], S[], - ops, true, xIdx, yIdx ) - - return Op -end - -### Normal Matrix Code ### - -mutable struct DiagNormalOp{T,V} <: AbstractLinearOperator{T} - nrow :: Int - ncol :: Int - symmetric :: Bool - hermitian :: Bool - prod! :: Function - tprod! :: Nothing - ctprod! :: Nothing - nprod :: Int - ntprod :: Int - nctprod :: Int - args5 :: Bool - use_prod5! :: Bool - allocated5 :: Bool - Mv5 :: Vector{T} - Mtu5 :: Vector{T} - normalOps::V - idx::Vector{Int64} - y::Vector{T} -end - -LinearOperators.storage_type(op::DiagNormalOp) = typeof(op.Mv5) - -function DiagNormalOp(normalOps, N, idx, y::Vector{T}) where {T} - - function produ!(y, normalOps, idx, x) - @floop for i=1:length(normalOps) - mul!(view(y,idx[i]:idx[i+1]-1), normalOps[i], view(x,idx[i]:idx[i+1]-1)) - end - return y - end - - return DiagNormalOp(N, N, false, false - , (res,x) -> produ!(res, normalOps, idx, x) - , nothing - , nothing - , 0, 0, 0, false, false, false, T[], T[] - , normalOps, idx, y) +function Base.copy(S::LinearOperator{T}) where T + deepcopy(S) end -function LinearOperatorCollection.normalOperator(S::DiagOp, W=opEye(eltype(S), size(S,1))) - weights = W*ones(S.nrow) - - T = promote_type(eltype(S), eltype(W)) - - if S.equalOps - # this optimization is only allowed if all ops are the same - - # we promote the weights to be of the same type as T, which will be required - # when creating the temporary vector in normalOperator in a later stage - opInner = normalOperator(S.ops[1], WeightingOp(T; weights=T.(weights[S.yIdx[1]:S.yIdx[2]-1].^2))) - op = DiagNormalOp([copy(opInner) for i=1:length(S.ops)], size(S,2), S.xIdx, zeros(T, S.ncol) ) - else - op = DiagNormalOp([normalOperator(S.ops[i], WeightingOp(T; weights=T.(weights[S.yIdx[i]:S.yIdx[i+1]-1].^2))) - for i in 1:length(S.ops)], size(S,2), S.xIdx, zeros(T, S.ncol) ) - end - - return op -end +fftParams(::Type{<:AbstractArray}) = (;:flags => FFTW.MEASURE) +# https://github.com/JuliaLang/julia/issues/35543 +stripParameters(arrT::Type{<:AbstractArray}) = Base.typename(arrT).wrapper end # module \ No newline at end of file diff --git a/MRIOperators/src/MapSliceOp.jl b/MRIOperators/src/MapSliceOp.jl index 8dcfcfba..25772747 100644 --- a/MRIOperators/src/MapSliceOp.jl +++ b/MRIOperators/src/MapSliceOp.jl @@ -1,10 +1,10 @@ export MapSliceOp -function mapSliceForeward(A, x::Vector, size1::Tuple, dim::Int) +function mapSliceForeward(A, x::AbstractVector, size1::Tuple, dim::Int) return vec( mapslices(u->A*u, reshape(x, size1), dims=dim) ) end -function mapSliceBackward(A, y::Vector, size2::Tuple, dim::Int) +function mapSliceBackward(A, y::AbstractVector, size2::Tuple, dim::Int) return vec( mapslices(v->adjoint(A)*v, reshape(y, size2), dims=dim) ) end @@ -19,9 +19,10 @@ generates an operator that applies `trafo` to each slice of dimension `dim` of a * `size2` - size of the resulting Array of applying trafo to all slices * (`T=ComplexF64`) - type of the transformation """ -function MapSliceOp(trafo, dim::Int64, size1::Tuple, size2::Tuple; T=ComplexF64) +function MapSliceOp(trafo, dim::Int64, size1::Tuple, size2::Tuple; T=ComplexF64, S = LinearOperators.storage_type(trafo)) return LinearOperator(prod(size2), prod(size1), false, false , (res,x) -> ( res .= mapSliceForeward(trafo, x, size1, dim) ) , nothing - , (res,y) -> ( res .= mapSliceBackward(trafo, y, size2, dim) ) ) + , (res,y) -> ( res .= mapSliceBackward(trafo, y, size2, dim) ) + , S = S) end diff --git a/MRIOperators/src/NFFTOp.jl b/MRIOperators/src/NFFTOp.jl deleted file mode 100644 index 30609d31..00000000 --- a/MRIOperators/src/NFFTOp.jl +++ /dev/null @@ -1,4 +0,0 @@ - -function LinearOperatorCollection.NFFTOp( tr::Trajectory{T}; shape::Tuple, toeplitz=false, oversamplingFactor=1.25, kernelSize=3, kargs...) where T - return LinearOperatorCollection.NFFTOp( Complex{T}; shape, nodes=kspaceNodes(tr), toeplitz=toeplitz, oversamplingFactor=oversamplingFactor, kernelSize=kernelSize, kargs...) -end \ No newline at end of file diff --git a/MRIOperators/src/SensitivityOp.jl b/MRIOperators/src/SensitivityOp.jl index 651ca908..21224c63 100644 --- a/MRIOperators/src/SensitivityOp.jl +++ b/MRIOperators/src/SensitivityOp.jl @@ -41,7 +41,8 @@ function SensitivityOp(sensMaps::AbstractMatrix{T}, numContr=1) where T return LinearOperator{T}(numVox*numContr*numChan, numVox*numContr, false, false, (res,x) -> prod_smap!(res,sensMaps,x,numVox,numChan,numContr), nothing, - (res,x) -> ctprod_smap!(res,sensMapsC,x,numVox,numChan,numContr)) + (res,x) -> ctprod_smap!(res,sensMapsC,x,numVox,numChan,numContr), + S = typeof(similar(sensMaps, 0))) end """ diff --git a/MRIOperators/src/Shutter.jl b/MRIOperators/src/Shutter.jl index 42ed977f..2fdd27bc 100644 --- a/MRIOperators/src/Shutter.jl +++ b/MRIOperators/src/Shutter.jl @@ -1,8 +1,8 @@ export circularShutter!, circularShutterFreq!, getCircularMask -preserveType(orig::Array{T}, data::Array{U}) where {T<:Real, U<:Complex} = real(data) +preserveType(orig::AbstractArray{T}, data::AbstractArray{U}) where {T<:Real, U<:Complex} = real(data) -preserveType(orig::Array{T}, data::Array{U}) where {T, U} = data +preserveType(orig::AbstractArray{T}, data::AbstractArray{U}) where {T, U} = data function circularShutter!(I::AbstractMatrix, radiusFactor::Number=1.0) imgSize = size(I) @@ -40,7 +40,7 @@ function getCircularMask(shape::Tuple, radiusFactor::Number=1.0) return A end -function circularShutterFreq!(I::Matrix, radiusFactor::T=1.0) where T<:Number +function circularShutterFreq!(I::AbstractMatrix, radiusFactor::T=1.0) where T<:Number imgSize = size(I) fftI = fftshift(fft(I)) center = (imgSize[1]/2.0, imgSize[2]/2.0) @@ -56,7 +56,7 @@ function circularShutterFreq!(I::Matrix, radiusFactor::T=1.0) where T<:Number return preserveType(I, ifft(ifftshift(fftI))) end -function circularShutterFreq!(I::Array{T,3}, radiusFactor::T=1.0) where T<:Number +function circularShutterFreq!(I::AbstractArray{T,3}, radiusFactor::T=1.0) where T<:Number imgSize = size(I) fftI = fftshift(fft(I)) center = (imgSize[1]/2.0, imgSize[2]/2.0, imgSize[3]/2.0) diff --git a/MRIOperators/src/SparseOp.jl b/MRIOperators/src/SparseOp.jl index 51a8327b..d3c04ac1 100644 --- a/MRIOperators/src/SparseOp.jl +++ b/MRIOperators/src/SparseOp.jl @@ -10,7 +10,7 @@ generates the sparsifying transform (`<: AbstractLinearOperator`) given its name * `shape::NTuple{D,Int64}` - size of the Array to be transformed * (`kargs`) - additional keyword arguments """ -function SparseOp(T::Type,name::AbstractString, shape::NTuple{N,Int64}; kargs...) where N +function SparseOp(T::Type,name::AbstractString, shape::NTuple{N,Int64}; S = Vector{T}, kargs...) where N params = Dict(kargs) if name=="Wavelet" # if get(params, :multiEcho, false) @@ -18,12 +18,12 @@ function SparseOp(T::Type,name::AbstractString, shape::NTuple{N,Int64}; kargs... # else # return WaveletOp(params[:shape]) # end - return WaveletOp(T; shape) + return WaveletOp(T; shape, S = S) elseif name=="nothing" - return opEye(T,prod(shape)) + return opEye(T,prod(shape), S = S) else error("SparseOp $name is not yet defined.") end - return opEye(T, prod(shape)) + return opEye(T, prod(shape), S = S) end diff --git a/MRIOperators/src/SubspaceOp.jl b/MRIOperators/src/SubspaceOp.jl index f3a843c6..4b7aefac 100644 --- a/MRIOperators/src/SubspaceOp.jl +++ b/MRIOperators/src/SubspaceOp.jl @@ -1,6 +1,6 @@ export SubspaceOp -function prod_subspace!(y::Vector{T}, basis::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numContr, numBasis) where T +function prod_subspace!(y::AbstractVector{T}, basis::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numContr, numBasis) where T x_ = reshape(x,numVox,numBasis) y_ = reshape(y,numVox,numContr) @@ -8,7 +8,7 @@ function prod_subspace!(y::Vector{T}, basis::AbstractMatrix{T}, x::AbstractVecto return y end -function ctprod_subspace!(y::Vector{T},basis::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numContr, numBasis) where T +function ctprod_subspace!(y::AbstractVector{T},basis::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numContr, numBasis) where T x_ = reshape(x,numVox,numContr) y_ = reshape(y,numVox,numBasis) @@ -24,14 +24,15 @@ basis correspond to svd_object.V cropped to a certain level basis type should correspond to the type of the rawdata """ -function SubspaceOp(basis::AbstractMatrix{T},shape::NTuple{D,Int64},numContr) where {T,D} +function SubspaceOp(basis::AbstractMatrix{T},shape::NTuple{D,Int64},numContr, S = LinearOperators.storage_type(basis)) where {T,D} numVox = prod(shape) numBasis = size(basis,2) return LinearOperator{T}(numVox*numContr, numVox*numBasis, false, false, (res,x) -> prod_subspace!(res,basis,x,numVox,numContr,numBasis), nothing, - (res,x) -> ctprod_subspace!(res,basis,x,numVox,numContr,numBasis)) + (res,x) -> ctprod_subspace!(res,basis,x,numVox,numContr,numBasis), + S = S) end diff --git a/MRIOperators/test/gpu/cuda.jl b/MRIOperators/test/gpu/cuda.jl new file mode 100644 index 00000000..e7d0eb2c --- /dev/null +++ b/MRIOperators/test/gpu/cuda.jl @@ -0,0 +1,5 @@ +using CUDA + +arrayTypes = [CuArray] + +include(joinpath(@__DIR__(), "..", "runtests.jl")) \ No newline at end of file diff --git a/MRIOperators/test/gpu/rocm.jl b/MRIOperators/test/gpu/rocm.jl new file mode 100644 index 00000000..ebf32fb3 --- /dev/null +++ b/MRIOperators/test/gpu/rocm.jl @@ -0,0 +1,5 @@ +using AMDGPU + +arrayTypes = [ROCArray] + +include(joinpath(@__DIR__(), "..", "runtests.jl")) \ No newline at end of file diff --git a/MRIOperators/test/runtests.jl b/MRIOperators/test/runtests.jl index 05a54187..19b78dbd 100644 --- a/MRIOperators/test/runtests.jl +++ b/MRIOperators/test/runtests.jl @@ -1,5 +1,11 @@ -using Test, MRIBase, MRIOperators, MRISimulation, NFFT.FFTW -using LinearAlgebra, LinearOperatorCollection +using Test, MRIBase, MRIOperators, MRISimulation, MRIOperators.NFFT.FFTW +using LinearAlgebra, MRIOperators.LinearOperatorCollection +using JLArrays -include("testOperators.jl") \ No newline at end of file +areTypesDefined = @isdefined arrayTypes +arrayTypes = areTypesDefined ? arrayTypes : [Array, JLArray] + +@testset "MRIOperators" begin + include("testOperators.jl") +end \ No newline at end of file diff --git a/MRIOperators/test/testOperators.jl b/MRIOperators/test/testOperators.jl index 57097273..d5961d7d 100644 --- a/MRIOperators/test/testOperators.jl +++ b/MRIOperators/test/testOperators.jl @@ -1,6 +1,6 @@ -# test FourierOperators -function testFT(N=16) +# test FourierOperators, almost duplicate with NFFT tests in LinearOperatorCollection, except for ExplicitOp usage +function testFT(N=16; arrayType = Array) # random image x = zeros(ComplexF64,N,N) for i=1:N,j=1:N @@ -13,19 +13,20 @@ function testFT(N=16) F_adj = F' # Operators + xop = arrayType(vec(x)) tr = CartesianTrajectory(Float64,N,N) - F_nfft = NFFTOp(tr; shape=(N,N), symmetrize=false) - F_exp = ExplicitOp((N,N),tr,zeros(ComplexF64,N,N),symmetrize=false) + F_nfft = NFFTOp(ComplexF64; nodes = kspaceNodes(tr), shape=(N,N), symmetrize=false, S = typeof(xop)) + F_exp = ExplicitOp((N,N),tr,zeros(ComplexF64,N,N),symmetrize=false, S = typeof(xop)) # test agains FourierOperators y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N)) ) y_adj = vec( ifftshift(reshape(F_adj*vec(fftshift(x)),N,N)) ) - y_nfft = F_nfft*vec(x) - y_adj_nfft = adjoint(F_nfft) * vec(x) + y_nfft = Array(F_nfft * xop) + y_adj_nfft = Array(adjoint(F_nfft) * xop) - y_exp = F_exp*vec(x) - y_adj_exp = adjoint(F_exp) * vec(x) + y_exp = Array(F_exp * xop) + y_adj_exp = Array(adjoint(F_exp) * xop) @test y ≈ y_nfft rtol = 1e-2 @test y ≈ y_exp rtol = 1e-2 @@ -35,31 +36,32 @@ function testFT(N=16) # test AHA w/o Toeplitz F_nfft.toeplitz = false AHA = LinearOperatorCollection.normalOperator(F_nfft) - y_AHA_nfft = AHA * vec(x) + y_AHA_nfft = Array(AHA * xop) y_AHA = F' * F * vec(x) @test y_AHA ≈ y_AHA_nfft rtol = 1e-2 # test AHA with Toeplitz F_nfft.toeplitz = true AHA = LinearOperatorCollection.normalOperator(F_nfft) - y_AHA_nfft = AHA * vec(x) - y_AHA_nfft = adjoint(F_nfft) * F_nfft * vec(x) + y_AHA_nfft_1 = Array(AHA * xop) + y_AHA_nfft_2 = Array(adjoint(F_nfft) * F_nfft * xop) y_AHA = F' * F * vec(x) - @test y_AHA ≈ y_AHA_nfft rtol = 1e-2 + @test y_AHA_nfft_1 ≈ y_AHA_nfft_2 rtol = 1e-2 + @test y_AHA ≈ y_AHA_nfft_1 rtol = 1e-2 # test type stability; # TODO: Ensure type stability for Trajectory objects and test here nodes = Float32.(tr.nodes) - F_nfft = NFFTOp(ComplexF32; shape=(N,N), nodes,symmetrize=false) + F_nfft = NFFTOp(ComplexF32; shape=(N,N), nodes,symmetrize=false, S = typeof(ComplexF32.(xop))) - y_nfft = F_nfft * vec(ComplexF32.(x)) - y_adj_nfft = adjoint(F_nfft) * vec(ComplexF32.(x)) + y_nfft = F_nfft * ComplexF32.(xop) + y_adj_nfft = adjoint(F_nfft) * ComplexF32.(xop) @test Complex{eltype(nodes)} === eltype(y_nfft) @test Complex{eltype(nodes)} === eltype(y_adj_nfft) end -function testFT3d(N=12) +function testFT3d(N=12; arrayType = Array) # random image x = zeros(ComplexF64,N,N,N) for i=1:N,j=1:N,k=1:N @@ -72,19 +74,20 @@ function testFT3d(N=12) F_adj = F' # Operators + xop = arrayType(vec(x)) tr = CartesianTrajectory3D(Float64,N,N,numSlices=N) - F_nfft = NFFTOp(tr; shape=(N,N,N), symmetrize=false) - F_exp = ExplicitOp((N,N,N),tr,zeros(ComplexF64,N,N,N),symmetrize=false) + F_nfft = NFFTOp(ComplexF64; nodes = kspaceNodes(tr), shape=(N,N,N), symmetrize=false, S = typeof(xop)) + F_exp = ExplicitOp((N,N,N),tr,zeros(ComplexF64,N,N,N),symmetrize=false, S = typeof(xop)) # test agains FourierOperators y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N,N)) ) y_adj = vec( ifftshift(reshape(F_adj*vec(fftshift(x)),N,N,N)) ) - y_nfft = F_nfft*vec(x) - y_adj_nfft = adjoint(F_nfft) * vec(x) + y_nfft = Array(F_nfft * xop) + y_adj_nfft = Array(adjoint(F_nfft) * xop) - y_exp = F_exp*vec(x) - y_adj_exp = adjoint(F_exp) * vec(x) + y_exp = Array(F_exp * xop) + y_adj_exp = Array(adjoint(F_exp) * xop) @test y ≈ y_exp rtol = 1e-2 @test y ≈ y_nfft rtol = 1e-2 @@ -94,38 +97,40 @@ function testFT3d(N=12) # test AHA w/o Toeplitz F_nfft.toeplitz = false AHA = LinearOperatorCollection.normalOperator(F_nfft) - y_AHA_nfft = AHA * vec(x) + y_AHA_nfft = Array(AHA * xop) y_AHA = F' * F * vec(x) @test y_AHA ≈ y_AHA_nfft rtol = 1e-2 # test AHA with Toeplitz F_nfft.toeplitz = true AHA = LinearOperatorCollection.normalOperator(F_nfft) - y_AHA_nfft = AHA * vec(x) - y_AHA_nfft = adjoint(F_nfft) * F_nfft * vec(x) + y_AHA_nfft_1 = Array(AHA * xop) + y_AHA_nfft_2 = Array(adjoint(F_nfft) * F_nfft * xop) y_AHA = F' * F * vec(x) - @test y_AHA ≈ y_AHA_nfft rtol = 1e-2 + @test y_AHA_nfft_1 ≈ y_AHA_nfft_2 rtol = 1e-2 + @test y_AHA ≈ y_AHA_nfft_1 rtol = 1e-2 end -function testUndersampledFourierOp(N=16) +function testUndersampledFourierOp(N=16; arrayType = Array) x = rand(ComplexF64,N,N) tr = CartesianTrajectory(Float64, N÷2, N) + xop = arrayType(vec(x)) # FourierOperator - F_ft = fourierEncodingOp((N,N), tr, "fast") + F_ft = fourierEncodingOp((N,N), tr, "fast"; S = typeof(xop)) # Explicit operator - F = ExplicitOp((N,N),tr,zeros(ComplexF64,N,N)) + F = ExplicitOp((N,N),tr,zeros(ComplexF64,N,N); S = typeof(xop)) - y1 = F_ft*vec(x) - y2 = F*vec(x) - x1 = adjoint(F_ft)*y1 - x2 = adjoint(F)*y2 + y1 = Array(F_ft*xop) + y2 = Array(F*xop) + x1 = Array(adjoint(F_ft)*arrayType(y1)) + x2 = Array(adjoint(F)*arrayType(y2)) @test y1 ≈ y2 rtol = 1e-2 @test x1 ≈ x2 rtol = 1e-2 end ## test FieldmapNFFTOp -function testFieldmapFT(N=16) +function testFieldmapFT(N=16; arrayType = Array) # random image x = zeros(ComplexF64,N,N) for i=1:N,j=1:N @@ -143,34 +148,36 @@ function testFieldmapFT(N=16) F_adj = F' # Operators - F_nfft = FieldmapNFFTOp((N,N),tr,cmap,symmetrize=false) + xop = arrayType(vec(x)) + F_nfft = FieldmapNFFTOp((N,N),tr,cmap,symmetrize=false, S = typeof(xop)) # test agains FourierOperators y = F*vec(x) y_adj = F_adj*vec(x) - y_nfft = F_nfft*vec(x) - y_adj_nfft = adjoint(F_nfft) * vec(x) + y_nfft = Array(F_nfft*xop) + y_adj_nfft = Array(adjoint(F_nfft) * xop) @test (norm(y-y_nfft)/norm(y)) < 1e-2 @test (norm(y_adj-y_adj_nfft)/norm(y_adj)) < 1e-2 end -function testSparseOp(T::Type,shape) - W = SparseOp(Complex{T},"Wavelet", shape) - +function testSparseOp(T::Type,shape; arrayType = Array) x = zeros(Complex{T},shape) for i=1:shape[1],j=1:shape[2],k=1:shape[3] x[i,j,k] = rand() end - x=vec(x) - xapprox = W' * W * x + xop = arrayType(vec(x)) + W = SparseOp(Complex{T},"Wavelet", shape; S = typeof(xop)) + + + xapprox = reshape(Array(W' * W * xop), size(x)...) @test (norm(xapprox-x)/norm(x)) < 1e-3 end ## test FieldmapNFFTOp -function testCopySizes(N=16) +function testCopySizes(N=16; arrayType = Array) # random image x = zeros(ComplexF64,N,N) for i=1:N,j=1:N @@ -186,9 +193,9 @@ function testCopySizes(N=16) idx = CartesianIndices((N,N))[collect(1:N^2)] # Operators - - F_nfft = NFFTOp(tr; shape=(N,N), symmetrize=false) - F_fmap_nfft = FieldmapNFFTOp((N,N),tr,cmap,symmetrize=false) + xop = arrayType(vec(xop)) + F_nfft = NFFTOp(tr; shape=(N,N), symmetrize=false, S = typeof(xop)) + F_fmap_nfft = FieldmapNFFTOp((N,N),tr,cmap,symmetrize=false, S = typeof(xop)) # Copy the FieldmapNFFTOp operator and change the plans field of the new operator to empty F_fmap_nfft_copy = copy(F_fmap_nfft) @@ -198,17 +205,19 @@ function testCopySizes(N=16) end -function testOperators() - @testset "Linear Operator" begin - testFT() - testFT3d() - testFieldmapFT() - testUndersampledFourierOp() - testSparseOp(Float32,(128,80,1)) - testSparseOp(Float64,(128,80,1)) - testSparseOp(Float32,(128,128,80)) - testSparseOp(Float64,(128,128,80)) +function testOperators(arrayType = Array) + @testset "MRI Linear Operator: $arrayType" begin + arrayType == JLArray || @testset "FT" testFT(;arrayType) + arrayType == JLArray || @testset "FT3d" testFT3d(;arrayType) + @testset "FieldmapFT" testFieldmapFT(;arrayType) + @testset "Undersampled Fourier Op" testUndersampledFourierOp(;arrayType) + @testset "Sparse Op F32 2D" testSparseOp(Float32,(128,80,1);arrayType) + @testset "Sparse Op F64 2D" testSparseOp(Float64,(128,80,1);arrayType) + @testset "Sparse Op F32 3D" testSparseOp(Float32,(128,128,80);arrayType) + @testset "Sparse Op F64 3D" testSparseOp(Float64,(128,128,80);arrayType) end end -testOperators() +for arrayType in arrayTypes + testOperators(arrayType) +end \ No newline at end of file diff --git a/MRISimulation/Project.toml b/MRISimulation/Project.toml index 52c5d833..29ab55de 100644 --- a/MRISimulation/Project.toml +++ b/MRISimulation/Project.toml @@ -1,7 +1,7 @@ name = "MRISimulation" uuid = "8988da37-ea20-4fa6-9af7-8a6f6f9a8970" author = ["Tobias Knopp "] -version = "0.1.2" +version = "0.1.3" [deps] Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -15,7 +15,7 @@ MRIOperators = "fb1137e3-90a6-46ce-a672-6e1e53d120f2" [compat] julia = "1.6" MRIBase = "0.3, 0.4" -MRIOperators = "0.2" +MRIOperators = "0.3" Reexport = "1" ImageUtils = "0.2.5" StatsBase = "0.33, 0.34" diff --git a/Project.toml b/Project.toml index 3c520aa6..b28e0503 100644 --- a/Project.toml +++ b/Project.toml @@ -16,17 +16,21 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" RegularizedLeastSquares = "1e9c538a-f78c-5de5-8ffb-0b6dbe892d23" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +[weakdeps] +GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" + [compat] AxisArrays = "0.4.6" FFTW = "1.0" FLoops = "0.2" +GPUArrays = "8, 9, 10" ImageUtils = "0.2.8" MRIBase = "0.3, 0.4" -MRIOperators = "0.2" +MRIOperators = "0.3" PrecompileTools = "1" ProgressMeter = "1.2" Reexport = "0.2, 1" -RegularizedLeastSquares = "0.11" +RegularizedLeastSquares = "0.16" Unitful = "1.2" julia = "1.6" @@ -40,3 +44,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "ImageUtils", "MRISimulation", "MRISampling", "Scratch", "FFTW"] + +[extensions] +MRIRecoGPUArraysExt = "GPUArrays" \ No newline at end of file diff --git a/ext/MRIRecoGPUArraysExt/MRIRecoGPUArraysExt.jl b/ext/MRIRecoGPUArraysExt/MRIRecoGPUArraysExt.jl new file mode 100644 index 00000000..5156297a --- /dev/null +++ b/ext/MRIRecoGPUArraysExt/MRIRecoGPUArraysExt.jl @@ -0,0 +1,8 @@ +module MRIRecoGPUArraysExt + +using MRIReco, GPUArrays, MRIReco.FLoops + +MRIReco.executor(::Type{<:AbstractGPUArray}) = SequentialEx() +MRIReco.copyOpsFn(::Type{<:AbstractGPUArray}) = identity + +end \ No newline at end of file diff --git a/src/Reconstruction/DirectReconstruction.jl b/src/Reconstruction/DirectReconstruction.jl index 6e6d8829..74460cb6 100644 --- a/src/Reconstruction/DirectReconstruction.jl +++ b/src/Reconstruction/DirectReconstruction.jl @@ -13,7 +13,9 @@ input: function reconstruction_direct(acqData::AcquisitionData{T} , reconSize::NTuple{D,Int64} , weights::Vector{Vector{Complex{T}}} - , correctionMap::Array{Complex{T}}=Complex{T}[]) where {D,T} + , correctionMap::Array{Complex{T}}=Complex{T}[] + , arrayType = Array + , S = Vector{Complex{T}}) where {D,T} encDims = ndims(trajectory(acqData)) if encDims!=D @@ -26,17 +28,17 @@ function reconstruction_direct(acqData::AcquisitionData{T} p = Progress(numSl*numChan*numContr*numRep, dt=1, desc="Direct Reconstruction...") for i = 1:numSl - F = encodingOps_simple(acqData, reconSize, slice=i, correctionMap=correctionMap) + F = encodingOps_simple(acqData, reconSize, slice=i, correctionMap=correctionMap, S = S) for k = 1:numContr for j = 1:numChan for l = 1:numRep - kdata = kData(acqData,k,j,i,rep=l) .* (weights[k].^2) + kdata = arrayType(kData(acqData,k,j,i,rep=l)) .* arrayType((weights[k].^2)) I = adjoint(F[k]) * kdata if isCircular( trajectory(acqData, k) ) circularShutter!(reshape(I, reconSize), 1.0) end - Ireco[:,i,k,j,l] = I + Ireco[:,i,k,j,l] = Array(I) next!(p) end diff --git a/src/Reconstruction/IterativeReconstruction.jl b/src/Reconstruction/IterativeReconstruction.jl index 84bc725c..4c137ddc 100644 --- a/src/Reconstruction/IterativeReconstruction.jl +++ b/src/Reconstruction/IterativeReconstruction.jl @@ -18,10 +18,11 @@ function reconstruction_simple( acqData::AcquisitionData{T} ; reconSize::NTuple{D,Int64} , reg::Vector{<:AbstractRegularization} , sparseTrafo - , weights::Vector{Vector{Complex{T}}} + , weights::Vector{vecTc} , solver::Type{<:AbstractLinearSolver} , encodingOps=nothing - , params...) where {D, T <: AbstractFloat} + , arrayType = Array + , params...) where {D, T <: AbstractFloat, vecTc <: AbstractVector{Complex{T}}} encDims = ndims(trajectory(acqData)) if encDims!=length(reconSize) @@ -30,7 +31,7 @@ function reconstruction_simple( acqData::AcquisitionData{T} numContr, numChan, numSl, numRep = numContrasts(acqData), numChannels(acqData), numSlices(acqData), numRepetitions(acqData) - encParams = getEncodingOperatorParams(;params...) + encParams = getEncodingOperatorParams(;arrayType, params...) # set sparse trafo in reg temp = [] @@ -44,9 +45,9 @@ function reconstruction_simple( acqData::AcquisitionData{T} end reg = identity.(temp) - # reconstruction + # reconstructiond Ireco = zeros(Complex{T}, prod(reconSize), numSl, numContr, numChan, numRep) - #@floop + #@floop executor(arrayType) for l = 1:numRep, k = 1:numSl if encodingOps!=nothing F = encodingOps[:,k] @@ -56,18 +57,17 @@ function reconstruction_simple( acqData::AcquisitionData{T} for j = 1:numContr W = WeightingOp(Complex{T}; weights=weights[j]) for i = 1:numChan - kdata = kData(acqData,j,i,k,rep=l).* weights[j] - EFull = ∘(W, F[j])#, isWeighting=true) - EFullᴴEFull = normalOperator(EFull) - solv = createLinearSolver(solver, EFull; AᴴA=EFullᴴEFull, reg=reg, params...) + kdata = arrayType(kData(acqData,j,i,k,rep=l)).* weights[j] + EFull = ProdOp(W, F[j])#, isWeighting=true) + EFullᴴEFull = normalOperator(EFull; normalOpParams(arrayType)...) + solv = createLinearSolver(solver, EFull; AHA=EFullᴴEFull, reg=reg, params...) - I = solve(solv, kdata, startVector=get(params,:startVector,Complex{T}[]), - solverInfo=get(params,:solverInfo,nothing)) + I = solve!(solv, kdata, x0=get(params,:startVector,0)) if isCircular( trajectory(acqData, j) ) circularShutter!(reshape(I, reconSize), 1.0) end - Ireco[:,k,j,i,l] = I + Ireco[:,k,j,i,l] = Array(I) end end end @@ -94,10 +94,11 @@ function reconstruction_multiEcho(acqData::AcquisitionData{T} ; reconSize::NTuple{D,Int64} , reg::Vector{<:AbstractRegularization} , sparseTrafo - , weights::Vector{Vector{Complex{T}}} + , weights::Vector{vecTc} , solver::Type{<:AbstractLinearSolver} , encodingOps=nothing - , params...) where {D , T <: AbstractFloat} + , arrayType = Array + , params...) where {D , T <: AbstractFloat, vecTc <: AbstractVector{Complex{T}}} encDims = ndims(trajectory(acqData)) if encDims!=length(reconSize) @@ -105,7 +106,7 @@ function reconstruction_multiEcho(acqData::AcquisitionData{T} end numContr, numChan, numSl, numRep = numContrasts(acqData), numChannels(acqData), numSlices(acqData), numRepetitions(acqData) - encParams = getEncodingOperatorParams(;params...) + encParams = getEncodingOperatorParams(;arrayType, params...) # set sparse trafo in reg temp = [] @@ -123,19 +124,19 @@ function reconstruction_multiEcho(acqData::AcquisitionData{T} # reconstruction Ireco = zeros(Complex{T}, prod(reconSize)*numContr, numChan, numSl, numRep) - @floop for l = 1:numRep, i = 1:numSl + @floop executor(arrayType) for l = 1:numRep, i = 1:numSl if encodingOps != nothing F = encodingOps[i] else F = encodingOp_multiEcho(acqData, reconSize, slice=i; encParams...) end for j = 1:numChan - kdata = multiEchoData(acqData, j, i,rep=l) .* vcat(weights...) + kdata = arrayType(multiEchoData(acqData, j, i,rep=l)) .* vcat(weights...) EFull = ∘(W, F)#, isWeighting=true) - EFullᴴEFull = normalOperator(EFull) - solv = createLinearSolver(solver, EFull; AᴴA=EFullᴴEFull, reg=reg, params...) + EFullᴴEFull = normalOperator(EFull; normalOpParams(arrayType)...) + solv = createLinearSolver(solver, EFull; AHA=EFullᴴEFull, reg=reg, params...) - Ireco[:,j,i,l] = solve(solv,kdata; params...) + Ireco[:,j,i,l] = Array(solve!(solv,kdata)) # TODO circular shutter end end @@ -164,7 +165,7 @@ are reconstructed independently. * `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData * `L_inv::Array{Complex{<:AbstractFloat}}` - noise decorrelation matrix * `solver::Type{<:AbstractLinearSolver}` - name of the solver to use -* `senseMaps::Array{Complex{<:AbstractFloat}}` - coil sensitivities +* `senseMaps::AbstractArray{Complex{<:AbstractFloat}}` - coil sensitivities * (`normalize::Bool=false`) - adjust regularization parameter according to the size of k-space data * (`params::Dict{Symbol,Any}`) - Dict with additional parameters """ @@ -172,12 +173,13 @@ function reconstruction_multiCoil(acqData::AcquisitionData{T} ; reconSize::NTuple{D,Int64} , reg::Vector{<:AbstractRegularization} , sparseTrafo - , weights::Vector{Vector{Complex{T}}} - , L_inv::Union{LowerTriangular{Complex{T}, Matrix{Complex{T}}}, Nothing} + , weights::Vector{vecTc} + , L_inv::Union{LowerTriangular{Complex{T}, <:AbstractMatrix{Complex{T}}}, Nothing} , solver::Type{<:AbstractLinearSolver} - , senseMaps::Array{Complex{T}} + , senseMaps::AbstractArray{Complex{T}} , encodingOps=nothing - , params...) where {D , T} + , arrayType = Array + , params...) where {D , T, vecTc <: AbstractVector{Complex{T}}} encDims = ndims(trajectory(acqData)) if encDims!=length(reconSize) @@ -185,7 +187,7 @@ function reconstruction_multiCoil(acqData::AcquisitionData{T} end numContr, numChan, numSl, numRep = numContrasts(acqData), numChannels(acqData), numSlices(acqData), numRepetitions(acqData) - encParams = getEncodingOperatorParams(;params...) + encParams = getEncodingOperatorParams(;arrayType, params...) # noise decorrelation senseMapsUnCorr = decorrelateSenseMaps(L_inv, senseMaps, numChan) @@ -206,7 +208,7 @@ function reconstruction_multiCoil(acqData::AcquisitionData{T} # solve optimization problem Ireco = zeros(Complex{T}, prod(reconSize), numSl, numContr, numRep) let reg = reg # Fix @floop warning due to conditional/multiple assignment to reg - @floop for l = 1:numRep, k = 1:numSl + @floop executor(arrayType) for l = 1:numRep, k = 1:numSl if encodingOps != nothing E = encodingOps[:,k] else @@ -215,20 +217,20 @@ function reconstruction_multiCoil(acqData::AcquisitionData{T} for j = 1:numContr W = WeightingOp(Complex{T}; weights=weights[j], rep=numChan) - kdata = multiCoilData(acqData, j, k, rep=l) .* repeat(weights[j], numChan) + kdata = arrayType((multiCoilData(acqData, j, k, rep=l))) .* repeat(weights[j], numChan) if !isnothing(L_inv) kdata = vec(reshape(kdata, :, numChan) * L_inv') end - EFull = ∘(W, E[j], isWeighting=true) - EFullᴴEFull = normalOperator(EFull) - solv = createLinearSolver(solver, EFull; AᴴA=EFullᴴEFull, reg=reg, params...) - I = solve(solv, kdata; params...) + EFull = ∘(W, E[j]) + EFullᴴEFull = normalOperator(EFull; normalOpParams(arrayType)...) + solv = createLinearSolver(solver, EFull; AHA=EFullᴴEFull, reg=reg, params...) + I = solve!(solv, kdata) if isCircular( trajectory(acqData, j) ) circularShutter!(reshape(I, reconSize), 1.0) end - Ireco[:,k,j,l] = I + Ireco[:,k,j,l] = Array(I) end end end @@ -249,7 +251,7 @@ Different slices are reconstructed independently. * `sparseTrafo::AbstractLinearOperator` - sparsifying transformation * `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData * `solver::Type{<:AbstractLinearSolver}` - name of the solver to use -* `senseMaps::Array{Complex{<:AbstractFloat}}` - coil sensitivities +* `senseMaps::AbstractArray{Complex{<:AbstractFloat}}` - coil sensitivities * (`normalize::Bool=false`) - adjust regularization parameter according to the size of k-space data * (`params::Dict{Symbol,Any}`) - Dict with additional parameters """ @@ -257,12 +259,13 @@ function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData{T} ; reconSize::NTuple{D,Int64} , reg::Vector{<:AbstractRegularization} , sparseTrafo - , weights::Vector{Vector{Complex{T}}} + , weights::Vector{vecTc} , L_inv::Union{LowerTriangular{Complex{T}, Matrix{Complex{T}}}, Nothing} , solver::Type{<:AbstractLinearSolver} - , senseMaps::Array{Complex{T}} + , senseMaps::AbstractArray{Complex{T}} , encodingOps=nothing - , params...) where {D, T} + , arrayType = Array + , params...) where {D, T, vecTc <: AbstractVector{Complex{T}}} encDims = ndims(trajectory(acqData)) if encDims!=length(reconSize) @@ -270,7 +273,7 @@ function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData{T} end numContr, numChan, numSl, numRep = numContrasts(acqData), numChannels(acqData), numSlices(acqData), numRepetitions(acqData) - encParams = getEncodingOperatorParams(;params...) + encParams = getEncodingOperatorParams(;arrayType, params...) # noise decorrelation senseMapsUnCorr = decorrelateSenseMaps(L_inv, senseMaps, numChan) @@ -290,23 +293,23 @@ function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData{T} W = WeightingOp(Complex{T}; weights=vcat(weights...), rep=numChan ) Ireco = zeros(Complex{T}, prod(reconSize)*numContr, numSl, numRep) - @floop for l = 1:numRep, i = 1:numSl + @floop executor(arrayType) for l = 1:numRep, i = 1:numSl if encodingOps != nothing E = encodingOps[i] else E = encodingOp_multiEcho_parallel(acqData, reconSize, senseMapsUnCorr; slice=i, encParams...) end - kdata = multiCoilMultiEchoData(acqData, i) .* repeat(vcat(weights...), numChan) + kdata = arrayType(multiCoilMultiEchoData(acqData, i)) .* repeat(vcat(weights...), numChan) if !isnothing(L_inv) - kdata = vec(reshape(kdata, :, numChan) * L_inv') + kdata = arrayType(vec(reshape(kdata, :, numChan) * L_inv')) end EFull = ∘(W, E)#, isWeighting=true) - EFullᴴEFull = normalOperator(EFull) - solv = createLinearSolver(solver, EFull; AᴴA=EFullᴴEFull, reg=reg, params...) + EFullᴴEFull = normalOperator(EFull; normalOpParams(arrayType)...) + solv = createLinearSolver(solver, EFull; AHA=EFullᴴEFull, reg=reg, params...) - Ireco[:,i,l] = solve(solv, kdata; params...) + Ireco[:,i,l] = Array(solve!(solv, kdata)) end @@ -335,7 +338,7 @@ Different slices are reconstructed independently. * `sparseTrafo::AbstractLinearOperator` - sparsifying transformation * `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData * `solver::Type{<:AbstractLinearSolver}` - name of the solver to use -* `senseMaps::Array{Complex{<:AbstractFloat}}` - coil sensitivities +* `senseMaps::AbstractArray{Complex{<:AbstractFloat}}` - coil sensitivities * (`normalize::Bool=false`) - adjust regularization parameter according to the size of k-space data * (`params::Dict{Symbol,Any}`) - Dict with additional parameters """ @@ -343,11 +346,12 @@ function reconstruction_multiCoilMultiEcho_subspace(acqData::AcquisitionData{T} ; reconSize::NTuple{D,Int64} , reg::Vector{<:AbstractRegularization} , sparseTrafo - , weights::Vector{Vector{Complex{T}}} + , weights::Vector{vecTc} , solver::Type{<:AbstractLinearSolver} - , senseMaps::Array{Complex{T}} + , senseMaps::AbstractArray{Complex{T}} , encodingOps=nothing - , params...) where {D, T} + , arrayType = Array + , params...) where {D, T, vecTc <: AbstractVector{Complex{T}}} encDims = ndims(trajectory(acqData)) if encDims!=length(reconSize) @@ -357,7 +361,7 @@ function reconstruction_multiCoilMultiEcho_subspace(acqData::AcquisitionData{T} numContr, numChan, numSl, numRep = numContrasts(acqData), numChannels(acqData), numSlices(acqData), numRepetitions(acqData) numBasis = size(params[:basis],2) - encParams = getEncodingOperatorParams(;params...) + encParams = getEncodingOperatorParams(;arrayType, params...) # set sparse trafo in reg temp = [] @@ -374,7 +378,7 @@ function reconstruction_multiCoilMultiEcho_subspace(acqData::AcquisitionData{T} W = WeightingOp(Complex{T}; weights=vcat(weights...), rep=numChan ) Ireco = zeros(Complex{T}, prod(reconSize)*numBasis, numSl, numRep) - @floop for l = 1:numRep, i = 1:numSl + @floop executor(arrayType) for l = 1:numRep, i = 1:numSl if encodingOps != nothing E = encodingOps[i] else @@ -383,13 +387,13 @@ function reconstruction_multiCoilMultiEcho_subspace(acqData::AcquisitionData{T} E = ∘(E,subOp) end - kdata = multiCoilMultiEchoData(acqData, i) .* repeat(vcat(weights...), numChan) + kdata = arrayType(multiCoilMultiEchoData(acqData, i)) .* repeat(vcat(weights...), numChan) EFull = ∘(W, E)#, isWeighting=true) - EFullᴴEFull = normalOperator(EFull) - solv = createLinearSolver(solver, EFull; AᴴA=EFullᴴEFull, reg=reg, params...) + EFullᴴEFull = normalOperator(EFull; normalOpParams(arrayType)...) + solv = createLinearSolver(solver, EFull; AHA=EFullᴴEFull, reg=reg, params...) - Ireco[:,i,l] = solve(solv, kdata; params...) + Ireco[:,i,l] = Array(solve!(solv, kdata)) end diff --git a/src/Reconstruction/RecoParameters.jl b/src/Reconstruction/RecoParameters.jl index c0e6f3b4..e2df65a0 100644 --- a/src/Reconstruction/RecoParameters.jl +++ b/src/Reconstruction/RecoParameters.jl @@ -3,12 +3,13 @@ export defaultRecoParams function defaultRecoParams() params = Dict{Symbol,Any}() params[:reco] = "direct" - params[:sparseTrafoName] = "Wavelet" params[:reg] = L1Regularization(0.0) params[:normalizeReg] = NoNormalization() params[:solver] = ADMM - params[:ρ] = 5.e-2 + params[:rho] = 5.e-2 params[:iterations] = 30 + params[:arrayType] = Array + params[:kwargWarning] = false return params end @@ -19,7 +20,10 @@ function setupDirectReco(acqData::AcquisitionData{T}, recoParams::Dict) where T # field map cmap = get(recoParams, :cmap, Complex{T}[]) - return reconSize, weights, cmap + arrayType = get(recoParams, :arrayType, Array) + S = typeof(arrayType{Complex{T}}(undef, 0)) + + return reconSize, weights, cmap, arrayType, S end @@ -80,22 +84,26 @@ function setupIterativeReco!(acqData::AcquisitionData{T}, recoParams::Dict) wher reconSize = recoParams[:reconSize] end + # Array type + arrayType = get(recoParams, :arrayType, Array) + S = typeof(arrayType{Complex{T}}(undef, 0)) + # density weights densityWeighting = get(recoParams,:densityWeighting,true) if densityWeighting - weights = samplingDensity(acqData,reconSize) + weights = map(weight -> S(weight), samplingDensity(acqData,reconSize)) else numContr = numContrasts(acqData) - weights = Array{Vector{Complex{T}}}(undef,numContr) + weights = Array{S}(undef,numContr) for contr=1:numContr numNodes = size(acqData.kdata[contr],1) - weights[contr] = [1.0/sqrt(prod(reconSize)) for node=1:numNodes] + weights[contr] = S([1.0/sqrt(prod(reconSize)) for node=1:numNodes]) end end # bare regularization (without sparsifying transform) reg = get(recoParams,:reg,L1Regularization(zero(T))) - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] # sparsifying transform if haskey(recoParams, :sparseTrafo) @@ -105,7 +113,7 @@ function setupIterativeReco!(acqData::AcquisitionData{T}, recoParams::Dict) wher end # Construct SparseOp for each string instance - sparseTrafos = map(x-> x isa String ? SparseOp(Complex{T}, x, reconSize; recoParams...) : x, sparseTrafos) + sparseTrafos = map(x-> x isa String ? SparseOp(Complex{T}, x, reconSize; S = S, recoParams...) : x, sparseTrafos) # Fill up SparseOps for remaining reg terms with nothing temp = Union{Nothing, eltype(sparseTrafos)}[nothing for i = 1:length(reg)] @@ -127,40 +135,54 @@ function setupIterativeReco!(acqData::AcquisitionData{T}, recoParams::Dict) wher solver = get(recoParams, :solver, FISTA) # sensitivity maps - senseMaps = get(recoParams, :senseMaps, Complex{T}[]) + senseMaps = arrayType(get(recoParams, :senseMaps, S())) if red3d && !isempty(senseMaps) # make sure the dimensions match the trajectory dimensions senseMaps = permutedims(senseMaps,[2,3,1,4]) end # noise data acquisition [samples, coils] - noiseData = get(recoParams, :noiseData, Complex{T}[]) + noiseData = arrayType(get(recoParams, :noiseData, S())) if isempty(noiseData) L_inv = nothing else psi = convert(typeof(noiseData), covariance(noiseData)) - L = cholesky(psi, check = true) + # Check if approx. hermitian and avoid check in cholesky + # GPU arrays can have float rounding differences + @assert isapprox(adjoint(psi), psi) + L = cholesky(psi; check = false) L_inv = inv(L.L) #noise decorrelation matrix end - recoParams[:reconSize] = reconSize recoParams[:weights] = weights recoParams[:L_inv] = L_inv recoParams[:sparseTrafo] = sparseTrafo recoParams[:reg] = reg - recoParams[:normalize] = normalize - recoParams[:encOps] = encOps + recoParams[:normalizeReg] = normalize + recoParams[:encodingOps] = encOps recoParams[:solver] = solver recoParams[:senseMaps] = senseMaps + recoParams[:arrayType] = arrayType + recoParams[:S] = S return recoParams end -function getEncodingOperatorParams(; kargs...) - encKeys = [:correctionMap, :method, :toeplitz, :oversamplingFactor, :kernelSize, :K, :K_tol] - return Dict([key=>kargs[key] for key in intersect(keys(kargs),encKeys)]) +function getEncodingOperatorParams(; arrayType = Array, kargs...) + encKeys = [:correctionMap, :method, :toeplitz, :oversamplingFactor, :kernelSize, :K, :K_tol, :S] + dict = Dict{Symbol, Any}([key=>kargs[key] for key in intersect(keys(kargs),encKeys)]) + push!(dict, :copyOpsFn => copyOpsFn(arrayType)) + return dict end # convenience methods volumeSize(reconSize::NTuple{2,Int}, numSlice::Int) = (reconSize..., numSlice) volumeSize(reconSize::NTuple{3,Int}, numSlice::Int) = reconSize + +executor(::Type{<:AbstractArray}) = nothing +copyOpsFn(::Type{<:AbstractArray}) = copy +normalOpParams(::Type{aT}) where aT <: AbstractArray = (; :copyOpsFn => copyOpsFn(aT), MRIOperators.fftParams(aT)...) + +executor(f::Function) = executor(f{Complex{Float32}}(undef, 0)) +copyOpsFn(f::Function) = copyOpsFn(f{Complex{Float32}}(undef, 0)) +normalOpParams(f::Function) = normalOpParams(f{Complex{Float32}}(undef, 0)) \ No newline at end of file diff --git a/src/Reconstruction/Reconstruction.jl b/src/Reconstruction/Reconstruction.jl index fdae4495..111e4b8b 100644 --- a/src/Reconstruction/Reconstruction.jl +++ b/src/Reconstruction/Reconstruction.jl @@ -36,8 +36,8 @@ function reconstruction(acqData::AcquisitionData, recoParams::Dict) # iterative reco if recoParams[:reco] == "direct" - reconSize, weights, cmap = setupDirectReco(acqData, recoParams) - return reconstruction_direct(acqData, reconSize[1:encodingDims], weights, cmap) + reconSize, weights, cmap, arrayType, S = setupDirectReco(acqData, recoParams) + return reconstruction_direct(acqData, reconSize[1:encodingDims], weights, cmap, arrayType, S) else setupIterativeReco!(acqData, recoParams) recoParams[:reconSize] = recoParams[:reconSize][1:encodingDims] diff --git a/src/Tools/NoiseDeCorr.jl b/src/Tools/NoiseDeCorr.jl index ce6c362c..72907f55 100644 --- a/src/Tools/NoiseDeCorr.jl +++ b/src/Tools/NoiseDeCorr.jl @@ -6,8 +6,8 @@ export decorrelateSenseMaps, covariance multiplies the senseMaps by the noise uncorrelation matrix (L_inv). """ -function decorrelateSenseMaps(L_inv::Union{LowerTriangular{Complex{T}, Matrix{Complex{T}}}, Nothing}, - senseMaps::Array{Complex{T},4}, +function decorrelateSenseMaps(L_inv::Union{LowerTriangular{Complex{T}, <:AbstractMatrix{Complex{T}}}, Nothing}, + senseMaps::AbstractArray{Complex{T},4}, numChan::Int64) where {T} if isnothing(L_inv) senseMapsUnCorr = senseMaps @@ -27,7 +27,7 @@ The `noiseData` array should be number of noise samples by number of coils. computes the covariance of the noise acquisition. """ -function covariance(noiseData::Array{Complex{T},2}) where {T} +function covariance(noiseData::AbstractArray{Complex{T},2}) where {T} N = size(noiseData, 1) cov = (1/(N-1)) .* (noiseData' * noiseData) diff --git a/src/Tools/Regridding.jl b/src/Tools/Regridding.jl index d3304485..4d8032cd 100644 --- a/src/Tools/Regridding.jl +++ b/src/Tools/Regridding.jl @@ -28,7 +28,7 @@ function regrid(acqData::AcquisitionData{T,2}, kspaceSize::NTuple{2,Int64}; solver = RegularizedLeastSquares.CGNR(W*E[j], iterations=cgnr_iter) for i = 1:numChan kdata = kData(acqData,j,i,k).* dcf[j] - img .= solve(solver, kdata) + img .= solve!(solver, kdata) kdata_cart[j,k,1][:,i] .= F*img end end diff --git a/test/gpu/cuda.jl b/test/gpu/cuda.jl new file mode 100644 index 00000000..e7d0eb2c --- /dev/null +++ b/test/gpu/cuda.jl @@ -0,0 +1,5 @@ +using CUDA + +arrayTypes = [CuArray] + +include(joinpath(@__DIR__(), "..", "runtests.jl")) \ No newline at end of file diff --git a/test/gpu/rocm.jl b/test/gpu/rocm.jl new file mode 100644 index 00000000..ebf32fb3 --- /dev/null +++ b/test/gpu/rocm.jl @@ -0,0 +1,5 @@ +using AMDGPU + +arrayTypes = [ROCArray] + +include(joinpath(@__DIR__(), "..", "runtests.jl")) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 742122ce..bb904890 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,10 @@ using Scratch const tmpdir = @get_scratch!("tmp") @info "If you want to check the output of the tests, please head to $tmpdir." -include("testIO.jl") -include("testReconstruction.jl") -include("testSpecificApplications.jl") +areTypesDefined = @isdefined arrayTypes +arrayTypes = areTypesDefined ? arrayTypes : [Array] # , JLArray] +@testset "MRIReco" begin + include("testIO.jl") + include("testReconstruction.jl") + include("testSpecificApplications.jl") +end \ No newline at end of file diff --git a/test/testReconstruction.jl b/test/testReconstruction.jl index 7a63474c..43f94fe6 100644 --- a/test/testReconstruction.jl +++ b/test/testReconstruction.jl @@ -1,5 +1,5 @@ # test gridding reco -function testGriddingReco(N=32) +function testGriddingReco(N=32;arrayType = Array) # image x = shepp_logan(N) @@ -15,13 +15,14 @@ function testGriddingReco(N=32) #reco params[:reco] = "direct" params[:reconSize] = (N,N) + params[:arrayType] = arrayType x_approx = reconstruction(acqData, params) @test (norm(vec(x)-vec(x_approx))/norm(vec(x))) < 1e-2 end # test gridding reco -function testConvertKspace(N=32) +function testConvertKspace(N=32;arrayType = Array) # image x = shepp_logan(N) @@ -37,6 +38,7 @@ function testConvertKspace(N=32) #reco params[:reco] = "direct" params[:reconSize] = (N,N) + params[:arrayType] = arrayType x_approx = reconstruction(acqData, params) kspace = kDataCart(acqData) @@ -57,6 +59,7 @@ function testConvertKspace(N=32) acqCS = AcquisitionData(kspace_cs) params[:reco] = "direct" params[:reconSize] = (N,N) + params[:arrayType] = arrayType x_cs = reconstruction(acqCS, params) x_cs2 = ifftshift(ifft(ifftshift(kspace_cs))) @@ -64,7 +67,7 @@ function testConvertKspace(N=32) @test (norm(vec(abs.(x_cs))-vec(abs.(x_cs2[:,:,1,1,1,1])))/norm(vec(abs.(x_cs2)))) < 1e-10 end -function testConvertKspace3D(N=32) +function testConvertKspace3D(N=32;arrayType = Array) # image x = shepp_logan(N) x = repeat(x,1,1,N) @@ -82,6 +85,7 @@ function testConvertKspace3D(N=32) #reco params[:reco] = "direct" params[:reconSize] = (N,N,N) + params[:arrayType] = arrayType x_approx = reconstruction(acqData, params) kspace = kDataCart(acqData) @@ -110,7 +114,7 @@ function testConvertKspace3D(N=32) @test (norm(vec(abs.(x_cs))-vec(abs.(x_cs2[:,:,:,1,1,1])))/norm(vec(abs.(x_cs2)))) < 1e-10 end -function testGriddingReco3d(N=32) +function testGriddingReco3d(N=32;arrayType = Array) sh = ComplexF64.(shepp_logan(N)) x = cat(sh,0.9*sh,0.8*sh,0.7*sh,0.6*sh,0.5*sh,0.4*sh,0.3*sh,dims=3) @@ -127,13 +131,14 @@ function testGriddingReco3d(N=32) # reco params[:reco] = "direct" params[:reconSize] = (N,N,8) + params[:arrayType] = arrayType x_approx = vec(reconstruction(acqData, params)) @test (norm(vec(x)-vec(x_approx))/norm(vec(x))) < 1e-2 end # test CS reco -function testCSReco(N=32,redFac=1.1;sampling="poisson") +function testCSReco(N=32,redFac=1.1;sampling="poisson",arrayType = Array) # image x = shepp_logan(N) @@ -156,13 +161,14 @@ function testCSReco(N=32,redFac=1.1;sampling="poisson") params[:reg] = L1Regularization(1.e-3) # regularization params[:solver] = ADMM # solver params[:iterations] = 1000 - params[:ρ] = 1.0e-1 + params[:rho] = 1.0e-1 + params[:arrayType] = arrayType x_approx = reconstruction(acqData, params) @test (norm(vec(x)-vec(x_approx))/norm(vec(x))) < 1e-1 end -function testCSRecoMultCoil(N=32;type = ComplexF64) +function testCSRecoMultCoil(N=32;type = ComplexF64,arrayType = Array) # image x = shepp_logan(N) smaps = birdcageSensitivity(32,2,3.0) @@ -190,9 +196,10 @@ function testCSRecoMultCoil(N=32;type = ComplexF64) params[:reg] = TVRegularization(2.e-3, shape = (N,N)) # regularization params[:solver] = ADMM # solver params[:iterations] = 100 - params[:ρ] = 1.0e-1 + params[:rho] = 1.0e-1 params[:absTol] = 1.e-5 params[:relTol] = 1.e-4 + params[:arrayType] = arrayType x_approx = reshape( reconstruction(acqData, params), 32,32,2) @@ -206,7 +213,7 @@ function testCSRecoMultCoil(N=32;type = ComplexF64) end # test CSSense Reco -function testCSSenseReco(N=32,redFac=1.1) +function testCSSenseReco(N=32,redFac=1.1;arrayType = Array) # image x = shepp_logan(N) @@ -235,15 +242,16 @@ function testCSSenseReco(N=32,redFac=1.1) params[:reg] = L1Regularization(1.e-3) # regularization params[:solver] = ADMM params[:iterations] = 1000 - params[:ρ] = 1.0e-1 + params[:rho] = 1.0e-1 params[:absTol] = 1.e-5 params[:relTol] = 1.e-4 + params[:arrayType] = arrayType x_approx = vec(reconstruction(acqData, params)) @test (norm(vec(x)-x_approx)/norm(vec(x))) < 1e-1 end -function testCSRecoMultiCoilCGNR(;N=32,redFac=1.1,type = ComplexF32) +function testCSRecoMultiCoilCGNR(;N=32,redFac=1.1,type = ComplexF32,arrayType = Array) x = shepp_logan(N) # coil sensitivites @@ -274,9 +282,10 @@ function testCSRecoMultiCoilCGNR(;N=32,redFac=1.1,type = ComplexF32) params[:reg] = L2Regularization(1.e-3) # regularization params[:solver] = CGNR params[:iterations] = 1 - params[:ρ] = 1.0e-1 + params[:rho] = 1.0e-1 params[:absTol] = 1.e-5 params[:relTol] = 1.e-4 + params[:arrayType] = arrayType x_approx = vec(reconstruction(acqData,params)) # This test is expected to have a higher error since CGNR with CS will not work @@ -284,7 +293,7 @@ function testCSRecoMultiCoilCGNR(;N=32,redFac=1.1,type = ComplexF32) @test (norm(vec(x)-x_approx)/norm(vec(x))) < 2e-1 end -function testOffresonanceReco(N = 128; accelMethod="nfft") +function testOffresonanceReco(N = 128; accelMethod="nfft",arrayType = Array) I = shepp_logan(N) I = circularShutterFreq!(I,1) @@ -312,13 +321,14 @@ function testOffresonanceReco(N = 128; accelMethod="nfft") params[:reconSize] = (N,N) params[:correctionMap] = cmap params[:method] = accelMethod + params[:arrayType] = arrayType Ireco = reconstruction(acqData, params) @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 1e-1 end -function testSENSEReco(N = 64, T = ComplexF64) +function testSENSEReco(N = 64, T = ComplexF64;arrayType = Array) numCoils = 8 I = T.(shepp_logan(N)) I = circularShutterFreq!(I,1) @@ -346,13 +356,14 @@ function testSENSEReco(N = 64, T = ComplexF64) params[:iterations] = 50 params[:solver] = CGNR params[:senseMaps] = coilsens + params[:arrayType] = arrayType Ireco = reconstruction(acqData, params) @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 1e-1 end -function testSENSEnoiseUnCorr(N = 64, T = ComplexF64) +function testSENSEnoiseUnCorr(N = 64, T = ComplexF64;arrayType = Array) numCoils = 8 R = 4 Img = T.(shepp_logan(N)) @@ -386,6 +397,7 @@ function testSENSEnoiseUnCorr(N = 64, T = ComplexF64) params[:reg] = L2Regularization(0.0) params[:iterations] = 150 params[:solver] = CGNR + params[:arrayType] = arrayType Ireco = reconstruction(acqData, params) @@ -395,7 +407,7 @@ function testSENSEnoiseUnCorr(N = 64, T = ComplexF64) @test (norm(vec(Img)-vec(Ireco))/norm(vec(Img)-vec(IrecoUnCorr))) > 1 end -function testOffresonanceSENSEReco(N = 64) +function testOffresonanceSENSEReco(N = 64;arrayType = Array) numCoils = 8 I = shepp_logan(N) I = circularShutterFreq!(I,1) @@ -426,13 +438,14 @@ function testOffresonanceSENSEReco(N = 64) params[:solver] = ADMM params[:senseMaps] = coilsens params[:correctionMap] = cmap + params[:arrayType] = arrayType Ireco = reconstruction(acqData, params) @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 1.6e-1 end -function testSenseMultiEcho(N=32,T = ComplexF32) +function testSenseMultiEcho(N=32,T = ComplexF32;arrayType = Array) # image x = T.(shepp_logan(N)) rmap = 20.0*ones(N,N) @@ -460,6 +473,7 @@ function testSenseMultiEcho(N=32,T = ComplexF32) params[:iterations] = 1 params[:solver] = CGNR params[:senseMaps] = reshape(coilsens,N,N,1,8) + params[:arrayType] = arrayType x_approx = reshape(reconstruction(acqData,params),N,N,:) @@ -478,7 +492,7 @@ function testSenseMultiEcho(N=32,T = ComplexF32) @test relErrorEcho1 < 1e-6 end -function testSenseMultiEchoMeasCoils(N=32;T = ComplexF32) +function testSenseMultiEchoMeasCoils(N=32;T = ComplexF32,arrayType = Array) # image x = T.(shepp_logan(N)) rmap = 20.0*ones(N,N) @@ -506,6 +520,7 @@ function testSenseMultiEchoMeasCoils(N=32;T = ComplexF32) params[:iterations] = 1 params[:solver] = CGNR params[:senseMaps] = reshape(coilsens,N,N,1,8) + params[:arrayType] = arrayType x_approx = reshape(reconstruction(acqData,params),N,N,:) @@ -517,6 +532,7 @@ function testSenseMultiEchoMeasCoils(N=32;T = ComplexF32) params[:iterations] = 1 params[:solver] = CGNR params[:senseMaps] = reshape(coilsens,N,N,1,8) + params[:arrayType] = arrayType x_approx2= reshape(reconstruction(acqData,params),N,N,:) @@ -524,7 +540,7 @@ function testSenseMultiEchoMeasCoils(N=32;T = ComplexF32) @test relErrorEcho1 < 1e-6 end -function testRecoMultiEcho(N=32) +function testRecoMultiEcho(N=32;arrayType = Array) # image x = ComplexF64.(shepp_logan(N)) rmap = 20.0*ones(N,N) @@ -559,6 +575,7 @@ function testRecoMultiEcho(N=32) params[:reg] = L2Regularization(1.e-3) params[:iterations] = 1 params[:solver] = CGNR + params[:arrayType] = arrayType x_approx2 = reshape(reconstruction(acqData,params),N,N,:) @@ -566,7 +583,7 @@ function testRecoMultiEcho(N=32) @test relErrorEcho3 < 1e-3 end -function testCSReco3d(N=128) +function testCSReco3d(N=128;arrayType = Array) sh = ComplexF64.(shepp_logan(N)) I = cat(sh,0.9*sh,0.8*sh,0.7*sh,0.6*sh,0.5*sh,0.4*sh,0.3*sh,dims=3) I = permutedims(I,[3,1,2]) @@ -590,16 +607,17 @@ function testCSReco3d(N=128) params[:reg] = L1Regularization(1.e-3) #"TV" # regularization params[:solver] = ADMM # solver params[:iterations] = 1000 - params[:ρ] = 1.0e-1 + params[:rho] = 1.0e-1 params[:absTol] = 1.e-4 params[:relTol] = 1.e-4 + params[:arrayType] = arrayType Ireco = reconstruction(acqData, params) @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 1e-1 end -function testCSSenseReco3d(N=128) +function testCSSenseReco3d(N=128;arrayType = Array) sh = ComplexF64.(shepp_logan(N)) I = cat(sh,0.9*sh,0.8*sh,0.7*sh,0.6*sh,0.5*sh,0.4*sh,0.3*sh,dims=3) I = permutedims(I,[3,1,2]) @@ -629,16 +647,17 @@ function testCSSenseReco3d(N=128) params[:reg] = L1Regularization(1.e-3) #"TV" # regularization params[:solver] = ADMM # solver params[:iterations] = 1000 - params[:ρ] = 1.0e-1 + params[:rho] = 1.0e-1 params[:absTol] = 1.e-5 params[:relTol] = 1.e-4 + params[:arrayType] = arrayType Ireco = reconstruction(acqData, params) @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 1e-1 end -function testRegridding(N=64) +function testRegridding(N=64;arrayType = Array) # image x = shepp_logan(N) @@ -665,6 +684,7 @@ function testRegridding(N=64) params[:solver] = CGNR params[:reg] = L2Regularization(0.0) params[:iterations] = 3 + params[:arrayType] = arrayType x_reg = collect(reshape(reconstruction(acqDataReg, params),N,N)) circularShutter!(x_reg) @@ -680,53 +700,55 @@ function testRegridding(N=64) @test (norm(vec(x_rad).-vec(x_reg))/norm(vec(x_rad))) < 1e-2 end -function testReco(N=32) - @testset "Reconstructions" begin +function testReco(N=32;arrayType = Array) + @testset "Reconstructions: $arrayType" begin @testset "MultiEcho" begin - testRecoMultiEcho() - testSenseMultiEcho(N) - testSenseMultiEchoMeasCoils(N) + @testset testRecoMultiEcho(;arrayType) + @testset testSenseMultiEcho(N;arrayType) + @testset testSenseMultiEchoMeasCoils(N;arrayType) end @testset "Convert k-space" begin - testConvertKspace() - testConvertKspace3D() + @testset "K-Space" testConvertKspace(;arrayType) + @testset "K-Space 3D" testConvertKspace3D(;arrayType) end @testset "Gridding" begin - testGriddingReco() - testGriddingReco3d() - testRegridding() + @testset "Gridding" testGriddingReco(;arrayType) + @testset "Gridding 3D" testGriddingReco3d(;arrayType) + @testset "Regridding" testRegridding(;arrayType) end @testset "CS" begin sampling = ["random", "poisson", "vdPoisson"] # "lines" - for samp in sampling - testCSReco(sampling=samp) + @testset "Sampling" for samp in sampling + testCSReco(sampling=samp;arrayType) end - testCSRecoMultCoil(type = ComplexF64) - testCSRecoMultCoil(type = ComplexF32) - testCSSenseReco() - testCSReco3d() - testCSSenseReco3d() - testCSRecoMultiCoilCGNR(type = ComplexF64) - testCSRecoMultiCoilCGNR(type = ComplexF32) + @testset "MultiCoil F64" testCSRecoMultCoil(type = ComplexF64;arrayType) + @testset "MultiCoil F32" testCSRecoMultCoil(type = ComplexF32;arrayType) + @testset "Sense" testCSSenseReco(;arrayType) + @testset "3D" testCSReco3d(;arrayType) + @testset "Sense 3D" testCSSenseReco3d(;arrayType) + @testset "Multi Coil CGNR" testCSRecoMultiCoilCGNR(type = ComplexF64;arrayType) + @testset "Multi Coil CGNR" testCSRecoMultiCoilCGNR(type = ComplexF32;arrayType) end @testset "off-resonance" begin accelMethods = ["nfft", "hist", "leastsquare"] for a in accelMethods - !Sys.iswindows() && testOffresonanceReco(accelMethod=a) - !Sys.iswindows() && testOffresonanceSENSEReco() + !Sys.iswindows() && testOffresonanceReco(accelMethod=a;arrayType) + !Sys.iswindows() && testOffresonanceSENSEReco(;arrayType) end end @testset "SENSE" begin - testSENSEReco(64, ComplexF32) - testSENSEReco(64, ComplexF64) - testSENSEnoiseUnCorr(64, ComplexF64) + @testset "SENSE F32" testSENSEReco(64, ComplexF32;arrayType) + @testset "SENSE F64" testSENSEReco(64, ComplexF64;arrayType) + @testset "SENSE F64 Uncorr" testSENSEnoiseUnCorr(64, ComplexF64;arrayType) end end end -testReco() +for arrayType in arrayTypes + testReco(;arrayType) +end \ No newline at end of file