Examples typos

examples/qmri_sg_challenge_2024_t1.ipynb

@@ -5,7 +5,7 @@
    "id": "0f82262f",
    "metadata": {},
    "source": [
-    "# QMRI Challenge ISMRM 2024 - T1 mapping"
+    "# QMRI Challenge ISMRM 2024 - $T_1$ mapping"
    ]
   },
   {
@@ -40,7 +40,7 @@
    "source": [
     "### Overview\n",
     "The dataset consists of images obtained at 10 different inversion times using a turbo spin echo sequence. Each\n",
-    "inversion time is saved in a separate DICOM file. In order to obtain a T1 map, we are going to:\n",
+    "inversion time is saved in a separate DICOM file. In order to obtain a $T_1$ map, we are going to:\n",
     "- download the data from Zenodo\n",
     "- read in the DICOM files (one for each inversion time) and combine them in an IData object\n",
     "- define a signal model and data loss (mean-squared error) function\n",
@@ -105,7 +105,7 @@
     "fig, axes = plt.subplots(1, 3, squeeze=False)\n",
     "for idx, ax in enumerate(axes.flatten()):\n",
     "    ax.imshow(torch.abs(idata_multi_ti.data[idx, 0, 0, :, :]))\n",
-    "    ax.set_title(f'TI = {idata_multi_ti.header.ti[idx]:.0f}ms')"
+    "    ax.set_title(f'TI = {idata_multi_ti.header.ti[idx]:.3f}s')"
    ]
   },
   {
@@ -116,9 +116,9 @@
     "### Signal model and loss function\n",
     "We use the model $q$\n",
     "\n",
-    "$q(TI) = M_0 (1 - e^{-TI/T1})$\n",
+    "$q(TI) = M_0 (1 - e^{-TI/T_1})$\n",
     "\n",
-    "with the equilibrium magnetization $M_0$, the inversion time $TI$, and $T1$. We have to keep in mind that the DICOM\n",
+    "with the equilibrium magnetization $M_0$, the inversion time $TI$, and $T_1$. We have to keep in mind that the DICOM\n",
     "images only contain the magnitude of the signal. Therefore, we need $|q(TI)|$:"
    ]
   },
@@ -162,7 +162,7 @@
    "source": [
     "Now we can simply combine the two into a functional to solve\n",
     "\n",
-    "$ \\min_{M_0, T1} || |q(M_0, T1, TI)| - x||_2^2$"
+    "$ \\min_{M_0, T_1} || |q(M_0, T_1, TI)| - x||_2^2$"
    ]
   },
   {
@@ -187,11 +187,11 @@
     "To increase our chances of reaching the global minimum, we can ensure that our starting\n",
     "values are already close to the global minimum. We need a good starting point for each pixel.\n",
     "\n",
-    "One option to get a good starting point is to calculate the signal curves for a range of T1 values and then check\n",
+    "One option to get a good starting point is to calculate the signal curves for a range of $T_1$ values and then check\n",
     "for each pixel which of these signal curves fits best. This is similar to what is done for MR Fingerprinting. So we\n",
     "are going to:\n",
-    "- define a list of realistic T1 values (we call this a dictionary of T1 values)\n",
-    "- calculate the signal curves corresponding to each of these T1 values\n",
+    "- define a list of realistic $T_1$ values (we call this a dictionary of $T_1$ values)\n",
+    "- calculate the signal curves corresponding to each of these $T_1$ values\n",
     "- compare the signal curves to the signals of each voxel (we use the maximum of the dot-product as a metric of how\n",
     "well the signals fit to each other)"
    ]
@@ -203,8 +203,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Define 100 T1 values between 100 and 3000 ms\n",
-    "t1_dictionary = torch.linspace(100, 3000, 100)\n",
+    "# Define 100 T1 values between 0.1 and 3.0 s\n",
+    "t1_dictionary = torch.linspace(0.1, 3.0, 100)\n",
     "\n",
     "# Calculate the signal corresponding to each of these T1 values. We set M0 to 1, but this is arbitrary because M0 is\n",
     "# just a scaling factor and we are going to normalize the signal curves.\n",
@@ -227,8 +227,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# The image with the longest inversion time is a good approximation of the equilibrium magnetization\n",
-    "m0_start = torch.abs(idata_multi_ti.data[torch.argmax(idata_multi_ti.header.ti), ...])"
+    "# The maximum absolute value observed is a good approximation for m0\n",
+    "m0_start = torch.amax(torch.abs(idata_multi_ti.data), 0)"
    ]
   },
   {
@@ -242,10 +242,10 @@
     "fig, axes = plt.subplots(1, 2, figsize=(8, 2), squeeze=False)\n",
     "colorbar_ax = [make_axes_locatable(ax).append_axes('right', size='5%', pad=0.05) for ax in axes[0, :]]\n",
     "im = axes[0, 0].imshow(m0_start[0, 0, ...])\n",
-    "axes[0, 0].set_title('M0 start values')\n",
+    "axes[0, 0].set_title('$M_0$ start values')\n",
     "fig.colorbar(im, cax=colorbar_ax[0])\n",
-    "im = axes[0, 1].imshow(t1_start[0, 0, ...], vmin=0, vmax=2500)\n",
-    "axes[0, 1].set_title('T1 start values')\n",
+    "im = axes[0, 1].imshow(t1_start[0, 0, ...], vmin=0, vmax=2.5)\n",
+    "axes[0, 1].set_title('$T_1$ start values')\n",
     "fig.colorbar(im, cax=colorbar_ax[1])"
    ]
   },
@@ -266,7 +266,7 @@
    "source": [
     "# Hyperparameters for optimizer\n",
     "max_iter = 2000\n",
-    "lr = 1e0\n",
+    "lr = 1e-1\n",
     "\n",
     "# Run optimization\n",
     "params_result = adam(functional, [m0_start, t1_start], max_iter=max_iter, lr=lr)\n",
@@ -283,7 +283,7 @@
     "### Visualize the final results\n",
     "To get an impression of how well the fit has worked, we are going to calculate the relative error between\n",
     "\n",
-    "$E_{relative} = \\sum_{TI}\\frac{|(q(M_0, T1, TI) - x)|}{|x|}$\n",
+    "$E_{relative} = \\sum_{TI}\\frac{|(q(M_0, T_1, TI) - x)|}{|x|}$\n",
     "\n",
     "on a voxel-by-voxel basis"
    ]
@@ -304,10 +304,10 @@
     "fig, axes = plt.subplots(1, 3, figsize=(10, 2), squeeze=False)\n",
     "colorbar_ax = [make_axes_locatable(ax).append_axes('right', size='5%', pad=0.05) for ax in axes[0, :]]\n",
     "im = axes[0, 0].imshow(m0[0, 0, ...])\n",
-    "axes[0, 0].set_title('M0')\n",
+    "axes[0, 0].set_title('$M_0$')\n",
     "fig.colorbar(im, cax=colorbar_ax[0])\n",
-    "im = axes[0, 1].imshow(t1[0, 0, ...], vmin=0, vmax=2500)\n",
-    "axes[0, 1].set_title('T1')\n",
+    "im = axes[0, 1].imshow(t1[0, 0, ...], vmin=0, vmax=2.5)\n",
+    "axes[0, 1].set_title('$T_1$')\n",
     "fig.colorbar(im, cax=colorbar_ax[1])\n",
     "im = axes[0, 2].imshow(relative_absolute_error[0, 0, ...], vmin=0, vmax=1.0)\n",
     "axes[0, 2].set_title('Relative error')\n",

examples/qmri_sg_challenge_2024_t1.py

@@ -1,5 +1,5 @@
 # %% [markdown]
-# # QMRI Challenge ISMRM 2024 - T1 mapping
+# # QMRI Challenge ISMRM 2024 - $T_1$ mapping
 
 # %%
 # Imports
@@ -22,7 +22,7 @@
 # %% [markdown]
 # ### Overview
 # The dataset consists of images obtained at 10 different inversion times using a turbo spin echo sequence. Each
-# inversion time is saved in a separate DICOM file. In order to obtain a T1 map, we are going to:
+# inversion time is saved in a separate DICOM file. In order to obtain a $T_1$ map, we are going to:
 # - download the data from Zenodo
 # - read in the DICOM files (one for each inversion time) and combine them in an IData object
 # - define a signal model and data loss (mean-squared error) function
@@ -53,15 +53,15 @@
 fig, axes = plt.subplots(1, 3, squeeze=False)
 for idx, ax in enumerate(axes.flatten()):
     ax.imshow(torch.abs(idata_multi_ti.data[idx, 0, 0, :, :]))
-    ax.set_title(f'TI = {idata_multi_ti.header.ti[idx]:.0f}ms')
+    ax.set_title(f'TI = {idata_multi_ti.header.ti[idx]:.3f}s')
 
 # %% [markdown]
 # ### Signal model and loss function
 # We use the model $q$
 #
-# $q(TI) = M_0 (1 - e^{-TI/T1})$
+# $q(TI) = M_0 (1 - e^{-TI/T_1})$
 #
-# with the equilibrium magnetization $M_0$, the inversion time $TI$, and $T1$. We have to keep in mind that the DICOM
+# with the equilibrium magnetization $M_0$, the inversion time $TI$, and $T_1$. We have to keep in mind that the DICOM
 # images only contain the magnitude of the signal. Therefore, we need $|q(TI)|$:
 
 # %%
@@ -76,7 +76,7 @@
 # %% [markdown]
 # Now we can simply combine the two into a functional to solve
 #
-# $ \min_{M_0, T1} || |q(M_0, T1, TI)| - x||_2^2$
+# $ \min_{M_0, T_1} || |q(M_0, T_1, TI)| - x||_2^2$
 # %%
 functional = mse @ model
 
@@ -88,17 +88,17 @@
 # To increase our chances of reaching the global minimum, we can ensure that our starting
 # values are already close to the global minimum. We need a good starting point for each pixel.
 #
-# One option to get a good starting point is to calculate the signal curves for a range of T1 values and then check
+# One option to get a good starting point is to calculate the signal curves for a range of $T_1$ values and then check
 # for each pixel which of these signal curves fits best. This is similar to what is done for MR Fingerprinting. So we
 # are going to:
-# - define a list of realistic T1 values (we call this a dictionary of T1 values)
-# - calculate the signal curves corresponding to each of these T1 values
+# - define a list of realistic $T_1$ values (we call this a dictionary of $T_1$ values)
+# - calculate the signal curves corresponding to each of these $T_1$ values
 # - compare the signal curves to the signals of each voxel (we use the maximum of the dot-product as a metric of how
 # well the signals fit to each other)
 
 # %%
-# Define 100 T1 values between 100 and 3000 ms
-t1_dictionary = torch.linspace(100, 3000, 100)
+# Define 100 T1 values between 0.1 and 3.0 s
+t1_dictionary = torch.linspace(0.1, 3.0, 100)
 
 # Calculate the signal corresponding to each of these T1 values. We set M0 to 1, but this is arbitrary because M0 is
 # just a scaling factor and we are going to normalize the signal curves.
@@ -114,18 +114,18 @@
 t1_start = rearrange(t1_dictionary[idx_best_match], '(y x)->1 1 y x', y=n_y, x=n_x)
 
 # %%
-# The image with the longest inversion time is a good approximation of the equilibrium magnetization
-m0_start = torch.abs(idata_multi_ti.data[torch.argmax(idata_multi_ti.header.ti), ...])
+# The maximum absolute value observed is a good approximation for m0
+m0_start = torch.amax(torch.abs(idata_multi_ti.data), 0)
 
 # %%
 # Visualize the starting values
 fig, axes = plt.subplots(1, 2, figsize=(8, 2), squeeze=False)
 colorbar_ax = [make_axes_locatable(ax).append_axes('right', size='5%', pad=0.05) for ax in axes[0, :]]
 im = axes[0, 0].imshow(m0_start[0, 0, ...])
-axes[0, 0].set_title('M0 start values')
+axes[0, 0].set_title('$M_0$ start values')
 fig.colorbar(im, cax=colorbar_ax[0])
-im = axes[0, 1].imshow(t1_start[0, 0, ...], vmin=0, vmax=2500)
-axes[0, 1].set_title('T1 start values')
+im = axes[0, 1].imshow(t1_start[0, 0, ...], vmin=0, vmax=2.5)
+axes[0, 1].set_title('$T_1$ start values')
 fig.colorbar(im, cax=colorbar_ax[1])
 
 # %% [markdown]
@@ -134,7 +134,7 @@
 # %%
 # Hyperparameters for optimizer
 max_iter = 2000
-lr = 1e0
+lr = 1e-1
 
 # Run optimization
 params_result = adam(functional, [m0_start, t1_start], max_iter=max_iter, lr=lr)
@@ -146,7 +146,7 @@
 # ### Visualize the final results
 # To get an impression of how well the fit has worked, we are going to calculate the relative error between
 #
-# $E_{relative} = \sum_{TI}\frac{|(q(M_0, T1, TI) - x)|}{|x|}$
+# $E_{relative} = \sum_{TI}\frac{|(q(M_0, T_1, TI) - x)|}{|x|}$
 #
 # on a voxel-by-voxel basis
 
@@ -158,10 +158,10 @@
 fig, axes = plt.subplots(1, 3, figsize=(10, 2), squeeze=False)
 colorbar_ax = [make_axes_locatable(ax).append_axes('right', size='5%', pad=0.05) for ax in axes[0, :]]
 im = axes[0, 0].imshow(m0[0, 0, ...])
-axes[0, 0].set_title('M0')
+axes[0, 0].set_title('$M_0$')
 fig.colorbar(im, cax=colorbar_ax[0])
-im = axes[0, 1].imshow(t1[0, 0, ...], vmin=0, vmax=2500)
-axes[0, 1].set_title('T1')
+im = axes[0, 1].imshow(t1[0, 0, ...], vmin=0, vmax=2.5)
+axes[0, 1].set_title('$T_1$')
 fig.colorbar(im, cax=colorbar_ax[1])
 im = axes[0, 2].imshow(relative_absolute_error[0, 0, ...], vmin=0, vmax=1.0)
 axes[0, 2].set_title('Relative error')

examples/qmri_sg_challenge_2024_t2_star.ipynb

@@ -5,7 +5,7 @@
    "id": "5efa18e9",
    "metadata": {},
    "source": [
-    "# QMRI Challenge ISMRM 2024 - T2* mapping"
+    "# QMRI Challenge ISMRM 2024 - $T_2^*$ mapping"
    ]
   },
   {
@@ -39,7 +39,7 @@
    "source": [
     "### Overview\n",
     "The dataset consists of gradient echo images obtained at 11 different echo times, each saved in a separate DICOM file.\n",
-    "In order to obtain a T2* map, we are going to:\n",
+    "In order to obtain a $T_2^*$ map, we are going to:\n",
     "- download the data from Zenodo\n",
     "- read in the DICOM files (one for each echo time) and combine them in an IData object\n",
     "- define a signal model (mono-exponential decay) and data loss (mean-squared error) function\n",
@@ -100,6 +100,8 @@
    "source": [
     "te_dicom_files = data_folder.glob('**/*.dcm')\n",
     "idata_multi_te = IData.from_dicom_files(te_dicom_files)\n",
+    "# scaling the signal down to make the optimization easier\n",
+    "idata_multi_te.data[...] = idata_multi_te.data / 1500\n",
     "\n",
     "# Move the data to the GPU\n",
     "if flag_use_cuda:\n",
@@ -120,7 +122,7 @@
     "fig, axes = plt.subplots(1, 3, squeeze=False)\n",
     "for idx, ax in enumerate(axes.flatten()):\n",
     "    ax.imshow(torch.abs(idata_multi_te.data[idx, 0, 0, :, :]).cpu())\n",
-    "    ax.set_title(f'TE = {idata_multi_te.header.te[idx]:.0f}ms')"
+    "    ax.set_title(f'TE = {idata_multi_te.header.te[idx]:.3f}s')"
    ]
   },
   {
@@ -131,9 +133,9 @@
     "### Signal model and loss function\n",
     "We use the model $q$\n",
     "\n",
-    "$q(TE) = M_0 e^{-TE/T2^*}$\n",
+    "$q(TE) = M_0 e^{-TE/T_2^*}$\n",
     "\n",
-    "with the equilibrium magnetization $M_0$, the echo time $TE$, and $T2^*$"
+    "with the equilibrium magnetization $M_0$, the echo time $TE$, and $T_2^*$"
    ]
   },
   {
@@ -176,7 +178,7 @@
    "source": [
     "Now we can simply combine the two into a functional which will then solve\n",
     "\n",
-    "$ \\min_{M_0, T2^*} ||q(M_0, T2^*, TE) - x||_2^2$"
+    "$ \\min_{M_0, T_2^*} ||q(M_0, T_2^*, TE) - x||_2^2$"
    ]
   },
   {
@@ -207,11 +209,11 @@
     "# The shortest echo time is a good approximation of the equilibrium magnetization\n",
     "m0_start = torch.abs(idata_multi_te.data[torch.argmin(idata_multi_te.header.te), ...])\n",
     "# 20 ms as a starting value for T2*\n",
-    "t2star_start = torch.ones(m0_start.shape, dtype=torch.float32, device=m0_start.device) * 20\n",
+    "t2star_start = torch.ones(m0_start.shape, dtype=torch.float32, device=m0_start.device) * 20e-3\n",
     "\n",
     "# Hyperparameters for optimizer\n",
     "max_iter = 20000\n",
-    "lr = 1e0\n",
+    "lr = 1e-3\n",
     "\n",
     "if flag_use_cuda:\n",
     "    functional.cuda()\n",
@@ -235,7 +237,7 @@
     "### Visualize the final results\n",
     "To get an impression of how well the fit has worked, we are going to calculate the relative error between\n",
     "\n",
-    "$E_{relative} = \\sum_{TE}\\frac{|(q(M_0, T2^*, TE) - x)|}{|x|}$\n",
+    "$E_{relative} = \\sum_{TE}\\frac{|(q(M_0, T_2^*, TE) - x)|}{|x|}$\n",
     "\n",
     "on a voxel-by-voxel basis."
    ]
@@ -257,11 +259,11 @@
     "colorbar_ax = [make_axes_locatable(ax).append_axes('right', size='5%', pad=0.05) for ax in axes[0, :]]\n",
     "\n",
     "im = axes[0, 0].imshow(m0[0, 0, ...].cpu())\n",
-    "axes[0, 0].set_title('M0')\n",
+    "axes[0, 0].set_title('$M_0$')\n",
     "fig.colorbar(im, cax=colorbar_ax[0])\n",
     "\n",
-    "im = axes[0, 1].imshow(t2star[0, 0, ...].cpu(), vmin=0, vmax=500)\n",
-    "axes[0, 1].set_title('T2*')\n",
+    "im = axes[0, 1].imshow(t2star[0, 0, ...].cpu(), vmin=0, vmax=5)\n",
+    "axes[0, 1].set_title('$T_2^*$')\n",
     "fig.colorbar(im, cax=colorbar_ax[1])\n",
     "\n",
     "im = axes[0, 2].imshow(relative_absolute_error[0, 0, ...].cpu(), vmin=0, vmax=0.1)\n",