Skip to content

Commit

Permalink
Merge branch 'master' into stack-random
Browse files Browse the repository at this point in the history
  • Loading branch information
@paquiteau
paquiteau authored Dec 16, 2024
2 parents d4edc5d + 9c41feb commit f6aad60
Show file tree
Hide file tree
Showing 26 changed files with 2,091 additions and 524 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/test-ci.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,7 @@ jobs:
shell: bash
run: |
python -m pip install --upgrade pip
python -m pip install -e .[test,dev]
python -m pip install -e .[extra,test,dev]
python -m pip install finufft pooch brainweb-dl torch fastmri
- name: Install GPU related interfaces
Expand Down
26 changes: 13 additions & 13 deletions docs/nufft.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,14 +174,14 @@ To achieve a clinically feasible scan time, each frame or contrast is acquired w
.. math::
\tilde{\boldsymbol{y}} = \begin{bmatrix}
\mathcal{F}_\Omega_1 & 0 & \cdots & 0 \\
0 & \mathcal{F}_\Omega_2 & \cdots & 0 \\
\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
0 & 0 & \cdots & \mathcal{F}_{\Omega_T}
\end{bmatrix}
\boldsymbol{x} + \boldsymbol{n}
where :math:`\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
where :math:`\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 :math:`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:

Expand All @@ -195,17 +195,17 @@ training dataset or an ensemble of simulated Bloch responses. The signal model c
.. math::
\tilde{\boldsymbol{y}} = \begin{bmatrix}
\mathcal{F}_\Omega_1 & 0 & \cdots & 0 \\
0 & \mathcal{F}_\Omega_2 & \cdots & 0 \\
\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
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 \\
\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
0 & 0 & \cdots & \mathcal{F}_{\Omega_T}
\end{bmatrix}
\Phi \boldsymbol{\alpha} + \boldsymbol{n}
Expand All @@ -214,15 +214,15 @@ the projection operator :math:`\boldsymbol{\Phi}` commutes with the Fourier tran

.. math::
\tilde{\boldsymbol{y}} = \Phi \mathcal{F}_Omega(\boldsymbol{\alpha}) + \boldsymbol{n}
\tilde{\boldsymbol{y}} = \Phi \mathcal{F}_\Omega(\boldsymbol{\alpha}) + \boldsymbol{n}
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
======================================


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

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
Expand Down Expand Up @@ -259,4 +259,4 @@ 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] 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.
.. [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.
44 changes: 22 additions & 22 deletions examples/example_2D_trajectories.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
# External
import matplotlib.pyplot as plt
import numpy as np
from utils import show_argument, show_trajectory
from utils import show_trajectories, show_trajectory

# Internal
import mrinufft as mn
Expand Down Expand Up @@ -75,7 +75,7 @@

arguments = [8, 16, 32, 64]
function = lambda x: mn.initialize_2D_radial(x, Ns, tilt=tilt, in_out=in_out)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -89,7 +89,7 @@

arguments = [8, 16, 32, 64]
function = lambda x: mn.initialize_2D_radial(Nc, x, tilt=tilt, in_out=in_out)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -105,7 +105,7 @@

arguments = ["uniform", "golden", "mri-golden", np.pi / 17]
function = lambda x: mn.initialize_2D_radial(Nc, Ns, tilt=x, in_out=in_out)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -128,7 +128,7 @@

arguments = [True, False]
function = lambda x: mn.initialize_2D_radial(Nc, Ns, tilt=tilt, in_out=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -170,7 +170,7 @@
function = lambda x: mn.initialize_2D_spiral(
Nc, Ns, tilt=tilt, nb_revolutions=x, in_out=in_out
)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -194,7 +194,7 @@

arguments = ["galilean", "archimedes", "fermat", 1 / 4]
function = lambda x: mn.initialize_2D_spiral(Nc, Ns, tilt=tilt, spiral=x, in_out=in_out)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -221,7 +221,7 @@
Ns,
patch_center=x,
)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -267,7 +267,7 @@
Ns,
spiral_reduction=x,
)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -293,7 +293,7 @@
Ns,
patch_center=x,
)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -333,7 +333,7 @@
function = lambda x: mn.initialize_2D_cones(
Nc, Ns, tilt=tilt, in_out=in_out, nb_zigzags=x
)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -348,7 +348,7 @@

arguments = [0.2, 1, 2, 3]
function = lambda x: mn.initialize_2D_cones(Nc, Ns, tilt=tilt, in_out=in_out, width=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -413,7 +413,7 @@

arguments = [2, 3, 4, 6]
function = lambda x: mn.initialize_2D_propeller(Nc, Ns, nb_strips=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -450,7 +450,7 @@

arguments = [Nc, int(2 * Nc / 3), int(Nc / 3)]
function = lambda x: mn.initialize_2D_rings(Nc=x, Ns=Ns, nb_rings=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)

# %%
#
Expand All @@ -461,7 +461,7 @@

arguments = [Nc, int(4 * Nc / 3), 2 * Nc]
function = lambda x: mn.initialize_2D_rings(Nc=x, Ns=Ns, nb_rings=Nc)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -501,7 +501,7 @@

arguments = [0, 1, 5, 10]
function = lambda x: mn.initialize_2D_rosette(Nc, Ns, in_out=in_out, coprime_index=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -541,7 +541,7 @@
function = lambda x: mn.initialize_2D_polar_lissajous(
Nc, Ns, in_out=in_out, coprime_index=x
)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -564,7 +564,7 @@
function = lambda x: mn.initialize_2D_polar_lissajous(
Nc, Ns, in_out=in_out, nb_segments=x
)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -588,7 +588,7 @@
function = lambda x: mn.initialize_2D_polar_lissajous(
Nc, Ns, in_out=io, coprime_index=cpi, nb_segments=x
)
show_argument(
show_trajectories(
function, arguments, one_shot=one_shot, subfig_size=subfigure_size
)

Expand Down Expand Up @@ -625,7 +625,7 @@

arguments = [1, 2.5, 5, 10]
function = lambda x: mn.initialize_2D_waves(Nc, Ns, nb_zigzags=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand All @@ -644,7 +644,7 @@

arguments = [0, 1, 1.5, 3]
function = lambda x: mn.initialize_2D_waves(Nc, Ns, width=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down Expand Up @@ -679,7 +679,7 @@

arguments = [1, 1.5, 2, 3]
function = lambda x: mn.initialize_2D_lissajous(Nc, Ns, density=x)
show_argument(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)
show_trajectories(function, arguments, one_shot=one_shot, subfig_size=subfigure_size)


# %%
Expand Down
Loading

0 comments on commit f6aad60

Please sign in to comment.