STIR_demonstration_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the necessary modules\n",
    "from sirf.STIR import (ImageData, AcquisitionData,\n",
    "                       SPECTUBMatrix, AcquisitionModelUsingMatrix,\n",
    "                       MessageRedirector)\n",
    "from src.simind import *\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "import pandas as pd\n",
    "import os\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Redirect STIR messages\n",
    "msg = MessageRedirector()\n",
    "\n",
    "# Hardcoded arguments that were originally passed via argparse\n",
    "total_activity = 258.423  # in MBq\n",
    "time_per_projection = 43  # in seconds\n",
    "photon_multiplier = 0.001\n",
    "photopeak_energy = 208  # keV\n",
    "window_lower = 187.56  # keV\n",
    "window_upper = 229.24  # keV\n",
    "source_type = 'lu177'\n",
    "collimator = 'G8-MEGP'\n",
    "kev_per_channel = 10.0\n",
    "max_energy = 498.3\n",
    "mu_map_path = 'data/Lu177/registered_CTAC.hv'\n",
    "image_path = 'data/Lu177/osem_reconstruction_postfilter_555.hv'\n",
    "measured_data_path = 'data/Lu177/SPECTCT_NEMA_128_EM001_DS_en_1_Lu177_EM.hdr'\n",
    "measured_additive = '/home/sam/working/STIR_users_MIC2023/data/Lu177/STIR_TEW.hs'\n",
    "output_dir = 'simind_output'\n",
    "output_prefix = 'output'\n",
    "input_smc_file_path = 'input/input.smc'\n",
    "scoring_routine = 1\n",
    "collimator_routine = 0\n",
    "photon_direction = 3\n",
    "crystal_thickness = 7.25  # mm\n",
    "crystal_half_length_radius = 393.6 / 2  # mm\n",
    "crystal_half_width = 511.7 / 2  # mm\n",
    "flag_11 = True\n",
    "\n",
    "# Calculating the number of energy spectra channels\n",
    "num_energy_spectra_channels = max_energy // kev_per_channel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the acquisition model function\n",
    "def get_acquisition_model(measured_data, additive_data, image, mu_map_stir):\n",
    "    acq_matrix = SPECTUBMatrix()\n",
    "    acq_matrix.set_attenuation_image(mu_map_stir)\n",
    "    acq_matrix.set_keep_all_views_in_cache(True)\n",
    "    acq_matrix.set_resolution_model(1.81534, 0.02148, False)\n",
    "    \n",
    "    acq_model = AcquisitionModelUsingMatrix(acq_matrix)\n",
    "    acq_model.set_up(measured_data, image)\n",
    "    \n",
    "    try:\n",
    "        acq_model.set_additive(additive_data)\n",
    "    except Exception as e:\n",
    "        print(e)\n",
    "        print(\"Could not set additive data\")\n",
    "        \n",
    "    return acq_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load image, mu_map, measured data, and measured additive data\n",
    "image = ImageData(image_path)\n",
    "mu_map = ImageData(mu_map_path)\n",
    "measured_data = AcquisitionData(measured_data_path)\n",
    "measured_additive = AcquisitionData(measured_additive)\n",
    "\n",
    "# Adjusting mu_map for STIR\n",
    "mu_map_stir = mu_map.clone()\n",
    "mu_map_stir.fill(np.flip(mu_map.as_array(), axis=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change to the correct directory\n",
    "os.chdir(\"/home/sam/working/STIR_users_MIC2023\")\n",
    "\n",
    "# Set up the SIMIND simulator\n",
    "simulator = SimindSimulator(template_smc_file_path=input_smc_file_path,\n",
    "                            output_dir=output_dir, output_prefix=output_prefix,\n",
    "                            source=image, mu_map=mu_map, template_sinogram=measured_data)\n",
    "\n",
    "# Configure the simulator parameters\n",
    "simulator.add_comment(\"Demonstration of SIMIND simulation\")\n",
    "simulator.set_windows(window_lower, window_upper, 0)\n",
    "simulator.add_index(\"photon_energy\", photopeak_energy)\n",
    "simulator.add_index(\"scoring_routine\", scoring_routine)\n",
    "simulator.add_index(\"collimator_routine\", collimator_routine)\n",
    "simulator.add_index(\"photon_direction\", photon_direction)\n",
    "simulator.add_index(\"source_activity\", total_activity * time_per_projection)\n",
    "simulator.add_index(\"crystal_thickness\", crystal_thickness / 10)  # cm\n",
    "simulator.add_index(\"crystal_half_length_radius\", crystal_half_length_radius / 10)  # cm\n",
    "simulator.add_index(\"crystal_half_width\", crystal_half_width / 10)  # cm\n",
    "simulator.config.set_flag(11, flag_11)\n",
    "simulator.add_index(\"step_size_photon_path_simulation\", min(*image.voxel_sizes()) / 10)  # cm\n",
    "simulator.add_index(\"energy_resolution\", 9.5)  # percent\n",
    "simulator.add_index(\"intrinsic_resolution\", 0.31)  # cm\n",
    "\n",
    "# Set runtime switches\n",
    "simulator.add_runtime_switch(\"CC\", collimator)\n",
    "simulator.add_runtime_switch(\"NN\", photon_multiplier)\n",
    "simulator.add_runtime_switch(\"FI\", source_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the simulation\n",
    "simulator.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve simulation outputs\n",
    "simind_total = simulator.get_output_total()\n",
    "simind_scatter = simulator.get_output_scatter()\n",
    "simind_true = simind_total - simind_scatter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the acquisition model\n",
    "acq_model = get_acquisition_model(measured_data, measured_additive, image, mu_map_stir)\n",
    "stir_forward_projection = acq_model.forward(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display counts\n",
    "print(f\"simind total counts: {simind_total.sum()}\")\n",
    "print(f\"simind true counts: {simind_true.sum()}\")\n",
    "print(f\"simind scatter counts: {simind_scatter.sum()}\")\n",
    "print(\"\\n\")\n",
    "print(f\"measured total counts: {measured_data.sum()}\")\n",
    "print(f\"stir true counts: {stir_forward_projection.sum()}\")\n",
    "print(f\"measured additive counts: {measured_additive.sum()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the data to be plotted\n",
    "data_list = [\n",
    "    ((simind_total), \"simind total\"),\n",
    "    ((simind_true), \"simind true\"),\n",
    "    ((simind_scatter), \"simind scatter\"),\n",
    "    ((measured_data), \"measured\"),\n",
    "    ((stir_forward_projection), \"stir forward\"),\n",
    "    ((measured_additive), \"measured additive\")\n",
    "]\n",
    "\n",
    "data_list = [(data.as_array(), title) for data, title in data_list]\n",
    "\n",
    "# Plot the axial slice\n",
    "axial_slice = 66\n",
    "vmax = max([data[0][axial_slice].max() for data, _ in data_list])\n",
    "\n",
    "# Define consistent font size and colormap\n",
    "font_size = 14\n",
    "colormap = 'viridis'\n",
    "\n",
    "# Create a figure and a GridSpec with 3 rows\n",
    "fig = plt.figure(figsize=(len(data_list) * 4, 7 * 2))\n",
    "gs = GridSpec(3, len(data_list), height_ratios=[2, 0.15, 3])  # Adjusted GridSpec for clarity\n",
    "\n",
    "# Create image subplots\n",
    "ax_images = [fig.add_subplot(gs[0, i]) for i in range(len(data_list))]\n",
    "\n",
    "for i, (data, title) in enumerate(data_list):\n",
    "    im = ax_images[i].imshow(data[0, axial_slice], vmin=0, vmax=vmax, cmap=colormap)\n",
    "    ax_images[i].set_title(f\"{title}: {np.trunc(data.sum())} \", fontsize=font_size)\n",
    "    ax_images[i].axis('off')\n",
    "\n",
    "# Create a colorbar\n",
    "cbar_ax = fig.add_subplot(gs[1, :])\n",
    "fig.colorbar(im, cax=cbar_ax, orientation='horizontal', pad=0.02)\n",
    "cbar_ax.set_xlabel('Counts', fontsize=font_size)\n",
    "cbar_ax.xaxis.set_label_position('top')\n",
    "\n",
    "# Plot the line plot\n",
    "ax_line = fig.add_subplot(gs[2, :])\n",
    "colours = plt.cm.viridis(np.linspace(0, 1, len(data_list)))\n",
    "for i, (data, title) in enumerate(data_list):\n",
    "    ax_line.plot(data[0, axial_slice][60], linewidth=2, color=colours[i], linestyle='-', label=title)\n",
    "\n",
    "# Enhance the line plot appearance\n",
    "ax_line.set_xlabel('Projection angle', fontsize=font_size)\n",
    "ax_line.set_ylabel('Intensity', fontsize=font_size)\n",
    "ax_line.set_title(f'Profile Through Sinogram', fontsize=font_size + 2)\n",
    "ax_line.grid(True, which='both', linestyle='--', linewidth=0.5)\n",
    "ax_line.legend(loc='upper left', fontsize=font_size)\n",
    "ax_line.set_xlim(0, 128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Adjust layout and save the figure\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(output_dir, \"comparison_axial.png\"))\n",
    "plt.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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": 2
}