MRIFieldmaps

Module MRIFieldmaps exports methods for fieldmap estimation.

Echotime{T} = Union{AbstractVector{<:T}, NTuple{N,<:T} where N}

The echo times can be a vector or a Tuple.

RealU

A data type that is just Number but is to be thought of as Union{Real, Unitful} without needing a dependence on the Unitful package.

finit = b0init(ydata, echotime; kwargs...)

Classic B0 field map estimation based on the phase difference of complex images at two different echo times. If sensitivity maps (smap) are provided, complex coil combination is done first. This code works with images of arbitrary dimensions (2D, 3D, etc.), and with multiple coils.

In the usual case where echotime has units of seconds, the returned B0 fieldmap will have units of Hz.

If df is nonempty (which always holds for water-fat case), then perform discrete maximum-likelihood estimation using fdict.

In

• ydata (dims..., nc, ne) ne sets of complex images for nc ≥ 1 coils
• echotime (ne) vector of ne ≥ 2 echo times (only first 2 are used)

Options

• smap (dims..., nc) complex coil maps; default nothing
• threshold set finit values where |y1| < threshold * max(|y1|) to the mean of the "good" values where |y1| ≥ threshold * max(|y1|). default: 0.1

Options for water-fat case:

• df Δf values in water-fat imaging (def: []) units Hz, e.g., [440] at 3T
• relamp relative amplitudes in multi-species water-fat (def: [])
• fband frequency bandwidth for fdict; defaultfloor(1 / minimum(echo time spacing))`
• nf number of discrete frequencies to try; default 1+floor(fband) so ≈1Hz spacing
• fdict "dictionary" of discrete frequency values to try; default LinRange(-1/2,1/2,nf) * fband

Out

• finit initial B0 fieldmap estimate in Hz

(fhat, times, out) = b0map(ydata, echotime; kwargs...)

Field map estimation from multiple (ne ≥ 2) echo-time images, using preconditioned nonlinear CG (NCG) with a monotonic line search. This code works with images of arbitrary dimensions (2D, 3D, etc.), and with multiple coils.

Caution: single coil data must be reshaped to size (dims..., 1, ne).

The cost function for the single-coil case is: cost(w) = ∑_{j=1}^{#voxel} ∑_{m=1}^ne ∑_{n=1}^ne |y_{mj} y_{nj}| wj (1 - cos(w_j * (t_m - t_n) + ∠y_{ni} - ∠y_{mj})) + R(w) where t_n denotes the echo time of the nth scan and R(w) = 0.5 * | C * w |^2 is a quadratic roughness regularizer based on 1st-order or 2nd-order finite differences. See the documentation for the general multi-coil case.

The initial field map finit and output fhat are field maps in Hz, but internally the code works with ω = 2π f (rad/s).

In

• ydata (dims..., nc, ne) ne sets of complex images for nc ≥ 1 coils
• echotime::Echotime (ne ≥ 2) echo time offsets (in seconds)

Options

• finit (dims) initial fieldmap estimate (in Hz); default from b0init()
• smap (dims..., nc) complex coil maps; default nothing
• mask (dims...) logical reconstruction mask; default: trues(size(finit))
• ninner # of inner iterations for monotonic line search default: 3 inner iterations
• b0init_args::NamedTuple = (;) options for b0init, such as threshold
• niter # of outer iterations (def: 30)
• order order of the finite-difference matrix (def: 2)
• l2b log2 of regularization parameter (def: -6)
• gamma_type CG direction:
  • :PR = Polak-Ribiere (default)
  • :FR = Fletcher-Reeves
• precon Preconditioner:
  • :I (nothing)
  • :diag
  • :chol may require too much memory
  • :ichol (default)
• reset # of iterations before resetting direction (def: Inf)
• df Δf values in water-fat imaging (def: [0]) units Hz, e.g., [440] at 3T
• relamp relative amplitudes in multi-peak water-fat (def: [1])
• lldl_args::NamedTuple options for lldl, default: (;memory=2)
• track::Bool if true then track cost and save all iterations (def: false)
• chat::Bool = true # @info updates each iteration
• chat_iter::Int = 10 # print progress report every few iterations.

Out

• fhat final fieldmap estimate in Hz
• times (niter+1) wall time for each iteration
• out::NamedTuple that contains:
  • (xw, xf) (dims) water / fat images if !iszero(df)
  • finit (dims) initial fieldmap
• if track == true then out also contains:
  • costs (niter+1) (nonconvex) cost for each iteration
  • fhats (dims, niter+1) fieldmap estimates every iteration

The algorithm is based on the paper: C Y Lin, J A Fessler, "Efficient Regularized Field Map Estimation in 3D MRI", IEEE TCI 2020 https://doi.org/10.1109/TCI.2020.3031082 https://arxiv.org/abs/2005.08661

Please cite that paper if you use this code.

b0map()

This version expects masked column-like inputs. For expert use only.

In

• finit (np) initial estimate in Hz (np is # of pixels in mask)
• zdata (np, ne) ne sets of coil-combined measurements
• sos (np) sum-of-squares of coil maps
• echotime (ne) vector of ne echo time offsets
• mask (N) logical reconstruction mask

b0model(fmap, xw, echotime; kwargs...)

Compute complex images for B0 field mapping. This function is used mainly for simulation, and for code testing.

Model: images[j,c,l] = smap[j,c] exp(ı 2π fmap[j] t_l) exp(-relax[j] t_l) x[j,l] where x[j,l] = xw[j] + xf[j] * sum_{p=0}^P α_p exp(ı 2π Δf_p t_l).

Field map estimation from multiple (ne ≥ 2) echo-time images,

In

• fmap (dims) fieldmap (in Hz)
• xw (dims) water magnetization component
• echotime (ne) vector or Tuple of ne echo time offsets (in sec)

Options

• smap (dims..., nc) complex coil maps, default ones(size(fmap)..., 1)
• xw (dims) fat magnetization component, default zeros(size(fmap))
• df Δf values in water-fat imaging (def: [0]) units Hz, e.g., [440] at 3T
• relamp relative amplitudes in multi-peak water-fat (def: [1])
• relax (dims) R2* relaxation in same units as fmap, i.e., 1/s

Out

• ydata (dims..., nc, ne) ne sets of complex images for nc coils

(ydata, scalefactor) = b0scale(ydata, echotime; dmax)

Scale complex images ydata to account for R2* effects and for magnitude variations using median(di) where

• ri = sum_j sum_k |y_{ij} y_{ik}|^2 (t_k - t_j)^2
• di = ri / sum_k |y_{ik}|^2

Only values where di > dmax * maximum(di) affect scalefactor, so it is fine to pass unmasked images here.

This normalization simplifies regularization parameter selection for regularized B0 fieldmap estimation. See eqn (9) and (15) of Funai & Fessler, Oct. 2008, IEEE T-MI, https://doi.org/10.1109/TMI.2008.923956

In

• ydata (dims..., ne) scan images for ne different echo times
• echotime (ne) echo times (units of sec if fieldmap is in Hz)

Option

• dmax::Real threshold for relative di value (default 0.1)

Out

• ydata (dims..., ne) scaled scan images
• scalefactor = sqrt(median(rj))

zdata, sos = coil_combine(ydata [, smap])

For estimating a B0 field map from complete multi-coil image data, it suffices to first do complex coil combination, while tracking the sum-of-squares sos for proper weighting. When sensitivity maps are provided, often sos is all 1's and 0's.

Uses a phase contrast coil combination approach (reference below) when sensitivity maps are not provided.

In

• ydata (dims..., nc, ne) ne ≥ 1 sets of complex images for nc ≥ 1 coils
• smap (dims..., nc) complex coil maps (optional)

Out

• zdata (dims..., ne) complex coil combination: sum_c smap[c]' * ydata[c] ./ sos or sum_c ydata[c,1]' * ydata[c] ./ sos (if smap not provided or isnothing(smap))
• sos (dims...) sum-of-squares: sum_c |smap[c]|^2 or sqrt.(sum_c |ydata[c,1]|^2 (if smap not provided or isnothing(smap))

See equation [13] in M A Bernstein et al., "Reconstructions of Phase Contrast, Phased Array Multicoil Data", MRM 1994. https://doi.org/10.1002/mrm.1910320308

pc = phasecontrast(ydata)
fhat = phasecontrast(ydata, echotime)

Compute the phase contrast between two multicoil data sets. If echotime is provided, return a field map by converting from radians to Hz (if echotime is in seconds, as is typical).

In

• ydata (dims..., nc, ne) ne sets of complex images for nc ≥ 1 coils
• echotime::Echotime (ne = 2) echo time offsets

Out

• pc (dims...) phase contrast (or field map, if echotime is provided): ∠(sum_c ydata[c,1]' * ydata[c,2])

See equation [13] in M A Bernstein et al., "Reconstructions of Phase Contrast, Phased Array Multicoil Data", MRM 1994. https://doi.org/10.1002/mrm.1910320308

spdiff(dims::Dims; order=1, kwargs...)

Return Vector of length(dims) sparse finite-difference matrices of order order, one for each dimension.

Typically one will vcat the vector output to make a sparse finite-difference matrix suitable for the vec of a multi-dimensional array.

The kwargs are passed to spdiff1 (for order = 1) or spdiff2 (for order = 2). These functions are called once for each dimension. The options are ending and T.

Examples:

• spdiff((4,5,6))[1] == kron(I(6*5), spdiff1(4))
• spdiff((4,5,6); order=1)[2] == kron(I(6), spdiff1(5), I(4))
• spdiff((4,5,6); order=2)[3] == kron(spdiff2(6), I(4*5))

spdiff1(n::Int; ending::Symbol = :remove, T::Type = Int32)

Sparse n × n matrix for 1st-order finite differences.

Option

• T element type, default Int32 to save memory
• ending
  • :remove (default) remove first difference
  • :zero keep first row, akin to zero boundary conditions

spdiff2(n::Int; ending::Symbol = :remove, T::Type = Int32)

Sparse n × n matrix for 2nd-order finite differences.

Option

• T element type, default Int32 to save memory
• ending
  • :remove (default) remove first and last finite difference
  • :zero keep first and last rows, akin to zero boundary conditions
  • :first use 1st-order finite differences for first and last rows