Example for WASABI

examples/_toc.yml

@@ -14,4 +14,7 @@ chapters:
   - file: t2star_multi_echo_flash
   - file: t1rho_se_single_line
   - file: t1_t2_spiral_cmrf
+- file: system.md
+  sections:
+  - file: b0_b1_wasabi
 

examples/b0_b1_wasabi.ipynb

@@ -0,0 +1,275 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# B0 and B1 Mapping using WASABI\n",
+    "\n",
+    "Simultaneous mapping of water shift (B0) and B1 (WASABI) using using an off-resonant preparation pulse. This induces Rabi oscillations which can be seen as sinc-like oscillations in the frequency-offset dimension. B0 can then be estimated by its symmetry axis and B1 by its oscillation frequency."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1",
+   "metadata": {},
+   "source": [
+    "### Imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import tempfile\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import MRzeroCore as mr0\n",
+    "import numpy as np\n",
+    "import torch\n",
+    "from einops import rearrange\n",
+    "from mrpro.algorithms.reconstruction import DirectReconstruction\n",
+    "from mrpro.data import KData\n",
+    "from mrpro.data import SpatialDimension\n",
+    "from mrpro.data.traj_calculators import KTrajectoryCartesian\n",
+    "from mrpro.operators import DictionaryMatchOp\n",
+    "from mrpro.operators.models import WASABI\n",
+    "\n",
+    "from mrseq.scripts.b0_b1_wasabi import main as create_seq\n",
+    "from mrseq.utils import sys_defaults"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3",
+   "metadata": {},
+   "source": [
+    "### Settings\n",
+    "We are going to use a numerical phantom with a matrix size of 128 x 128. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "image_matrix_size = [128, 128]\n",
+    "\n",
+    "frequency_offsets = np.linspace(-240, 240, 31)\n",
+    "norm_offset = -35e3\n",
+    "\n",
+    "tmp = tempfile.TemporaryDirectory()\n",
+    "fname_mrd = Path(tmp.name) / 'b0_b1_wasabi.mrd'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5",
+   "metadata": {},
+   "source": [
+    "### Create the digital phantom\n",
+    "\n",
+    "We use the standard Brainweb phantom from [MRzero](https://github.com/MRsources/MRzero-Core)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "phantom = mr0.util.load_phantom(image_matrix_size)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "### Create the WASABI sequence\n",
+    "\n",
+    "To create the WASABI sequence, we use the previously imported [b0_b1_wasabi script](../src/mrseq/scripts/b0_b1_wasabi.py).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sequence, fname_seq = create_seq(\n",
+    "    system=sys_defaults,\n",
+    "    test_report=False,\n",
+    "    timing_check=False,\n",
+    "    fov_xy=float(phantom.size.numpy()[0]),\n",
+    "    n_readout=image_matrix_size[0],\n",
+    "    n_phase_encoding=image_matrix_size[0],\n",
+    "    frequency_offsets=frequency_offsets,\n",
+    "    norm_offset=norm_offset,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9",
+   "metadata": {},
+   "source": [
+    "### Simulate the sequence\n",
+    "Now, we pass the sequence and the phantom to the MRzero simulation and save the simulated signal as an (ISMR)MRD file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n",
+    "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-1)\n",
+    "mr0.sig_to_mrd(fname_mrd, signal, sequence)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11",
+   "metadata": {},
+   "source": [
+    "### Reconstruct the images at different inversion times\n",
+    "\n",
+    "We use [MRpro](https://github.com/PTB-MR/MRpro) for the image reconstruction."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "kdata = KData.from_file(fname_mrd, trajectory=KTrajectoryCartesian())\n",
+    "kdata.header.encoding_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0] * 2)\n",
+    "kdata.header.recon_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0])\n",
+    "recon = DirectReconstruction(kdata, csm=None)\n",
+    "idata = recon(kdata)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13",
+   "metadata": {},
+   "source": [
+    "We can now plot the images at different offset frequencies."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "idat = idata.rss().abs().numpy().squeeze()\n",
+    "offsets = np.concatenate(([norm_offset], frequency_offsets))\n",
+    "fig, ax = plt.subplots(4, idat.shape[0] // 4, figsize=(3 * idat.shape[0] // 4, 3 * 4))\n",
+    "ax = ax.flatten()\n",
+    "for i in range(idat.shape[0]):\n",
+    "    ax[i].imshow(idat[i, :, :], cmap='gray')\n",
+    "    ax[i].set_title(f'Offset = {int(offsets[i])} Hz')\n",
+    "    ax[i].set_xticks([])\n",
+    "    ax[i].set_yticks([])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15",
+   "metadata": {},
+   "source": [
+    "### Estimate B0 and B1\n",
+    "We use a dictionary matching approach to estimate the T1 maps. Afterward, we compare them to the input and ensure they match."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "16",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dictionary = DictionaryMatchOp(WASABI(offsets=torch.as_tensor(offsets, dtype=torch.float32)))\n",
+    "dictionary.append(\n",
+    "    torch.linspace(-100.0, 100.0, 100),\n",
+    "    torch.linspace(0.1, 1.5, 100)[None, :],\n",
+    "    torch.linspace(-100.0, 100.0, 100)[None, None, :],\n",
+    "    torch.linspace(-100.0, 100.0, 100)[None, None, None, :],\n",
+    ")\n",
+    "b0_match, rb1_match, c_match, a_match = dictionary(idata.data)\n",
+    "\n",
+    "b0_input = np.roll(rearrange(phantom.B0.numpy().squeeze()[::-1, ::-1], 'x y -> y x'), shift=(1, 1), axis=(0, 1))\n",
+    "obj_mask = np.zeros_like(b0_input)\n",
+    "obj_mask[b0_input > 0] = 1\n",
+    "b0_measured = b0_match.numpy().squeeze() * obj_mask\n",
+    "\n",
+    "rb1_input = np.abs(\n",
+    "    np.roll(rearrange(phantom.B1.numpy().squeeze()[::-1, ::-1], 'x y -> y x'), shift=(1, 1), axis=(0, 1))\n",
+    ")\n",
+    "rb1_measured = rb1_match.abs().numpy().squeeze() * obj_mask\n",
+    "\n",
+    "fig, ax = plt.subplots(1, 4)\n",
+    "ax[0].imshow(b0_match.squeeze().numpy() * obj_mask)\n",
+    "ax[1].imshow(rb1_match.abs().squeeze().numpy() * obj_mask)\n",
+    "ax[2].imshow(c_match.squeeze().numpy() * obj_mask)\n",
+    "ax[3].imshow(a_match.squeeze().numpy() * obj_mask)\n",
+    "\n",
+    "fig, ax = plt.subplots(2, 3, figsize=(15, 6))\n",
+    "for cax in ax.flatten():\n",
+    "    cax.set_xticks([])\n",
+    "    cax.set_yticks([])\n",
+    "\n",
+    "im = ax[0, 0].imshow(b0_input, vmin=-40, vmax=40, cmap='bwr')\n",
+    "fig.colorbar(im, ax=ax[0, 0], label='Input B0 (Hz)')\n",
+    "\n",
+    "im = ax[0, 1].imshow(b0_measured, vmin=-1, vmax=1, cmap='bwr')\n",
+    "fig.colorbar(im, ax=ax[0, 1], label='Measured B0 (Hz)')\n",
+    "\n",
+    "im = ax[0, 2].imshow(b0_measured - b0_input, vmin=-1.8, vmax=1.8, cmap='bwr')\n",
+    "fig.colorbar(im, ax=ax[0, 2], label='Difference B0 (Hz)')\n",
+    "\n",
+    "relative_error = np.sum(np.abs(b0_input - b0_measured)) / np.sum(np.abs(b0_input))\n",
+    "print(f'Relative error {relative_error}')\n",
+    "# assert relative_error < 0.08\n",
+    "\n",
+    "im = ax[1, 0].imshow(rb1_input, vmin=0, vmax=1.2, cmap='magma')\n",
+    "fig.colorbar(im, ax=ax[1, 0], label='Input relative B1 ')\n",
+    "\n",
+    "im = ax[1, 1].imshow(rb1_measured, vmin=-1, vmax=1, cmap='magma')\n",
+    "fig.colorbar(im, ax=ax[1, 1], label='Measured relative B1')\n",
+    "\n",
+    "im = ax[1, 2].imshow(rb1_measured - rb1_input, vmin=-1.8, vmax=1.8, cmap='bwr')\n",
+    "fig.colorbar(im, ax=ax[1, 2], label='Difference relative B1')\n",
+    "\n",
+    "relative_error = np.sum(np.abs(rb1_input - rb1_measured)) / np.sum(np.abs(rb1_input))\n",
+    "print(f'Relative error {relative_error}')\n",
+    "# assert relative_error < 0.08"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "mrseq_mrzero",
+   "language": "python",
+   "name": "python3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

examples/system.md

@@ -0,0 +1,5 @@
+# System Characterization
+
+Sequences to characterize the MR system. Available examples are:
+
+- [B0 and B1 mapping using WASABI](b0_b1_wasabi.ipynb)