Skip to content

Commit

Permalink
Merge pull request #149 from mind-inria/debug_doc
Browse files Browse the repository at this point in the history
Debug doc
  • Loading branch information
philouc authored Jun 28, 2024
2 parents 6d44381 + 0334ba0 commit ba54f49
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 25 deletions.
46 changes: 23 additions & 23 deletions docs/nufft.rst
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
The NUFFT Operator
====================

This document gives a general overview of the Non-Uniform Fast fourier transform (NUFFT), and its application in MRI.
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
Expand All @@ -18,15 +18,15 @@ The 1D-NUDFT [1]_ is defined as:

.. math::
X_k = \sum_{n=0}^{N-1} x(p_n) e^{ -2\pi i p_n \nu_k}
X_k = \sum_{n=0}^{N-1} x(p_n) e^{ -2\imath \pi p_n \nu_k}
where :math:`X_k` is the frequency point at :math:`\nu_k`.

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

There exists 3 types of NUDFT:
There exist 3 types of NUDFT:

* Type 1: :math:`p_n = n/N` and :math:`\nu_k` are non uniformly spaced
* Type 1: :math:`p_n = n/N` and :math:`\nu_k` are non-uniformly spaced
* Type 2: :math:`p_n` are non uniformly spaced and :math:`\nu_k = k/M`
* Type 3: :math:`p_n` and :math:`\nu_k` are non uniformly spaced
* If :math:`p_n=n/N` and :math:`\nu_k=k/M` then the NUDFT is simply the Discrete Fourier Transform.
Expand All @@ -38,28 +38,28 @@ 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.
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 required data, using a non Cartesian (i.e. non uniformly distributed in *every* direction) sampling scheme offer great opportunity.
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:

.. math::
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
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:

- :math:`x(\boldsymbol{u})` is the spatially varying image contrast acquired.
- :math:`y_1, \dots, y_M` are the sampled points at frequency locations :math:`\Omega=\lbrace \boldsymbol{\nu}_1, \dots, \boldsymbol{\nu}_M \in [-1/2, 1/2]^d\rbrace`.
Typically images (:math:`d=2`) or volumes (:math:`d=3`) are acquired.
- :math:`n_i` is a zero-mean complex valued Gaussian Noise, modeling the "thermal noise" of the scanner.
- :math:`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 :math:`N`:

.. math::
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
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 :math:`\boldsymbol{u}_j` are the :math:`N` spatial locations of image voxels.
This is stated using the operator notation:
Expand Down Expand Up @@ -94,12 +94,12 @@ As the sampling locations :math:`\Omega` are non-uniform and the image locations

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

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).
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 :math:`L` coils is:

Expand All @@ -123,33 +123,33 @@ Where :math:`S_1, \dots, S_L` are the sensitivity maps of each coil. Such sensit
..
TODO Add ref to SENSE and CG-Sense
Off-Resonance Correction Model
Off-resonance correction model
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The constant magnetic field applied in a MRI machine :math:`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 :math:`\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:
The constant magnetic field applied in a MRI machine :math:`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 :math:`\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:

.. math::
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}
y(t_i) = \int_{\mathbb{R}^d} x(\boldsymbol{u}) e^{-2\imath\pi \boldsymbol{u} \cdot\boldsymbol{\nu_i} + \Delta\omega(\boldsymbol{u}) t_i} d\boldsymbol{u}
where :math:`t_i` is the time at which the frequency :math:`\nu_i` is acquired. Similarly at the reconstruction we have

.. math::
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}
x(\boldsymbol{u_n}) = \sum_{m}^M y(t_m) e^{2\imath\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.
With these mixed-domain field pertubations, the Fourier model does not hold anymore and the 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 :math:`b_{m, \ell}` and :math:`c_{\ell, n}`:

.. math::
e^{i\Delta\omega(\boldsymbol{u_n}) t_m} = \sum_{\ell=1}^L b_{m, \ell}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 :math:`L \ll M, N` regular Fourier transforms are performed to approximate the non-Fourier transform.

.. math::
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}}
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]_:

Expand All @@ -171,15 +171,15 @@ 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 :math:`O(MN)` complexity to `O(M\log(N))`). Thus, the main steps of the *Non-Uniform Fast Fourier Transform* are for the type 1:
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 `O(M\log(N))`). Thus, the main steps of the *Non-Uniform Fast Fourier Transform* are for the type 1:

1. Spreading/Interpolation of the non-uniform point to an oversampled Cartesian grid (typically with twice the resolution of the final image)
1. Spreading/Interpolation of the non-uniform points to an oversampled Cartesian grid (typically with twice the resolution of the final image)
2. Perform the (I)FFT on this image
3. 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.
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.

Type 2 transform performs those steps in reversed order.
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.
Expand All @@ -200,7 +200,7 @@ Apart from MRI, The NUFFT operator is also used for:
- Astronomical Imaging
- ...

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

References
==========
Expand Down
4 changes: 2 additions & 2 deletions docs/trajectory_gradspec.rst
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
The trajectory binary file specification
=========================================

The k-space trajectories are transformed to a binary file which is processed by scanner through the arbitrary gradients sequences from NeuroSpin.
The k-space trajectories are transformed to a binary file which is processed by a given scanner through arbitrary gradient sequences. At Neurospin, we have specifically tested MR systems from Siemens-Healthineers vendor.
This file mainly specifies an arbitrary gradient profile which is played on scanner at gradient raster time rate (10e-6 seconds).

The binary file format is specified as follows:
Expand Down Expand Up @@ -40,7 +40,7 @@ The binary file format is specified as follows:
+----------------+-------+---------+---------+------------------------------------------------------------------------+


:mod:`mrinufft.trajectories.io` module helps to convert a trajectory as numpy array to a binary file and vice versa.
:mod:`mrinufft.trajectories.io` module helps convert a trajectory as numpy array to a binary file and vice versa.

All the trajectory FLOAT's are specified with `float32` always.

Expand Down

0 comments on commit ba54f49

Please sign in to comment.