From e0d793dacc923547ba4f9a715ba4a7ce50d9da30 Mon Sep 17 00:00:00 2001 From: Sebastian Flassbeck Date: Wed, 24 Jan 2024 14:25:52 -0500 Subject: [PATCH 1/9] Overloaded calcFBP function to allow for easier non-MRF fbps --- src/CoilMaps.jl | 67 ++++++++++++++++++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 20 deletions(-) diff --git a/src/CoilMaps.jl b/src/CoilMaps.jl index 174abe3..c38f8c2 100644 --- a/src/CoilMaps.jl +++ b/src/CoilMaps.jl @@ -1,32 +1,59 @@ -function calcFilteredBackProjection(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}, Ncoils::Int; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D, verbose = false) where {N,T} +function applyDensityCompensation!(dataU::AbstractArray{Complex{T}}, trj::Matrix{T}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D) where {T} + if density_compensation == :radial_3D + dataU .*= transpose(sum(abs2, trj,dims=1)) + elseif density_compensation == :radial_2D + dataU .*= transpose(sqrt.(sum(abs2, trj,dims=1))) + elseif density_compensation == :none + # do nothing here + elseif isa(density_compensation, Symbol) + error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector{T}}`") + else + @simd for i ∈ CartesianIndices(dataU) + dataU .*= density_compensation + end + end +end + +function calcFilteredBackProjection(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D, verbose = false) where {N,T} p = plan_nfft(reduce(hcat,trj), img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE) + Ncoils = size(data,3) xbp = Array{Complex{T}}(undef, img_shape..., Ncoils) dataU = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt img_idx = CartesianIndices(img_shape) t = @elapsed for icoil ∈ axes(data,3) - if density_compensation == :radial_3D - @simd for i ∈ CartesianIndices(dataU) - dataU[i] = data[i,icoil] * conj(U[i[2],1]) * sum(abs2, @view trj[i[2]][:,i[1]]) - end - elseif density_compensation == :radial_2D - @simd for i ∈ CartesianIndices(dataU) - dataU[i] = data[i,icoil] * conj(U[i[2],1]) * sqrt(sum(abs2, @view trj[i[2]][:,i[1]])) + dataU = data[:,:,icoil] .* U[:,1]' + for it in axes(data,2) + if typeof(density_compensation) <:AbstractVector + @views applyDensityCompensation!(dataU[:,it],trj[it]; density_compensation = density_compensation[it]) + else + @views applyDensityCompensation!(dataU[:,it],trj[it]; density_compensation) end - elseif density_compensation == :none - # no density compensation; premultiply data with inverse of sampling density before calling function - @simd for i ∈ CartesianIndices(dataU) - dataU[i] = data[i,icoil] * conj(U[i[2],1]) - end - elseif isa(density_compensation, Symbol) - error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector{T}}`") + end + @views mul!(xbp[img_idx,icoil], adjoint(p), vec(dataU)) + end + verbose && println("BP for coils maps: $t s") + return xbp +end + + +function calcFilteredBackProjection(data::AbstractArray{Complex{T}}, trj::Matrix{T}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D, verbose = false) where {N,T} + p = plan_nfft(trj, img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE) + Ncoils = size(data,2) + + xbp = Array{Complex{T}}(undef, img_shape..., Ncoils) + + img_idx = CartesianIndices(img_shape) + t = @elapsed for icoil ∈ axes(data,2) + data_temp = data[:,icoil] # size = Ncycles*Nr + + if typeof(density_compensation) <:AbstractVector + @views applyDensityCompensation!(data_temp,trj; density_compensation = density_compensation[it]) else - @simd for i ∈ CartesianIndices(dataU) - dataU[i] = data[i,icoil] * conj(U[i[2],1]) * density_compensation[i[2]][i[1]] - end + @views applyDensityCompensation!(data_temp,trj; density_compensation) end - @views mul!(xbp[img_idx,icoil], adjoint(p), vec(dataU)) + @views mul!(xbp[img_idx,icoil], adjoint(p), vec(data_temp)) end verbose && println("BP for coils maps: $t s") return xbp @@ -38,7 +65,7 @@ function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:A Ndims = length(img_shape) imdims = ntuple(i->i, Ndims) - xbp = calcFilteredBackProjection(data, trj, U, img_shape, Ncoils; density_compensation, verbose) + xbp = calcFilteredBackProjection(data, trj, U, img_shape; density_compensation, verbose) img_idx = CartesianIndices(img_shape) kbp = fftshift(xbp, imdims) From f9a0f38fa16f7fb8fbd6eedf4920ebec8c7d81d8 Mon Sep 17 00:00:00 2001 From: Sebastian Flassbeck Date: Wed, 24 Jan 2024 20:46:12 -0500 Subject: [PATCH 2/9] fixed some bug and cleaned some more --- src/CoilMaps.jl | 42 +++++++++++++++++++++--------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/CoilMaps.jl b/src/CoilMaps.jl index c38f8c2..073c020 100644 --- a/src/CoilMaps.jl +++ b/src/CoilMaps.jl @@ -1,15 +1,23 @@ -function applyDensityCompensation!(dataU::AbstractArray{Complex{T}}, trj::Matrix{T}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D) where {T} +function applyDensityCompensation!(data::AbstractArray{Complex{T}}, trj; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D) where {T} + if typeof(trj) <: AbstractMatrix + trj = [trj] + end + if density_compensation == :radial_3D - dataU .*= transpose(sum(abs2, trj,dims=1)) + for it in axes(data, 2) + data[:,it] .*= transpose(sum(abs2, trj[it],dims=1)) + end elseif density_compensation == :radial_2D - dataU .*= transpose(sqrt.(sum(abs2, trj,dims=1))) + for it in axes(data, 2) + data[:,it] .*= transpose(sqrt.(sum(abs2, trj[it],dims=1))) + end elseif density_compensation == :none # do nothing here elseif isa(density_compensation, Symbol) error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector{T}}`") else - @simd for i ∈ CartesianIndices(dataU) - dataU .*= density_compensation + for it in axes(data, 2) + data[:,it] .*= density_compensation[it] end end end @@ -19,18 +27,14 @@ function calcFilteredBackProjection(data::AbstractArray{Complex{T},3}, trj::Abst Ncoils = size(data,3) xbp = Array{Complex{T}}(undef, img_shape..., Ncoils) - dataU = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt + data_temp = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt img_idx = CartesianIndices(img_shape) t = @elapsed for icoil ∈ axes(data,3) - dataU = data[:,:,icoil] .* U[:,1]' - for it in axes(data,2) - if typeof(density_compensation) <:AbstractVector - @views applyDensityCompensation!(dataU[:,it],trj[it]; density_compensation = density_compensation[it]) - else - @views applyDensityCompensation!(dataU[:,it],trj[it]; density_compensation) - end - end - @views mul!(xbp[img_idx,icoil], adjoint(p), vec(dataU)) + data_temp = data[:,:,icoil] .* U[:,1]' + + applyDensityCompensation!(data_temp,trj; density_compensation) + + @views mul!(xbp[img_idx,icoil], adjoint(p), vec(data_temp)) end verbose && println("BP for coils maps: $t s") return xbp @@ -47,12 +51,8 @@ function calcFilteredBackProjection(data::AbstractArray{Complex{T}}, trj::Matrix t = @elapsed for icoil ∈ axes(data,2) data_temp = data[:,icoil] # size = Ncycles*Nr - if typeof(density_compensation) <:AbstractVector - @views applyDensityCompensation!(data_temp,trj; density_compensation = density_compensation[it]) - else - @views applyDensityCompensation!(data_temp,trj; density_compensation) - end - + applyDensityCompensation!(data_temp,trj; density_compensation) + @views mul!(xbp[img_idx,icoil], adjoint(p), vec(data_temp)) end verbose && println("BP for coils maps: $t s") From db70dfa63ba5da16386bcc42b1f2ab25007d25c1 Mon Sep 17 00:00:00 2001 From: Sebastian Flassbeck Date: Mon, 29 Jan 2024 16:26:31 -0500 Subject: [PATCH 3/9] merged the two functions again. --- src/CoilMaps.jl | 36 ++++++++++-------------------------- 1 file changed, 10 insertions(+), 26 deletions(-) diff --git a/src/CoilMaps.jl b/src/CoilMaps.jl index 073c020..8014b84 100644 --- a/src/CoilMaps.jl +++ b/src/CoilMaps.jl @@ -1,8 +1,4 @@ function applyDensityCompensation!(data::AbstractArray{Complex{T}}, trj; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D) where {T} - if typeof(trj) <: AbstractMatrix - trj = [trj] - end - if density_compensation == :radial_3D for it in axes(data, 2) data[:,it] .*= transpose(sum(abs2, trj[it],dims=1)) @@ -22,7 +18,14 @@ function applyDensityCompensation!(data::AbstractArray{Complex{T}}, trj; density end end -function calcFilteredBackProjection(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D, verbose = false) where {N,T} + +function calculateBackProjection(data, trj, img_shape; U = ones(1,1), density_compensation=:radial_3D, verbose = false) + T = typeof(real(data[1])) + if typeof(trj) <: AbstractMatrix + trj = [trj] + data = reshape(data,size(data,1),1,size(data,2)) + end + p = plan_nfft(reduce(hcat,trj), img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE) Ncoils = size(data,3) xbp = Array{Complex{T}}(undef, img_shape..., Ncoils) @@ -30,7 +33,7 @@ function calcFilteredBackProjection(data::AbstractArray{Complex{T},3}, trj::Abst data_temp = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt img_idx = CartesianIndices(img_shape) t = @elapsed for icoil ∈ axes(data,3) - data_temp = data[:,:,icoil] .* U[:,1]' + data_temp .= data[:,:,icoil] .* U[:,1]' applyDensityCompensation!(data_temp,trj; density_compensation) @@ -41,31 +44,12 @@ function calcFilteredBackProjection(data::AbstractArray{Complex{T},3}, trj::Abst end -function calcFilteredBackProjection(data::AbstractArray{Complex{T}}, trj::Matrix{T}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D, verbose = false) where {N,T} - p = plan_nfft(trj, img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE) - Ncoils = size(data,2) - - xbp = Array{Complex{T}}(undef, img_shape..., Ncoils) - - img_idx = CartesianIndices(img_shape) - t = @elapsed for icoil ∈ axes(data,2) - data_temp = data[:,icoil] # size = Ncycles*Nr - - applyDensityCompensation!(data_temp,trj; density_compensation) - - @views mul!(xbp[img_idx,icoil], adjoint(p), vec(data_temp)) - end - verbose && println("BP for coils maps: $t s") - return xbp -end - - function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D, kernel_size = ntuple(_->6, N), calib_size = ntuple(_->24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose = false) where {N,T} Ncoils = size(data,3) Ndims = length(img_shape) imdims = ntuple(i->i, Ndims) - xbp = calcFilteredBackProjection(data, trj, U, img_shape; density_compensation, verbose) + xbp = calculateBackProjection(data, trj, img_shape; U, density_compensation, verbose) img_idx = CartesianIndices(img_shape) kbp = fftshift(xbp, imdims) From 896972a2bd129af89c8c2923fac2929516a8ac5f Mon Sep 17 00:00:00 2001 From: Jakob Asslaender Date: Mon, 29 Jan 2024 19:11:34 -0500 Subject: [PATCH 4/9] some fixes --- src/CoilMaps.jl | 67 ++++++++++++++++++++++++------------------------- 1 file changed, 33 insertions(+), 34 deletions(-) diff --git a/src/CoilMaps.jl b/src/CoilMaps.jl index 8014b84..e3b6ed7 100644 --- a/src/CoilMaps.jl +++ b/src/CoilMaps.jl @@ -1,53 +1,52 @@ -function applyDensityCompensation!(data::AbstractArray{Complex{T}}, trj; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D) where {T} - if density_compensation == :radial_3D - for it in axes(data, 2) - data[:,it] .*= transpose(sum(abs2, trj[it],dims=1)) - end - elseif density_compensation == :radial_2D - for it in axes(data, 2) - data[:,it] .*= transpose(sqrt.(sum(abs2, trj[it],dims=1))) - end - elseif density_compensation == :none - # do nothing here - elseif isa(density_compensation, Symbol) - error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector{T}}`") - else - for it in axes(data, 2) - data[:,it] .*= density_compensation[it] +function applyDensityCompensation!(data, trj; density_compensation=:radial_3D) + for it in axes(data, 2) + if density_compensation == :radial_3D + data[:, it] .*= transpose(sum(abs2, trj[it], dims=1)) + elseif density_compensation == :radial_2D + data[:, it] .*= transpose(sqrt.(sum(abs2, trj[it], dims=1))) + elseif density_compensation == :none + # do nothing here + elseif isa(density_compensation, AbstractVector{<:AbstractVector}) + data[:, it] .*= density_compensation[it] + else + error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector}`") end end end -function calculateBackProjection(data, trj, img_shape; U = ones(1,1), density_compensation=:radial_3D, verbose = false) - T = typeof(real(data[1])) +function calculateBackProjection(data::AbstractArray{T}, trj, img_shape; U=ones(T, 1, 1), density_compensation=:radial_3D, verbose=false) where T if typeof(trj) <: AbstractMatrix trj = [trj] - data = reshape(data,size(data,1),1,size(data,2)) end - - p = plan_nfft(reduce(hcat,trj), img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE) - Ncoils = size(data,3) - xbp = Array{Complex{T}}(undef, img_shape..., Ncoils) - data_temp = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt + if ndims(data) == 2 + data = reshape(data, size(data, 1), 1, size(data, 2)) + end + + Ncoils = size(data, 3) + + p = plan_nfft(reduce(hcat, trj), img_shape; precompute=TENSOR, blocking=true, fftflags=FFTW.MEASURE) + xbp = Array{T}(undef, img_shape..., Ncoils) + + data_temp = similar(@view data[:, :, 1]) # size = Ncycles*Nr x Nt img_idx = CartesianIndices(img_shape) - t = @elapsed for icoil ∈ axes(data,3) - data_temp .= data[:,:,icoil] .* U[:,1]' + t = @elapsed for icoil ∈ axes(data, 3) + data_temp .= data[:, :, icoil] .* U[:, 1]' - applyDensityCompensation!(data_temp,trj; density_compensation) + applyDensityCompensation!(data_temp, trj; density_compensation) - @views mul!(xbp[img_idx,icoil], adjoint(p), vec(data_temp)) + @views mul!(xbp[img_idx, icoil], adjoint(p), vec(data_temp)) end - verbose && println("BP for coils maps: $t s") + verbose && println("Backprojection: $t s") return xbp end -function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol, <:AbstractVector{<:AbstractVector{T}}}=:radial_3D, kernel_size = ntuple(_->6, N), calib_size = ntuple(_->24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose = false) where {N,T} - Ncoils = size(data,3) +function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol,<:AbstractVector{<:AbstractVector{T}}}=:radial_3D, kernel_size=ntuple(_ -> 6, N), calib_size=ntuple(_ -> 24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose=false) where {N,T} + Ncoils = size(data, 3) Ndims = length(img_shape) - imdims = ntuple(i->i, Ndims) + imdims = ntuple(i -> i, Ndims) xbp = calculateBackProjection(data, trj, img_shape; U, density_compensation, verbose) @@ -57,13 +56,13 @@ function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:A kbp = fftshift(kbp, imdims) m = CartesianIndices(calib_size) .+ CartesianIndex((img_shape .- calib_size) .÷ 2) - kbp = kbp[m,:] + kbp = kbp[m, :] t = @elapsed begin cmaps = espirit(kbp, img_shape, kernel_size, eigThresh_1=eigThresh_1, eigThresh_2=eigThresh_2, nmaps=nmaps) end verbose && println("espirit: $t s") - cmaps = [cmaps[img_idx,ic,1] for ic=1:Ncoils] + cmaps = [cmaps[img_idx, ic, 1] for ic = 1:Ncoils] return cmaps end From 9dc7203fb8115af2382db6a5419071e9543e2d3e Mon Sep 17 00:00:00 2001 From: Jakob Asslaender Date: Mon, 29 Jan 2024 19:17:24 -0500 Subject: [PATCH 5/9] merge backprojection (draft) --- src/BackProjection.jl | 43 ++++++++++++++++++++++++++++++++----- src/CoilMaps.jl | 49 ++----------------------------------------- 2 files changed, 40 insertions(+), 52 deletions(-) diff --git a/src/BackProjection.jl b/src/BackProjection.jl index 3c1bc17..2396c1c 100644 --- a/src/BackProjection.jl +++ b/src/BackProjection.jl @@ -1,4 +1,19 @@ -function calculateBackProjection(data::AbstractArray{T}, trj, U, cmaps; verbose = false) where {T} + +function calculateBackProjection(data::AbstractArray{T}, trj, img_shape::NTuple{N,Int}; U=ones(T, 1, 1), density_compensation=:radial_3D, verbose=false) where {N,T} + cmaps = [ones(T, img_shape)] + xbp = calculateBackProjection(data, trj, cmaps; U, density_compensation, verbose) + return xbp +end + +function calculateBackProjection(data::AbstractArray{T}, trj, cmaps::AbstractVector{<:AbstractArray}; U=ones(T, 1, 1), density_compensation=:radial_3D, verbose=false) where T + if typeof(trj) <: AbstractMatrix + trj = [trj] + end + + if ndims(data) == 2 + data = reshape(data, size(data, 1), 1, size(data, 2)) + end + test_dimension(data, trj, U, cmaps) _, Ncoef = size(U) @@ -8,14 +23,16 @@ function calculateBackProjection(data::AbstractArray{T}, trj, U, cmaps; verbose xbp = zeros(T, img_shape..., Ncoef) xtmp = Array{T}(undef, img_shape) - dataU = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt + data_temp = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt img_idx = CartesianIndices(img_shape) for icoef ∈ axes(U,2) t = @elapsed for icoil ∈ eachindex(cmaps) - @simd for i ∈ CartesianIndices(dataU) - @inbounds dataU[i] = data[i,icoil] * conj(U[i[2],icoef]) + @simd for i ∈ CartesianIndices(data_temp) + @inbounds data_temp[i] = data[i,icoil] * conj(U[i[2],icoef]) end - mul!(xtmp, adjoint(p), vec(dataU)) + applyDensityCompensation!(data_temp, trj; density_compensation) + + mul!(xtmp, adjoint(p), vec(data_temp)) xbp[img_idx,icoef] .+= conj.(cmaps[icoil]) .* xtmp end verbose && println("coefficient = $icoef: t = $t s"); flush(stdout) @@ -23,6 +40,22 @@ function calculateBackProjection(data::AbstractArray{T}, trj, U, cmaps; verbose return xbp end +function applyDensityCompensation!(data, trj; density_compensation=:radial_3D) + for it in axes(data, 2) + if density_compensation == :radial_3D + data[:, it] .*= transpose(sum(abs2, trj[it], dims=1)) + elseif density_compensation == :radial_2D + data[:, it] .*= transpose(sqrt.(sum(abs2, trj[it], dims=1))) + elseif density_compensation == :none + # do nothing here + elseif isa(density_compensation, AbstractVector{<:AbstractVector}) + data[:, it] .*= density_compensation[it] + else + error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector}`") + end + end +end + function test_dimension(data, trj, U, cmaps) Nt, _ = size(U) img_shape = size(cmaps[1]) diff --git a/src/CoilMaps.jl b/src/CoilMaps.jl index e3b6ed7..28ef4bb 100644 --- a/src/CoilMaps.jl +++ b/src/CoilMaps.jl @@ -1,54 +1,9 @@ -function applyDensityCompensation!(data, trj; density_compensation=:radial_3D) - for it in axes(data, 2) - if density_compensation == :radial_3D - data[:, it] .*= transpose(sum(abs2, trj[it], dims=1)) - elseif density_compensation == :radial_2D - data[:, it] .*= transpose(sqrt.(sum(abs2, trj[it], dims=1))) - elseif density_compensation == :none - # do nothing here - elseif isa(density_compensation, AbstractVector{<:AbstractVector}) - data[:, it] .*= density_compensation[it] - else - error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector}`") - end - end -end - - -function calculateBackProjection(data::AbstractArray{T}, trj, img_shape; U=ones(T, 1, 1), density_compensation=:radial_3D, verbose=false) where T - if typeof(trj) <: AbstractMatrix - trj = [trj] - end - - if ndims(data) == 2 - data = reshape(data, size(data, 1), 1, size(data, 2)) - end - - Ncoils = size(data, 3) - - p = plan_nfft(reduce(hcat, trj), img_shape; precompute=TENSOR, blocking=true, fftflags=FFTW.MEASURE) - xbp = Array{T}(undef, img_shape..., Ncoils) - - data_temp = similar(@view data[:, :, 1]) # size = Ncycles*Nr x Nt - img_idx = CartesianIndices(img_shape) - t = @elapsed for icoil ∈ axes(data, 3) - data_temp .= data[:, :, icoil] .* U[:, 1]' - - applyDensityCompensation!(data_temp, trj; density_compensation) - - @views mul!(xbp[img_idx, icoil], adjoint(p), vec(data_temp)) - end - verbose && println("Backprojection: $t s") - return xbp -end - - -function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation::Union{Symbol,<:AbstractVector{<:AbstractVector{T}}}=:radial_3D, kernel_size=ntuple(_ -> 6, N), calib_size=ntuple(_ -> 24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose=false) where {N,T} +function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation=:radial_3D, kernel_size=ntuple(_ -> 6, N), calib_size=ntuple(_ -> 24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose=false) where {N,T} Ncoils = size(data, 3) Ndims = length(img_shape) imdims = ntuple(i -> i, Ndims) - xbp = calculateBackProjection(data, trj, img_shape; U, density_compensation, verbose) + xbp = calculateBackProjection(data, trj, img_shape; U=U[:,1], density_compensation, verbose) img_idx = CartesianIndices(img_shape) kbp = fftshift(xbp, imdims) From 8f420cbc7cc6b0a4bbbcdfc7f3ab5ca83e975035 Mon Sep 17 00:00:00 2001 From: Jakob Asslaender Date: Mon, 29 Jan 2024 20:46:40 -0500 Subject: [PATCH 6/9] bugfixes --- src/BackProjection.jl | 46 +++++++++++++++++++++++++++++------- src/CoilMaps.jl | 1 + src/NFFTNormalOpBasisFunc.jl | 6 ++--- 3 files changed, 42 insertions(+), 11 deletions(-) diff --git a/src/BackProjection.jl b/src/BackProjection.jl index 2396c1c..0d31bc2 100644 --- a/src/BackProjection.jl +++ b/src/BackProjection.jl @@ -1,11 +1,40 @@ +function calculateBackProjection(data::AbstractArray{T}, trj, img_shape::NTuple{N,Int}; U = N==3 ? I(size(data,2)) : I(1), density_compensation=:none, verbose=false) where {N,T} + if typeof(trj) <: AbstractMatrix + trj = [trj] + end -function calculateBackProjection(data::AbstractArray{T}, trj, img_shape::NTuple{N,Int}; U=ones(T, 1, 1), density_compensation=:radial_3D, verbose=false) where {N,T} - cmaps = [ones(T, img_shape)] - xbp = calculateBackProjection(data, trj, cmaps; U, density_compensation, verbose) + if ndims(data) == 2 + data = reshape(data, size(data, 1), 1, size(data, 2)) + end + Ncoils = size(data, 3) + Ncoef = size(U,2) + + p = plan_nfft(reduce(hcat, trj), img_shape; precompute=TENSOR, blocking=true, fftflags=FFTW.MEASURE) + xbp = Array{T}(undef, img_shape..., Ncoef, Ncoils) + + data_temp = similar(@view data[:, :, 1]) # size = Ncycles*Nr x Nt + img_idx = CartesianIndices(img_shape) + verbose && println("calculating backprojection..."); flush(stdout) + for icoef ∈ axes(U,2) + t = @elapsed for icoil ∈ axes(data, 3) + @simd for i ∈ CartesianIndices(data_temp) + @inbounds data_temp[i] = data[i,icoil] * conj(U[i[2],icoef]) + end + applyDensityCompensation!(data_temp, trj; density_compensation) + @views mul!(xbp[img_idx, icoef, icoil], adjoint(p), vec(data_temp)) + end + verbose && println("coefficient = $icoef: t = $t s"); flush(stdout) + end return xbp end -function calculateBackProjection(data::AbstractArray{T}, trj, cmaps::AbstractVector{<:AbstractArray}; U=ones(T, 1, 1), density_compensation=:radial_3D, verbose=false) where T +function calculateBackProjection(data::AbstractArray{T}, trj, U, cmaps::AbstractVector{<:AbstractArray{T}}; density_compensation=:none, verbose=false) where T + @warn "calculateBackProjection(data, trj, U, cmaps) has been deprecated – call calculateBackProjection(data, trj, cmaps; U=U) with U as a keyword argument instead." maxlog=1 + return calculateBackProjection(data, trj, cmaps; U, density_compensation, verbose) +end + + +function calculateBackProjection(data::AbstractArray{T,N}, trj, cmaps::AbstractVector{<:AbstractArray{T}}; U = N==3 ? I(size(data,2)) : I(1), density_compensation=:none, verbose=false) where {N,T} if typeof(trj) <: AbstractMatrix trj = [trj] end @@ -16,7 +45,7 @@ function calculateBackProjection(data::AbstractArray{T}, trj, cmaps::AbstractVec test_dimension(data, trj, U, cmaps) - _, Ncoef = size(U) + Ncoef = size(U,2) img_shape = size(cmaps[1]) p = plan_nfft(reduce(hcat,trj), img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE) @@ -25,6 +54,7 @@ function calculateBackProjection(data::AbstractArray{T}, trj, cmaps::AbstractVec data_temp = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt img_idx = CartesianIndices(img_shape) + verbose && println("calculating backprojection..."); flush(stdout) for icoef ∈ axes(U,2) t = @elapsed for icoil ∈ eachindex(cmaps) @simd for i ∈ CartesianIndices(data_temp) @@ -57,9 +87,9 @@ function applyDensityCompensation!(data, trj; density_compensation=:radial_3D) end function test_dimension(data, trj, U, cmaps) - Nt, _ = size(U) - img_shape = size(cmaps[1]) - Ncoils = length(cmaps) + Nt = size(U,1) + img_shape = size(cmaps)[1:end-1] + Ncoils = size(cmaps)[end] Nt != size(data, 2) && ArgumentError( "The second dimension of data ($(size(data, 2))) and the first one of U ($Nt) do not match. Both should be number of time points.", diff --git a/src/CoilMaps.jl b/src/CoilMaps.jl index 28ef4bb..35cdb3c 100644 --- a/src/CoilMaps.jl +++ b/src/CoilMaps.jl @@ -4,6 +4,7 @@ function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:A imdims = ntuple(i -> i, Ndims) xbp = calculateBackProjection(data, trj, img_shape; U=U[:,1], density_compensation, verbose) + xbp = dropdims(xbp, dims=ndims(xbp)-1) img_idx = CartesianIndices(img_shape) kbp = fftshift(xbp, imdims) diff --git a/src/NFFTNormalOpBasisFunc.jl b/src/NFFTNormalOpBasisFunc.jl index c2bc485..9442716 100644 --- a/src/NFFTNormalOpBasisFunc.jl +++ b/src/NFFTNormalOpBasisFunc.jl @@ -61,7 +61,7 @@ function NFFTNormalOpBasisFunc( img_shape, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Tc}; - cmaps = (1,), + cmaps=[ones(T, img_shape)], verbose = false, Λ_kmask_indcs = calculateToeplitzKernelBasis(2 .* img_shape, trj, U; verbose = verbose), num_fft_threads = round(Int, Threads.nthreads()/size(U, 2)) @@ -154,13 +154,13 @@ function NFFTNormalOpBasisFuncLO( img_shape, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Tc}; - cmaps = (1,), + cmaps=[ones(T, img_shape)], verbose = false, Λ_kmask_indcs = calculateToeplitzKernelBasis(2 .* img_shape, trj, U; verbose = verbose), num_fft_threads = round(Int, Threads.nthreads()/size(U, 2)) ) where {T, Tc <: Union{T, Complex{T}}} - S = NFFTNormalOpBasisFunc(img_shape, trj, U; cmaps = cmaps, Λ_kmask_indcs = Λ_kmask_indcs, num_fft_threads = num_fft_threads) + S = NFFTNormalOpBasisFunc(img_shape, trj, U; cmaps, Λ_kmask_indcs, num_fft_threads) return NFFTNormalOpBasisFuncLO(S) end From 79610bbc70a5b922411897b0c312b137df4e5ca5 Mon Sep 17 00:00:00 2001 From: Sebastian Flassbeck Date: Fri, 16 Feb 2024 16:58:02 -0500 Subject: [PATCH 7/9] Made call of calcCoilMaps consistent with call of calculateBackProjection --- src/CoilMaps.jl | 2 +- src/deprecated.jl | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) create mode 100644 src/deprecated.jl diff --git a/src/CoilMaps.jl b/src/CoilMaps.jl index 35cdb3c..882c738 100644 --- a/src/CoilMaps.jl +++ b/src/CoilMaps.jl @@ -1,4 +1,4 @@ -function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation=:radial_3D, kernel_size=ntuple(_ -> 6, N), calib_size=ntuple(_ -> 24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose=false) where {N,T} +function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, img_shape::NTuple{N,Int}; U = N==3 ? I(size(data,2)) : I(1), density_compensation=:radial_3D, kernel_size=ntuple(_ -> 6, N), calib_size=ntuple(_ -> 24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose=false) where {N,T} Ncoils = size(data, 3) Ndims = length(img_shape) imdims = ntuple(i -> i, Ndims) diff --git a/src/deprecated.jl b/src/deprecated.jl new file mode 100644 index 0000000..aa43a05 --- /dev/null +++ b/src/deprecated.jl @@ -0,0 +1,10 @@ +function calculateBackProjection(data::AbstractArray{T}, trj, U, cmaps::AbstractVector{<:AbstractArray{T}}; density_compensation=:none, verbose=false) where T + @warn "calculateBackProjection(data, trj, U, cmaps) has been deprecated – call calculateBackProjection(data, trj, cmaps; U=U) with U as a keyword argument instead." maxlog=1 + return calculateBackProjection(data, trj, cmaps; U, density_compensation, verbose) +end + + +function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:AbstractMatrix{T}}, U::AbstractMatrix{Complex{T}}, img_shape::NTuple{N,Int}; density_compensation=:radial_3D, kernel_size=ntuple(_ -> 6, N), calib_size=ntuple(_ -> 24, N), eigThresh_1=0.01, eigThresh_2=0.9, nmaps=1, verbose=false) where {N,T} + @warn "calcCoilMaps(data, trj, U, img_shape) has been deprecated – call calcCoilMaps(data, trj, img_shape; U=U) with U as a keyword argument instead." maxlog=1 + return calcCoilMaps(data, trj, img_shape; U, density_compensation, kernel_size, calib_size, eigThresh_1, eigThresh_2, nmaps, verbose) +end From 429ae8a87fc595e3e213cc5d6423fa920d0ece8a Mon Sep 17 00:00:00 2001 From: Sebastian Flassbeck Date: Fri, 16 Feb 2024 17:01:49 -0500 Subject: [PATCH 8/9] added deprecated.jl to includes --- src/BackProjection.jl | 6 ------ src/MRFingerprintingRecon.jl | 1 + 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/src/BackProjection.jl b/src/BackProjection.jl index 0d31bc2..6c3f425 100644 --- a/src/BackProjection.jl +++ b/src/BackProjection.jl @@ -28,12 +28,6 @@ function calculateBackProjection(data::AbstractArray{T}, trj, img_shape::NTuple{ return xbp end -function calculateBackProjection(data::AbstractArray{T}, trj, U, cmaps::AbstractVector{<:AbstractArray{T}}; density_compensation=:none, verbose=false) where T - @warn "calculateBackProjection(data, trj, U, cmaps) has been deprecated – call calculateBackProjection(data, trj, cmaps; U=U) with U as a keyword argument instead." maxlog=1 - return calculateBackProjection(data, trj, cmaps; U, density_compensation, verbose) -end - - function calculateBackProjection(data::AbstractArray{T,N}, trj, cmaps::AbstractVector{<:AbstractArray{T}}; U = N==3 ? I(size(data,2)) : I(1), density_compensation=:none, verbose=false) where {N,T} if typeof(trj) <: AbstractMatrix trj = [trj] diff --git a/src/MRFingerprintingRecon.jl b/src/MRFingerprintingRecon.jl index 278d0e0..2ae5a34 100644 --- a/src/MRFingerprintingRecon.jl +++ b/src/MRFingerprintingRecon.jl @@ -21,5 +21,6 @@ include("NFFTNormalOpBasisFunc.jl") include("CoilMaps.jl") include("BackProjection.jl") include("Trajectories.jl") +include("deprecated.jl") end # module From 53f9f624727d68fab41f66ac17f2cc8c9c6180f0 Mon Sep 17 00:00:00 2001 From: Jakob Asslaender Date: Fri, 16 Feb 2024 18:06:13 -0500 Subject: [PATCH 9/9] bugfix in test --- test/reconstruct.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/reconstruct.jl b/test/reconstruct.jl index e958722..effbcdd 100644 --- a/test/reconstruct.jl +++ b/test/reconstruct.jl @@ -44,7 +44,9 @@ for it ∈ axes(data,2) end ## BackProjection -b = vec(calculateBackProjection(data, trj, U, [ones(T, Nx,Nx)], verbose = true)) +b = vec(calculateBackProjection(data, trj, (Nx, Nx); U, verbose = true)) +bcoil = vec(calculateBackProjection(data, trj, [ones(Complex{T}, Nx,Nx)]; U, verbose = true)) +@test bcoil ≈ b ## construct forward operator A = NFFTNormalOpBasisFuncLO((Nx,Nx), trj, U; verbose = true)