API

FFTNormalOp(img_shape, trj, U; cmaps)
FFTNormalOp(M, U; cmaps)
FFTNormalOp(Λ; cmaps, eltype_x)

Create normal operator of FFT operator. Differentiate between functions exploiting a pre-calculated kernel basis Λ and the functions which calculate Λ based on a passed trajectory trj or mask M.

Arguments

• img_shape::Tuple{Int}: Image dimensions
• traj::Vector{Matrix}: Trajectory
• U::Matrix: Basis coefficients of subspace
• cmaps::Matrix=(1,): Coil sensitivities
• sample_mask::AbstractArray{Bool}=trues(size(trj)[2:end]): Mask indicating which acquired k-space samples are included in the reconstruction
• M::Vector{Matrix}: Mask specific to FFTNormalOp(M, U; cmaps) method. Here, M represents a Cartesian binary mask that indicates in every time frame which phase-encoding lines were acquired. This means, for example, M is sized (Nx, Ny, Nt) for a 2D image. Explicit use of a trajectory with FFTNormalOp(img_shape, trj, U; cmaps) is generally recommended.
• Λ::Array{Complex{T},3}: Toeplitz kernel basis resulting from internal calculate_kernel_cartesian method.
• num_fft_threads::Int=round(Int, Threads.nthreads()/size(U, 2)) or `round(Int, Threads.nthreads()/size(Λ, 1)): Number of Threads for FFT
• eltype_x=eltype(Λ) define the type of x (in the product FFTNormalOp(Λ) * x). The default is the same eltype as Λ

References

1. Tamir JI, et al. “T2 shuffling: Sharp, multicontrast, volumetric fast spin-echo imaging”. Magn Reson Med. 77.1 (2017), pp. 180–195. https://doi.org/10.1002/mrm.26102
2. Assländer J, et al. “Low rank alternating direction method of multipliers reconstruction for MR fingerprinting”. Magn Reson Med 79.1 (2018), pp. 83–96. https://doi.org/10.1002/mrm.26639

NFFTNormalOp(img_shape, trj, U; cmaps, sample_mask, verbose, num_fft_threads)
NFFTNormalOp(img_shape, Λ, kmask_indcs; cmaps)

Create normal operator of NFFT operator. Differentiate between functions exploiting a pre-calculated Toeplitz kernel basis Λ and the function which calculates Λ based on a passed trajectory trj. When the basis functions U are real-valued, a real-only NUFFT is used to compute Λ, reducing the data volume and computation time for the spreading and interpolation steps by half.

Arguments

• img_shape::Tuple{Int}: Image dimensions
• trj::AbstractArray: Trajectory, use CuArray as input type to use CUDA code.
• U::AbstractMatrix: Basis coefficients of subspace
• cmaps::AbstractVector{Matrix}=(1,): Coil sensitivities, use AbstractVector{CuArray} as type for use with CUDA code.
• sample_mask::AbstractArray{Bool} = trues(size(trj)[2:end]): Mask indicating which acquired k-space samples are included in the reconstruction
• verbose::Boolean=false: Verbose level
• num_fft_threads::Int=round(Int, Threads.nthreads()/size(U, 2)) or round(Int, Threads.nthreads()/size(Λ, 1)): Number of threads for FFT
• Λ: Toeplitz kernel basis resulting from internal calculate_kernel_noncartesian method. Using a real basis U reduces its size by a factor of 2.

References

1. Wajer FTAW, and Pruessmann, KP. “Major Speedup of Reconstruction for Sensitivity Encoding with Arbitrary Trajectories”. In: Proc. Intl. Soc. Mag. Reson. Med 9 (2001).
2. Fessler JA, et al. "Toeplitz-Based Iterative Image Reconstruction for MRI With Correction for Magnetic Field Inhomogeneity". IEEE Trans. Signal Process., 53.9 (2006).
3. Mani M, et al. “Fast iterative algorithm for the reconstruction of multishot non-cartesian diffusion data”. Magn Reson Med. 74.4 (2015), pp. 1086–1094. https://doi.org/10.1002/mrm.25486
4. Uecker M, Zhang S, and Frahm J. “Nonlinear inverse reconstruction for real-time MRI of the human heart using undersampled radial FLASH". Magn Res Med. 63 (2010), pp. 1456–1462. https://doi.org/10.1002/mrm.22453

calculate_backprojection(data, trj, img_shape; U, sample_mask, density_compensation, verbose)
calculate_backprojection(data, trj, cmaps::AbstractVector{<:AbstractArray{Tc,N}}; U, sample_mask, density_compensation, verbose)
calculate_backprojection(data, trj, img_shape; U, sample_mask, density_compensation, verbose)
calculate_backprojection(data, trj, cmaps; U, sample_mask)

Calculate (filtered) backprojection.

Arguments

• data <: AbstractArray{Complex}: Complex dataset with axes (samples, time frames, channels). Time frames are reconstructed using the subspace defined in U. Use CuArray as input type to use CUDA GPU code.
• trj <: AbstractArray: Trajectory with sample coordinates corresponding to the dataset. For a Cartesian reconstruction, use T <: Int and define trj[idim,it,ik] ∈ (1, img_shape[idim]). If T <: Float, the NFFT is used. Use CuArray as input type to use CUDA code.

One of the following arguments needs to be supplied

• img_shape::NTuple{N,Int}: Shape of output image; in this case, the data is reconstructed per coil.
• cmaps::AbstractVector{<:AbstractArray{Tc}}: Coil sensitivities; when provided, the coils are combined into an array with no coil axis. Use AbstractVector{CuArray{Tc,N}} as type for use with CUDA code.

Optional Keyword Arguments

• U::Matrix=I(size(trj)[end]) or =I(1): Basis coefficients of subspace (only defined if data and trj have different timeframes)
• sample_mask::AbstractArray{Bool}=trues(size(trj)[2:end]): Mask indicating which acquired k-space samples are included in the reconstruction
• density_compensation=:none: Values of :radial_3D, :radial_2D, :none, or of type AbstractVector{<:AbstractVector}
• verbose::Boolean=false: Verbosity level

calculate_coil_maps(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::AbstractArray: Complex dataset with axes (samples, time frames, channels). Time frames are reconstructed using the subspace defined in U. Use CuArray as input type to use CUDA code.
• trj::AbstractArray: Trajectory with sample coordinates corresponding to the dataset Use CuArray as input type to use GPU code.
• img_shape::NTuple{N,Int}: Shape of output image

Keyword Arguments

• U::Matrix=ones(size(trj)[end]) : I(1): Basis coefficients of subspace
• sample_mask::AbstractArray{Bool}=trues(size(trj)[2:end]): Mask indicating which acquired k-space samples are included in the reconstruction
• density_compensation=:radial_3D: Values of :radial_3D, :radial_2D, :none, or of type AbstractArray
• 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
• Niter_cg: Number of CG iterations used for the reconstruction. The default is 100 and 5 for non-Cartesian and Cartesian trajectories, respectively. Using a number of iterations too large for a given dataset can introduce errors.
• verbose::Boolean=false: Verbosity level

Return

• cmaps::Vector{Array}: Coil sensitivities as Vector of arrays

References

1. Uecker M, Lai P, Murphy MJ, Virtue P, Elad M, Pauly JM, Vasanawala SS, and Lustig M. "ESPIRiT—an eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA". Magn. Reson. Med. 71 (2014), pp. 990-1001. https://doi.org/10.1002/mrm.24751

calculate_golden_means()

Function to calculate the 3D golden means [1].

References

1. Chan RW, Ramsay EA, Cunningham CH, and Plewes DB. "Temporal stability of adaptive 3D radial MRI using multidimensional golden means". Magn. Reson. Med. (2009), 61: 354-363. https://doi.org/10.1002/mrm.21837

grog_calib(data, trj, Nr)

Perform GROG kernel calibration based on the radial trajectory and acquired k-space data [1].

Arguments

• data::AbstractAbstractArray: Complex dataset with dimensions (samples per time frame [Nr], time frames, Rx channels)
• trj::AbstractArray: Trajectory with samples corresponding to the dataset, passed as AbstractArray with dimension (dims, samples per time frame, time frames)
• Nr::Int: Number of samples per readout

References

1. Seiberlich N, Breuer F, Blaimer M, Jakob P, and Griswold M. "Self-calibrating GRAPPA operator gridding for radial and spiral trajectories". Magn. Reson. Med. 59 (2008), pp. 930-935. https://doi.org/10.1002/mrm.21565

grog_gridding!(data, trj, lnG, Nr, img_shape)

Perform gridding of data based on pre-calculated GROG kernel.

Arguments

• data::AbstractVector{<:AbstractMatrix}: Complex dataset passed as AbstractVector of matrices
• trj::AbstractArray}: Trajectory with samples corresponding to the dataset passed as AbstractVector of matrices with Float32 entries
• lnG::AbstractArray{Matrix}: Natural logarithm of GROG kernel in all dimensions
• Nr::Int: Number of samples per read out
• img_shape::Tuple{<:Integer}: Image dimensions

Output

• trj::AbstractArray{<:Integer}: Cartesian trajectory with the elements trj[idim,ik,it] ∈ (1, img_shape[idim])

Dimensions:

• data: (samples, timesteps, coils, repetitions of sampling pattern)
• trj: (dims, samples, timesteps)
• lnG: (dims][Ncoils, Ncoils)

radial_grog!(data, trj, Nr, img_shape)

Perform GROG kernel calibration and gridding of data in-place [1]. The trajectory is returned with integer values.

Arguments

• data::AbstractArray}: Complex dataset with dimensions (samples per time frame [Nr], time frames, Rx channels)
• trj::AbstractArray: Trajectory with samples corresponding to the dataset, passed as AbstractArray with dimension (dims, samples per time frame, time frames)
• Nr::Int: Number of samples per read out
• img_shape::Tuple{Int}: Image dimensions

Output

• trj::AbstractArray{<:Integer}}: Cartesian trajectory with the elements trj[idim,ik,it] ∈ (1, img_shape[idim])

References

1. Seiberlich N, Breuer F, Blaimer M, Jakob P, and Griswold M. "Self-calibrating GRAPPA operator gridding for radial and spiral trajectories". Magn. Reson. Med. 59 (2008), pp. 930-935. https://doi.org/10.1002/mrm.21565

traj_2d_radial_goldenratio(Nr, Ncyc, Nt; thetaRot, phiRot, delay, N)

Function to calculate a 2D radial trajectory with golden-angle spacing between subsequent readouts [1]. The use of tiny golden angles [2] is supported by modifying N.

Arguments

• Nr::Int: Number of read out samples
• Ncyc::Int: Number of cycles
• Nt::Int: Number of time steps in the trajectory
• thetaRot::Float = 0: Fixed rotation angle along theta
• phiRot::Float = 0: Fixed rotation angle along phi
• delay::Tuple{Float, Float, Float} = (0, 0, 0): Gradient delays in (HF, AP, LR)
• N::Int = 1: Number of tiny golden angle

References

1. Winkelmann S, Schaeffter T, Koehler T, Eggers H, Doessel O. "An optimal radial profile order based on the Golden Ratio for time-resolved MRI". IEEE TMI 26:68-76 (2007)
2. Wundrak S, Paul J, Ulrici J, Hell E, Geibel MA, Bernhardt P, Rottbauer W, Rasche V. "Golden ratio sparse MRI using tiny golden angles". Magn. Reson. Med. 75:2372-2378 (2016)

traj_cartesian(T, Nx, Ny, Nz, Nt)

Generate a 3D Cartesian trajectory.

Arguments

• Nx::Int: Number of readout samples
• Ny::Int: Number of phase encoding lines
• Nz::Int: Number of phase encoding lines (third dimension)
• Nt::Int: Number of time steps in the trajectory

Optional Keyword Argument

• T::Type=Int: Type of output trajectory. If T <: Float, trajectory is defined ∈ (-0.5, 0.5). If T <: Int, trajectory consists of values ∈ (1, N) instead.

traj_kooshball(Nr, theta, phi; thetaRot, phiRot, delay)

Function to calculate a 3D radial kooshball trajectory with custom sets of projection angles.

Arguments

• Nr::Int: Number of read out samples
• theta::Array{Float,2}: Array with dimensions: Ncyc, Nt defining the angles theta for each cycle and time step.
• phi::Array{Float,2}: Array with dimensions: Ncyc, Nt defining the angles phi for each cycle and time step.
• thetaRot::Float = 0: Fixed rotation angle along theta
• phiRot::Float = 0: Fixed rotation angle along phi
• delay::Tuple{Float, Float, Float} = (0, 0, 0): Gradient delays in (HF, AP, LR)

traj_kooshball_goldenratio(Nr, Ncyc, Nt; thetaRot, phiRot, delay)

Function to calculate a 3D radial kooshball trajectory with a golden-means angular spacing of k-space readouts [1].

Arguments

• Nr::Int: Number of read out samples
• Ncyc::Int: Number of cycles
• Nt::Int: Number of time steps in the trajectory
• thetaRot::Float = 0: Fixed rotation angle along theta
• phiRot::Float = 0: Fixed rotation angle along phi
• delay::Tuple{Float, Float, Float}= (0, 0, 0): Gradient delays in (HF, AP, LR)

References

1. Chan RW, Ramsay EA, Cunningham CH, and Plewes DB. "Temporal stability of adaptive 3D radial MRI using multidimensional golden means". Magn. Reson. Med. 61 (2009) pp. 354-363. https://doi.org/10.1002/mrm.21837