docs: add more information on NUFFT with ORC.

Author: paquiteau

CONTRIBUTING.md

@@ -149,6 +149,8 @@ To build the documentation locally, you can run the following commands :
     # you can specify a specfic backend to run the example with an environment variable:
     export MRINUFFT_BACKEND=finufft
     python -m sphinx-build docs docs_build 
+    # If you don't care about examples you can use
+    python -m sphinx-build docs docs_build -D plot_gallery=0
     # view the documentation in your browser
     python -m http.server --directory docs_build 8080
     # open localhost:8080 in your browser

docs/mrinufft_convention.rst

@@ -1,5 +1,5 @@
 .. include:: <isonum.txt>
-
+.. _mri-nufft-interface:
 ===============================
 MRI-NUFFT Interfaces Convention
 ===============================

docs/nufft.rst

@@ -1,5 +1,6 @@
 .. include:: <isonum.txt>
 .. _NUFFT:
+
 ====================
  The NUFFT Operator
 ====================
@@ -96,6 +97,10 @@ 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 :ref:`nufft-interface` for the standard, and :class:`FourierOperatorBase <mrinufft.operators.base.FourierOperatorBase>`
+
+
 Parallel Imaging Model
 ~~~~~~~~~~~~~~~~~~~~~~
 
@@ -124,6 +129,7 @@ Where :math:`S_1, \dots, S_L` are the sensitivity maps of each coil. Such sensit
     TODO Add ref to SENSE and CG-Sense
 
 .. _nufft-orc:
+
 Off-resonance correction model
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -152,13 +158,23 @@ Yielding the following model, where :math:`L \ll M, N` regular 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 :math:`B=(b_{m, \ell}) \in \mathbb{C}^{M\times L}` and :math:`C=(c_\ell, n) \in \mathbb{C}^{L\times N}` can be (optimally) estimated for any given :math:`L` by solving the following matrix factorisation problem [3]_:
+The coefficients :math:`B=(b_{m, \ell}) \in \mathbb{C}^{M\times L}` and :math:`C=(c_\ell, n) \in \mathbb{C}^{L\times N}` can be (optimally) estimated for any given :math:`L` by solving the following matrix factorisation problem [4]_:
 
 .. math::
 
    \hat{B}, \hat{C} = \arg\min_{B,C} \| E- BC\|_{fro}^2
 
-Where :math:`E_mn = e^{i\Delta\omega_0(u_n)t_m}`.
+Where :math:`E_{mn} = e^{i\Delta\omega_0(u_n)t_m}`.
+
+.. note::
+   
+   The estimation of the B and C methods are provided in the :py:mod:`mrinufft.extras.field_map` module.
+   Other methods like MTI [5]_ and MFI [6]_ are also available.
+
+.. tip::
+
+   You can use the method :func:`.with_off_resonance_correction <mrinufft.operators.base.FourierOperatorBase.with_off_resonance_correction>` to augment an existing operator with off-resonance correction capability.
+
 
 
 .. TODO Add Reference to the Code doing this.
@@ -219,8 +235,10 @@ the projection operator :math:`\boldsymbol{\Phi}` commutes with the Fourier tran
 
 that is, computation now involves :math:`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.
 
-The Non Uniform Fast Fourier Transform
-======================================
+.. _nufft-algo:
+
+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 :math:`O(MN)` complexity to :math:`O(M\log(N))`). Thus, the main steps of the *Non-Uniform Fast Fourier Transform* are for the type 1:
@@ -229,20 +247,46 @@ In order to lower the computational cost of the Non-Uniform Fourier Transform, t
 2. Perform the (I)FFT on this image
 3. Downsampling to the final grid, and apply some bias correction.
 
-This package exposes 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.
+This package exposes interfaces to the main NUFFT libraries available (See :mod:`mrinufft.operators.interfaces`). 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 transforms perform those steps in reversed order.
 
-.. TODO Add Reference to all the NUFFT methods article
-   Maybe to Fessler never-going-to-be-published book.
-
 
 Density Compensation
 ====================
 
 
+In non-uniform sampling, such as radial or spiral MRI, the acquired k-space samples :math:`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 :math:`w_m`—called the density compensation function (DCF)—is applied to the measured data :math:`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:
+
+.. math::
+
+    x_n = \sum_{m=1}^M y_m \, w_m \, e^{-i 2\pi k_m \cdot x_n}
+
+where:
+
+- :math:`x_n` is the reconstructed pixel value at position :math:`x_n`,
+- :math:`y_m` are the measured non-uniform k-space data,
+- :math:`w_m` is the density compensation weight for sample :math:`m`,
+- :math:`k_m` is the k-space sampling location.
+
+The choice of :math:`w_m` depends on the trajectory:
+
+- **Analytical weights**: For simple trajectories (e.g., radial), :math:`w_m` may have closed forms.
+- **Voronoi weights**: :math:`w_m` corresponds to the area/volume of Voronoi cells around each sample :math:`k_m`.
+- **Iterative methods**: For arbitrary trajectories, :math:`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 :class:`FourierOperatorBase <mrinufft.operators.base.FourierOperatorBase>` and the ``compute_density`` API for more details. Several geometry-based and NUFFT-based DCF methods are available in the :mod:`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 :func:`pipe <mrinufft.operators.interfaces.cufinufft.MRICufinufft.pipe>` method)
 
+   
 Other Application
 =================
 Apart from MRI, The NUFFT operator is also used for:
@@ -261,3 +305,5 @@ References
 .. [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] 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

src/mrinufft/extras/field_map.py

@@ -1,6 +1,12 @@
-"""Field map generator module."""
+"""Field map Module.
+
+This module provides methods to generate dummy B0map as well as estimation of the
+bilinearization of the off-resonance model (See :ref:`nufft-orc` for a detailed
+explanation).
+"""
 
 import numpy as np
+from collections.abc import Callable
 
 from numpy.typing import ArrayLike, NDArray
 
@@ -149,6 +155,10 @@ def _make_sphere(shape, frac_radius=0.3):
     time points; nt = len(readout_time).
 E: NDArray
     [nbins, nt] exponential off-resonance phase matrix at input histogram bins.
+
+See Also
+--------
+    :ref:`nufft-orc`
 """,
 )
 
@@ -180,7 +190,7 @@ def get_complex_fieldmap_rad(
 
     See Also
     --------
-    :ref:`_nufft-orc`
+    :ref:`nufft-orc`
 
     """
     xp = get_array_module(b0_map)
@@ -242,6 +252,14 @@ def _create_histogram(
 class C_lazy:
     """A lazy version of the C interpolator.
 
+    Instead of storing C as a ``(L, Nx, Ny, Nz)`` array, we store the quantize
+    version of shape L, N_bins, and the indexes mapping each voxel of the
+    complex-valued field-map to the bins values.
+
+    When accessing the L-th interpolator, the full ``(Nx, Ny, Nz)`` associated array is
+    generated.
+
+
     Parameters
     ----------
     C_sr: NDArray of shape (L, Nbinsreal, Nbinsimag)
@@ -282,14 +300,17 @@ def __getitem__(self, i: int) -> NDArray:
 
 @with_numpy_cupy
 def _full_C(
-    field_map: NDArray,
+    field_map: NDArray[np.complex64],
     C_small: NDArray,
     n_bins: tuple[int, int],
     lazy: bool = False,
 ):
     """
     Generate a full spatial interpolator C from a quantized version.
 
+    If `lazy=True` uses the an array like structure to generate the
+    interpolators on the fly.
+
     Parameters
     ----------
     field_map: NDArray
@@ -300,7 +321,12 @@ def _full_C(
         Number of bins for the real and imaginary part
     lazy: bool, default False
         If True, returns a lazy version of the C matrix.
-    Return
+
+    Returns
+    -------
+    NDArray if lazy = False
+
+    C_Lazy if lazy = True
     """
     xp = get_array_module(field_map)
     fr = field_map.real.ravel()
@@ -329,7 +355,9 @@ def _full_C(
         return xp.ascontiguousarray(C_big.reshape(-1, *field_map.shape))
 
 
-def _get_svds(xp, partial_svd):
+def _get_svds(
+    xp, partial_svd
+) -> Callable[[NDArray, int], tuple[NDArray, NDArray, NDArray]]:
     if partial_svd:
         if xp.__name__ == "cupy":
             from cupyx.scipy.sparse.linalg import svds
@@ -355,22 +383,29 @@ def compute_svd_coefficients(
     lazy: bool = False,
     partial_svd: bool = True,
 ):
-    """
-    Compute off-resonance correction coefficients using tSVD.
+    r"""
+    Compute off-resonance correction coefficients using an SVD.
+
+    Solves :math:`\arg\min_{B,C} = \|E - BC \|_{F}^2`
+
+    In practise it utilizes truncated SVD to extract L leading singular vectors from
+    the weighted exponentiated phase matrix sampled at field map histogram centers.
 
-    Utilizes truncated SVD to extract L leading singular vectors from the weighted
-    exponentiated phase matrix sampled at field map histogram centers. Outputs tSVD
-    basis (B), interpolation matrix (C), and exponential phase matrix (E).
 
     Parameters
     ----------
     ${base_params}
-
     partial_svd: bool, default True
         If True, only compute the L components required for the estimation, not
         the full SVD.
 
     ${returns}
+
+    References
+    ----------
+    Sutton BP, Noll DC, Fessler JA. Fast, iterative image reconstruction for MRI
+    in the presence of field inhomogeneities. IEEE Trans Med Imaging. 2003
+    Feb;22(2):178-88. doi: 10.1109/tmi.2002.808360. PMID: 12715994.
     """
     field_map = get_complex_fieldmap_rad(field_map)
     xp = get_array_module(field_map)
@@ -416,6 +451,13 @@ def compute_mti_coefficients(
     ${base_params}
 
     ${returns}
+
+    References
+    ----------
+    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
     """
     field_map = get_complex_fieldmap_rad(field_map)
     xp = get_array_module(field_map)
@@ -483,6 +525,13 @@ def compute_mfi_coefficients(
     ${base_params}
 
     ${returns}
+
+    References
+    ----------
+    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
+
     """
     # Format the input and apply the weight option
 

src/mrinufft/operators/base.py

@@ -206,7 +206,7 @@ def _safe_squeeze(self, arr):
         return arr
 
     @abstractmethod
-    def op(self, data):
+    def op(self, data: NDArray) -> NDArray:
         """Compute operator transform.
 
         Parameters
@@ -222,7 +222,7 @@ def op(self, data):
         pass
 
     @abstractmethod
-    def adj_op(self, coeffs):
+    def adj_op(self, coeffs: NDArray) -> NDArray:
         """Compute adjoint operator transform.
 
         Parameters
@@ -237,7 +237,7 @@ def adj_op(self, coeffs):
         """
         pass
 
-    def data_consistency(self, image_data, obs_data):
+    def data_consistency(self, image_data: NDArray, obs_data: NDArray) -> NDArray:
         """Compute the gradient data consistency.
 
         This is the naive implementation using adj_op(op(x)-y).

src/mrinufft/operators/off_resonance.py

@@ -54,17 +54,22 @@ class MRIFourierCorrected(FourierOperatorBase):
         - If ``{"name":name, **kwargs}`` use an existing
           method in `extra_field_map.py` parameterize by kwargs.
         - If``tuple[NDArray, NDArray]`` use this directly as the decomposition
-         (B and C)
+          (B and C)
 
     Notes
     -----
     The total field map used to calculate the field coefficients is
     ``field_map = R2*_map + 1j * B0_map``. If R2* is not provided,
     the field is purely imaginary: ``field_map = 1j * B0_map``.
 
+    You can also  use the method :py:func:`.with_off_resonance_correction
+    <mrinufft.operators.base.FourierOperatorBase.with_off_resonance_correction>`
+    to augment an existing operator with off-resonance correction capability.
+
+
     See Also
     --------
-    :ref:`_nufft-orc`
+    :ref:`nufft-orc`
 
     """