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 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	Stream handle (unused on CPU-only builds, NULL).
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.
batchPipelineCtx	Lazily created batched pipeline context (see forwardBatched/adjointBatched).
cudaPipelineCtx	Batch size the batched context was created for (0 = not created).
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.
forwardBatched	Batched forward NUFFT: `batch` stacked images -> stacked k-space. Equivalent to applying operator to each slice of `d` independently, but on the Metal pipeline all slices are encoded into a constant number of command buffers with a single batched FFT (used by SENSE to fold all coils into one GPU round trip). Falls back to a per-slice loop when no batched GPU pipeline is available. `d` : Stacked input images, length batch * n1 (slice b occupies [bn1, (b+1)n1)). `batch` : Number of stacked slices. Return* : Stacked k-space vectors, length batch * n2.
adjointBatched	Batched adjoint NUFFT: `batch` stacked k-space vectors -> stacked images. Equivalent to applying operator/ to each slice of `d` independently; see forwardBatched for the batching semantics. `d` : Stacked input k-space, length batch * n2. `batch` : Number of stacked slices. Return : Stacked output images, length batch * n1.
forwardBatched	Batched forward NUFFT with gathered input: like forwardBatched, but the input slices live in separate vectors (e.g. persistent page-aligned scratch buffers). On the Metal pipeline the gather rides on the pipeline's host-to-GPU copy, so no extra stacking pass is needed. `slices` : One image vector of length n1 per batch slice. Return : Stacked k-space vectors, length slices.size() * n2.
adjointBatched	Batched adjoint NUFFT with scattered output: like adjointBatched, but each output image is written into the caller-provided vector outSlices[b] (each of length n1, e.g. persistent page-aligned scratch buffers usable for zero-copy Metal weighting afterwards). `d` : Stacked input k-space, length batch * n2. `batch` : Number of stacked slices; outSlices.size() must match. `outSlices` : Receives one image vector of length n1 per slice.
adjointBatched	Batched adjoint NUFFT with gathered input and scattered output: input k-space slices and output image slices both live in separate vectors (e.g. persistent page-aligned scratch buffers). `slices` : One k-space vector of length n2 per batch slice. `outSlices` : Receives one image vector of length n1 per slice; outSlices.size() must equal slices.size().
Gnufft	Copy constructor for thread-safe parallel use.
gridDim	Oversampled grid dimension for one axis: round(os * N). Round (not ceil) so that effective ratios reproduce their grid exactly despite float representation error.
gridDimEven	Grid dimension rounded up to even (required for x/y axes).
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).
getBatchPipeline	Return (lazily creating) a batched pipeline context for `batch` slices, or nullptr when the batched GPU path is unavailable (non-float type, no single-slice pipeline, batch <= 1, or allocation failure).
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).

batchPipelineCtx ¶

mutable MetalNufftPipelineContext* batchPipelineCtx

Lazily created batched pipeline context (see forwardBatched/adjointBatched).

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).

cudaPipelineCtx ¶

forge::cuda::CudaNufftPipeline* cudaPipelineCtx

Batch size the batched context was created for (0 = not created).
CUDA NUFFT pipeline context; nullptr when T1 != float.

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 ¶

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

Stream handle (unused on CPU-only builds, NULL).

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.

Gnufft(const Gnufft& other)

Copy constructor for thread-safe parallel use.
Shares read-only data (kx/ky/kz, LUT, trajectory) with the original,
but allocates independent scratch buffers for thread-safe operation.
The copy does NOT own kx/ky/kz/LUT and will not free them.

adjointBatched ¶

forgeCol<forgeComplex<T1>> adjointBatched(const forgeCol<forgeComplex<T1>>& d, uword batch) const

Batched adjoint NUFFT: `batch` stacked k-space vectors -> stacked images.

Equivalent to applying operator/ to each slice of d independently; see forwardBatched for the batching semantics.

d : Stacked input k-space, length batch * n2.

batch : Number of stacked slices.

Return : Stacked output images, length batch * n1.

void adjointBatched(const forgeCol<forgeComplex<T1>>& d, uword batch, std::vector<forgeCol<forgeComplex<T1>>>& outSlices) const

Batched adjoint NUFFT with scattered output: like adjointBatched, but
each output image is written into the caller-provided vector
outSlices[b] (each of length n1, e.g. persistent page-aligned scratch
buffers usable for zero-copy Metal weighting afterwards).

d : Stacked input k-space, length batch * n2.

batch : Number of stacked slices; outSlices.size() must match.

outSlices : Receives one image vector of length n1 per slice.

void adjointBatched(const std::vector<forgeCol<forgeComplex<T1>>>& slices, std::vector<forgeCol<forgeComplex<T1>>>& outSlices) const

Batched adjoint NUFFT with gathered input and scattered output: input
k-space slices and output image slices both live in separate vectors
(e.g. persistent page-aligned scratch buffers).

slices : One k-space vector of length n2 per batch slice.

outSlices : Receives one image vector of length n1 per slice; outSlices.size() must equal slices.size().

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.

forwardBatched ¶

forgeCol<forgeComplex<T1>> forwardBatched(const forgeCol<forgeComplex<T1>>& d, uword batch) const

Batched forward NUFFT: `batch` stacked images -> stacked k-space.

Equivalent to applying operator* to each slice of d independently, but on the Metal pipeline all slices are encoded into a constant number of command buffers with a single batched FFT (used by SENSE to fold all coils into one GPU round trip). Falls back to a per-slice loop when no batched GPU pipeline is available.

d : Stacked input images, length batch * n1 (slice b occupies [bn1, (b+1)n1)).

batch : Number of stacked slices.

Return : Stacked k-space vectors, length batch * n2.

forgeCol<forgeComplex<T1>> forwardBatched(const std::vector<forgeCol<forgeComplex<T1>>>& slices) const

Batched forward NUFFT with gathered input: like forwardBatched, but the
input slices live in separate vectors (e.g. persistent page-aligned
scratch buffers). On the Metal pipeline the gather rides on the
pipeline's host-to-GPU copy, so no extra stacking pass is needed.

slices : One image vector of length n1 per batch slice.

Return : Stacked k-space vectors, length slices.size() * n2.

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.

getBatchPipeline ¶

MetalNufftPipelineContext* getBatchPipeline(uword batch) const

Return (lazily creating) a batched pipeline context for `batch` slices,
or nullptr when the batched GPU path is unavailable (non-float type,
no single-slice pipeline, batch <= 1, or allocation failure).

gridDim ¶

static int gridDim(T1 os, uword N)

Oversampled grid dimension for one axis: round(os * N).
Round (not ceil) so that effective ratios reproduce their grid exactly
despite float representation error.

gridDimEven ¶

static int gridDimEven(T1 os, uword N)

Grid dimension rounded up to even (required for x/y axes).

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¶

batchPipelineCtx¶

beta¶

cudaPipelineCtx¶

gridData¶

gridNumElems¶

gridOS¶

imagXformedData¶

imagXformedImg¶

imageNumElems¶

ix¶

iy¶

iz¶

kernelWidth¶

kx¶

ky¶

kz¶

metalCtx¶

pGridData¶

pSamples¶

pipelineCtx¶

plan¶

realXformedData¶

realXformedImg¶

samples¶

stream¶

Operator Details¶

operator*¶

operator/¶

Function Details¶

Gnufft¶

adjointBatched¶

adjointSpatialInterp¶

forwardBatched¶

forwardSpatialInterp¶

getBatchPipeline¶

gridDim¶

gridDimEven¶

metalAdjointImpl¶

metalForwardImpl¶

~Gnufft¶

XformedData ¶

XformedImg ¶

batchPipelineCtx ¶

beta ¶

cudaPipelineCtx ¶

gridData ¶

gridNumElems ¶

gridOS ¶

imagXformedData ¶

imagXformedImg ¶

imageNumElems ¶

ix ¶

iy ¶

iz ¶

kernelWidth ¶

kx ¶

ky ¶

kz ¶

metalCtx ¶

pGridData ¶

pSamples ¶

pipelineCtx ¶

plan ¶

realXformedData ¶

realXformedImg ¶

samples ¶

stream ¶

Gnufft ¶

adjointBatched ¶

adjointSpatialInterp ¶

forwardBatched ¶

forwardSpatialInterp ¶

getBatchPipeline ¶

gridDim ¶

gridDimEven ¶

metalAdjointImpl ¶

metalForwardImpl ¶

~Gnufft ¶