Changedata2vector (#38)

* Changed BackProjection to work with data in vector format. Changed CoilMaps accordingly. Main changes are in calculateBackProjection, which only works if data is a vector of vectors or vector of matrices (no data type check).

* Changed call to calculateBackProjection to have U = ones(Nt) to work with data as vector - changed call to cmaps too (altough it is not actually used like that in the recon)

* Resolved slow calcBackProjection: now using a long vector and increasing indices, adapted applydensitycomp to work with the format (basically back to original), corrected bug with variable N not specified, other small changes

* small change variable name

* Move GROG code to vectorized dataformat

Not benchmarked against previous version!

* fix GROG to accomodate variable number of k-space points in different time frames

* Bug fixes for data removal in calculateBackProjection

* fix backprojection default U; add a method that takes data of a single time frame as a matrix w/o the timeframe vector wrapper

* fix gridded backprojection

* fix 2D trajectory

* set GROG gridded trj to integer and return trajectory instead of modifying it in place

* Update documentation

* Allow for data vectors of floats

* Minor interface changes + documentation

* test beautifying

* remove @inbounds

* added @views in calculateBackProjection

* Fix Nrep bug (U is a 3D array if Nrep > 1)

src/BackProjection.jl

@@ -1,130 +1,108 @@
 """
  calculateBackProjection(data, trj, img_shape; U, density_compensation, verbose)
- calculateBackProjection(data, trj, cmaps; U, density_compensation, verbose)
+ calculateBackProjection(data, trj, cmaps::AbstractVector{<:AbstractArray{cT,N}}; U, density_compensation, verbose)
+ calculateBackProjection(data, trj, cmaps_img_shape; U, density_compensation, verbose)
+ calculateBackProjection(data, trj, cmaps; U)
 
-Calculate backprojection
+Calculate (filtered) backprojection
 
 # Arguments
-- `data::AbstractArray{T}`: Basis coefficients of subspace
-- `trj::Vector{Matrix{Float32}}`: Trajectory
+- `data <: Union{AbstractVector{<:AbstractMatrix{cT}},AbstractMatrix{cT}}`: Complex dataset either as AbstractVector of matrices or single matrix. The optional outer matrix defines different time frames that are reconstruced in the subspace defined in U.
+- `trj <: Union{:AbstractVector{<:AbstractMatrix{T}},AbstractMatrix{T}}`: Trajectory with samples corresponding to the dataset either as AbstractVector of matrices or single matrix.
 - `img_shape::NTuple{N,Int}`: Shape of image
-- `U::Matrix{ComplexF32}`: Basis coefficients of subspace
-- `density_compensation`: Values of `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector}``
-- `verbose::Boolean`: Verbosity level
-- `cmaps::::AbstractVector{<:AbstractArray{T}}`: Coil sensitivities
 
-"""
-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
+# Optional Keyword Arguments
+- `cmaps::::AbstractVector{<:AbstractArray{T}}`: Coil sensitivities as AbstractVector of arrays
+- `cmaps_img_shape`: Either equal `img_shape` or `cmaps`
+- `U::Matrix` = I(length(data)) or = I(1): Basis coefficients of subspace (only defined if data and trj have different timeframes)
+- `density_compensation`=:`none`: Values of `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector}`
+- `verbose::Boolean`=`false`: Verbosity level
 
- if ndims(data) == 2
- data = reshape(data, size(data, 1), 1, size(data, 2))
- end
- Ncoils = size(data, 3)
+# Notes
+- The type of the elements of the trajectory define if a gridded backprojection (eltype(trj[1]) or eltype(trj) <: Int) or a non-uniform (else) is performed.
+"""
+function calculateBackProjection(data::AbstractVector{<:AbstractArray{cT}}, trj::AbstractVector{<:AbstractMatrix{T}}, img_shape::NTuple{N,Int}; U=I(length(data)), density_compensation=:none, verbose=false) where {T <: Real, cT <: Complex{T},N}
  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)
+ trj_v = reduce(hcat, trj)
+ p = plan_nfft(trj_v, img_shape; precompute=TENSOR, blocking=true, fftflags=FFTW.MEASURE)
+
+ Ncoil = size(data[1], 2)
+ xbp = Array{cT}(undef, img_shape..., Ncoef, Ncoil)
+
+ trj_l = [size(trj[it],2) for it in eachindex(trj)]
+ data_temp = Vector{cT}(undef,sum(trj_l))
 
- 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])
+ for icoef ∈ axes(U, 2)
+ t = @elapsed for icoil ∈ axes(data[1], 2)
+ @simd for it in eachindex(data)
+ idx1 = sum(trj_l[1:it-1]) + 1
+ idx2 = sum(trj_l[1:it])
+ @views data_temp[idx1:idx2] .= data[it][:,icoil] .* conj(U[it,icoef])
  end
- applyDensityCompensation!(data_temp, trj; density_compensation)
- @views mul!(xbp[img_idx, icoef, icoil], adjoint(p), vec(data_temp))
+ applyDensityCompensation!(data_temp, trj_v; density_compensation)
+
+ @views mul!(xbp[img_idx, icoef, icoil], adjoint(p), data_temp)
  end
  verbose && println("coefficient = $icoef: t = $t s"); flush(stdout)
  end
  return xbp
 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
-
- if ndims(data) == 2
- data = reshape(data, size(data, 1), 1, size(data, 2))
- end
-
+function calculateBackProjection(data::AbstractVector{<:AbstractMatrix{cT}}, trj::AbstractVector{<:AbstractMatrix{T}}, cmaps::AbstractVector{<:AbstractArray{cT,N}}; U=I(length(data)), density_compensation=:none, verbose=false) where {T <: Real, cT <: Complex{T}, N}
  test_dimension(data, trj, U, cmaps)
 
- Ncoef = size(U,2)
+ trj_v = reduce(hcat, trj)
+ Ncoef = size(U, 2)
  img_shape = size(cmaps[1])
 
- p = plan_nfft(reduce(hcat,trj), img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE)
- xbp = zeros(T, img_shape..., Ncoef)
- xtmp = Array{T}(undef, img_shape)
+ p = plan_nfft(trj_v, img_shape; precompute=TENSOR, blocking = true, fftflags = FFTW.MEASURE)
+ xbp = zeros(cT, img_shape..., Ncoef)
+ xtmp = Array{cT}(undef, img_shape)
 
- data_temp = similar(@view data[:,:,1]) # size = Ncycles*Nr x Nt
+ trj_l = [size(trj[it],2) for it in eachindex(trj)]
+ data_temp = Vector{cT}(undef,sum(trj_l))
  img_idx = CartesianIndices(img_shape)
  verbose && println("calculating backprojection..."); flush(stdout)
- for icoef ∈ axes(U,2)
+ for icoef ∈ axes(U, 2)
  t = @elapsed for icoil ∈ eachindex(cmaps)
- @simd for i ∈ CartesianIndices(data_temp)
- @inbounds data_temp[i] = data[i,icoil] * conj(U[i[2],icoef])
+ @simd for it in eachindex(data)
+ idx1 = sum(trj_l[1:it-1]) + 1
+ idx2 = sum(trj_l[1:it])
+ @views data_temp[idx1:idx2] .= data[it][:,icoil] .* conj(U[it,icoef])
  end
- applyDensityCompensation!(data_temp, trj; density_compensation)
-
- mul!(xtmp, adjoint(p), vec(data_temp))
+ applyDensityCompensation!(data_temp, trj_v; density_compensation)
+ mul!(xtmp, adjoint(p), data_temp)
  xbp[img_idx,icoef] .+= conj.(cmaps[icoil]) .* xtmp
  end
  verbose && println("coefficient = $icoef: t = $t s"); flush(stdout)
  end
  return xbp
 end
 
-"""
- calculateBackProjection_gridded(data, trj, U, cmaps)
-
-Calculate gridded backprojection
+function calculateBackProjection(data::AbstractArray{cT}, trj::AbstractMatrix{T}, cmaps_img_shape; density_compensation=:none, verbose=false) where {T <: Real, cT <: Complex{T}}
+ return calculateBackProjection([data], [trj], cmaps_img_shape; U=I(1), density_compensation, verbose)
+end
 
-# Arguments
-- `data::Matrix{ComplexF32}`: Basis coefficients of subspace
-- `trj::Vector{Matrix{Float32}}`: Trajectory
-- `U::Matrix{ComplexF32}`: Basis coefficients of subspace
-- `cmaps::Matrix{ComplexF32}`: Coil sensitivities
-
-# Note
-In case of repeated sampling (Nrep > 1), a joint basis reconstruction is required.
-Therefore, the basis needs to have a temporal dimension of Nt⋅Nrep with Nt as time dimension defined by the trajectory.
-"""
-function calculateBackProjection_gridded(data, trj, U, cmaps)
- Ncoil = length(cmaps)
+# Method for GROG gridded data / trajectory
+function calculateBackProjection(data::AbstractVector{<:AbstractArray}, trj::AbstractVector{<:AbstractMatrix{<:Integer}}, cmaps; U=I(length(data)))
  Ncoeff = size(U, 2)
  img_shape = size(cmaps[1])
  img_idx = CartesianIndices(img_shape)
 
- Nt = length(trj)
- Nrep = size(data, 4)
-
- if (1 != Nrep) # Avoid error during reshape that joins rep and t dim
- data = permutedims(data, (1,2,4,3))
- @assert Nt*Nrep == size(U, 1) "Mismatch between data and basis"
- else
- @assert Nt == size(U, 1) "Mismatch between trajectory and basis"
- end
- data = reshape(data, :, Nt*Nrep, Ncoil)
-
- dataU = similar(data, img_shape..., Ncoeff)
- xbp = zeros(eltype(data), img_shape..., Ncoeff)
+ dataU = similar(data[1], img_shape..., Ncoeff)
+ xbp = zeros(eltype(data[1]), img_shape..., Ncoeff)
 
  Threads.@threads for icoef ∈ axes(U, 2)
- for icoil ∈ axes(data, 3)
+ for icoil ∈ axes(data[1], 2)
  dataU[img_idx, icoef] .= 0
 
- for i ∈ CartesianIndices(@view data[:, :, 1, 1])
- t_idx = mod(i[2] + Nt - 1, Nt) + 1 # "mod" to incorporate repeated sampling pattern, "mod(i[2]+Nt-1,Nt)+1" to compensate for one indexing
- k_idx = ntuple(j -> mod1(Int(trj[t_idx][j, i[1]]) - img_shape[j] ÷ 2, img_shape[j]), length(img_shape)) # incorporates ifftshift
+ for it ∈ eachindex(data), is ∈ axes(data[it], 1), irep ∈ axes(data[it], 3)
+ k_idx = ntuple(j -> mod1(Int(trj[it][j, is]) - img_shape[j] ÷ 2, img_shape[j]), length(img_shape)) # incorporates ifftshift
  k_idx = CartesianIndex(k_idx)
-
- @views dataU[k_idx, icoef] += data[i[1], i[2], icoil] * conj(U[i[2], icoef])
+ dataU[k_idx, icoef] += data[it][is, icoil, irep] * conj(U[it, icoef, irep])
  end
 
  @views ifft!(dataU[img_idx, icoef])
@@ -139,25 +117,23 @@ 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
+ if density_compensation == :radial_3D
+ data .*= transpose(sum(abs2, trj, dims=1))
+ elseif density_compensation == :radial_2D
+ data .*= transpose(sqrt.(sum(abs2, trj, dims=1)))
+ elseif density_compensation == :none
+ # do nothing here
+ elseif isa(density_compensation, AbstractVector{<:AbstractVector})
+ data .*= reduce(hcat, density_compensation)
+ else
+ error("`density_compensation` can only be `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector}`")
  end
 end
 
 function test_dimension(data, trj, U, cmaps)
- Nt = size(U,1)
+ Nt = size(U, 1)
  img_shape = size(cmaps)[1:end-1]
- Ncoils = size(cmaps)[end]
+ Ncoil = 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.",
@@ -174,8 +150,8 @@ function test_dimension(data, trj, U, cmaps)
  "`trj` has the length $(length(trj)) and the 2ⁿᵈ dimension of data is $(size(data,2)). They should match and reflect the number of time points.",
  )
 
- Ncoils != size(data, 3) && ArgumentError(
- "The last dimension of `cmaps` is $(Ncoils) and the 3ⁿᵈ dimension of data is $(size(data,3)). They should match and reflect the number of coils.",
+ Ncoil != size(data, 3) && ArgumentError(
+ "The last dimension of `cmaps` is $(Ncoil) and the 3ⁿᵈ dimension of data is $(size(data,3)). They should match and reflect the number of coils.",
  )
 
 end

src/CoilMaps.jl

@@ -1,5 +1,31 @@
-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)
+"""
+ calcCoilMaps(data, trj, img_shape; U, density_compensation, kernel_size, calib_size, eigThresh_1, eigThresh_2, nmaps, verbose)
+
+Estimate coil sensitivity maps using ESPIRiT [1].
+
+# Arguments
+- `data::AbstractVector{<:AbstractMatrix{Complex{T}}}`: Complex dataset either as AbstractVector of matrices or single matrix. The optional outer vector defines different time frames that are combined using the subspace defined in `U`
+- `trj::AbstractVector{<:AbstractMatrix{T}}`: Trajectory with samples corresponding to the dataset either as AbstractVector of matrices or single matrix.
+- `img_shape::NTuple{N,Int}`: Shape of image
+
+# Keyword Arguments
+- `U::Matrix` = N==3 ? ones(size(data,1)) : I(1): Basis coefficients of subspace (only defined if `data` and `trj` are vectors of matrices)
+- `density_compensation`=:`radial_3D`: Values of `:radial_3D`, `:radial_2D`, `:none`, or of type `AbstractVector{<:AbstractVector}`
+- `kernel_size`=`ntuple(_ -> 6, N)`: Kernel size
+- `calib_size`=`ntuple(_ -> 24, N)`: Size of calibration region
+- `eigThresh_1`=0.01: Threshold of first eigenvalue
+- `eigThresh_2`=0.9: Threshold of second eigenvalue
+- `nmaps`=1: Number of estimated maps
+- `verbose::Boolean`=`false`: Verbosity level
+
+# return
+- `cmaps::::Vector{<:Array{Complex{T}}}`: Coil sensitivities as Vector of arrays
+
+# References
+[1] Uecker, M., Lai, P., Murphy, M.J., Virtue, P., Elad, M., Pauly, J.M., Vasanawala, S.S. and Lustig, M. (2014), ESPIRiT—an eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA. Magn. Reson. Med., 71: 990-1001. https://doi.org/10.1002/mrm.24751
+"""
+function calcCoilMaps(data::AbstractVector{<:AbstractMatrix{Complex{T}}}, trj::AbstractVector{<:AbstractMatrix{T}}, img_shape::NTuple{N,Int}; U = ones(length(data)), 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}
+ Ncoil = size(data[1], 2)
  Ndims = length(img_shape)
  imdims = ntuple(i -> i, Ndims)
 
@@ -19,6 +45,10 @@ function calcCoilMaps(data::AbstractArray{Complex{T},3}, trj::AbstractVector{<:A
  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:Ncoil]
  return cmaps
 end
+
+function calcCoilMaps(data::AbstractMatrix{Complex{T}}, trj::AbstractMatrix{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}
+ calcCoilMaps([data], [trj], img_shape; U = I(1), density_compensation, kernel_size, calib_size, eigThresh_1, eigThresh_2, nmaps, verbose)
+end

src/FFTNormalOpBasisFunc.jl

@@ -13,10 +13,11 @@ Differentiate between functions exploiting a pre-calculated kernel basis `Λ` an
 # Arguments
 - `img_shape::Tuple{Int}`: Image dimensions
 - `traj::Vector{Matrix{Float32}}`: Trajectory
-- `U::Matrix{ComplexF32}`: Basis coefficients of subspace
+- `U::Matrix{ComplexF32}`=(1,): Basis coefficients of subspace
 - `cmaps::Matrix{ComplexF32}`: Coil sensitivities
 - `M::Vector{Matrix{Float32}}`: Mask
 - `Λ::Array{Complex{T},3}`: Toeplitz kernel basis
+- `num_fft_threads::Int` = `round(Int, Threads.nthreads()/size(U, 2))` or `round(Int, Threads.nthreads()/size(Λ, 1)): Number of Threads for FFT
 """
 function FFTNormalOp(img_shape, trj, U; cmaps=(1,), num_fft_threads = round(Int, Threads.nthreads()/size(U, 2)))
  Λ = calculateKernelBasis(img_shape, trj, U)
@@ -74,24 +75,13 @@ end
 
 function calculateKernelBasis(img_shape, trj, U)
  Ncoeff = size(U, 2)
- Nt = length(trj) # number of time points
- Nrep = size(U, 1) / Nt
-
- @assert isinteger(Nrep) && (Nrep != 0) "Mismatch between trajectory and basis"
-
  Λ = zeros(eltype(U), Ncoeff, Ncoeff, img_shape...)
-
- for it ∈ 1:size(U, 1)
-
- t_idx = mod(it + Nt - 1, Nt) + 1 # "mod" to incorporate repeated sampling pattern, "mod(i[2]+Nt-1,Nt)+1" to compensate for one indexing
-
- for ix ∈ axes(trj[t_idx], 2)
-
- k_idx = ntuple(j -> mod1(Int(trj[t_idx][j, ix]) - img_shape[j] ÷ 2, img_shape[j]), length(img_shape)) # incorporates ifftshift
+ for it ∈ axes(U, 1)
+ for ix ∈ axes(trj[it], 2)
+ k_idx = ntuple(j -> mod1(Int(trj[it][j, ix]) - img_shape[j] ÷ 2, img_shape[j]), length(img_shape)) # incorporates ifftshift
  k_idx = CartesianIndex(k_idx)
-
- for ic ∈ CartesianIndices((Ncoeff, Ncoeff))
- Λ[ic[1], ic[2], k_idx] += conj(U[it, ic[1]]) * U[it, ic[2]]
+ for ic ∈ CartesianIndices((Ncoeff, Ncoeff)), irep ∈ axes(U, 3)
+ Λ[ic[1], ic[2], k_idx] += conj(U[it, ic[1], irep]) * U[it, ic[2], irep]
  end
  end
  end