MRIFieldmaps overview

This page summarizes the Julia package MRIFieldmaps.

This page comes from a single Julia file: 01-overview.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: 01-overview.ipynb, or open it in binder here: 01-overview.ipynb.

Setup

Packages needed here. Use Pkg.add as illustrated here when using a package for the first time.

using MRIFieldmaps: spdiff
using MIRTjim: jim, prompt; jim(:prompt, true)
using InteractiveUtils: versioninfo

The following line is helpful when running this file as a script; this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

Overview

This package provides algorithms for estimating fieldmaps in MRI. Currently it supports B0 map estimation. A future extension could support B1+ map estimation.

This page just discusses the regularization; other examples illustrate specific fieldmap estimators.

Regularization

Some methods in this package use sparse matrices to perform finite-difference operations. Those operations could be performed with a LinearMap based on diff, but the fastest algorithms here use preconditioning based on incomplete Cholesky factorization of a Hessian matrix, and such factorization is supported for sparse matrices, but probably not for general linear maps.

It may be helpful to visualize these sparse matrices.

Here are 1D 1st-order and 2nd-order finite differences.

D1 = spdiff((8,); order=1)[1]
D2 = spdiff((8,); order=2)[1]
clim = (-1, 2)
color = :cividis
jim(
 jim(Matrix(D1)', "1st-order"; color, clim, prompt=false),
 jim(Matrix(D2)', "2nd-order"; color, clim, prompt=false),
)

There are other boundary conditions available:

D1z = spdiff((8,); order=1, ending=:zero)[1]
D2z = spdiff((8,); order=2, ending=:zero)[1]
D2f = spdiff((8,); order=2, ending=:first)[1]
jim(
 jim(Matrix(D1z)', "1st-order, zero ends"; color, clim, prompt=false),
 jim(Matrix(D2z)', "2nd-order, zero ends"; color, clim, prompt=false),
 jim(Matrix(D2f)', "2nd-order, 1st ends"; color, clim, prompt=false),
 ; layout = (3,1),
)

For multi-dimensional arrays, regularizers need finite-differences along each dimension, and these are constructed with Kronecker products and are applied to the vec of an array.

dims = (9,8)
d1 = spdiff(dims; order=1)
d2 = spdiff(dims; order=2)
jim(
 jim(Matrix(d1[1])', "1st-order, 1st dim"; color, clim, prompt=false),
 jim(Matrix(d1[2])', "1st-order, 2nd dim"; color, clim, prompt=false),
 jim(Matrix(d2[1])', "2nd-order, 1st dim"; color, clim, prompt=false),
 jim(Matrix(d2[2])', "2nd-order, 2nd dim"; color, clim, prompt=false),
)

Here is an illustration of applying these finite-difference matrices to a simple test phantom. Note the use of vec and reshape for display.

dims = (40,30)
x = LinRange(-1, 1, dims[1])
y = LinRange(-1, 1, dims[2])
phantom = @. abs(x) + abs(y') < 0.5
pl = Array{Any}(undef, 2, 2)
for order in 1:2, d in 1:2
    sp = spdiff(dims; order)
    dif = reshape(sp[d] * vec(phantom), dims)
    pl[d,order] = jim(x, y, dif, "dim$d differences\norder=$order")
end
pp = jim(x, y, phantom, "Test image")
jim([[pp pp]; pl]...)

Support mask

Often we want to estimate a fieldmap over some spatial support "mask" that is smaller than the entire image, e.g., only in voxels where the signal is sufficiently large. See the ImageGeoms.jl documentation about the related embed and maskit operations.

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

12-element Vector{SubString{String}}:
 "Julia Version 1.10.3"
 "Commit 0b4590a5507 (2024-04-30 10:59 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()

Status `~/work/MRIFieldmaps.jl/MRIFieldmaps.jl/docs/Project.toml`
  [e30172f5] Documenter v1.4.1
  [98b081ad] Literate v2.18.0
  [23992714] MAT v0.10.6
  [7035ae7a] MIRT v0.17.0
  [170b2178] MIRTjim v0.24.0
  [5af58b90] MRIFieldmaps v0.0.4 `~/work/MRIFieldmaps.jl/MRIFieldmaps.jl`
  [91a5bcdd] Plots v1.40.4
  [2913bbd2] StatsBase v0.34.3
  [1986cc42] Unitful v1.20.0
  [42071c24] UnitfulRecipes v1.6.1
  [b77e0a4c] InteractiveUtils
  [9a3f8284] Random

This page was generated using Literate.jl.