The NUFFT Operator

This document gives a general overview of the Non-Uniform Fast fourier transform (NUFFT), and its application to MRI.

The Non-Uniform Discrete Fourier Transform

The Non-Uniform Discrete Fourier Transform (NUDFT) is a generalization of the Discrete Fourier Transform (DFT) to non-uniform sampling.

For a signal \(x\) sampled at location \(p_0, p_1, \ldots, p_{N-1}\) we want to get the frequency points (non uniformly spaced) at \(\nu_0, \nu_1, \ldots, \nu_{M-1}\)

The 1D-NUDFT [1] is defined as:

\[X_k = \sum_{n=0}^{N-1} x(p_n) e^{ -2\imath \pi p_n \nu_k}\]

where \(X_k\) is the frequency point at \(\nu_k\).

The multidimensional cases are derived by using vectorized location e.g \(\boldsymbol{p}_n\) and \(\boldsymbol{\nu_k}\) as for the classical DFT.

There exist 3 types of NUDFT:

• Type 1: \(p_n = n/N\) and \(\nu_k\) are non-uniformly spaced

• Type 2: \(p_n\) are non uniformly spaced and \(\nu_k = k/M\)

• Type 3: \(p_n\) and \(\nu_k\) are non uniformly spaced

• If \(p_n=n/N\) and \(\nu_k=k/M\) then the NUDFT is simply the Discrete Fourier Transform.

The naive implementation of the NUDFT is \(O(NM)\). This becomes problematic for large N and M. The NUFFT is a fast algorithm to compute the NUDFT. The NUFFT is a generalization of the Fast Fourier Transform (FFT) to non-uniform sampling. The underlying principles of the NUFFT algorithm are described briefly down in the document

Application in MRI

In Magnetic Resonance Imaging (MRI) the raw data is acquired in the k-space, ideally corresponding to the Fourier domain. Traditional sampling schemes of the k-space usually consist of acquired lines in a specific direction, in a Cartesian fashion.

In order to accelerate the acquisition of k-space data, one may use a non-Cartesian (i.e. non-uniformly distributed in every direction) sampling scheme as the latter offers increased sampling efficiency, i.e. broader k-space coverage in a given time period.

The acquisition model is usually described as:

\[y(\boldsymbol{\nu}_i) = \int_{\mathbb{R}^d} x(\boldsymbol{u}) e^{-2\imath\pi \boldsymbol{u} \cdot \boldsymbol{\nu_i}} d\boldsymbol{u} + n_i, \quad i=1,\dots,M\]

Where:

• \(x(\boldsymbol{u})\) is the spatially varying image contrast acquired.

• \(y_1, \dots, y_M\) are the sampled points at frequency locations \(\Omega=\lbrace \boldsymbol{\nu}_1, \dots, \boldsymbol{\nu}_M \in [-1/2, 1/2]^d\rbrace\). Typically images (\(d=2\)) or volumes (\(d=3\)) are acquired.

• \(n_i\) is a zero-mean complex-valued Gaussian Noise, modeling the “thermal noise” of the scanner.

In practice the equation above is discretized, and the integral is replaced by a sum, using a finite number of samples \(N\):

\[y(\boldsymbol{\nu}_i) = \sum_{j=1}^N x(\boldsymbol{u}_j) e^{-2\imath\pi\boldsymbol{u}_j\cdot\boldsymbol{\nu_i}} + n_i, \quad i=1,\dots,M\]

Where \(\boldsymbol{u}_j\) are the \(N\) spatial locations of image voxels. This is stated using the operator notation:

\[\boldsymbol{y} = \mathcal{F}_\Omega (\boldsymbol{x}) + \boldsymbol{n}\]

As the sampling locations \(\Omega\) are non-uniform and the image locations \(\boldsymbol{u}_j\) are uniform, \(\mathcal{F}_\Omega\) is a NUDFT operator, and the equation above describe a Type 2 NUDFT. Similarly the adjoint operator is:

Correspondance Table between NUFFT and MRI acquisition model.

NUFFT Type	Operation	MRI Transform	Operator
Type 1	Adjoint	Kspace → Image	\(\mathcal{F}_\Omega^*\)
Type 2	Forward	Image → Kspace	\(\mathcal{F}_\Omega\)

Attention

In order to reconstruct \(x\) from \(y\), one has to solve the inverse problem, stated usually as

\[\hat{x} = \arg\min_x \frac{1}{2} \|\mathcal{F}_\Omega(\boldsymbol{x}) - \boldsymbol{y}\|_2^2 + g(\boldsymbol{x}).\]

This package focuses solely on computing \(\mathcal{F}_\Omega\boldsymbol{x}\) or \(\mathcal{F}_\Omega^**\boldsymbol{y}\). solving this problem is not addressed here, but you can check pysap-mri for this purpose.

Extension of the Acquisition model

The MRI acquisition model can be extended in two main ways. First by taking into account Parallel Imaging, where multiple coils are receiving data, each with a dedicated sensitivity profile.

Tip

MRI-NUFFT provides the FourierOperator interface to implement all the physical model described below. See MRI-NUFFT Interfaces Convention for the standard, and FourierOperatorBase

Parallel Imaging Model

In MRI the acquired signal can be received by multiple antennas (“coils”). Each coil possesses a specific sensitivity profile (i.e. each sees the object differently due to its physical layout).

The acquisition model for parallel imaging with \(L\) coils is:

\[y_{i,\ell} = \int_{\mathbb{R}^d} S_\ell(\boldsymbol{u})x(\boldsymbol{u}) e^{-2i\pi \boldsymbol{u} \cdot \boldsymbol{\nu_i}} d\boldsymbol{u} + n_{i,\ell}\]

Or using the operator notation:

\[\begin{split}\tilde{\boldsymbol{y}} = \begin{bmatrix} \mathcal{F}_\Omega S_1 \\ \vdots \\ \mathcal{F}_\Omega S_L \\ \end{bmatrix} \boldsymbol{x} + \boldsymbol{n}_\ell = \mathcal{F}_\Omega S \otimes \boldsymbol{x} + \tilde{\boldsymbol{n}}\end{split}\]

Where \(S_1, \dots, S_L\) are the sensitivity maps of each coil. Such sensitivity maps can be acquired separetely by acquiring low frequency of the kspace, or estimated from the data.

Off-resonance correction model

The constant magnetic field applied in a MRI machine \(B0\) (with a typical intensity 1.5, 3 or 7 Tesla) is inherently disturbed at tissue interfaces with owing to different magnetic susceptibilities (such as air-tissue interfaces in the nose and ear canals). Those field perturbations introduce a spatially varying phase shift in the frequencies acquired (noted \(\Delta\omega_0\)), making the acquisition model deviating from the convenient Fourier model. Fortunately, this inhomogeneity map can be acquired separatly or estimated (but this goes beyond the scope of this package) and integrated in the model as:

\[y(t_m) = \int_{\mathbb{R}^d} x(\boldsymbol{u}) e^{-2\imath\pi \boldsymbol{u} \cdot\boldsymbol{\nu_m} + \Delta\omega(\boldsymbol{u}) t_m} d\boldsymbol{u}\]

where \(t_i\) is the time at which the frequency \(\nu_m\) is acquired. Similarly at the reconstruction we have

\[x(\boldsymbol{u_n}) = \sum_{m}^M y(t_m) e^{2\imath\pi \boldsymbol{u} \cdot \boldsymbol{\nu_m}} e^{i\Delta\omega(\boldsymbol{u_n}) t_m}\]

With these mixed-domain field pertubations, the Fourier model does not hold anymore and the (NU)FFT algorithm cannot be used any longer. The main approach (initially proposed by Noll et al. [2]) is to approximate the mixed-domain exponential term by splitting it into single-domain weights \(b_{m, \ell}\) and \(c_{\ell, n}\):

\[e^{i\Delta\omega(\boldsymbol{u_n}) t_m} = \sum_{\ell=1}^L b_{m, \ell}\, c_{\ell, n}\]

Yielding the following model, where \(L \ll M, N\) regular Fourier transforms are performed to approximate the non-Fourier transform.

\[x(\boldsymbol{u_n}) = \sum_{\ell=1}^L c_{\ell, n} \sum_{m}^M y(t_m) b_{m, \ell} e^{2\imath\pi \boldsymbol{u} \cdot \boldsymbol{\nu_i}}\]

The coefficients \(B=(b_{m, \ell}) \in \mathbb{C}^{M\times L}\) and \(C=(c_\ell, n) \in \mathbb{C}^{L\times N}\) can be (optimally) estimated for any given \(L\) by solving the following matrix factorisation problem [4]:

\[\hat{B}, \hat{C} = \arg\min_{B,C} \| E- BC\|_{fro}^2\]

Where \(E_{mn} = e^{i\Delta\omega_0(u_n)t_m}\).

Note

The estimation of the B and C methods are provided in the mrinufft.extras.field_map module. Other methods like MTI [5] and MFI [6] are also available.

Tip

You can use the method .with_off_resonance_correction to augment an existing operator with off-resonance correction capability.

Subspace Projection Model

In several MRI applications, such as dynamic or quantitative MRI, a single acquisition provides a stack of two- or three-dimensional images, each representing a single time frame (for dynamic MRI) or a single contrast (for quantitative MRI). To achieve a clinically feasible scan time, each frame or contrast is acquired with a different aggressively undersampled k-space trajectory. In this context, the single-coil acquisition model becomes:

\[\begin{split}\tilde{\boldsymbol{y}} = \begin{bmatrix} \mathcal{F}_{\Omega_1} & 0 & \cdots & 0 \\ 0 & \mathcal{F}_{\Omega_2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathcal{F}_{\Omega_T} \end{bmatrix} \boldsymbol{x} + \boldsymbol{n}\end{split}\]

where \(\mathcal{F}_{\Omega_1}, \dots, \mathcal{F}_{\Omega_T}\) are the Fourier operators corresponding to each individual frame. Some applications (e.g., MR Fingerprinting [3]) may consists of thousands of total frames \(T\), leading to repeated Fourier Transform operations and high computational burden. However, the 1D signal series arising from similar voxels, e.g., with similar relaxation properties, are typically highly correlated. For this reason, the image series can be represented as:

\[\boldsymbol{x} = \Phi\Phi^H \boldsymbol{x}\]

where \(\Phi\) is an orthonormal basis spanning a low dimensional subspace whose rank \(K \ll T\) which can be obtained performing a Singular Value Decomposition of a low resolution fully sampled training dataset or an ensemble of simulated Bloch responses. The signal model can be then written as:

\[\begin{split}\tilde{\boldsymbol{y}} = \begin{bmatrix} \mathcal{F}_{\Omega_1} & 0 & \cdots & 0 \\ 0 & \mathcal{F}_{\Omega_2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathcal{F}_{\Omega_T} \end{bmatrix} \Phi \Phi^H \boldsymbol{x} + \boldsymbol{n} = \begin{bmatrix} \mathcal{F}_{\Omega_1} & 0 & \cdots & 0 \\ 0 & \mathcal{F}_{\Omega_2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathcal{F}_{\Omega_T} \end{bmatrix} \Phi \boldsymbol{\alpha} + \boldsymbol{n}\end{split}\]

where \(\boldsymbol{\alpha} = \Phi^H \boldsymbol{x}\) are the spatial coefficients representing the image series. Since the elements of \(\Phi^H\) do not depend on the specific k-space frequency points, the projection operator \(\boldsymbol{\Phi}\) commutes with the Fourier transform, and the signal equation finally becomes:

\[\tilde{\boldsymbol{y}} = \Phi \mathcal{F}_\Omega(\boldsymbol{\alpha}) + \boldsymbol{n}\]

that is, computation now involves \(K \ll T\) Fourier Transform operations, each with the same sampling trajectory, which can be computed by levaraging efficient NUFFT implementations for conventional static MRI.

Stacked NUFFT

If the k-space trajectory consists of a stacked of equally (or a subsampling of) spaced 2D planes of the 3D k-space, the NUFFT operator can be optimized by performing A 2D NUFFT on each plane, followed by a 1D FFT along the third dimension, resulting in a 2.5D NUFFT operator, lowering the computational cost and memory footprint.

Note

You can use the stacked nufft operator by using a stacked-* backend, and provide a 3D stacked trajectory. See mrinufft.operators.stacked for more details.

The Non Uniform Fast Fourier Transform in practice

In order to lower the computational cost of the Non-Uniform Fourier Transform, the main idea is to move back to a regular grid where an FFT would be performed (going from a typical \(O(MN)\) complexity to \(O(M\log(N))\)).

The three steps of the type 1 NUFFT (non-uniform to uniform) are:

1. Spreading

Each non-uniform source point \(x(p_n)\) is convolved with a compactly supported spreading kernel \(\phi\) and accumulated onto an oversampled Cartesian grid of size \(\sigma N\) (with oversampling factor \(\sigma \approx 2\)):

\[\hat{g}_k = \sum_{n=0}^{N-1} x(p_n)\, \phi\!\left(\frac{k}{\sigma} - p_n\right), \quad k = 0, \ldots, \sigma N - 1\]

The kernel \(\phi\) must be concentrated in the spatial domain (to limit the width of the convolution, i.e. the number of grid points touched per source point) while having a well-behaved Fourier transform (for the deapodization step). Two common choices are:

Kaiser-Bessel (KB) kernel — used by libraries such as NFFT, MIRT, and BART:

\[\phi_{\mathrm{KB}}(z) = \frac{I_0\!\left(\beta\sqrt{1-z^2}\right)}{I_0(\beta)}\, \mathbf{1}_{|z| \leq 1}\]

where \(I_0\) is the zeroth-order modified Bessel function and \(\beta\) is a shape parameter. Its Fourier transform is known analytically, which makes the deapodization step cheap.

Exponential-of-semicircle (ES) kernel — used by FINUFFT [8]:

\[\phi_{\mathrm{ES}}(z) = e^{\beta(\sqrt{1-z^2}-1)}\, \mathbf{1}_{|z| \leq 1}\]

The ES kernel is simpler and faster to evaluate than KB and achieves essentially the same error-convergence rate, but has no known analytic Fourier transform; its deapodization correction is therefore obtained by numerical quadrature (negligible extra cost).

In both cases \(\beta\) and the half-width \(W\) (number of grid cells touched per source point) govern the accuracy–cost tradeoff.

2. FFT

A standard FFT is applied to the oversampled grid \(\hat{g}\), yielding the oversampled frequency representation. The oversampled grid is then cropped (or subsampled for type 2) to obtain the \(N\) desired output frequencies.

3. Deapodization

The spreading step is equivalent to a convolution with \(\phi\), which in the frequency domain corresponds to a multiplication by \(\hat{\phi}\). This introduces a systematic bias that is corrected by dividing each output by the Fourier transform of the kernel, evaluated at the output frequencies:

\[X_k = \frac{\hat{G}_k}{\hat{\phi}(k/N)}, \quad k = 0, \ldots, N-1\]

The accuracy of the NUFFT depends primarily on the kernel half-width \(W\) (larger means fewer approximation errors but more work per point) and the oversampling factor \(\sigma\).

Type 2 transforms (uniform to non-uniform, i.e. the MRI forward model) perform these steps in reversed order: deapodization, IFFT, interpolation at the non-uniform locations.

This package exposes interfaces to the main NUFFT libraries available (See mrinufft.operators.interfaces). The choice of spreading kernel and correction are the main theoretical differences between backends.

Density Compensation

In non-uniform sampling, such as radial or spiral MRI, the acquired k-space samples \(k_m\) are not equally spaced. As a result, each sample does not contribute equally to the final image. To account for the non-uniform sampling density, a set of weights \(w_m\)—called the density compensation function (DCF)—is applied to the measured data \(y_m\).

In the adjoint NUFFT (type 2), which maps from non-uniform k-space onto the image grid, the operation can be mathematically written as:

\[x_n = \sum_{m=1}^M y_m \, w_m \, e^{-i 2\pi \nu_m \cdot x_n}\]

where:

• \(x_n\) is the reconstructed pixel value at position \(x_n\),

• \(y_m\) are the measured non-uniform k-space data,

• \(w_m\) is the density compensation weight for sample \(m\),

• \(\nu_m\) is the k-space sampling location.

The choice of \(w_m\) depends on the trajectory:

• Analytical weights: For simple trajectories (e.g., radial), \(w_m\) may have closed forms.

• Voronoi weights: \(w_m\) corresponds to the area/volume of Voronoi cells around each sample \(\nu_m\).

• Iterative methods: For arbitrary trajectories, \(w_m\) can be estimated via iterative algorithms that minimize reconstruction artifacts.

The DCF is typically applied before the NNUFFT ensuring each k-space measurement contributes proportionally to its neighborhood. Proper density compensation is crucial for artifact-free, quantitatively accurate image reconstruction.

Note

In mri-nufft, density compensation can be specified when initializing the NUFFT operator (via the density argument) as either a precomputed array, a method name (e.g., 'voronoi', 'pipe'), or by providing your own function. See FourierOperatorBase and the compute_density API for more details. Several geometry-based and NUFFT-based DCF methods are available in the mrinufft.density module.

Tip

For consistent scaling, density compensation weights should be normalized, so that the total signal energy is preserved across different trajectories and density choices. If you supply your own weights (See the for instance the normalization done for the pipe method)

Toeplitz NUFFT

In most application of the NUFFT operator, we will frequently need to apply the Gram Operator \(T=A^HWA\) where \(A\) is the discretized NUFFT operator, and \(W\) is a diagonal matrix of weights \(\boldsymbol{w}\) (e.g density compensation weights). Direct application of this operator is computationally expensive as it requires two NUFFT operations.

However, the operator \(T\) is a Toeplitz operator, and can be efficiently applied using only one NUFFT operation and one inverse NUFFT operation, by precomputing the convolution kernel \(K\) associated to \(T\). For a given image \(x\), the operation is:

\[\begin{split}(\boldsymbol{T}\boldsymbol{x})_n &= (\boldsymbol{A^HWAx})_n = \sum_{m=1}^M w_m e^{2i\pi \boldsymbol{\nu}_m \cdot \boldsymbol{u}_n} (\boldsymbol{Ax})_m\\ &= \sum_{m=1}^M w_m e^{2i\pi \boldsymbol{\nu}_m \cdot \boldsymbol{u}_n} \sum_{n'=1}^N e^{-2i\pi \boldsymbol{\nu}_m \cdot \boldsymbol{u}_{n'}} \boldsymbol{x}_{n'}\\ &= \sum_{n'=1}^N \sum_{m=1}^M w_m e^{2i\pi \boldsymbol{\nu}_m \cdot (\boldsymbol{u}_n-\boldsymbol{u}_{n'})} x_{n'} \\ &= \sum_{n'=1}^N K[n-n'] x_{n'} \\ &= K \ast x_n\end{split}\]

A few things can be said about this formulation:

• \(K\) is the point-spread-function (PSF) associated to a given k-space sampling trajectory.

• To capture all possible correlation (lags) between voxel in a grid of size \(N\) the resulting convolution kernel must have a size of \(2N\).

• for real-valued \(\boldsymbol{W}\) the Kernel should have hermitian symetry, since \(\boldsymbol{A^HWA}\) is an hermitian operator.

• \(T\) has Toeplitz structure in 1D and a BTTB (Block Toeplitz with Toeplitz Block) in 2D and 3D.

• \(T\) is not directly diagonalizable by a Fourier Transform, but Circulant matrix are, so to apply the Toeplitz kernel efficiently, we will zero-pad the image into a 2N grid, and we have:

\[Tx = \operatorname{iFFT}( \hat{K} \cdot FFT(x))\]

With \(\hat{K}\) the Fourier transform of the kernel.

Note

• Multi-coil operators (e.g with coil sensitivities) wraps around the Gram Operator easily,but don’t have a global Toeplitz structure, however they have a “toeplitz-like structure” that can be used [9]. It is not clear yet if this is more efficient than the classic SENSE operator, and is not implemented in MRI-NUFFT yet.

• For off-resonance correction, see [3], this is not implemented as well.

Computing the Toeplitz Kernel

The Toeplitz Kernel can be written as:

\[K[n] = \sum_{m=1}^{M} w_m e^{2i\pi \boldsymbol{\nu}_m r_n} \quad\text{with } n \in [-N,N]\]

This is the Adjoint NUFFT (type-1) of the input weights \(\boldsymbol{w}\), but evaluated on the extended lag grid \(n \in [-N, N]\) of size \(2N\). Computing this adjoint directly would require an oversized NUFFT plan (twice the resolution per dimension, i.e. \(2^d\) more memory: \(2^3=8\times\) in 3D, on top of the NUFFT’s own internal gridding oversampling). We avoid this entirely: the kernel is assembled from adjoint NUFFTs computed on the native \(N\)-sized grid, so the memory footprint is that of a single adjoint.

Two observations make this possible:

Fourier-shift modulation. A single adjoint NUFFT on the \(N\)-grid yields the kernel only on the inner lags \(n \in [-N/2, N/2)\). Modulating the weights by a linear phase before the adjoint shifts the window of computed lags: by the shift theorem,

\[\sum_{m=1}^{M} \big(w_m e^{i \boldsymbol{\omega}_m \cdot \boldsymbol{s}}\big)\, e^{i \boldsymbol{\omega}_m r_n} = \sum_{m=1}^{M} w_m\, e^{i \boldsymbol{\omega}_m (r_n + \boldsymbol{s})} = K[n + \boldsymbol{s}]\]

where \(\boldsymbol{\omega}_m = 2\pi\boldsymbol{\nu}_m\) are the sampling locations in radians. Choosing the integer shift \(\boldsymbol{s} = \pm N/2\) along each axis slides the window onto the outer lags, recovering the full \(2N\) support at \(N\)-resolution.

Hermitian symmetry. For real-valued weights \(\boldsymbol{w}\), the operator \(\boldsymbol{A^HWA}\) is Hermitian, so \(K[-n] = K[n]^*\). This halves the work: only the half-space \(n_0 \geq 0\) of the first axis is computed by modulation (the remaining sign combinations of the other axes), and the opposite half is filled by conjugate symmetry.

The kernel is therefore built from \(2^{d-1}\) adjoint NUFFTs (2 in 2D, 4 in 3D), each on the native \(N\)-grid, covering every combination of inner/outer lags. Because the shifts land the zero lag at index \([0, \dots, 0]\) by construction, no post-hoc re-centering is needed. This scheme is exact for any trajectory but requires even grid sizes. See the detail of the implementation in mrinufft.operators.toeplitz module.

Autodifferentiation

TBA

Autodifferentiation

Several MRI applications require differentiating through the NUFFT operator, most notably:

• Trajectory optimisation: learning the k-space sampling locations \(\boldsymbol{\nu}_m\) to minimise a reconstruction metric.

• Deep-learning reconstruction: end-to-end training of networks whose forward model includes the NUFFT.

• Field-map estimation: jointly estimating the image and the off-resonance map \(\Delta\omega\).

NUFFT backends (finufft, cufinufft, …) do not expose native autograd, so mri-nufft implements custom torch.autograd.Function wrappers in mrinufft.operators.autodiff. Three gradient modes are supported.

Gradient with respect to data

Because \(\mathcal{F}_\Omega\) is a linear operator, the vector-Jacobian product (VJP) needed by reverse-mode autodiff reduces to the adjoint (or forward) operator:

\[\frac{\partial L}{\partial \boldsymbol{x}} = \mathcal{F}_\Omega^*\!\left(\frac{\partial L}{\partial \boldsymbol{y}}\right), \qquad \frac{\partial L}{\partial \boldsymbol{y}} = \mathcal{F}_\Omega\!\left(\frac{\partial L}{\partial \boldsymbol{x}}\right)\]

No additional computation beyond the paired operator is needed.

Gradient with respect to the trajectory

The partial derivative of a single k-space sample \(y_m\) with respect to a trajectory coordinate \(\nu_{m,d}\) is:

\[\frac{\partial y_m}{\partial \nu_{m,d}} = \sum_{n=1}^N (-2i\pi\, u_{n,d})\, x_n\, e^{-2i\pi\, \boldsymbol{u}_n \cdot \boldsymbol{\nu}_m} = -2i\pi\, \bigl[\mathcal{F}_\Omega(\boldsymbol{x} \odot \boldsymbol{r}_d)\bigr]_m\]

where \(\boldsymbol{r}_d = (u_{1,d}, \ldots, u_{N,d})^T\) is the vector of pixel coordinates along dimension \(d\). The VJP with respect to the full trajectory is then:

\[\frac{\partial L}{\partial \nu_{m,d}} = -i\, \overline{\delta y_m}\cdot\bigl[\mathcal{F}_\Omega(\boldsymbol{x} \odot \boldsymbol{r}_d)\bigr]_m\]

where \(\delta \boldsymbol{y} = \partial L / \partial \boldsymbol{y}\) is the upstream gradient. Computing this requires \(d\) additional forward NUFFT evaluations — one per spatial dimension — applied to the image weighted by the pixel coordinates. This efficient approximation is described in [7].

Similarly, for the adjoint operator \(x_n = [\mathcal{F}_\Omega^* \boldsymbol{y}]_n\), the VJP is:

\[\frac{\partial L}{\partial \nu_{m,d}} = i\, y_m\cdot\bigl[\mathcal{F}_\Omega(\overline{\delta \boldsymbol{x}} \odot \boldsymbol{r}_d)\bigr]_m\]

Gradient with respect to the field map

For the off-resonance corrected model (see Off-resonance correction model), with:

\[y_m = \int x(\boldsymbol{u})\,e^{-2i\pi\boldsymbol{u}\cdot\boldsymbol{\nu}_m + \Delta\omega(\boldsymbol{u})\,t_m}\,d\boldsymbol{u}\]

the gradient of a real loss \(L\) with respect to the (real-valued) field map \(\Delta\omega(\boldsymbol{u}_n)\) is:

\[\frac{\partial L}{\partial \Delta\omega(\boldsymbol{u}_n)} = x_n \,\bigl[\mathcal{F}_{\Omega,t}^*(\delta\boldsymbol{y})\bigr]_n\]

where \(\mathcal{F}_{\Omega,t}^*\) is the corrected adjoint operator with the readout time \(t_m\) absorbed into the weights.

Tip

Use MRINufftAutoGrad to wrap any existing NUFFT operator for PyTorch autodiff. The flags wrt_data, wrt_traj, and wrt_field_map individually enable each gradient mode.

from mrinufft import get_operator
from mrinufft.operators.autodiff import MRINufftAutoGrad

nufft = get_operator("finufft")(samples, shape)
# gradient w.r.t. image data only (default)
autograd_nufft = MRINufftAutoGrad(nufft, wrt_data=True)
# gradient w.r.t. k-space trajectory
autograd_nufft = MRINufftAutoGrad(nufft, wrt_data=True, wrt_traj=True)

Other Application

Apart from MRI, The NUFFT operator is also used for:

• Electron tomography

• Probability Density Function estimation

• Astronomical Imaging

• …

These applications are not covered by this package, do it yourself!

References

[1]

https://en.m.wikipedia.org/wiki/Non-uniform_discrete_Fourier_transform

[2]

Noll, D. C., Meyer, C. H., Pauly, J. M., Nishimura, D. G., Macovski, A., “A homogeneity correction method for magnetic resonance imaging with time-varying gradients”, IEEE Transaction on Medical Imaging (1991), pp. 629-637.

[3] (1,2)

Fessler, J. A., Lee, S., Olafsson, V. T., Shi, H. R., Noll, D. C., “Toeplitz-based iterative image reconstruction for MRI with correction for magnetic field inhomogeneity”, IEEE Transactions on Signal Processing 53.9 (2005), pp. 3393–3402.

[4]

D.F. McGivney et al., “SVD Compression for Magnetic Resonance Fingerprinting in the Time Domain,” IEEE Transactions on Medical Imaging (2014), pp. 2311-2322.

[5]

D.C. Noll, C. H. Meyer, J. M. Pauly, D. G. Nishimura and A. Macovski, “A homogeneity correction method for magnetic resonance imaging with time-varying gradients,” in IEEE Transactions on Medical Imaging, vol. 10, no. 4, pp. 629-637, Dec. 1991, doi: 10.1109/42.108599

[6]

Man, L.-C., Pauly, J.M. and Macovski, A. (1997), Multifrequency interpolation for fast off-resonance correction. Magn. Reson. Med., 37: 785-792. https://doi.org/10.1002/mrm.1910370523

[7]

Wang, G., Fessler, J. A., “Efficient approximation of Jacobian matrices involving a non-uniform fast Fourier transform (NUFFT).” IEEE Transactions on Computational Imaging 9 (2023), pp. 43–54.

[8]

Barnett, A. H., Magland, J., af Klinteberg, L., “A parallel non-uniform fast Fourier transform library based on an exponential of semicircle kernel.” SIAM Journal on Scientific Computing 41.5 (2019), pp. C479–C504.

[9]

Chen, J., Christodoulou, A.G., Fan, Z. (2026). Direct Inversion Formula of the Multi-coil MR Operator Under Arbitrary Trajectories. In: Gee, J.C., et al. Medical Image Computing and Computer Assisted Intervention – MICCAI 2025. MICCAI 2025. Lecture Notes in Computer Science, vol 15963. Springer, Cham. https://doi.org/10.1007/978-3-032-04965-0_14