Gnufft¶

template <typename T1> class Gnufft

Non-uniform Fast Fourier Transform operator using Kaiser-Bessel gridding.

Implements an efficient NUFFT via convolution with a Kaiser-Bessel (KB) kernel onto an oversampled Cartesian grid, followed by an FFT. Supports both CPU and GPU execution (via OpenACC or Apple Metal). An LUT (look-up table) is used for fast kernel evaluation.

Forward algorithm (image → k-space):

graph LR
    A["$$\mathbf{x}$$ (image)"] --> B["$$\mathbf{x} / \hat{\kappa}$$ (deapodize)"]
    B --> C["Zero-pad to $$N \cdot \text{gridos}$$"]
    C --> D["FFT"]
    D --> E["$$\kappa$$ interpolation"]
    E --> F["$$\mathbf{d}$$ (k-space)"]

Adjoint algorithm (k-space → image):

graph RL
    G["$$\mathbf{d}$$ (k-space)"] --> H["$$\kappa$$ gridding"]
    H --> I["IFFT"]
    I --> J["Crop to $$N$$"]
    J --> K["$$/ \hat{\kappa}$$ (deapodize)"]
    K --> L["$$\hat{\mathbf{x}}$$ (image)"]

Density compensation: Non-uniform k-space sampling requires density compensation weights (DCF) to correct for variable sampling density. The Pipe iterative algorithm computes these weights using the gridding/degridding pair:

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

where \(G^H G\) is the forward-then-adjoint gridding operation. This is implemented in griddingSupport.h via calc3DDensityCompensation().

See : Pipe & Menon, "Sampling Density Compensation in MRI," MRM, 1999. https://doi.org/10.1002/(SICI)1522-2594(199901)41:1<179::AID-MRM25>3.0.CO;2-V

The Kaiser-Bessel kernel is:

\[\kappa(u) = \frac{1}{W} I_0\!\left(\beta\sqrt{1 - (2u/W)^2}\right)\]

for \(|u| \leq W/2\), zero otherwise. \(W\) is the kernel width (kernelWidth), \(\beta\) is the shape parameter, and \(I_0\) is the modified Bessel function of the first kind.

Parameter guidance: - gridos (grid oversampling factor): typically 2.0. Higher values reduce aliasing artifacts at the cost of memory and computation. - beta is computed automatically from gridos and kernelWidth to minimize the worst-case interpolation error.

k-space coordinates must lie within \([-N/2, N/2)\) in each dimension.

Use G * x for the forward transform and G / d for the adjoint.

See : Gdft for the exact (but slower) field-corrected DFT.

See : gridding.h for the low-level gridding and degridding routines.

See : Jackson et al., "Selection of a Convolution Function for Fourier Inversion Using Gridding," IEEE TMI, 1991. https://doi.org/10.1109/42.75611

See : Pipe & Menon, "Sampling Density Compensation in MRI," MRM, 1999. https://doi.org/10.1002/(SICI)1522-2594(199901)41:1<179::AID-MRM25>3.0.CO;2-V

See : Beatty et al., "Rapid Gridding Reconstruction with a Minimal Oversampling Ratio," IEEE TMI, 2005. https://doi.org/10.1109/TMI.2004.842452

Gnufft<float> G(dataLength, 2.0f, Nx, Ny, Nz, kx, ky, kz, ix, iy, iz);
Col<cx_float> kdata = G * image;  // forward: gridding + FFT
Col<cx_float> recon = G / kdata;  // adjoint: FFT + degridding

T1 : Floating-point precision type (float or double).

Variables¶

Name	Description
kx	Number of image pixels (n1 = Nx * Ny * Nz).
ky	k-space y-coordinates (length n2, range [-Ny/2, Ny/2)).
kz	k-space z-coordinates (length n2, range [-Nz/2, Nz/2)).
ix	Image-space x-coordinates (length n1, stored for reference).
iy	Image-space y-coordinates (length n1, stored for reference).
iz	Image-space z-coordinates (length n1, stored for reference).
gridOS	Grid oversampling factor.
beta	Kaiser-Bessel kernel look-up table.
kernelWidth	Kaiser-Bessel kernel half-width in grid units.
stream	GPU stream handle (OpenACC/CUDA, NULL on CPU builds).
plan	cuFFT plan handle for GPU FFT on the oversampled grid.
plan	FFT plan handle (unused on CPU-only builds).
metalCtx	Metal gridding context (legacy); nullptr when T1 != float or Metal unavailable.
pipelineCtx	Full GPU NUFFT pipeline context (macOS 14+); nullptr if unavailable.
imageNumElems	Total number of elements in the (non-oversampled) image grid.
gridNumElems	Total number of elements in the oversampled grid.
pGridData	Real-valued grid buffers (image-size and oversampled) for host and device.
pSamples	Real-valued interleaved k-space sample buffer (length 2*n2).
gridData	Complex grid buffers (image-size and oversampled) for host and device.
samples	Complex k-space sample buffer (length n2).
XformedData	Scratch buffer for transformed k-space data (forward path).
XformedImg	Scratch buffer for transformed image data (adjoint path).
realXformedData	Real part of transformed k-space data (split-complex layout).
imagXformedData	Imaginary part of transformed k-space data (split-complex layout).
realXformedImg	Real part of transformed image data (split-complex layout).
imagXformedImg	Imaginary part of transformed image data (split-complex layout).

Operators¶

Name	Description
operator*	Forward NUFFT: image -> k-space via KB gridding. Applies deapodization, zero-pads to the oversampled grid, FFT-shifts, performs an FFT, and interpolates at the non-uniform k-space locations. `d` : Input image vector of length `n1.` Return : Output k-space vector of length `n2.`
operator/	Adjoint NUFFT: k-space -> image via KB gridding. Spreads k-space samples onto the oversampled grid using the KB kernel, performs an inverse FFT, crops to image size, and applies deapodization. `d` : Input k-space vector of length `n2.` Return : Output image vector of length `n1.`
operator*	Forward NUFFT (forgeCol overload for Metal path).
operator/	Adjoint NUFFT (forgeCol overload for Metal path).

Functions¶

Name	Description
Gnufft	Default constructor.
Gnufft	Construct a Gnufft operator with the given trajectory and image dimensions.
~Gnufft	Destructor.
metalForwardImpl	Metal forward pipeline: deapodize → zero-pad → FFT → gridding. Writes result to pSamples (n2 complex elements).
metalAdjointImpl	Metal adjoint pipeline: gridding → IFFT → crop → deapodize. Writes result to pGridData (n1 complex elements).
forwardSpatialInterp	Forward spatial interpolation only (no deapodization or FFT).
adjointSpatialInterp	Adjoint spatial interpolation only (no deapodization or FFT).

Variable Details¶

XformedData ¶

mutable Col<CxT1> XformedData

Scratch buffer for transformed k-space data (forward path).

XformedImg ¶

mutable Col<CxT1> XformedImg

Scratch buffer for transformed image data (adjoint path).

beta ¶

T1 beta

Kaiser-Bessel kernel look-up table.
Number of entries in the KB kernel LUT.
KB kernel shape parameter (computed from gridOS and kernelWidth).

gridData ¶

complex<T1>*gridData, *gridData_d, *gridData_os, *gridData_os_d

Complex grid buffers (image-size and oversampled) for host and device.

gridNumElems ¶

uword gridNumElems

Total number of elements in the oversampled grid.

gridOS ¶

T1 gridOS

Grid oversampling factor.

imagXformedData ¶

mutable Col<T1> imagXformedData

Imaginary part of transformed k-space data (split-complex layout).

imagXformedImg ¶

mutable Col<T1> imagXformedImg

Imaginary part of transformed image data (split-complex layout).

imageNumElems ¶

uword imageNumElems

Total number of elements in the (non-oversampled) image grid.

ix ¶

Col<T1> ix

Image-space x-coordinates (length n1, stored for reference).

iy ¶

Col<T1> iy

Image-space y-coordinates (length n1, stored for reference).

iz ¶

Col<T1> iz

Image-space z-coordinates (length n1, stored for reference).

kernelWidth ¶

T1 kernelWidth

Kaiser-Bessel kernel half-width in grid units.

kx ¶

T1* kx

Number of image pixels (n1 = Nx * Ny * Nz).
Number of k-space samples.
Image x-dimension in pixels.
Image y-dimension in pixels.
Image z-dimension in pixels.
k-space x-coordinates (length n2, range [-Nx/2, Nx/2)).

ky ¶

T1* ky

k-space y-coordinates (length n2, range [-Ny/2, Ny/2)).

kz ¶

T1* kz

k-space z-coordinates (length n2, range [-Nz/2, Nz/2)).

metalCtx ¶

MetalGriddingContext* metalCtx

Metal gridding context (legacy); nullptr when T1 != float or Metal unavailable.

pGridData ¶

T1 *pGridData, *pGridData_d, *pGridData_os, *pGridData_os_d

Real-valued grid buffers (image-size and oversampled) for host and device.

pSamples ¶

T1* pSamples

Real-valued interleaved k-space sample buffer (length 2*n2).

pipelineCtx ¶

MetalNufftPipelineContext* pipelineCtx

Full GPU NUFFT pipeline context (macOS 14+); nullptr if unavailable.

plan ¶

cufftHandle plan

cuFFT plan handle for GPU FFT on the oversampled grid.

plan ¶

void* plan

FFT plan handle (unused on CPU-only builds).

realXformedData ¶

mutable Col<T1> realXformedData

Real part of transformed k-space data (split-complex layout).

realXformedImg ¶

mutable Col<T1> realXformedImg

Real part of transformed image data (split-complex layout).

samples ¶

complex<T1>* samples

Complex k-space sample buffer (length n2).

stream ¶

void* stream

GPU stream handle (OpenACC/CUDA, NULL on CPU builds).

Operator Details¶

operator*¶

Col<CxT1> operator*(const Col<CxT1>& d) const

Forward NUFFT: image -> k-space via KB gridding.

Applies deapodization, zero-pads to the oversampled grid, FFT-shifts, performs an FFT, and interpolates at the non-uniform k-space locations.

d : Input image vector of length n1.

Return : Output k-space vector of length n2.

forgeCol<forgeComplex<T1>> operator*(const forgeCol<forgeComplex<T1>>& d) const

Forward NUFFT (forgeCol overload for Metal path).

d : Input image vector of length n1.

Return : Output k-space vector of length n2.

operator/¶

Col<CxT1> operator/(const Col<CxT1>& d) const

Adjoint NUFFT: k-space -> image via KB gridding.

Spreads k-space samples onto the oversampled grid using the KB kernel, performs an inverse FFT, crops to image size, and applies deapodization.

d : Input k-space vector of length n2.

Return : Output image vector of length n1.

forgeCol<forgeComplex<T1>> operator/(const forgeCol<forgeComplex<T1>>& d) const

Adjoint NUFFT (forgeCol overload for Metal path).

d : Input k-space vector of length n2.

Return : Output image vector of length n1.

Function Details¶

Gnufft ¶

Gnufft()

Default constructor. Produces an uninitialized operator.

Gnufft(uword dataLength, T1 gridos, uword nx, uword ny, uword nz, const Col<T1>& k1, const Col<T1>& k2, const Col<T1>& k3, const Col<T1>& i1, const Col<T1>& i2, const Col<T1>& i3)

Construct a Gnufft operator with the given trajectory and image dimensions.

dataLength : Number of k-space samples (n2, output of forward transform).

gridos : Grid oversampling factor (typically 2.0).

nx : Image x-dimension in pixels.

ny : Image y-dimension in pixels.

nz : Image z-dimension in pixels (use 1 for 2-D).

k1 : k-space x-coordinates, length dataLength (units: cycles, range [-nx/2, nx/2)).

k2 : k-space y-coordinates, length dataLength.

k3 : k-space z-coordinates, length dataLength.

i1 : Image-space x-coordinates, length nxnynz (not used in gridding, stored for reference).

i2 : Image-space y-coordinates, length nxnynz.

i3 : Image-space z-coordinates, length nxnynz.

adjointSpatialInterp ¶

Col<CxT1> adjointSpatialInterp(const Col<CxT1>& d) const

Adjoint spatial interpolation only (no deapodization or FFT).

d : Input k-space vector of length n2.

Return : Output image vector of length n1.

forwardSpatialInterp ¶

Col<CxT1> forwardSpatialInterp(const Col<CxT1>& d) const

Forward spatial interpolation only (no deapodization or FFT).

d : Input image vector of length n1.

Return : Spatially interpolated k-space vector of length n2.

metalAdjointImpl ¶

void metalAdjointImpl(const T1* dataPtr) const

Metal adjoint pipeline: gridding → IFFT → crop → deapodize.
Writes result to pGridData (n1 complex elements).

metalForwardImpl ¶

void metalForwardImpl(const T1* dataPtr) const

Metal forward pipeline: deapodize → zero-pad → FFT → gridding.
Writes result to pSamples (n2 complex elements).

~Gnufft ¶

~Gnufft()

Destructor. Frees the KB kernel LUT and all internal grid buffers.

Gnufft¶

Variables¶

Operators¶

Functions¶

Variable Details¶

XformedData¶

XformedImg¶

beta¶

gridData¶

gridNumElems¶

gridOS¶

imagXformedData¶

imagXformedImg¶

imageNumElems¶

ix¶

iy¶

iz¶

kernelWidth¶

kx¶

ky¶

kz¶

metalCtx¶

pGridData¶

pSamples¶

pipelineCtx¶

plan¶

plan¶

realXformedData¶

realXformedImg¶

samples¶

stream¶

Operator Details¶

operator*¶

operator/¶

Function Details¶

Gnufft¶

adjointSpatialInterp¶

forwardSpatialInterp¶

metalAdjointImpl¶

metalForwardImpl¶

~Gnufft¶

XformedData ¶

XformedImg ¶

beta ¶

gridData ¶

gridNumElems ¶

gridOS ¶

imagXformedData ¶

imagXformedImg ¶

imageNumElems ¶

ix ¶

iy ¶

iz ¶

kernelWidth ¶

kx ¶

ky ¶

kz ¶

metalCtx ¶

pGridData ¶

pSamples ¶

pipelineCtx ¶

plan ¶

plan ¶

realXformedData ¶

realXformedImg ¶

samples ¶

stream ¶

Gnufft ¶

adjointSpatialInterp ¶

forwardSpatialInterp ¶

metalAdjointImpl ¶

metalForwardImpl ¶

~Gnufft ¶