notebooks/PET/reconstruct_measured_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reconstruct phantom data\n",
    "This exercise shows how to handle data from the Siemens mMR. It shows how to get from listmode data to sinograms, get a randoms estimate, and reconstruct using normalisation, randoms, attenuation and scatter.\n",
    "\n",
    "It is recommended you complete the [OSEM_reconstruction notebook](OSEM_reconstruction.ipynb) first. Even better would be to look at the OSEM part of the [ML_reconstruct notebook](ML_reconstruct.ipynb) as well.\n",
    "\n",
    "This exercise uses data from a phantom acquisition at UCL on a Siemens mMR. The phantom is the NEMA phantom (essentially a torso-shaped perspex box, with some spherical inserts). You will need to download that data. Please use the read INSTALL.md or DocForParticipants.md for details. \n",
    "\n",
    "The script should work for other data of course, but you will need to adapt filenames.\n",
    "\n",
    "You can also adjust it to use other reconstruction algorithms than OSEM with very little changes.\n",
    "\n",
    "Note that we currently don't show here how to extract the data from the console. Please\n",
    "[check our wiki for more information](https://github.com/CCPPETMR/SIRF/wiki/PET-raw-data)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Authors: Kris Thielemans and Evgueni Ovtchinnikov  \n",
    "First version: 8th of September 2016  \n",
    "vSecond Version: 17th of May 2018-  \n",
    "Third version: June 2021\n",
    "\n",
    "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  \n",
    "Copyright 2015 - 2017, 2021 Rutherford Appleton Laboratory STFC.  \n",
    "Copyright 2015 - 2018, 2020-2021 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": "markdown",
   "metadata": {},
   "source": [
    "# Initial set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% make sure figures appears inline and animations works\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run `download_data.sh -p` if you didn't yet!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "bash ../../scripts/download_data.sh -p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the working directory for the notebook\n",
    "import notebook_setup\n",
    "from sirf_exercises import cd_to_working_dir\n",
    "cd_to_working_dir('PET', 'reconstruct_measured_data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "from sirf.Utilities import show_2D_array, examples_data_path\n",
    "from sirf.STIR import *\n",
    "from sirf_exercises import exercises_data_path\n",
    "\n",
    "# Find the path to the data directory\n",
    "data_path = exercises_data_path('PET', 'mMR', 'NEMA_IQ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert the Siemens Normalisation into STIR format, some keys slightly differ\n",
    "# $data_path/20170809_NEMA[_MUMAP]_UCL.n.hdr are the extracted headers from the scanner\n",
    "# norm.n.hdr and umap.v.hdr are the output converted files\n",
    "#!convertSiemensInterfileToSTIR.sh $data_path/20170809_NEMA_UCL.n.hdr norm.n.hdr\n",
    "!cp $data_path/20170809_NEMA_UCL.n.hdr norm.n.hdr\n",
    "!convertSiemensInterfileToSTIR.sh $data_path/20170809_NEMA_MUMAP_UCL.v.hdr umap.v.hdr\n",
    "\n",
    "# These files will sometimes have lines terminated with CR instead of CRLF, fix this\n",
    "# This just means, if we see CR<character> and <character> isn't LF, replace it with CRLF<character>\n",
    "!sed -i.bak \"s/\\r\\([^\\n]\\)/\\r\\n\\1/g\" norm.n.hdr\n",
    "!sed -i.bak \"s/\\r\\([^\\n]\\)/\\r\\n\\1/g\" umap.v.hdr\n",
    "\n",
    "# Now add absolute data path to the header file\n",
    "# This command prepends the data path to the data file so that the header in our working folder points to the data\n",
    "# You won't need to do this for your own data if the data file is in the same directory.\n",
    "!sed -i.bak2 -e \"s#\\(!name of data file:=\\)#\\\\1{data_path}/#\" umap.v.hdr\n",
    "!sed -i.bak2 -e \"s#\\(!name of data file:=\\)#\\\\1{data_path}/#\" norm.n.hdr\n",
    "\n",
    "# Advanced: if you'd like to have a look at what changed in the umap, uncomment below\n",
    "# Lines starting with < are the original Siemens\n",
    "# and lines starting with > are the STIR converted\n",
    "# !diff $data_path/20170809_NEMA_MUMAP_UCL.v.hdr umap.v.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the above failed, it probably means that `convertSiemensInterfileToSTIR.sh` wasn't installed. The following line might help."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !cp $SIRF_PATH/../STIR/scripts/IO/convertSiemensInterfileToSTIR.sh $SIRF_INSTALL_PATH/bin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% set filenames \n",
    "# input files\n",
    "list_file = os.path.join(data_path, '20170809_NEMA_60min_UCL.l.hdr')\n",
    "norm_file = 'norm.n.hdr'\n",
    "attn_file = 'umap.v.hdr'\n",
    "# output filename prefixes\n",
    "sino_file = 'sino'\n",
    "%ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# redirect STIR messages to some files\n",
    "# you can check these if things go wrong\n",
    "_ = MessageRedirector('info.txt', 'warnings.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating sinograms from listmode data\n",
    "Modern PET scanners can store data in listmode format. This is essentially a long list of all events detected by the scanner. We are interested here in the *prompts* (the coincidence events) and the *delayed events* (which form an estimate of the *accidental coincidences* in the prompts.\n",
    "\n",
    "We show how to histogram the prompts into a sinogram etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First create a template for the sinogram\n",
    "This template is used to specify the sizes of the output sinogram.\n",
    "\n",
    "It is often the case in PET that we use sinograms with \"larger\" bins, i.e. combine data from several detector pairs into a single bin. This reduces size of the final sinogram, and decreases computation time. The terminology here is somewhat complicated, but *span* uses \"axial compression\" (higher span means smaller data size), *max_ring_diff* specifies the maximum ring difference to store, and *view_mash_factor* can be used to reduce the number of views (or azimutal angles). You could check the [STIR glossary](http://stir.sourceforge.net/documentation/STIR-glossary.pdf) for more detail.\n",
    "\n",
    "Siemens uses span=11, max_ring_diff=60 and view_mash_factor=1 for the mMR. Here we will use a smaller data size to reduce computation time for the exercise. Feel free to change these numbers (if you know what you are doing...). (Note that the list mode data stores data only up to ring difference 60, even though the scanner has 64 rings)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "template_acq_data = AcquisitionData('Siemens_mMR', span=11, max_ring_diff=15, view_mash_factor=2)\n",
    "template_acq_data.write('template.hs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create listmode-to-sinograms converter object\n",
    "lm2sino = ListmodeToSinograms()\n",
    "\n",
    "# set input, output and template files\n",
    "lm2sino.set_input(list_file)\n",
    "lm2sino.set_output_prefix(sino_file)\n",
    "lm2sino.set_template('template.hs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set timing interval (in secs) since start of acquisition\n",
    "# (the listmode file provided is for 1 hour).\n",
    "# you can vary this to see the effect on noise. Increasing it will mean somewhat longer\n",
    "# processing time in the following steps (but not in the reconstruction).\n",
    "lm2sino.set_time_interval(0, 600)  # 0 - 600 is the first 10 minutes\n",
    "# set up the converter\n",
    "lm2sino.set_up()\n",
    "# create the prompts sinogram\n",
    "lm2sino.process()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the content of the directory. there should be a `sino*.hs`, `'.s` pair.\n",
    "# The `.hs` file is an Interfile header pointing to the binary data.\n",
    "%ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check the prompts sinograms\n",
    "The 3D PET data returned by `as_array` are organised by 2D sinogram. The exact order of the sinograms\n",
    "is complicated for 3D PET, but SIRF orders them by *segment* (roughly: average ring difference). The first\n",
    "segment corresponds to \"segment 0\", i.e. detector pairs which are (roughly) in the same \n",
    "detector ring. For a clinical scanner with `N` rings, there will be `2N-1` (2D) sinograms in segment 0. See also information in the [image_creation_and_simulation notebook](image_creation_and_simulation.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get access to the sinograms\n",
    "acq_data = lm2sino.get_output()\n",
    "# copy the acquisition data into a Python array\n",
    "acq_array = acq_data.as_array()[0,:,:,:]  # first index is for ToF, which we don't have here\n",
    "# how many counts total?\n",
    "print('num prompts: %d' % acq_array.sum())\n",
    "# print the data sizes. \n",
    "print('acquisition data dimensions: %dx%dx%d' % acq_array.shape)\n",
    "# use a slice number for display that is appropriate for the NEMA phantom\n",
    "# showing a \"middle\" sinogram in segment 0.\n",
    "z = 71\n",
    "show_2D_array('Acquisition data', acq_array[z,:,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also show a horizontal profile through this sinogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(acq_array[z,0,:],'.')\n",
    "plt.xlabel('radial distance')\n",
    "plt.ylabel('counts')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, there seems to be lots of things \"wrong\" with this sinogram. In the next few sections, we will show how to incorporate various effects into the acquisition model. We will also reconstruct the data with these increasingly accurate models for the acquisition as illustration.\n",
    "\n",
    "For simplicity, we will use OSEM and use only a few sub-iterations for speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create initial image estimate of dimensions and voxel sizes\n",
    "# compatible with the scanner geometry (included in the AcquisitionData\n",
    "# object acq_data) and initialize each voxel to 1.0\n",
    "nxny = (127, 127)\n",
    "initial_image = acq_data.create_uniform_image(1.0, nxny)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reconstruct the prompt data without any extra modelling\n",
    "The simplest model that you can use for PET data is line integrals. So let's do that here and see how well it does.\n",
    "\n",
    "In SIRF, we currently have 2 options:\n",
    "- the ray tracing matrix multiplication which we have used in other demos\n",
    "- Alternatively, you could use the [parallelproj acquisition model](https://github.com/gschramm/parallelproj/) which is faster if you have SIRF/STIR with GPU capability. However, note that our wrapper of parallelproj currently is inefficient for a large number of subsets, so you might want to reduce the number of subsets below when using this choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  by ray tracing\n",
    "acq_model = AcquisitionModelUsingRayTracingMatrix()\n",
    "acq_model.set_num_tangential_LORs(10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using parallelproj\n",
    "# acq_model = AcquisitionModelUsingParallelproj()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now tell the acquisition model more about the geometry of the data it will need to handle\n",
    "acq_model.set_up(acq_data, initial_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define objective function to be maximized as\n",
    "# Poisson logarithmic likelihood (with linear model for mean)\n",
    "obj_fun = make_Poisson_loglikelihood(acq_data)\n",
    "obj_fun.set_acquisition_model(acq_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the reconstruction object\n",
    "recon = OSMAPOSLReconstructor()\n",
    "recon.set_objective_function(obj_fun)\n",
    "\n",
    "# Choose a number of subsets.\n",
    "# For the mMR, best performance requires to not use a multiple of 9 as there are gaps\n",
    "# in the sinograms, resulting in unbalanced subsets (which isn't ideal for OSEM).\n",
    "num_subsets = 21\n",
    "# Feel free to increase these.\n",
    "# (Clinical reconstructions use around 60 subiterations, e.g. 21 subsets, 3 full iterations)\n",
    "num_subiterations = 12\n",
    "recon.set_num_subsets(num_subsets)\n",
    "recon.set_num_subiterations(num_subiterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = initial_image.clone()\n",
    "recon.set_up(image)\n",
    "# set the initial image estimate\n",
    "recon.set_current_estimate(image)\n",
    "# reconstruct\n",
    "recon.process()\n",
    "# show reconstructed image\n",
    "image_array = recon.get_output().as_array()\n",
    "show_2D_array('Reconstructed image', image_array[z,:,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously this image doesn't look great. We have used an acquisition model that just uses line integrals. Clearly not good enough!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add detector sensitivity modelling\n",
    "Probably the obvious feature of the sinogram above are the diagonal lines. There are due to 2 effects:\n",
    "- Each crystal pair will have different detection efficiency.\n",
    "- The Siemens mMR (and other scanners) has some \"gaps\" between detector blocks, which are accomodated in the sinograms by creating a \"virtual\" crystal between each block. These will never detect any counts.\n",
    "\n",
    "We need to take this into account\n",
    "in our acquisition model. The scanner provides a *normalisation file* to do this (the terminology\n",
    "originates from the days that we were \"normalising\" by  dividing by the detected counts \n",
    "by the sensitivities. Note that this cannot be done with the \"virtual\" crystals as it would lead to 0/0).\n",
    "\n",
    "In SIRF, you can incorporate this effect in the acquisition model by using an `AcquisitionSensitivityModel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create it from the supplied file\n",
    "asm_norm = AcquisitionSensitivityModel(norm_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this to create sinogram data with the detection efficiencies. Let's do that here and display."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "asm_norm.set_up(acq_data)\n",
    "det_efficiencies=acq_data.get_uniform_copy(1)\n",
    "asm_norm.unnormalise(det_efficiencies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(det_efficiencies.as_array()[0,z,:,:],clim=None)\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(det_efficiencies.as_array()[0,z,0,:],'.')\n",
    "plt.xlabel('radial distance')\n",
    "plt.ylabel('efficiency')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the scale of these detection efficiencies is a bit arbitrary. At this point in time (SIRF 3.4 with STIR 5.1), a global calibration factor is not yet taken into account for instance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see the same diagonal patterns here. Therefore including this into our acquisition model will likely give much better results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add it to the acquisition model\n",
    "acq_model.set_acquisition_sensitivity(asm_norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acq_model.set_up(acq_data, initial_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update the objective function\n",
    "obj_fun.set_acquisition_model(acq_model)\n",
    "recon.set_objective_function(obj_fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reconstruct\n",
    "recon.set_up(initial_image)\n",
    "recon.set_current_estimate(initial_image)\n",
    "recon.process()\n",
    "# show reconstructed image\n",
    "image_array = recon.get_output().as_array()\n",
    "show_2D_array('Reconstructed image', image_array[z,:,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_image.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add attenuation modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read attenuation image\n",
    "attn_image = ImageData(attn_file)\n",
    "z = 71\n",
    "attn_image.show(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "attn_acq_model = AcquisitionModelUsingRayTracingMatrix()\n",
    "asm_attn = AcquisitionSensitivityModel(attn_image, attn_acq_model)\n",
    "# converting attenuation into attenuation factors (see previous exercise)\n",
    "asm_attn.set_up(acq_data)\n",
    "attn_factors = acq_data.get_uniform_copy(1)\n",
    "print('applying attenuation (please wait, may take a while)...')\n",
    "asm_attn.unnormalise(attn_factors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use these in the final attenuation model\n",
    "asm_attn = AcquisitionSensitivityModel(attn_factors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's display the attenuation factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(attn_factors.as_array()[0,z,:,:],clim=None)\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(attn_factors.as_array()[0,z,0,:],'.')\n",
    "plt.xlabel('radial distance')\n",
    "plt.ylabel('attenuation factor')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have two acquisition_sensitivity_models: for detection sensitivity and for\n",
    "count loss due to attenuation. We combine them by \"chaining\" them together (which will\n",
    "model the multiplication of both sensitivities)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# chain attenuation and normalisation\n",
    "asm = AcquisitionSensitivityModel(asm_norm, asm_attn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update the acquisition model etc\n",
    "acq_model.set_acquisition_sensitivity(asm)\n",
    "obj_fun.set_acquisition_model(acq_model)\n",
    "recon.set_objective_function(obj_fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reconstruct\n",
    "recon.set_up(initial_image)\n",
    "recon.set_current_estimate(initial_image)\n",
    "recon.process()\n",
    "# show reconstructed image\n",
    "image_array = recon.get_output().as_array()\n",
    "show_2D_array('Reconstructed image', image_array[z,:,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add a background term for modelling the randoms\n",
    "PET data includes \"accidental coincidences\" (often called \"randoms\"). These occur when annihilation phantoms of 2 different annihilations are detected within the coincidence window. This gives a global background to the data. So we need to model this!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate the *randoms* background\n",
    "Siemens stores *delayed coincidences*. These form a very noisy estimate of the\n",
    "background due to accidental coincidences in the data. However, that estimate is too noisy\n",
    "to be used in iterative image reconstruction. Note that the acquisition model should give you an estimate of the \"mean\" of the data (i.e. is noiseless).\n",
    "\n",
    "SIRF uses an algorithm from STIR that gives a much less noisy estimate. The help message \n",
    "gives some information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(lm2sino)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the randoms estimate\n",
    "# This will take a while\n",
    "randoms = lm2sino.estimate_randoms()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the randoms-estimate\n",
    "A (2D) sinogram of the randoms has a similar pattern diagonal lines. This is related to the\n",
    "detector efficiencies of course, but we cannot get into that here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(randoms.as_array()[0,z,:,:],clim=None)\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(randoms.as_array()[0,z,0,:],'.')\n",
    "plt.xlabel('radial distance')\n",
    "plt.ylabel('estimated randoms')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Include the randoms estimate into the acquisition model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acq_model.set_background_term(randoms)\n",
    "obj_fun.set_acquisition_model(acq_model)\n",
    "recon.set_objective_function(obj_fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "recon.set_current_estimate(initial_image)\n",
    "recon.set_up(initial_image)\n",
    "recon.process()\n",
    "# show reconstructed image\n",
    "image_array = recon.get_output().as_array()\n",
    "show_2D_array('Reconstructed image', image_array[z,:,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scatter\n",
    "Finally, some of the detected counts will be from coincidences where one (or both) of the annihilation photons have scattered. Although the scanner tries to reject these by using energy windowing, there will still be a non-negligible fraction. This is due to the energy resolution of current PET scanners, and trade-offs made between detecting scattered counts and rejecting unscattered counts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most common way to estimate scatter in PET is \"model-based\". Essentially it is an iterative loop between image estimation and modelling the scatter based on the current image estimate (and the attenuation image). We can run the STIR implementation from inside SIRF as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "se = ScatterEstimator()\n",
    "\n",
    "se.set_input(acq_data)\n",
    "se.set_attenuation_image(attn_image)\n",
    "se.set_randoms(randoms)\n",
    "se.set_asm(asm_norm)\n",
    "# Unfortunately, the ScatterEstimator currently needs attenuation \"correction\" factors, which\n",
    "# is what we need to multiply by to correct for attenuation, while we computed the attenuation\n",
    "# factors above.\n",
    "# Fortunately, these are simply the inverse.\n",
    "acf_factors = attn_factors.get_uniform_copy()\n",
    "acf_factors.fill(1/attn_factors.as_array())\n",
    "# I could also have used the following (but it would take more time)\n",
    "#asm_attn.normalise(acf_factors)\n",
    "se.set_attenuation_correction_factors(acf_factors)\n",
    "\n",
    "# set the number of iterations used for the scatter algorithm.\n",
    "# The default is 5, but 3 is often enough, so we will use that here to reduce computation time.\n",
    "se.set_num_iterations(3)\n",
    "# optionally let it write intermediate scatter estimates to file\n",
    "se.set_output_prefix('scatter_estimate')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# go and compute it! (might take a minute or 2)\n",
    "se.set_up()\n",
    "se.process()\n",
    "scatter_estimate = se.get_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's display it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(scatter_estimate.as_array()[0,z,:,:],clim=None)\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(scatter_estimate.as_array()[0,z,0,:],'.')\n",
    "plt.xlabel('radial distance')\n",
    "plt.ylabel('estimated scatter')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the scatter estimate is a fairly smooth background, where the detection efficiencies are again superimposed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## reconstruct including scatter and all other terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acq_model.set_background_term(randoms+scatter_estimate)\n",
    "acq_model.set_up(acq_data, initial_image)\n",
    "obj_fun.set_acquisition_model(acq_model)\n",
    "recon.set_objective_function(obj_fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "recon.set_up(initial_image)\n",
    "recon.set_current_estimate(initial_image)\n",
    "recon.process()\n",
    "# show reconstructed image\n",
    "image_array = recon.get_output().as_array()\n",
    "show_2D_array('Reconstructed image', image_array[z,:,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary of acquisition modelling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We needed\n",
    "- detection efficiency, attenuation (both are multiplicative effects and incorporated via an `AcquisitionSensitivityModel`)\n",
    "- randoms and scatter (both are additive effects and incorporated by adding a background term)\n",
    "\n",
    "Let's see how well the two latter terms fit the measure data. This is easiest to check \"outside\" the body, i.e. where we did actually not expect any counts at all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(acq_data.as_array()[0,z,0,:],'.',markersize=2,label='prompts')\n",
    "plt.plot(randoms.as_array()[0,z,0,:],'og',markersize=4,label='randoms estimate')\n",
    "plt.plot((randoms+scatter_estimate).as_array()[0,z,0,:],'xr',markersize=2,label='randoms+scatter estimate')\n",
    "plt.xlabel('radial distance')\n",
    "plt.ylabel('counts')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This doesn't look so good, but of course, there is a lot of noise in the prompt data (and you can see that the detected counts are 0,1,2,...).\n",
    "\n",
    "So let's sum over all sinograms to reduce noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(numpy.sum(acq_data.as_array()[0,:,0,:],axis=0),'.',markersize=2,label='prompts')\n",
    "plt.plot(numpy.sum(randoms.as_array()[0,:,0,:],axis=0),'og',markersize=4,label='randoms estimate')\n",
    "plt.plot(numpy.sum((randoms+scatter_estimate).as_array()[0,:,0,:],axis=0),'xr',markersize=2,label='randoms+scatter estimate')\n",
    "plt.xlabel('radial distance')\n",
    "plt.ylabel('counts')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What now?\n",
    "Here are some suggestions for things to try:\n",
    "- now that you have a good acquisition model, you probably wnat to increase the number of subiterations a bit to get a better quality image.\n",
    "- change the duration used when going from listmode data to sinograms (longer duration, more counts, less noise), and run at `span=11` and `view_mashing=1` (which is what the mMR does clinically).\n",
    "- use the final `acq_model` to forward project the reconstructed image. Does it fit the data?\n",
    "- write a function that takes the listmode data, normalisation file and attenuation image, computes `acq_data` and the acquisition model. Once you have this, you can use your `OSEM` function from a previous notebook (or any other reconstruction method). You can then create a Python script that does everything from start to finish."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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": 4
}