NUFFT Overview

This example illustrates how to use Nonuniform FFT (NUFFT) for image reconstruction in MRI using the Julia language.

This page comes from a single Julia file: 1-nufft.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: 1-nufft.ipynb, or open it in binder here: 1-nufft.ipynb.

Some MRI scans use non-Cartesian sampling patterns (like radial and spiral k-space trajectories, among others), often in the interest of acquisition speed.

Image reconstruction from fully sampled Cartesian k-space data typically uses simple inverse FFT operations, whereas non-Cartesian sampling requires more complicated methods.

The examples here illustrate both non-iterative (aka "gridding") and iterative methods for non-Cartesian MRI reconstruction. For simplicity the examples consider the case of single-coil data and ignore the effects of B0 field inhomogeneity.

First we add the Julia packages that are need for these examples. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "ImagePhantoms"
        "Unitful"
        "Plots"
        "LaTeXStrings"
        "MIRTjim"
        "MIRT"
        "InteractiveUtils"
    ])
end

Now tell this Julia session to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using ImagePhantoms: shepp_logan, SheppLoganEmis, spectrum, phantom #, Gauss2
using Unitful: mm # Allows use of physical units (mm here)
using Plots; default(label="", markerstrokecolor=:auto)
using LaTeXStrings # for LaTeX in plot labels, e.g., L"\alpha_n"
using MIRTjim: jim, prompt # jiffy image display
using MIRT: Anufft, diffl_map, ncg
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 image is displayed.

isinteractive() && jim(:prompt, true);

Radial k-space sampling

We focus on radial sampling as a simple representative non-Cartesian case. Consider imaging a 256mm × 256mm field of FOV with the goal of reconstructing a 128 × 128 pixel image. The following radial and angular k-space sampling is reasonable.

N = 128
FOV = 256mm # physical units!
Δx = FOV / N # pixel size
kmax = 1 / 2Δx
kr = ((-N÷2):(N÷2)) / (N÷2) * kmax # radial sampling in k-space
Nϕ = 3N÷2 # theoretically should be about π/2 * N
kϕ = (0:Nϕ-1)/Nϕ * π # angular samples
νx = kr * cos.(kϕ)' # N × Nϕ k-space sampling in cycles/mm
νy = kr * sin.(kϕ)'
plot(νx, νy,
    xlabel=L"\nu_x", ylabel=L"\nu_y",
    aspect_ratio = 1,
    title = "Radial k-space sampling",
)

isinteractive() && prompt();

For the NUFFT routines considered here, the sampling arguments must be "Ω" values: digital frequencies have pseudo-units of radians / pixel.

Ωx = (2π * Δx) * νx # N × Nϕ grid of k-space sample locations
Ωy = (2π * Δx) * νy # in pseudo-units of radians / sample

scatter(Ωx, Ωy,
    xlabel=L"\Omega_x", ylabel=L"\Omega_y",
    xticks=((-1:1)*π, ["-π", "0", "π"]),
    yticks=((-1:1)*π, ["-π", "0", "π"]),
    xlims=(-π,π) .* 1.1,
    ylims=(-π,π) .* 1.1,
    aspect_ratio = 1, markersize = 0.5,
    title = "Radial k-space sampling",
)