notebooks/SPECT/SPECT_acquisition_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "962178b0",
   "metadata": {},
   "source": [
    "# Acquisition Models for SPECT\n",
    "This demonstration shows how to set-up and use SIRF acquisition models for SPECT. You should have tried the `introduction` notebook first. The current notebook briefly repeats some items without explanation.\n",
    "\n",
    "This demo is a jupyter notebook, i.e. intended to be run step by step.\n",
    "You could export it as a Python file and run it one go, but that might\n",
    "make little sense as the figures are not labelled.\n",
    "\n",
    "Forward projection demo: creates an image, projects it to simulate\n",
    "acquisition data and backprojects\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "441e4cb1",
   "metadata": {},
   "source": [
    "Authors: Daniel Deidda, Sam Porter, Kris Thielemans\n",
    "\n",
    "First version: 13th of May 2022 \n",
    "\n",
    "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  \n",
    "Copyright 2022 National Physical Laboratory.  \n",
    "Copyright 2015 - 2019, 2022 University College London.  \n",
    "\n",
    "This is software developed for the Collaborative Computational\n",
    "Project in Synergistic Reconstruction for Biomedical Imaging\n",
    "(http://www.ccpsynerbi.ac.uk/).\n",
    "\n",
    "SPDX-License-Identifier: Apache-2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f80fd2f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the working directory for the notebook\n",
    "import notebook_setup\n",
    "from sirf_exercises import cd_to_working_dir\n",
    "from sirf.Utilities import examples_data_path\n",
    "from sirf.STIR import show_2D_array\n",
    "import sirf\n",
    "# Initial imports etc\n",
    "import numpy\n",
    "import os\n",
    "import sys\n",
    "import shutil\n",
    "import brainweb\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "# import engine module\n",
    "from sirf.STIR import MessageRedirector\n",
    "msg_red = MessageRedirector('info.txt', 'warnings.txt', 'errors.txt')\n",
    "import sirf.STIR as spect\n",
    "cd_to_working_dir('SPECT', 'SPECT_acquisition_model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "539c95b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_sample_image(image):\n",
    "    '''fill the image with some simple geometric shapes.'''\n",
    "    image.fill(0)\n",
    "    # create a shape\n",
    "    shape = spect.EllipticCylinder()\n",
    "    shape.set_length(400)\n",
    "    shape.set_radii((100, 40))\n",
    "    shape.set_origin((0, 60, 10))\n",
    "\n",
    "    # add the shape to the image\n",
    "    image.add_shape(shape, scale = 1)\n",
    "\n",
    "    # add another shape\n",
    "    shape.set_radii((30, 30))\n",
    "    shape.set_origin((60, -30, 10))\n",
    "    image.add_shape(shape, scale = 1.5)\n",
    "\n",
    "    # add another shape\n",
    "    shape.set_origin((-60, -30, 10))\n",
    "    image.add_shape(shape, scale = 0.75)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ec9dcbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "templ_sino = spect.AcquisitionData(os.path.join(examples_data_path('SPECT'), '','template_sinogram.hs'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "850fd91f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create image with suitable sizes\n",
    "image = templ_sino.create_uniform_image()\n",
    "create_sample_image(image)\n",
    "image.write(\"simulated_image.hv\")\n",
    "# z-pixel coordinate of the xy-cross-section to show\n",
    "z = image.dimensions()[0]//2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fdf9db4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show the phantom image\n",
    "image_array = image.as_array()\n",
    "show_2D_array('Phantom image', image_array[z,:,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8fd2b50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select acquisition model that implements the geometric\n",
    "# forward projection by a ray tracing matrix multiplication\n",
    "acq_model_matrix = spect.SPECTUBMatrix();\n",
    "acq_model = spect.AcquisitionModelUsingMatrix(acq_model_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e472daf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# require same number slices and equal z-sampling for projection data & image\n",
    "image = image.zoom_image(zooms=(0.5, 1.0, 1.0))\n",
    "print('projecting image...')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e54f20c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# project the image to obtain simulated acquisition data\n",
    "    # data from raw_data_file is used as a template\n",
    "acq_model.set_up(templ_sino, image)\n",
    "simulated_data = templ_sino.get_uniform_copy()\n",
    "acq_model.forward(image, 0, 1, simulated_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b7d94e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show simulated acquisition data\n",
    "simulated_data_as_array = simulated_data.as_array()\n",
    "middle_slice=simulated_data_as_array.shape[0]//2\n",
    "show_2D_array('Forward projection', simulated_data_as_array[0, middle_slice,:,:])\n",
    "# write data\n",
    "simulated_data.write(\"simulated_data.hs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b306045",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('backprojecting the forward projection...')\n",
    "# backproject the computed forward projection\n",
    "back_projected_image = acq_model.backward(simulated_data, 0, 1)\n",
    "\n",
    "back_projected_image_as_array = back_projected_image.as_array()\n",
    "show_2D_array('Backprojection', back_projected_image_as_array[z,:,:])\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "19b3ca16cdfc0eb05e0a6512c0e4ce5cf00e8dfec5d97ffb38b4ca3f11434ee7"
  },
  "kernelspec": {
   "display_name": "Python 3.9.5 64-bit ('base': conda)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}