[feature] add subspace projection model.

docs/nufft.rst

@@ -166,6 +166,57 @@ Where :math:`E_mn = e^i\Delta\omega_0(u_n)t_m`.
 
 .. _nufft-algo:
 
+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:
+
+.. math::
+
+   \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}
+
+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:
+
+.. math::
+
+   \boldsymbol{x} = \Phi\Phi^H \boldsymbol{x}
+
+where :math:`\Phi` is an orthonormal basis spanning a low dimensional subspace whose rank :math:`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:
+
+.. math::
+
+   \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}
+
+where :math:`\boldsymbol{\alpha} = \Phi^H \boldsymbol{x}` are the spatial coefficients representing the image series. Since the elements of :math:`\Phi^H` do not depend on the specific k-space frequency points,
+the projection operator :math:`\boldsymbol{\Phi}` commutes with the Fourier transform, and the signal equation finally becomes:
+
+.. math::
+
+   \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
 ======================================
@@ -206,5 +257,6 @@ 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
+.. [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.

examples/example_subspace.py

@@ -0,0 +1,240 @@
+"""
+=======================
+Subspace NUFFT Operator
+=======================
+
+An example to show how to setup a subspace NUFFT operator.
+
+This example shows how to use the Subspace NUFFT operator to acquire data
+and reconstruct an image (or a volume for 3D imaging)
+while projecting onto a low rank spatiotemporal subspace.
+This can be useful e.g., for dynamic multi-contrast MRI.
+Hereafter a 2D spiral trajectory is used for demonstration.
+
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from mrinufft import display_2D_trajectory
+
+# %%
+# Data preparation
+# ================
+#
+# Image loading
+# -------------
+#
+# For realistic a 2D image we will use the BrainWeb dataset,
+# installable using ``pip install brainweb-dl``.
+
+from mrinufft.extras import get_brainweb_map
+
+M0, T1, T2 = get_brainweb_map(0)
+M0 = np.flip(M0, axis=(0, 1, 2))[90]
+T1 = np.flip(T1, axis=(0, 1, 2))[90]
+T2 = np.flip(T2, axis=(0, 1, 2))[90]
+
+fig1, ax1 = plt.subplots(1, 3)
+
+im0 = ax1[0].imshow(M0, cmap="gray")
+ax1[0].axis("off"), ax1[0].set_title("M0 [a.u.]")
+fig1.colorbar(im0, ax=ax1[0], fraction=0.046, pad=0.04)
+
+im1 = ax1[1].imshow(T1, cmap="magma")
+ax1[1].axis("off"), ax1[1].set_title("T1 [ms]")
+fig1.colorbar(im1, ax=ax1[1], fraction=0.046, pad=0.04)
+
+im2 = ax1[2].imshow(T2, cmap="viridis")
+ax1[2].axis("off"), ax1[2].set_title("T2 [ms]")
+fig1.colorbar(im2, ax=ax1[2], fraction=0.046, pad=0.04)
+
+plt.tight_layout()
+plt.show()
+
+# %%
+# Sequence Parameters
+# -------------------
+#
+# As an example, we simulate a simple monoexponential Spin Echo acquisition.
+# We assume that k-space center is sampled at each TE (as in spiral or radial imaging)
+# and a constant flip angle train is constant set to 180 degrees.
+# In this way, we obtain an image for each of the shots in the echo train.
+
+from mrinufft.extras import fse_simulation
+
+ETL = 48  # Echo Train Length
+ESP = 6.0  # Echo Spacing [ms]
+TE = np.arange(ETL, dtype=np.float32) * ESP  # [ms]
+TR = 3000.0  # [ms]
+
+# %%
+# Subspace Generation
+# -------------------
+#
+# Here, we generate temporal subspace basis
+# by performing the Singular Value Decomposition of an ensemble
+# of simulated signals corresponding to the same T1-T2 range
+# of our object.
+
+
+def make_grid(T1range, T2range, natoms=100):
+    """Prepare parameter grid for basis estimation."""
+    # Create linear grid
+    T1grid = np.linspace(T1range[0], T1range[1], num=natoms, dtype=np.float32)
+    T2grid = np.linspace(T2range[0], T2range[1], num=natoms, dtype=np.float32)
+
+    # Cartesian product
+    T1grid, T2grid = np.meshgrid(T1grid, T2grid)
+
+    return T1grid.ravel(), T2grid.ravel()
+
+
+def estimate_subspace_basis(train_data, ncoeff=4):
+    """Estimate subspace data via SVD of simulated training data."""
+    # Perform svd
+    _, _, basis = np.linalg.svd(train_data, full_matrices=False)
+
+    # Calculate basis (retain only ncoeff coefficients)
+    basis = basis[:ncoeff]
+
+    return basis
+
+
+# Get range from tissues
+T1range = (T1.min() + 1, T1.max())  # [ms]
+T2range = (T2.min() + 1, T2.max())  # [ms]
+
+# Prepare tissue grid
+T1grid, T2grid = make_grid(T1range, T2range)
+
+# Calculate training data
+train_data = fse_simulation(1.0, T1grid, T2grid, TE, TR).astype(np.float32)
+
+# Calculate basis
+basis = estimate_subspace_basis(train_data.T)
+
+fig2, ax2 = plt.subplots(1, 2)
+ax2[0].plot(TE, train_data[:, ::100]), ax2[0].set(
+    xlabel="TE [ms]", ylabel="signal [a.u.]"
+), ax2[0].set_title("training dataset")
+ax2[1].plot(TE, basis.T), ax2[1].set(xlabel="TE [ms]", ylabel="signal [a.u.]"), ax2[
+    1
+].set_title("subspace basis")
+
+plt.show()
+
+# %%
+# Here, we simulate Brainweb FSE data with the same
+# sequence parameters as those used for the subspace estimation.
+
+mri_data = fse_simulation(M0, T1, T2, TE, TR).astype(np.float32)
+mri_data = np.ascontiguousarray(mri_data)
+
+# Ground truth subspace coefficients
+ground_truth = (mri_data.T @ basis.T).T
+ground_truth = np.ascontiguousarray(ground_truth)
+ground_truth_display = np.concatenate(
+    [abs(coeff) / abs(coeff).max() for coeff in ground_truth], axis=-1
+)
+plt.figure()
+plt.imshow(ground_truth_display, cmap="gray"), plt.axis("off"), plt.title(
+    "ground truth subspace coefficients"
+)
+
+plt.show()
+
+# %%
+# Trajectory generation
+# ---------------------
+
+from mrinufft import initialize_2D_spiral
+from mrinufft.density import voronoi
+
+samples = initialize_2D_spiral(
+    Nc=ETL * 16, Ns=1200, nb_revolutions=10, tilt="mri-golden"
+)
+
+# assume trajectory is reordered as (ncontrasts, nshots_per_contrast, nsamples_per_shot, ndims)
+# with contrast axis sorted by ascending TEs
+samples = samples.reshape(ETL, 16, *samples.shape[1:])
+
+# compute density compensation
+density = voronoi(samples)
+density = density.reshape(ETL, 16, samples.shape[-2])
+
+display_2D_trajectory(samples.reshape(-1, *samples.shape[2:]))
+
+# %%
+# Operator setup
+# ==============
+
+from mrinufft import get_operator
+from mrinufft.operators import MRISubspace
+
+# Generate standard NUFFT operator
+nufft = get_operator("finufft")(
+    samples=samples,
+    shape=mri_data.shape[-2:],
+    density=density.ravel(),
+)
+
+# Generate subspace-projected NUFFT operator
+subspace_nufft = MRISubspace(nufft, subspace_basis=basis)
+
+# %%
+# Generate K-Space
+# ----------------
+#
+# We generate the k-space data using a non-projected operator.
+# This can be simply obtained by using an identity matrix
+# with of shape (ETL, ETL) (= number of contrasts) as a subspace basis.
+
+multicontrast_nufft = [
+    get_operator("finufft")(
+        samples=samples[n],
+        shape=mri_data.shape[-2:],
+        density=density[n].ravel(),
+    )
+    for n in range(basis.shape[-1])
+]
+kspace = np.stack([multicontrast_nufft[n].op(mri_data[n]) for n in range(ETL)], axis=0)
+
+# %%
+# Image reconstruction
+# -----------------------
+#
+# We now reconstruct both using the subspace expanded adjoint operator
+# and the zero-filling adjoint operator followed by projection in image space
+
+# Reconstruct with projection in image space
+zerofilled_data_adj = np.stack(
+    [multicontrast_nufft[n].adj_op(kspace[n]) for n in range(ETL)], axis=0
+)
+zerofilled_display = (zerofilled_data_adj.T @ basis.T).T
+zerofilled_display = np.concatenate(
+    [abs(coeff) / abs(coeff).max() for coeff in zerofilled_display], axis=-1
+)
+
+# Reconstruct with projection in k-space
+subspace_data_adj = subspace_nufft.adj_op(kspace)
+subspace_display = np.concatenate(
+    [abs(coeff) / abs(coeff).max() for coeff in subspace_data_adj], axis=-1
+)
+
+
+fig3, ax3 = plt.subplots(2, 1)
+ax3[0].imshow(zerofilled_display, cmap="gray")
+ax3[0].set_xticks([])
+ax3[0].set_yticks([])
+ax3[0].set(ylabel="kspace projection")
+ax3[0].set_title("reconstructed subspace coefficients")
+ax3[1].imshow(subspace_display, cmap="gray")
+ax3[1].set_xticks([])
+ax3[1].set_yticks([])
+ax3[1].set(ylabel="zerofill + projection")
+
+plt.show()
+
+# %%
+# The projected k-space is equivalent to the regular reconstruction followed by projection.

src/mrinufft/extras/__init__.py

@@ -1,13 +1,17 @@
-"""Sensitivity map estimation methods."""
+"""Additional supports routines."""
 
-from .data import make_b0map, make_t2smap
+from .data import get_brainweb_map
+from .field_map import make_b0map, make_t2smap
+from .sim import fse_simulation
 from .smaps import low_frequency
 from .utils import get_smaps
 
 
 __all__ = [
+    "fse_simulation",
+    "get_brainweb_map",
+    "get_smaps",
+    "low_frequency",
     "make_b0map",
     "make_t2smap",
-    "low_frequency",
-    "get_smaps",
 ]