The NUFFT Operator#

This document gives a general overview of the Non-Uniform Fast fourier transform (NUFFT), and its application in 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\pi i 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 exists 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 line in a specific direction, in a Cartesian fashion.

In order to accelerate the acquisition of required data, using a non Cartesian (i.e. non uniformly distributed in every direction) sampling scheme offer great opportunity.

The acquisition model is usually described as:

\[y(\boldsymbol{\nu}_i) = \int_{\mathbb{R}^d} x(\boldsymbol{u}) e^{-2i\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^{-2i\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 way. First by taking into account Parallel Imaging, where multiple coils are receiving data, each with a dedicated sensitivity profile.

Parallel Imaging Model#

In MRI the acquired signal can be received by multiple antenna (“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 interfaces with 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 deviate 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_i) = \int_{\mathbb{R}^d} x(\boldsymbol{u}) e^{-2i\pi \boldsymbol{u} \cdot\boldsymbol{\nu_i} + \Delta\omega(\boldsymbol{u}) t_i} d\boldsymbol{u}\]

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

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

With these mixed-domain field pertubations, the Fourier model does not hold anymore and the FFT algorithm can not be used. 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^{2i\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 [3]:

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

Where \(E_mn = e^i\Delta\omega_0(u_n)t_m\).

The Non Uniform Fast Fourier Transform#

In order to lower the computational cost of the Non-Uniform Fourier Transform, the main idea to move back to a regular grid where an FFT would be performed (going from a typical \(O(MN)\) complexity to O(Mlog(N))). Thus, the main steps of the Non-Uniform Fast Fourier Transform are for the type 1:

Spreading/Interpolation of the non-uniform point to an oversampled Cartesian grid (typically with twice the resolution of the final image)
Perform the (I)FFT on this image
Downsampling to the final grid, and apply some bias correction.

This package proposes interfaces to the main NUFFT libraries available. The choice of the spreading method (ie the interpolation kernel) in step 1. and the correction applied in step 3. are the main theoretical differences between the methods.

Type 2 transform performs those steps in reversed order.

Density Compensation#

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 !

The NUFFT Operator

Contents