This example illustrates how to use Nonuniform FFT (NUFFT) for image reconstruction in MRI using the Julia language.
This entire page was generated using a single Julia file: 1-nufft.jl. In any such Julia documentation, you can access the source code using the "Edit on GitHub" link in the top right.
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
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" "Plots" "Unitful" "UnitfulRecipes" "LaTeXStrings" "MIRTjim" "MIRT" "InteractiveUtils" ]) end
Now tell this Julia session to use the following packages for this example. Run
Pkg.add() in the preceding code block first, if needed.
using ImagePhantoms: shepp_logan, SheppLoganEmis, spectrum, phantom #, Gauss2 using Plots; default(label="", markerstrokecolor=:auto) using Unitful: mm # Allows use of physical units (mm here) using UnitfulRecipes # Ensures that plots show appropriate units when applicable 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 jl-file as a script; this way it will prompt user to hit a key after each image is displayed.
isinteractive() && jim(:prompt, true);
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", )