The NUFFT

Non-Cartesian samples don't fall on the grid an FFT needs. Gnufft bridges the gap by gridding: it convolves the samples onto an oversampled Cartesian grid with a small kernel, so the FFT can do the heavy lifting.

Forward and adjoint

graph LR
    subgraph "Forward (image → k-space)"
        A[Image] --> B[Deapodize]
        B --> C[Zero-pad]
        C --> D[Oversampled FFT]
        D --> E[KB interpolation]
        E --> F[k-space samples]
    end

graph RL
    subgraph "Adjoint (k-space → image)"
        G[k-space samples] --> H[KB gridding]
        H --> I[Inverse FFT]
        I --> J[Crop]
        J --> K[Deapodize]
        K --> L[Image]
    end

The interpolation kernel is a Kaiser-Bessel function of width 4 grid points, tabulated in a look-up table. Its shape parameter is set from the width and the grid oversampling factor to minimize the worst-case interpolation error (Beatty et al.). Interpolating with the kernel multiplies the image by the kernel's Fourier transform, so the image is divided by that function, the deapodization, to undo it.

The FFT shifts are folded into the neighbouring steps rather than run as separate passes. k-space coordinates must lie in [−N/2, N/2) along each axis.

Grid oversampling

A bigger grid makes the kernel's aliasing smaller but costs memory and FFT time. forgeSense uses an oversampling factor of 2. pcSenseTimeSeg uses 1.6, which was the fastest factor that kept the operator within 1% of an exact DFT in FORGE's qualification tests (GridOSQualificationTests.cpp). Because the kernel shape is chosen for the oversampling factor, a lower factor costs less accuracy than it would with a fixed kernel.

Single-kz-plane mode

In multiband and stack-of-spirals acquisitions, each shot of a 3D volume samples only one kz value. A 3D NUFFT of such a shot spends most of its time on an oversampled z grid that the data barely touches, and the z gridding error is large when Nz is small: in a synthetic Nz = 4 test, the 3D path was 46% off an exact DFT at kz = −2.

When all of an operator's samples share one kz, Gnufft instead runs an exact DFT along z, collapsing the Nz image planes onto that kz value, followed by a 2D NUFFT. This is exact in z and, for a 4-band shot, about an eighth of the FFT and gridding work. Operators with more than one kz value, such as EMPIRE's, are unchanged. FORGE_NUFFT_PLANE_MODE=0 turns it off for comparison.

Density compensation

The iterative reconstructions don't need density compensation: the solver accounts for the sampling density through the forward model. Direct gridding does need it, and gridCoilImages and the routines in IO/directRecon.h weight samples with calc3DDensityCompensation(), the fixed-point iteration of Pipe and Menon:

\[W_{n+1} = \frac{W_n}{C\,C^H W_n}\]

where C^H spreads samples onto the oversampled grid with the Kaiser-Bessel kernel and C interpolates back, with no FFT. gridCoilImages --dcf-iterations sets the number of iterations.

Further reading

• J. I. Jackson, C. H. Meyer, D. G. Nishimura and A. Macovski, "Selection of a convolution function for Fourier inversion using gridding," IEEE Transactions on Medical Imaging, 10(3), 1991.
• P. J. Beatty, D. G. Nishimura and J. M. Pauly, "Rapid gridding reconstruction with a minimal oversampling ratio," IEEE Transactions on Medical Imaging, 24(6), 2005.
• J. G. Pipe and P. Menon, "Sampling density compensation in MRI: rationale and an iterative numerical solution," Magnetic Resonance in Medicine, 41(1), 1999.