Merge pull request #227 from NicoleJurjew/PSMR_24_PET_introductory

Psmr24 pet introductory
SyneRBI · May 16, 2024 · 4581d44 · 4581d44
2 parents 71ce2bb + ab9718f
commit 4581d44
Show file tree

Hide file tree

Showing 13 changed files with 433 additions and 34 deletions.
diff --git a/notebooks/PET/DIY_OSEM.ipynb b/notebooks/PET/DIY_OSEM.ipynb
@@ -16,11 +16,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Authors: Kris Thielemans  \n",
+    "Authors: Kris Thielemans\n",
     "First version: June 2021\n",
     "\n",
+    "Solutions: Nicole Jurjew, May 2024\n",
+    "\n",
     "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  \n",
-    "Copyright 2021 University College London.\n",
+    "Copyright 2021, 2024 University College London.\n",
     "\n",
     "This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging (http://www.ccpsynerbi.ac.uk/).\n",
     "\n",
@@ -45,6 +47,10 @@
     "\n",
     "# Setup the working directory for the notebook\n",
     "import notebook_setup\n",
+    "import os\n",
+    "\n",
+    "# define the directory in which this notebook is saved\n",
+    "nb_dir = os.getcwd()\n",
     "from sirf_exercises import cd_to_working_dir\n",
     "cd_to_working_dir('PET', 'OSEM_reconstruction')"
    ]
@@ -58,7 +64,6 @@
     "#%% Initial imports etc\n",
     "import numpy\n",
     "import matplotlib.pyplot as plt\n",
-    "import os\n",
     "import sirf.STIR as pet\n",
     "from sirf.Utilities import examples_data_path\n",
     "from sirf_exercises import exercises_data_path\n",
@@ -144,11 +149,11 @@
     "$$ \\bar{y}=A x + b$$\n",
     "\n",
     "The MLEM update is\n",
-    "$$ x^{\\mathrm{new}} = \\frac{x}{A^t 1} A^t \\left(\\frac{y}{A x + b}\\right)$$\n",
+    "$$ x^{\\mathrm{new}} = \\frac{x}{A^T 1} A^T \\left(\\frac{y}{A x + b}\\right)$$\n",
     "\n",
-    "You hopefully recognise that the denominator of the factor on the right corresponds to the `forward` model applied ot the image $x$. Multiplication with the $A^t$ is the `backward` operation. So, we have used all the main operations already. We just need to do element-wise multiplication an division operation, but that's easy!\n",
+    "You hopefully recognise that the denominator of the factor on the right corresponds to the `forward` model applied ot the image $x$. Multiplication with the $A^T$ is the `backward` operation. So, we have used all the main operations already. We just need to do element-wise multiplication and division operations, but that's easy!\n",
     "\n",
-    "Let's first compute $A^t 1$, as this is a image that won't change over iterations. Note that the $1$ here represents a one-vector, i.e., an image filled with ones. It is often called the \"sensivity image\" as it is (proportional to) the probability that an event emitted in a voxel is detected by the scanner (without scattering). "
+    "Let's first compute $A^T 1$, as this is an image that won't change over iterations. Note that the $1$ here represents a one-vector, i.e., an image filled with ones. It is often called the \"sensivity image\" as it is (proportional to) the probability that an event emitted in a voxel is detected by the scanner (without scattering). "
    ]
   },
   {
@@ -233,6 +238,31 @@
     "    return estimated_image"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir_str = os.path.join(nb_dir, 'Solution_Snippets', 'DIY_OSEM_01_MLEM.py')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Uncomment the following line and run the cell to see the solution, to run the lines you'll need to run the cell a second time**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %load $dir_str"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -246,7 +276,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "estimated_image = MLEM(acquired_data, acq_model,image.get_uniform_copy(1), 40)"
+    "MLEM_image = MLEM(acquired_data, acq_model,image.get_uniform_copy(1), 40)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.subplot(1,1,1)\n",
+    "plt.imshow(MLEM_image.as_array()[im_slice,:,:])"
    ]
   },
   {
@@ -350,6 +391,51 @@
     "    return estimated_image"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir_str = os.path.join(nb_dir, 'Solution_Snippets', 'DIY_OSEM_02_OSEM.py')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Uncomment the following line and run the cell to see the solution, to run the lines you'll need to run the cell a second time**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %load $dir_str"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "OSEM_image = OSEM(acquired_data, acq_model,image.get_uniform_copy(1), 10)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.subplot(1,1,1)\n",
+    "plt.imshow(OSEM_image.as_array()[im_slice,:,:])"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

diff --git a/notebooks/PET/OSEM_reconstruction.ipynb b/notebooks/PET/OSEM_reconstruction.ipynb
@@ -17,8 +17,10 @@
     "Authors: Kris Thielemans and Evgueni Ovtchinnikov  \n",
     "First version: June 2021\n",
     "\n",
+    "Solutions: Nicole Jurjew, May 2024\n",
+    "\n",
     "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  \n",
-    "Copyright 2015 - 2018, 2021 University College London.\n",
+    "Copyright 2015 - 2018, 2021, 2024 University College London.\n",
     "\n",
     "This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging (http://www.ccpsynerbi.ac.uk/).\n",
     "\n",
@@ -91,9 +93,16 @@
     "import sirf.STIR as pet\n",
     "from sirf.Utilities import examples_data_path\n",
     "from sirf_exercises import exercises_data_path\n",
+    "from sirf_exercises import cd_to_working_dir\n",
     "\n",
     "# define the directory with input files for this notebook\n",
-    "data_path = os.path.join(examples_data_path('PET'), 'thorax_single_slice')"
+    "data_path = os.path.join(examples_data_path('PET'), 'thorax_single_slice')\n",
+    "\n",
+    "# define the directory in which this notebook is saved\n",
+    "nb_dir = os.getcwd()\n",
+    "\n",
+    "# change into working directory\n",
+    "cd_to_working_dir('PET', 'OSEM_reconstruction')\n"
    ]
   },
   {
@@ -522,6 +531,31 @@
     "    return acq_model"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir_str = os.path.join(nb_dir, 'Solution_Snippets', 'OSEM_reconstruction_01_create_acq_model.py')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Uncomment the following line and run the cell to see the solution, to run the lines you'll need to run the cell a second time**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %load $dir_str"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -549,6 +583,31 @@
     "    return recon.get_output()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir_str = os.path.join(nb_dir, 'Solution_Snippets', 'OSEM_reconstruction_02_OSEM.py')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Uncomment the following line and run the cell to see the solution, to run the lines you'll need to run the cell a second time**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %load $dir_str"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -566,6 +625,50 @@
     "my_reconstructed_image = OSEM(acquired_data, acq_model, image.get_uniform_copy(cmax), 30, 4)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "my_reconstructed_image.show(im_slice,title='recon img')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You will probably see artifacts in the image above, namely some very high activity images on the edges of the FOV. You can prevent that by adding a Cylinder processor, which essentially sets voxels outside the cylindrical FOV to 0. See below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir_str = os.path.join(nb_dir, 'Solution_Snippets', 'OSEM_reconstruction_02a_OSEMwCylProc.py')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %load $dir_str"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "my_reconstructed_image2 = OSEM(acquired_data, acq_model, image.get_uniform_copy(cmax), 30, 4)\n",
+    "my_reconstructed_image2.show(im_slice,title='recon img')"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -580,6 +683,24 @@
     "Hint: the easiest way would be to take the existing attenuation image, and use `get_uniform_copy(0)` to create an image where all $\\mu$-values are 0. Another (and more efficient) way would be to avoid creating the `AcquisitionSensitivityModel` at all."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dir_str = os.path.join(nb_dir, 'Solution_Snippets', 'OSEM_reconstruction_03_create_acq_model_noAtt.py')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %load $dir_str"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -593,15 +714,15 @@
    "source": [
     "In these exercises we have used attenuation images of the same size as the emission image. This is easier to code and for display etc, but it is not a requirement of SIRF.\n",
     "\n",
-    "In addition, we have simulated and reconstructed the data with the same `AcquisitionModel` (and preserved image sizes). This is also convenient, but not a requirement (as you've seen in the NAC exercise). In fact, do not write your next paper using this \"inverse crime\". The problem is too "
+    "In addition, we have simulated and reconstructed the data with the same `AcquisitionModel` (and preserved image sizes). This is also convenient, but not a requirement (as you've seen in the NAC exercise). In fact, do not write your next paper using this \"inverse crime\"."
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 2",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "python2"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {

diff --git a/notebooks/PET/Solution_Snippets/DIY_OSEM_01_MLEM.py b/notebooks/PET/Solution_Snippets/DIY_OSEM_01_MLEM.py
@@ -0,0 +1,9 @@
+def MLEM(acquired_data, acq_model, initial_image, num_iterations):
+    estimated_image = initial_image.clone()
+    sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1))
+    for i in range(num_iterations):
+        quotient = acquired_data/acq_model.forward(estimated_image)  # y / (Ax + b)
+        quotient.fill(numpy.nan_to_num(quotient.as_array()))
+        mult_update = acq_model.backward(quotient)/sensitivity         # A^t * quotient / A^t1
+        mult_update.fill(numpy.nan_to_num(mult_update.as_array()))
+    return estimated_image
diff --git a/notebooks/PET/Solution_Snippets/DIY_OSEM_02_OSEM.py b/notebooks/PET/Solution_Snippets/DIY_OSEM_02_OSEM.py
@@ -0,0 +1,21 @@
+def OSEM(acquired_data, acq_model, initial_image, num_iterations):
+    estimated_image=initial_image.clone()
+    # some stuff here - hint, this will be similar to your solution for MLEM
+    # but you will have to additionally iterate over your subsets 
+    for i in range(num_iterations):
+        for s in range(acq_model.num_subsets):
+
+            subs_sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1), subset_num=s)
+
+            quotient = acquired_data/acq_model.forward(estimated_image, subset_num=s)     # y / (Ax + b)
+            quotient.fill(numpy.nan_to_num(quotient.as_array()))
+
+            mult_update = acq_model.backward(quotient, subset_num=s)/subs_sensitivity     # A^t * quotient / A^t1
+            mult_update.fill(numpy.nan_to_num(mult_update.as_array()))
+
+            estimated_image *= mult_update                                                        # update (in place)
+
+            estimated_image.maximum(0)
+
+    #  some stuff here
+    return estimated_image