<jupyter_start><jupyter_text>Writing OSEM (or another reconstruction algorithm) yourself with SIRFThis notebook invites you to write MLEM and OSEM yourself using SIRF functionality, i.e. Do It Yourself OSEM!You should have completed the [OSEM_reconstruction notebook](OSEM_reconstruction.ipynb) first. The [ML_reconstruction notebook](ML_reconstruction.ipynb) might help as well.The notebook is currently set-up to use prepared data with a single slice of an XCAT phantom, with a low resolution scanner, such that all results can be obtained easily on a laptop. Of course, your code will have to work for any data. Authors: <NAME>  First version: June 2021CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  Copyright 2021 University College London.This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging (http://www.ccpsynerbi.ac.uk/).Licensed under the Apache License, Version 2.0 (the "License");you may not use this file except in compliance with the License.You may obtain a copy of the License at    http://www.apache.org/licenses/LICENSE-2.0Unless required by applicable law or agreed to in writing, softwaredistributed under the License is distributed on an "AS IS" BASIS,WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.See the License for the specific language governing permissions andlimitations under the License.  Initial set-up<jupyter_code>#%% make sure figures appears inline and animations works
%matplotlib notebook

# Setup the working directory for the notebook
import notebook_setup
from sirf_exercises import cd_to_working_dir
cd_to_working_dir('PET', 'OSEM_reconstruction')
#%% Initial imports etc
import numpy
import matplotlib.pyplot as plt
import os
import sirf.STIR as pet
from sirf.Utilities import examples_data_path
from sirf_exercises import exercises_data_path

# define the directory with input files for this notebook
data_path = os.path.join(examples_data_path('PET'), 'thorax_single_slice')
# set-up redirection of STIR messages to files
msg_red = pet.MessageRedirector('info.txt', 'warnings.txt', 'errors.txt')<jupyter_output><empty_output><jupyter_text>We will first create some simulated data from ground-truth imagessee previous notebooks for more information.<jupyter_code>#%% Read in images
image = pet.ImageData(os.path.join(data_path, 'emission.hv'))
attn_image = pet.ImageData(os.path.join(data_path, 'attenuation.hv'))
#%% save max for future displays
im_slice = image.dimensions()[0]//2
cmax = image.max()*.6
#%% create acquisition model
acq_model = pet.AcquisitionModelUsingRayTracingMatrix()
template = pet.AcquisitionData(os.path.join(data_path, 'template_sinogram.hs'))
acq_model.set_up(template, image)
#%% simulate data using forward projection
acquired_data = acq_model.forward(image)<jupyter_output><empty_output><jupyter_text>Maximum Likelihood Expectatin Maximisation (MLEM)Also called EMML. This is a standard algorithm, derived by using EM for the PET (or SPECT) problem. See the paper:<NAME>., and <NAME>. ‘Maximum Likelihood Reconstruction for Emission Tomography’. IEEE Transactions on Medical Imaging 1, no. 2 (1982): 113-122+.Our notation here: $x$ is the image, $y$ the measured data with $A$ the system matrix. This is different from the Shepp and Vardi paper, which uses $\lambda$ for the image, $n^*$ for the measured data, $p(b,d)$ for the elements of the system matrix, and it has no background.In our notation, the model for the mean of the data (i.e., modelling the expected measurement, given an image $x$) is$$ \bar{y}=A x + b$$The MLEM update is$$ x^{\mathrm{new}} = \frac{x}{A^t 1} A^t \left(\frac{y}{A x + b}\right)$$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!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).<jupyter_code>sensitivity = acq_model.backward(acquired_data.get_uniform_copy(1))<jupyter_output><empty_output><jupyter_text>Now we initialise the algorithm with a uniform image:<jupyter_code>estimated_image = image.get_uniform_copy(1)<jupyter_output><empty_output><jupyter_text>and we can do one MLEM iteration<jupyter_code>quotient = acquired_data/acq_model.forward(estimated_image)  # y / (Ax + b)
mult_update = acq_model.backward(quotient)/sensitivity       # A^t * quotient / A^t1
estimated_image *= mult_update                               # update (in place)<jupyter_output><empty_output><jupyter_text>And we can do some plots<jupyter_code>plt.figure()
plt.subplot(1,2,1)
plt.imshow(estimated_image.as_array()[im_slice,:,:])
plt.subplot(1,2,2)
plt.imshow(mult_update.as_array()[im_slice,:,:])<jupyter_output><empty_output><jupyter_text>Now you can of course duplicate some of these lines, or re-execute the above cells. However, it makes more sense to write a function to do this. Something like this:<jupyter_code>def MLEM(acquired_data, acq_model, initial_image, num_iterations):
    estimated_image = initial_image.clone()
    #  some stuff here
    return estimated_image<jupyter_output><empty_output><jupyter_text>And now you can run it!<jupyter_code>estimated_image = MLEM(acquired_data, acq_model,image.get_uniform_copy(1), 40)<jupyter_output><empty_output><jupyter_text>This was hopefully not too hard. Theree are a few problems though that you might encounter.- your image might display nice, but on closer investigation will probably contain NaNs (Not a Number). These come from divisions: 0/0 is not defined. They can occur in 2 places:    - division in the acquisition data term. This should in theory not occur if you start with a strictly positive image that is large enough to "cover" all of the projection data. Of course, in practice it could happen that your image is smaller than the Field of View (FOV), and you might need to handle this anyway.    - division of the images. This will occur wherever the sensitivty image is zero. This difficulty is not a surprise: if a voxel cannot be measured, its ML estimate is undefined.  We have the second problem of course, as by default the projector uses a circular FOV. You might wat to add a post-processing step that sets those values to zero. The STIR implementation of MLEM (`OSMAPOSL`) takes care of these corner cases, as well as any negatives in the data (arising when pre-correcting the data, as although this should not be done in theory, some people are interested in this anyway).  Ordered Subsets Expectation Maximisation (OSEM)Is discussed in previous notebooks, MLEM is great but slow. OSEM was introduced inHudson, H.M., and <NAME>. ‘Accelerated Image Reconstruction Using Ordered Subsets of Projection Data’. IEEE Transactions on Medical Imaging 13, no. 4 (December 1994): 601–9. https://doi.org/10.1109/42.363108.The update formula is exactly the same, except that at every update, only a subset of the data is used, i.e. restricting the data $y$, background $b$ and the matrix $A$ to a subset of all the rows. Clearly, for $N$ subsets, this reduces the number of computations for one image update with a factor $N$. While each update might be somewhat less accurate, it certainly works well in initial iterations. So how do we implement this in SIRF? Luckily, an `sirf.STIR` acquisition model can be told to use only a subset of the data. The class documentation will show you that you can set `num_subsets` and `subset_num`.(There is currently no way to change how the subsets are chosen, but only the number of subsets).Note: the call to `forward` and `backward` also supports extra arguments for specifying the subsets. However, these are deprecated and will be removed in a future release.Some interesting things to try:  - how does varying the `subset_num` affect the projection?  - how does varying the `num_subsets` affect the projection sparsity?  - what happens if you sum over all the subsets?<jupyter_code>acq_model.num_subsets = 4
acq_model.subset_num = 0  # for 4 subsets, use a number between 0 and 3
data = acq_model.forward(image)
data.show(im_slice)<jupyter_output><empty_output><jupyter_text>Unfortunately, SIRF currently has no way to restrict the data to a particular subset. There are 2 ways around that:- ignore it and do the divisions over all of the data anyway. This will lead to 1/0 etc, but as those elements of the data are not backprojected, it won't change the final image.- construct several "masks" as `AcquisitionData by forward projection an image full of ones for every subsets and doing some thresholding.    Clearly, the first option is easiest (although it does mean there is some overhead in computing extra additions/divisions). Let's see if it works ok!<jupyter_code>check = acq_model.backward(acquired_data/data)
check.show(im_slice)<jupyter_output><empty_output><jupyter_text>You should now be in a position to write your own OSEM algorithm. Don't forget that for a strict implementation of OSEM, you need to compute "subset sensitivities" for the division.<jupyter_code>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 
    return estimated_image<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07268632268632269, 'fraction_numerical': 0.007985257985257985, 'mean_word_length': 5.335278858625162, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 49, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '21225', 'license': 'Apache-2.0', 'n_tokens_mistral': 2862, 'n_tokens_neox': 2696, 'n_words': 1225, 'path': 'notebooks/PET/DIY_OSEM.ipynb', 'repo_name': 'mastergari/SIRF-Exercises', 'size': '9789'}