src/notebooks/BSREM_illustration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Demonstration of basic BSREM implementation with SIRF\n",
    "\n",
    "`BSREM1` and `BSREM2` are 2 different implementations of a (modified) BSREM algorithm. `BSREM1` uses the `sirf.STIR` objective function to compute gradients etc, while `BSREM2` computes these in terms of the `sirf.STIR.AcquisitionModel`. The actual algorithm is implemented in the base class `BSREMSkeleton`. Note that the implementations are in terms of \"lists of subset-data\" and corresponding \"list of acquisition models\" etc. (This is more efficient than using the subset functionality in the `sirf.STIR` acquisition models/objective function.)\n",
    "\n",
    "The source code for these implementations can be found in `src/Python`, or you could run `BSREM1??`\n",
    "\n",
    "This notebook is illustrates these algorithms with two datasets:\n",
    "- a simple 2D simulation (largely based on the `display_and_projection` notebook from the SIRF-Exercises)\n",
    "- the mMR NEMA IQ data (downloadable via the SIRF-exercises as well)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "Author: Kris Thielemans  \n",
    "\n",
    "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  \n",
    "Copyright 2015 - 2017 Rutherford Appleton Laboratory STFC.  \n",
    "Copyright 2015 - 2018, 2021, 2024 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": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% make sure figures appears inline and animations works\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Initial imports etc\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import sys\n",
    "import sirf.STIR as STIR\n",
    "from sirf.Utilities import examples_data_path\n",
    "#from sirf_exercises import exercises_data_path"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "Import functionality from the Python files in SIRF-Contribs.\n",
    "(Note that in most set-ups, this will be from the installed files.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sirf.contrib.partitioner import partitioner\n",
    "from sirf.contrib.BSREM.BSREM import BSREM1\n",
    "from sirf.contrib.BSREM.BSREM import BSREM2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Needed for get_subsets()\n",
    "STIR.AcquisitionData.set_storage_scheme('memory')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set-up redirection of STIR messages to files\n",
    "_ = STIR.MessageRedirector('info.txt', 'warnings.txt', 'errors.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fewer message from STIR and SIRF (set to higher value if you have problems)\n",
    "STIR.set_verbosity(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "# Some function definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% some handy function definitions\n",
    "def plot_2d_image(idx,vol,title,clims=None,cmap=\"viridis\"):\n",
    "    \"\"\"Customized version of subplot to plot 2D image\"\"\"\n",
    "    plt.subplot(*idx)\n",
    "    plt.imshow(vol,cmap=cmap)\n",
    "    if not clims is None:\n",
    "        plt.clim(clims)\n",
    "    plt.colorbar(shrink=.4)\n",
    "    plt.title(title)\n",
    "    plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_acq_model_and_obj_fun(acquired_data, additive_term, mult_factors, template_image):\n",
    "    '''\n",
    "    Create an acquisition model and objective function, corresponding to the given data.\n",
    "    '''\n",
    "    # We could construct this by hand here, but instead will just use `partitioner.data_partition`\n",
    "    # with 1 subset, which will then do the work for us.\n",
    "    num_subsets = 1\n",
    "    _, acq_models, obj_funs = partitioner.data_partition(acquired_data, additive_term, mult_factors, num_subsets, initial_image=template_image)\n",
    "    return (acq_models[0], obj_funs[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "def scale_initial_image(acquired_data, additive_term, mult_factor, template_image, obj_fun):\n",
    "    '''\n",
    "    Return a uniform image that has a reasonable \"scale\" (i.e. image values) for the data given.\n",
    "    \n",
    "    If there is an additive term, OSEM can be a bit slow to converge if the initial image has very wrong\n",
    "    image values. Here we find a scale such that the sum of the forward projection of the initial image is equal to the sum of the acquired data.\n",
    "    \n",
    "    WARNING: assumes that obj_fun has been set_up already\n",
    "    '''\n",
    "    data_sum = (acquired_data.sum() - (additive_term * mult_factors).sum())\n",
    "    ratio = data_sum / (obj_fun.get_subset_sensitivity(0).sum() * obj_fun.get_num_subsets())\n",
    "    return template_image.allocate(ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "def OSEM(obj_fun, initial_image, num_updates = 14, num_subsets = 2):\n",
    "    '''\n",
    "    run OSEM\n",
    "    \n",
    "    WARNING: this modified the `obj_fun` by setting its number of subsets. This is unfortunate of course.\n",
    "    '''\n",
    "    recon = STIR.OSMAPOSLReconstructor()\n",
    "    recon.set_objective_function(obj_fun)\n",
    "    recon.set_current_estimate(initial_image)\n",
    "    recon.set_num_subsets(num_subsets )\n",
    "    recon.set_num_subiterations(num_updates)\n",
    "    recon.set_up(initial_image)\n",
    "    recon.process()\n",
    "    return recon.get_output()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_RDP(penalty_strength, initial_image, kappa, max_scaling=1e-3):\n",
    "    '''\n",
    "    Construct a smoothed Relative Difference Prior (RDP)\n",
    "    \n",
    "    `initial_image` is used to determine a smoothing factor (epsilon).\n",
    "    `kappa` is used to pass voxel-dependent weights.\n",
    "    '''\n",
    "    prior = STIR.RelativeDifferencePrior()\n",
    "    # need to make it differentiable\n",
    "    epsilon = initial_image.max() * max_scaling\n",
    "    prior.set_epsilon(epsilon)\n",
    "    prior.set_penalisation_factor(penalty_strength)\n",
    "    prior.set_kappa(kappa)\n",
    "    prior.set_up(initial_image)\n",
    "    return prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_kappa_image(obj_fun, initial_image):\n",
    "    '''\n",
    "    Computes a \"kappa\" image for a prior as sqrt(H.1). This will attempt to give uniform \"perturbation response\".\n",
    "    See Yu-jung Tsai et al. TMI 2020 https://doi.org/10.1109/TMI.2019.2913889\n",
    "\n",
    "    WARNING: Assumes the objective function has been set-up already\n",
    "    '''\n",
    "    # This needs SIRF 3.7. If you don't have that yet, you should probably upgrade anyway!\n",
    "    Hessian_row_sum = obj_fun.multiply_with_Hessian(initial_image, initial_image.allocate(1))\n",
    "    return (-1*Hessian_row_sum).power(.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_prior_to_obj_funs(obj_funs, prior):\n",
    "    '''\n",
    "    Add prior evenly to every objective function in the obj_funs list.\n",
    "    \n",
    "    WARNING: modifies prior strength with 1/num_subsets (as currently needed for BSREM implementations)\n",
    "    WARNING: modifies elements of obj_funs\n",
    "    '''   \n",
    "    # evenly distribute prior over subsets\n",
    "    prior.set_penalisation_factor(prior.get_penalisation_factor() / len(obj_funs));\n",
    "    prior.set_up(initial_image)\n",
    "    for f in obj_funs:\n",
    "        f.set_prior(prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Try it with the thorax_single_slice data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "### First simulate some data\n",
    "(see the SIRF-exercises for more info, e.g. https://github.com/SyneRBI/SIRF-Exercises/blob/master/notebooks/PET/display_and_projection.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the directory with input files for this notebook\n",
    "data_path = os.path.join(examples_data_path('PET'), 'thorax_single_slice')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Read in images\n",
    "image = STIR.ImageData(os.path.join(data_path, 'emission.hv'))*0.05\n",
    "attn_image = STIR.ImageData(os.path.join(data_path, 'attenuation.hv'))\n",
    "template = STIR.AcquisitionData(os.path.join(data_path, 'template_sinogram.hs'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create attenuation\n",
    "acq_model_for_attn = STIR.AcquisitionModelUsingRayTracingMatrix()\n",
    "asm_attn = STIR.AcquisitionSensitivityModel(attn_image, acq_model_for_attn)\n",
    "asm_attn.set_up(template)\n",
    "attn_factors = asm_attn.forward(template.get_uniform_copy(1))\n",
    "asm_attn = STIR.AcquisitionSensitivityModel(attn_factors)\n",
    "# set it up\n",
    "asm_attn.set_up(template)\n",
    "# use \"standard\" terminology for multiplicative factor (Note: there is no norm here)\n",
    "mult_factors = attn_factors\n",
    "# fake background (randoms+scatter)\n",
    "background = template.get_uniform_copy(1)\n",
    "# find additive_term for the acq_model\n",
    "additive_term = background.clone()\n",
    "asm_attn.normalise(additive_term)\n",
    "\n",
    "# create acquisition model\n",
    "acq_model = STIR.AcquisitionModelUsingRayTracingMatrix()\n",
    "# we will increase the number of rays used for every Line-of-Response (LOR) as an example\n",
    "# (it is not required for this demo of course)\n",
    "acq_model.set_num_tangential_LORs(5)\n",
    "acq_model.set_acquisition_sensitivity(asm_attn)\n",
    "acq_model.set_additive_term(additive_term)\n",
    "# set-up\n",
    "acq_model.set_up(template,image)\n",
    "\n",
    "#%% simulate some data using forward projection\n",
    "acquired_data=acq_model.forward(image)\n",
    "\n",
    "print(acquired_data.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "### run initial OSEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% save max and central slice for future displays\n",
    "cmax = image.max()*.6\n",
    "im_slice = image.dimensions()[0] // 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_image = image.get_uniform_copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "acq_model, obj_fun = create_acq_model_and_obj_fun(acquired_data, additive_term, mult_factors, initial_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_image = scale_initial_image(acquired_data, additive_term, mult_factors, initial_image, obj_fun)\n",
    "OSEM_image = OSEM(obj_fun, initial_image, num_updates=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot_2d_image([1,1,1], OSEM_image.as_array()[im_slice,:,:], 'OSEM',[0,cmax])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29",
   "metadata": {},
   "source": [
    "### Construct the prior\n",
    "\n",
    "We will use a smoothed RDP with voxel-dependent weights (see function definitions above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa = compute_kappa_image(obj_fun, OSEM_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot_2d_image([1,1,1], kappa.as_array()[im_slice,:,:], 'kappa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = construct_RDP(1/36, OSEM_image, kappa)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33",
   "metadata": {},
   "source": [
    "### partition data (i.e. construct subsets) and adjust prior accordingly\n",
    "\n",
    "Use `partitioner.data_partition` to get a list of subset data, as well as corresponding acquisition models and objective functions. We'll add the prior evenly to each objective function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_subsets = 7\n",
    "data,acq_models, obj_funs = partitioner.data_partition(acquired_data,additive_term,mult_factors, num_subsets)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35",
   "metadata": {},
   "source": [
    "add_prior_to_obj_funs(obj_funs, prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36",
   "metadata": {},
   "source": [
    "Remember that this modified the penalty strength of the `prior` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "prior.get_penalisation_factor()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38",
   "metadata": {},
   "source": [
    "### compare 2 BSREM implementations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39",
   "metadata": {},
   "source": [
    "`BSREM1` etc are derived from `cil.Algorithm`, Check its help message for more information.\n",
    "\n",
    "Some notes on terminology:\n",
    "- CIL uses `loss` for the objective function value. This is somewhat confusing here, as we are maximising the objective function...\n",
    "- `cil.Algorithm` uses `iteration` for every call to `update`. This might not be what you expect for a subset algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40",
   "metadata": {},
   "source": [
    "Note: intentionally setting update_objective_intervals to be not a multiple of num_subsets such that we can see the oscillations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_iteration = 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42",
   "metadata": {},
   "outputs": [],
   "source": [
    "bsrem1 = BSREM1(data, obj_funs, initial=OSEM_image, initial_step_size=1, relaxation_eta=.05, update_objective_interval=5)\n",
    "bsrem1.run(iterations=max_iteration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "bsrem2=BSREM2(data, acq_models, prior, initial=OSEM_image, initial_step_size=1, relaxation_eta=.05, update_objective_interval=5)\n",
    "bsrem2.run(iterations=max_iteration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "tmp1=bsrem1.x\n",
    "tmp2=bsrem2.x\n",
    "plot_2d_image([1,2,1], tmp1.as_array()[im_slice,:,:], 'BSREM1 result',[0, cmax])\n",
    "plot_2d_image([1,2,2], (tmp1-tmp2).as_array()[im_slice,:,:], 'diff',[-cmax/100, cmax/100])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45",
   "metadata": {},
   "source": [
    "Plot objective function for both implementations. If all well, you should see only 1 curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(bsrem1.iterations, bsrem1.loss)\n",
    "plt.plot(bsrem2.iterations, bsrem2.loss);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close('all')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48",
   "metadata": {},
   "source": [
    "## Now use the NEMA IQ data acquired on the mMR"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49",
   "metadata": {},
   "source": [
    "The following names assume the data for the acquisition model has been written already by a preprocessing script. You might need to adjust this to your location and set-up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/home/kris/data/mMR/NEMA_IQ/tmpl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51",
   "metadata": {},
   "outputs": [],
   "source": [
    "acquired_data = STIR.AcquisitionData('prompts.hs')\n",
    "additive_term = STIR.AcquisitionData('additive_term.hs')\n",
    "mult_factors = STIR.AcquisitionData('mult_factors.hs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's use the attenuation image to get descent voxel-sizes etc\n",
    "# you might get a lot of warnings on \"unrecognized keyword\". Just ignore these...\n",
    "initial_image = STIR.ImageData('20170809_NEMA_MUMAP_UCL.hv')\n",
    "# crop it a bit to avoid wasting time\n",
    "initial_image=initial_image.zoom_image(zooms=(1,1,1), offsets_in_mm=(0,0,0), size=(-1,200,200))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run OSEM (see previous example) and construct kappa etc\n",
    "acq_model, obj_fun = create_acq_model_and_obj_fun(acquired_data, additive_term, mult_factors, initial_image)\n",
    "\n",
    "initial_image = scale_initial_image(acquired_data, additive_term, mult_factors, initial_image, obj_fun)\n",
    "OSEM_image = OSEM(obj_fun, initial_image, num_updates=14)\n",
    "\n",
    "kappa = compute_kappa_image(obj_fun, OSEM_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set to reasonable values for display\n",
    "cmax = OSEM_image.max()*1.1\n",
    "im_slice = 70 # for this acquisition, this might give a slcie through the spheres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot_2d_image([1,1,1], OSEM_image.as_array()[im_slice,:,:], 'OSEM',[0, OSEM_image.max()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot_2d_image([1,1,1], kappa.as_array()[im_slice,:,:], 'kappa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57",
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = construct_RDP(1/700, OSEM_image, kappa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot_2d_image([1,2,1], OSEM_image.as_array()[:,:,OSEM_image.dimensions()[1]//2 - 5], 'OSEM sagittal',[0,cmax])\n",
    "plot_2d_image([1,2,2], OSEM_image.as_array()[im_slice,:,:], 'OSEM transverse',[0,cmax])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_subsets = 7\n",
    "data,acq_models, obj_funs = partitioner.data_partition(acquired_data, additive_term, mult_factors, num_subsets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60",
   "metadata": {},
   "outputs": [],
   "source": [
    "add_prior_to_obj_funs(obj_funs, prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61",
   "metadata": {},
   "outputs": [],
   "source": [
    "bsrem1 = BSREM1(data, obj_funs, initial=OSEM_image, initial_step_size=.3, relaxation_eta=.05, update_objective_interval=10)\n",
    "bsrem1.run(iterations=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bsrem1.x.max())\n",
    "plt.figure()\n",
    "plot_2d_image([1,2,2], bsrem1.x.as_array()[im_slice,:,:], 'BSREM',[0, bsrem1.x.max()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot_2d_image([1,2,1], OSEM_image.as_array()[im_slice,:,:], 'initial OSEM',[0, cmax])\n",
    "plot_2d_image([1,2,2], bsrem1.x.as_array()[im_slice,:,:], 'BSREM',[0, cmax])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(bsrem1.iterations, bsrem1.loss);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close('all')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}