{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9fba7144",
   "metadata": {
    "papermill": {
     "duration": 0.010225,
     "end_time": "2024-11-19T22:42:34.546563",
     "exception": false,
     "start_time": "2024-11-19T22:42:34.536338",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Implementing New System Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "664a5b9e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:34.566821Z",
     "iopub.status.busy": "2024-11-19T22:42:34.566309Z",
     "iopub.status.idle": "2024-11-19T22:42:38.195912Z",
     "shell.execute_reply": "2024-11-19T22:42:38.195173Z"
    },
    "papermill": {
     "duration": 3.642095,
     "end_time": "2024-11-19T22:42:38.198020",
     "exception": false,
     "start_time": "2024-11-19T22:42:34.555925",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pytomography\n",
    "from pytomography.algorithms import OSEM, MLEM\n",
    "from pytomography.metadata import ObjectMeta, ProjMeta\n",
    "from pytomography.projectors import SystemMatrix\n",
    "from pytomography.likelihoods import NegativeMSELikelihood\n",
    "import matplotlib.pyplot as plt\n",
    "import torch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "608c20bb",
   "metadata": {
    "papermill": {
     "duration": 0.01003,
     "end_time": "2024-11-19T22:42:38.218561",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.208531",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "For now, we'll set `pytomography.device='cpu'` since all the objects we'll be creating/testing on will be on CPU, but this can be changed (provided all created objects are placed on GPU)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a5e573b5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.239785Z",
     "iopub.status.busy": "2024-11-19T22:42:38.239406Z",
     "iopub.status.idle": "2024-11-19T22:42:38.242308Z",
     "shell.execute_reply": "2024-11-19T22:42:38.241798Z"
    },
    "papermill": {
     "duration": 0.015537,
     "end_time": "2024-11-19T22:42:38.243908",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.228371",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pytomography.device = 'cpu'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8b3a6ea",
   "metadata": {
    "papermill": {
     "duration": 0.010315,
     "end_time": "2024-11-19T22:42:38.263998",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.253683",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "This tutorial demonstrates how to implement new system matrices in PyTomography that interface with all the available likelihoods/reconstruction algorithms. We will consider a very simple scenario of an imaging system where \n",
    "\n",
    "* data is acquired at two angles (0 and 90 degrees)\n",
    "* detector consists of an MxM grid that perfectly aligns with an object of shape $M \\times M \\times M$\n",
    "* each detector pixel has a different sensitivity $s_j$\n",
    "\n",
    "The value of each pixel is thus\n",
    "\n",
    "$$ p_j = s_j \\sum_{i \\in \\text{line}_j} x_i$$\n",
    "\n",
    "<img src=\"../images/example_scanner_sm.png\"  width=\"500\">\n",
    "\n",
    "\n",
    "\n",
    " We'll call this system the **Example Scanner**, (or *EXS* for short) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9d74de56",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.284976Z",
     "iopub.status.busy": "2024-11-19T22:42:38.284616Z",
     "iopub.status.idle": "2024-11-19T22:42:38.288191Z",
     "shell.execute_reply": "2024-11-19T22:42:38.287508Z"
    },
    "papermill": {
     "duration": 0.020966,
     "end_time": "2024-11-19T22:42:38.294306",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.273340",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "M = 3\n",
    "object_meta = ObjectMeta(dr=(1.5,1.5,1.5), shape=(M,M,M))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "622a1606",
   "metadata": {
    "papermill": {
     "duration": 0.009184,
     "end_time": "2024-11-19T22:42:38.314053",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.304869",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The objects we consider thus have dimensions of $3 \\times 3 \\times 3$ with voxel dimensions of $1.5 \\times 1.5 \\times 1.5$. The units of the voxel dimensions can be whatever (mm/cm): their units are defined by the system matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fd14578",
   "metadata": {
    "papermill": {
     "duration": 0.009194,
     "end_time": "2024-11-19T22:42:38.333767",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.324573",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Example 1: Non-ListMode"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4180e797",
   "metadata": {
    "papermill": {
     "duration": 0.009284,
     "end_time": "2024-11-19T22:42:38.353236",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.343952",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4364c1e4",
   "metadata": {
    "papermill": {
     "duration": 0.009267,
     "end_time": "2024-11-19T22:42:38.372683",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.363416",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We'll start by creating a projection metadata class which will interface with the system matrix\n",
    "* It should be a subclass of the `ProjMeta` class of PyTomography"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1a5f16fb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.394077Z",
     "iopub.status.busy": "2024-11-19T22:42:38.393595Z",
     "iopub.status.idle": "2024-11-19T22:42:38.397833Z",
     "shell.execute_reply": "2024-11-19T22:42:38.397193Z"
    },
    "papermill": {
     "duration": 0.016086,
     "end_time": "2024-11-19T22:42:38.399271",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.383185",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSProjMeta(ProjMeta):\n",
    "    def __init__(self, M, sensitivity_factor):\n",
    "        self.M = M\n",
    "        self.sensitivity_factor = sensitivity_factor\n",
    "        if (sensitivity_factor.shape[0]!=M)*(sensitivity_factor.shape[1]!=M):\n",
    "            raise ValueError(\"sensitivity_factor should have side dimensions M\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "66d71c60",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.422734Z",
     "iopub.status.busy": "2024-11-19T22:42:38.421980Z",
     "iopub.status.idle": "2024-11-19T22:42:38.428353Z",
     "shell.execute_reply": "2024-11-19T22:42:38.426653Z"
    },
    "papermill": {
     "duration": 0.019136,
     "end_time": "2024-11-19T22:42:38.429715",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.410579",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "M = object_meta.shape[0]\n",
    "# Note: the sensitivty factor is the same for each projection angle, a single detector is \"rotating\" between angle 0 and 90\n",
    "sensitivity_factor = torch.ones((M,M))+0.3*torch.rand((M,M))\n",
    "proj_meta = EXSProjMeta(M, sensitivity_factor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c35a1fb7",
   "metadata": {
    "papermill": {
     "duration": 0.010514,
     "end_time": "2024-11-19T22:42:38.450121",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.439607",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## System Matrix\n",
    "\n",
    "### Part 1: Understanding the Forward/Backward PRojections"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af64ab91",
   "metadata": {
    "papermill": {
     "duration": 0.009213,
     "end_time": "2024-11-19T22:42:38.468845",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.459632",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Before we begin to build the system matrix, lets understand the operations required to implement this system matrix. \n",
    "\n",
    "* Forward projection requires summing the object along its $x$-axis and $y$-axis respectively, and then concatenating together.\n",
    "* In PyTomography, our object is 3D, so it's shape is `[Lx,Ly,Lz]`; our projections have shape `[2,M,M]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "590d28dc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.488454Z",
     "iopub.status.busy": "2024-11-19T22:42:38.488065Z",
     "iopub.status.idle": "2024-11-19T22:42:38.496202Z",
     "shell.execute_reply": "2024-11-19T22:42:38.495585Z"
    },
    "papermill": {
     "duration": 0.019507,
     "end_time": "2024-11-19T22:42:38.497434",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.477927",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([2, 3, 3])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_object = torch.rand(object_meta.shape) \n",
    "# Sum object along x to get projection at 0 degrees\n",
    "sample_projection_0degrees = sample_object.sum(dim=0) * proj_meta.sensitivity_factor\n",
    "# Sum object along y to get projection at 90 degrees\n",
    "sample_projection_90degrees = sample_object.sum(dim=1) * proj_meta.sensitivity_factor\n",
    "# Concatenate to get the full set of projections\n",
    "sample_projections = torch.stack([sample_projection_0degrees, sample_projection_90degrees], dim=0)\n",
    "sample_projections.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f378ca",
   "metadata": {
    "papermill": {
     "duration": 0.009408,
     "end_time": "2024-11-19T22:42:38.516301",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.506893",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Let's plot the projections at each angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bb35936b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.536484Z",
     "iopub.status.busy": "2024-11-19T22:42:38.535931Z",
     "iopub.status.idle": "2024-11-19T22:42:38.736113Z",
     "shell.execute_reply": "2024-11-19T22:42:38.735462Z"
    },
    "papermill": {
     "duration": 0.212114,
     "end_time": "2024-11-19T22:42:38.737596",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.525482",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAACkCAYAAADWvuarAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAZmElEQVR4nO3da1QV59028Gs4bTSIxwCbBhDNC7FqjMEsJRYjjy0GjI9UmzT9kBifJg2NxiiLZdR80NYY22pTa6siDbzEWhujWB4PtJUaQVu1DRHTrHhIfINADTtIUjUacZ/+7wfrsIfNaSuwufdcv7VmLefeM3vuYV3r7z2nPZqICIiIqMcE+bsDRESBjoWWiKiHsdASEfUwFloioh7GQktE1MNYaImIehgLLRFRD2OhJSLqYSy0REQ9rE8X2g0bNkDTNIwZM6ZXt6tpGlauXNmt33nixAl885vfREREBAYNGoTZs2fjk08+6dZtUM8KlDyKCDZs2ID77rsPFosFVqsVP/zhD/Hvf/+7zeV/9atf6csmJibiRz/6ERwOR7f1xwz6dKEtKioCAHz44Yf4+9//7ufe3L4zZ85g6tSpsNvtePvtt1FUVISPPvoIaWlpuHjxor+7R10UKHnMy8vD4sWLMWvWLOzbtw9Lly7F9u3b8a1vfcurgK5evRovvfQSZs+ejT//+c944YUX8Nprr2H+/Pl+6r2ipI969913BYDMmDFDAMhzzz3Xa9sGICtWrOi273v88cdl2LBhcvnyZb3t/PnzEhoaKkuWLOm27VDPCZQ8/utf/5Lg4GB58cUXDe3bt28XAFJQUKC3NTU1SXh4uPzgBz8wLLt69WrRNE0+/PDDbumTGfTZEW1hYSEA4Cc/+QkefvhhvPXWW/jqq68My5w/fx6apmHdunV4/fXXkZiYiIiICKSmpuL48eNe3/mb3/wGSUlJsFgs+PrXv47t27fjmWeewfDhwzvtj81mw/PPP4977rkHYWFh+iGU0+nscD2n04l9+/Zhzpw5iIyM1NsTEhKQnp6OP/zhD134a5C/BUoejx8/DpfLhaysLEP7Y489BgAoKSnR2/70pz+hubkZ8+bNMyw7b948iAhKS0s77Sf9h78rfVu++uorGThwoDz00EMiIvLGG28IACkuLjYsV1NTIwBk+PDh8uijj0ppaamUlpbK2LFjZfDgwXLp0iV92S1btggAmTNnjuzbt09+97vfSVJSkiQkJEhCQoLhe9FqBNHQ0CBxcXGSkJAgW7Zskb/85S+yatUqsVgs8swzz3S4L2fOnBEAsnHjRq/P8vLyRNM0uX79uo9/IepNgZTHWyPXd955x9B+/fp10TRNrFar3rZ06VIBIFevXvX6nmHDhsn3vve9DrdFLfpkod26dasAkPz8fBER+fLLLyUiIkLS0tIMy90K9tixY8XpdOrt//jHPwSA/P73vxcREZfLJTExMTJx4kTD+rW1tRIaGtppsJ9//nmJiIiQ2tpaw3Lr1q0TAB0eQv3tb38z9MXTa6+9JgDk008/bf+PQX4XSHk8efKkAJBVq1YZ2g8ePCgAJCwsTG977rnnxGKxtPk9SUlJkpGR0e52yKhPnjooLCxEv3798OSTTwIAIiIi8Pjjj+PIkSP4+OOPvZafMWMGgoOD9fn7778fAFBbWwsAOHv2LGw2G5544gnDevHx8Zg8eXKn/dm3bx/S09MRGxsLp9OpT5mZmQCAysrKTr9D07Tb+oz8L5DyOG7cOEyZMgVr167Fzp07cenSJRw9ehQ5OTkIDg5GUJCxJDC33aPPFdpz587h8OHDmDFjBkQEly5dwqVLl/Cd73wHQMuVX09Dhw41zFssFgDA9evXAQCff/45ACA6Otpr3bbaWvvss8+wd+9ehIaGGqbRo0cDAJqamtpd91bfbvXB0xdffAFN0zBo0KBO+0D+EWh5BICdO3di8uTJeOKJJzB48GCkp6dj9uzZeOCBB/C1r33NsB/Nzc1e56KBm9kdMmRIp32lm0L83YHWioqKICLYtWsXdu3a5fX5m2++iVdffdUwYujMreB/9tlnXp/ZbLZO1x82bBjuv/9+rF69us3PY2Nj21135MiR6NevHz744AOvzz744APce++9CA8P77QP5B+BlkcAiIqKQllZGRobG2Gz2ZCQkIB+/fph06ZN+n8gADB27FgAN3M6ceJEQx+bmpp6/X5ilfWpQutyufDmm29i5MiReOONN7w+37dvH37+85/jj3/8o36VtCuSk5MRExODt99+G7m5uXp7XV0djh492mkwH3vsMZSVlWHkyJEYPHhw13cIQEhICGbOnIndu3fjZz/7GQYMGKBv+9ChQ1i8eLFP30e9JxDz6CkqKgpRUVEAbj6Mce3aNSxYsED//NFHH0V4eDiKi4sNhba4uBiapiE7O/u2t206/j1FbLR3714BID/96U/b/PzixYtisVgkOztbRFouPqxdu9ZrWbS6gOB5lXf//v36Vd74+HhJTEzscN1PP/1UEhIS5L777pNNmzbJwYMHZf/+/bJx40aZMWOG1NfXd7hfp0+floiICJkyZYqUlZXJ7t27ZcyYMRIbGyuNjY1d/OtQbwvUPBYUFEhBQYEcPHhQSkpK5NlnnxVN02TNmjVey7766quiaZosX75cKioqZO3atWKxWHr1PuJA0KcKbXZ2toSFhXVYfJ588kkJCQkRm83mU7BFbgbs3nvvlbCwMElKSpKioiKZNWuWjB8/vtN1L168KAsXLpTExEQJDQ2VIUOGSEpKirzyyitt3v7SWlVVlUybNk369+8vkZGRkp2dLefOnet0PfKfQM3jli1bZNSoUdK/f3/97onS0tJ2l//lL38pSUlJEhYWJvHx8bJixQqx2+0dboOMNBHzvgX30qVLSEpKQnZ2NgoKCvzdHTI55jFw9alztD3JZrNh9erVSE9Px9ChQ1FbW4tf/OIX+PLLL/HSSy/5u3tkMsyjuZim0FosFpw/fx4vvPACvvjiC/Tv3x+TJk1Cfn6+flsMUW9hHs3F1KcOiIh6Q597YIHUtGbNGjz00EMYMGAAoqKikJ2djbNnz3a6XmVlJVJSUhAeHo4RI0YgPz+/F3pLZuTPjLLQUreorKzE/Pnzcfz4cZSXl8PpdCIjIwPXrl1rd52amhpkZWUhLS0N1dXVWL58ORYuXGj4BSmi7uLPjPLUAfWIixcvIioqCpWVlZgyZUqby7z88svYs2cPTp8+rbfl5OTg/fffx7Fjx3qrq2RSvZlR01wMo441NzfDbrd7tYuI14+HWCwW/fn99ly+fBkAOnwe/tixY8jIyDC0TZ8+HYWFhXA4HAgNDe1q98kEVM5olwut25bU1UV71ciD8zpfqJf999f/6e8utOmX43/fZntzczMSEyJga3R5fRYREYGrV68a2lasWNHhO6xEBLm5ufjGN77R4fPwNpvN60dUoqOj4XQ60dTUBKvV2sHeeOurGd37VX9/d8HLO5dH+bsLbQrUjHJES7Db7bA1uvBR1T2IHNBy2v7Kl24kTfgX6uvrDW+H6GyksGDBAvzzn//EX//610633XokcutMFn+CjzypnlEWWtL1ixD0i2g5Ze/4T6AiIyMNIe7Iiy++iD179uDw4cO45557Olw2JibG69eqGhsbERIS4vVTg0SAuhnlXQekc8INh8fkhLvL64oIFixYgN27d+Odd95BYmJip+ukpqaivLzc0HbgwAFMmDCB52epTapmlIWWdA5xe01dNX/+fGzbtg3bt2/HgAEDYLPZYLPZ9B+7BoBly5bh6aef1udzcnJQW1uL3NxcnD59GkVFRSgsLEReXl637hcFDlUzykJLumYRr6mrNm/ejMuXL2Pq1KmwWq36tGPHDn2ZhoYG1NXV6fOJiYkoKytDRUUFHnjgAaxatQobNmzAnDlzunW/KHComlGeoyWdExoc0AzzXdWV27GLi4u92h555BGcOHGiy9shc1M1oyy0pHOIBodohnmivkTVjLLQks4hQXB4vBjZwWcGqY9RNaMstKS7ISEI8QjxDR8uNBD1BlUzykJLOrsEI9QjxHZFDsvIPFTNKAst6RwIggPBHvNEfYuqGWWhJZ1DguEQjxArcv6LzEPVjLLQks4hIbAbQqzGYRmZh6oZZaElXbOEIsgjxM3C51mob1E1oyy0pHNISKvDMjVGC2QeqmaUhZZ0TgQZQuyEIifAyDRUzSgLLekcEowQw2hBjRCTeaiaURZa0jW7Q6G5Qzzm1TgsI/NQNaMstKRzSggcEuIx78fOELVB1Yyy0JLOIcEIVvCwjMxD1Yyy0JLOIUGtQqzGc+RkHqpmlIWWdM5WT904FQkxmYeqGWWhJV2zOxTibnkP0g01MkwmompGWWhJ52x1/kuV0QKZh6oZVeP5NeoVDneQ1+SLw4cPY+bMmYiNjYWmaSgtLe1w+YqKCmia5jWdOXPmDvaCApmqGeWIlnROCTY8R+7raOHatWsYN24c5s2b59PL686ePYvIyEh9/u677/Zpu2QeqmaUhZZ0TncQgtweIXa7fFo/MzMTmZmZPm83KioKgwYN8nk9Mh9VM8pTB6S74Q7xmgDgypUrhunGjRvdut3x48fDarVi2rRpOHToULd+NwUWVTPKQks6F4LglJbJ9Z94xMXFYeDAgfq0Zs2abtme1WpFQUEBSkpKsHv3biQnJ2PatGk4fPhwt3w/BR5VM8pTB6RzuoOhGQ7Lbv67vr7ecH7KYrF0y/aSk5ORnJysz6empqK+vh7r1q3DlClTumUbFFhUzShHtKTzHCncmgAgMjLSMHVXiNsyadIkfPzxxz32/aQ2VTPKES3p7K5giMvj8UaPf/eW6upqWK3WXt8uqUHVjLLQks4lGjSPV4O4fPz1+qtXr+LcuXP6fE1NDU6ePIkhQ4YgPj4ey5Ytw4ULF7B161YAwPr16zF8+HCMHj0adrsd27ZtQ0lJCUpKSrpnhyjgqJpRFlrSOd1BgMcN4E4fbwavqqpCenq6Pp+bmwsAmDt3LoqLi9HQ0IC6ujr9c7vdjry8PFy4cAH9+vXD6NGjsX//fmRlZd3hnlCgUjWjLLSku9MQT506FdLBz9YVFxcb5pcsWYIlS5b4tA0yN1UzykJLOpcEtTos47VS6ltUzSgLLekcrS40OP1woYGoI6pmlIWWdG53EFweh2JuHw/LiHqaqhntcqEtuhLdk/24bf9v2v/1dxe8LG5I8XcXbosLGuBxFdcFNV58d8tD1Y/7uwttenf8Tn93wcuaj4f7uwu3RdWMckRLOlerCw0uRUYLZB6qZpSFlnRutwbN4/XNbkVe5UzmoWpGWWhJ53QFQVweowWXGqMFMg9VM8pCS7qbo4UgwzxRX6JqRlloSecWDZrHhQa3j483EvU0VTPKQkst3BrEc4SgyGiBTETRjLLQks7tCgI8znm5FTn/ReahakZZaEkn7puT5zxRX6JqRlloSSdiPCwTRc5/kXmomlEWWtJJq/Nfosj5LzIPVTPKQkstxPh4IxQZLZCJKJpRFlrSiUuDuDTDPFFfompGWWiphVsz3i6jyGEZmYiiGWWhJZ3mvjl5zhP1JapmVI2b0Kh33BoteE4+OHz4MGbOnInY2FhomobS0tJO16msrERKSgrCw8MxYsQI5Ofn32bnyRQUzSgLLbVwad6TD65du4Zx48bh17/+dZeWr6mpQVZWFtLS0lBdXY3ly5dj4cKFfAsutU/RjPLUAenu9LAsMzMTmZmZXV4+Pz8f8fHxWL9+PQBg1KhRqKqqwrp16zBnzhzfNk6moGpGOaIlnQZAE4/pP+1XrlwxTDdu3OiW7R07dgwZGRmGtunTp6OqqgoOh6NbtkGBRdWMstBSi3bOf8XFxWHgwIH6tGbNmm7ZnM1mQ3S08RVJ0dHRcDqdaGpq6pZtUIBRNKM8dUC69g7L6uvrERkZqbdbLJbu26ZmPMcmIm22EwHqZpSFlnSa6+bkOQ8AkZGRhhB3l5iYGNhsNkNbY2MjQkJCMHTo0G7fHqlP1Yzy1AG1uMNbZ3yVmpqK8vJyQ9uBAwcwYcIEhIaG9ui2SVGKZpSFlnS3Dss8J19cvXoVJ0+exMmTJwHcvDXm5MmTqKurAwAsW7YMTz/9tL58Tk4OamtrkZubi9OnT6OoqAiFhYXIy8vrrl2iAKNqRnnqgFq0Dq6PIa6qqkJ6ero+n5ubCwCYO3cuiouL0dDQoAcaABITE1FWVobFixdj48aNiI2NxYYNG3hrF7VP0Yyy0JLuTu9RnDp1qn6hoC3FxcVebY888ghOnDjh24bItFTNKAst6VR9jpzMQ9WMstBSC4HxUKz9//iJ/EPRjLLQkk7V0QKZh6oZZaElnaohJvNQNaMstKRr72Zwor5C1Yyy0JJO1dECmYeqGWWhJZ2qISbzUDWjLLTUwg3jFV1FQkwmomhGWWhJd+s3Pj3nifoSVTPKQks6zd3qQoMiowUyD1UzykJLOlXPf5F5qJpRFlrSqRpiMg9VM9rlQvu/nz3Qg924ff8T+Wd/d8HLL6zv+bsLt0XVEN/y7vid/u6CMg7ev93fXWhHx6+gUTWjHNGSTtUQk3momlEWWtJpLoHmEsM8UV+iakZZaEmn6miBzEPVjLLQkk7VEJN5qJpRvjOMdHf6PiYA2LRpExITExEeHo6UlBQcOXKk3WUrKiqgaZrXdObMmTvYCwpkqmaUI1rSaW5BkMc5L7fbt/NfO3bswKJFi7Bp0yZMnjwZW7ZsQWZmJk6dOoX4+Ph21zt79qzhVdF33323750nU1A1oxzRku7WT9B5Tr54/fXX8f3vfx/PPvssRo0ahfXr1yMuLg6bN2/ucL2oqCjExMToU3Bw8B3sBQUyVTPKQku6Ozkss9vteO+995CRkWFoz8jIwNGjRztcd/z48bBarZg2bRoOHTp0O10nk1A1ozx1QDrNLdA8DsVu/fvKlSuG5SwWCywWi6GtqakJLpcL0dHRhvbo6GjYbLY2t2e1WlFQUICUlBTcuHEDv/3tbzFt2jRUVFRgypQp3bFLFGBUzSgLLenau6IbFxdnWG7FihVYuXJl29+haYZ5EfFquyU5ORnJycn6fGpqKurr67Fu3ToWWmqTqhlloSWd5hJoQd43g9fX1xsuBLQeKQDAsGHDEBwc7DUyaGxs9BpBdGTSpEnYtm2br10nk1A1ozxHS7pbT914TgAQGRlpmNoKcVhYGFJSUlBeXm5oLy8vx8MPP9zlPlRXV8Nqtd7ZjlDAUjWjHNGSTpNWh2U+Pt2Ym5uLp556ChMmTEBqaioKCgpQV1eHnJwcAMCyZctw4cIFbN26FQCwfv16DB8+HKNHj4bdbse2bdtQUlKCkpKS7tolCjCqZpSFlnTtXWjoqu9+97v4/PPP8eMf/xgNDQ0YM2YMysrKkJCQAABoaGhAXV2dvrzdbkdeXh4uXLiAfv36YfTo0di/fz+ysrK6Z4co4KiaUU1EutTTmUde9OmLe8v//p++9zOJfVVQzEdttl+5cgUDBw5E2pQVCAkJ19udzmYcOfwjXL582XD+q69y25L83QVlXJcb/u5Cm+6y1rbZrnpGOaIl3Z2OFoh6mqoZZaElnaohJvNQNaMstKTTXAJN8751hqivUDWjLLTUwi03J895or5E0Yyy0JJOc7mhwW2YJ+pLVM0oCy3pNLdA87hJUZXzX2QeqmaUhZZauASAtJon6kMUzSgLLek0t7vVaEGNwzIyD1UzykJLLdytfhpJkRCTiSiaURZa0mlONzRR70IDmYeqGWWhpRYuN+BxRReKhJhMRNGMstBSC3EbD8VEjRCTiSiaURZaauFyAeLxtju3j2++I+ppimaUhZZaOF1AkHohJhNRNKMstNTCLTCc/1LkZnAyEUUzykJLLRQ9LCMTUTSjLLTUwuU2XlxQ5B5FMhFFM8pCSzpxuyAeowXPfxP1BapmlG/BpRZOp/fko02bNiExMRHh4eFISUnBkSNHOly+srISKSkpCA8Px4gRI5Cfn3+7vSczUDSjLLSkE5fLa/LFjh07sGjRIrzyyiuorq5GWloaMjMzDS+781RTU4OsrCykpaWhuroay5cvx8KFC/kWXGqXqhnlyxlNpLOXM/5X2OMI0UL1dqc48I59Z5dffDdx4kQ8+OCD2Lx5s942atQoZGdnY82aNV7Lv/zyy9izZw9Onz6tt+Xk5OD999/HsWPHfNk1AHw5oy9UfTmjqhnliJZ0dzJasNvteO+995CRkWFoz8jIwNGjR9tc59ixY17LT58+HVVVVXA4HL7vAAU8VTPKi2Gkc7iaIWgJrhM3g3TlyhXDchaLBRaLxdDW1NQEl8uF6OhoQ3t0dDRsNlub27PZbG0u73Q60dTUBKvVetv7QoFJ1Yx2udDuTftVVxclxYSFhSEmJgZ/tZV5fRYREYG4uDhD24oVK7By5co2v0vTNMO8iHi1dbZ8W+1d0d6pEfJ2l7874CPVM8oRLSE8PBw1NTWw2+1en7UVwtYjBQAYNmwYgoODvUYGjY2NXiOCW2JiYtpcPiQkBEOHDvV1NyiAqZ5RFloCcDPI4eHht71+WFgYUlJSUF5ejm9/+9t6e3l5OWbNmtXmOqmpqdi7d6+h7cCBA5gwYQJCQ0PbXIfMS+mMClE3eeuttyQ0NFQKCwvl1KlTsmjRIrnrrrvk/PnzIiKydOlSeeqpp/TlP/nkE+nfv78sXrxYTp06JYWFhRIaGiq7du3y1y5QgPNXRlloqVtt3LhREhISJCwsTB588EGprKzUP5s7d6488sgjhuUrKipk/PjxEhYWJsOHD5fNmzf3co/JbPyR0S7fR0tERLeH99ESEfUwFloioh7GQktE1MNYaImIehgLLRFRD2OhJSLqYSy0REQ9jIWWiKiHsdASEfUwFloioh7GQktE1MNYaImIetj/B8soic3TrBJ/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 400x150 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4,1.5))\n",
    "plt.subplot(121)\n",
    "plt.pcolormesh(sample_projections[0].cpu().T, vmin=0, vmax=2)\n",
    "plt.colorbar()\n",
    "plt.axis('off')\n",
    "plt.title('Angle 0')\n",
    "plt.subplot(122)\n",
    "plt.pcolormesh(sample_projections[1].cpu().T, vmin=0, vmax=2)\n",
    "plt.colorbar()\n",
    "plt.axis('off')\n",
    "plt.title('Angle 90')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "298ad3af",
   "metadata": {
    "papermill": {
     "duration": 0.009394,
     "end_time": "2024-11-19T22:42:38.757228",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.747834",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "* Back projection: The transpose of summing along X and Y is duplication along X and Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "539cdaa5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.781798Z",
     "iopub.status.busy": "2024-11-19T22:42:38.781305Z",
     "iopub.status.idle": "2024-11-19T22:42:38.786092Z",
     "shell.execute_reply": "2024-11-19T22:42:38.785534Z"
    },
    "papermill": {
     "duration": 0.016906,
     "end_time": "2024-11-19T22:42:38.787353",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.770447",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# First adjust projections by sensitivity factor\n",
    "sample_projections_angle_0_sensitivity_adjusted = sample_projections[0]*proj_meta.sensitivity_factor\n",
    "# Back project at angle 0 by duplication\n",
    "sample_object_BP_angle0 = sample_projections_angle_0_sensitivity_adjusted.unsqueeze(0).repeat(object_meta.shape[0],1,1)\n",
    "# Back project at angle 90 by duplication\n",
    "sample_projections_angle_90_sensitivity_adjusted = sample_projections[1]*proj_meta.sensitivity_factor\n",
    "sample_object_BP_angle90 = sample_projections_angle_90_sensitivity_adjusted.unsqueeze(1).repeat(1,object_meta.shape[0],1)\n",
    "# Back projected object is sum of each\n",
    "sample_object_BP = sample_object_BP_angle0 + sample_object_BP_angle90"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "189754b7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.811488Z",
     "iopub.status.busy": "2024-11-19T22:42:38.810939Z",
     "iopub.status.idle": "2024-11-19T22:42:38.854234Z",
     "shell.execute_reply": "2024-11-19T22:42:38.853705Z"
    },
    "papermill": {
     "duration": 0.057841,
     "end_time": "2024-11-19T22:42:38.855619",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.797778",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAK8AAADCCAYAAAA/+JLxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAK20lEQVR4nO3bf0hV9x/H8ddx93qn1zLS648ZGbWWOhYJRRGVZVFZsVpLZ4ymRa2ffwTZd8hoTZP6BsUoSu4NIgvaylGNwFFBOVvDDWN4o622LN0W4U3NPyoz7rX394/hXVev995Zdn3zfT1A6J77ued8jvfJuZ97JENEBEQKRYR7AkT9xXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+pFZZ4DcMI6ee7774b0HkUFhYGncOoUaOC7sftdsPhcGDSpEkYPnw4oqOjkZqaisWLF+PMmTPecU1NTTAMAxUVFd5tFRUVMAwDTU1NL/8EA+h5nlarFenp6SgpKcHjx499xvb8PVksFowbNw7bt29HZ2fnK53380zhOGhtba3P4x07dqC6uhqXLl3y2Z6RkTGg89i2bRvWrVvn97mKigo4HA689957QfezYsUKnD59Gps3b0ZJSQksFgvu3LmDc+fO4fz58wH3sXDhQtTW1iI5Obnf59Ffy5Ytw5YtWwAAjx49Qk1NDUpLS3Ht2jWcOnXKZ2xUVJT3/Wlvb8dXX32F0tJS3Lx5EydPnnzlcwcAyCBQUFAgVqs13NPwqq2tlcjISJkxY4a43e6AY+/cuSMA5LPPPvP7fFdXl/ffjY2NAkCOHDnyMqfbLwBk48aNvbavWLFCIiIi5MmTJ95tfb0/06dPFwBy9+7dAZ1rXwbtmvfBgwfYsGEDUlJSEBkZidGjR+PTTz/F06dPfcYZhoFNmzbB4XDgrbfegsViQUZGBk6cONGv4zY3N+P999+HzWZDZWUlTKbAH05tbW0A0OeVMyIi8K+4r2XDuXPnMHv2bMTGxiI6Ohrp6enYtWuXz5irV6/i3XffxfDhw/H6668jMzMTlZWVQc4wsNjYWBiGgddeey3o2ClTpgAA/vjjjxc6Zn+FZdkQTGdnJ2bNmoXbt2+jpKQE48ePx/fff49du3ahvr4eVVVVPuPPnj2L6upqlJaWwmq1ory8HMuXL4fJZMKyZctCPq7b7UZubi5aW1tRU1ODxMTEoK9JT0/HsGHDUFJSgoiICMydOzekdXIghw8fxpo1a5CVlQW73Y6EhAT8/vvvuH79undMdXU15s+fj8mTJ8NutyM2NhYnTpzABx98gI6ODhQWFgY9jojA4/EA+GfZcPToUeTn58NsNgd9fUNDAwDAZrP170RfVFiu9z30/Fiy2+0CQCorK33G7d69WwDIhQsXvNsASFRUlDQ3N3u3eTweSUtLkzfffPNfzWPDhg0CQOx2+796XVVVlcTHxwsAASBxcXGSm5srZ8+e9Rnnb9lw5MgRASCNjY0iIvLw4UMZOnSoTJs2TZ49e9bnMdPS0iQzM7PXsmbRokWSnJzss1zxp3uuPX9ycnLk0aNHPmO73x+32y1ut1taWlpk3759YhiGTJo0KYTf0MAYlPHm5eWJ1Wrt9ea5XC4BIJ988ol3GwBZtGhRr31u375dAMhff/0V0hy6I1q1apXf57u6urxvntvtFo/H4/N8R0eHnDlzRoqKimTGjBliNpt7rStDiff8+fMCQL788ss+53rr1i0BIHv27PGZk9vtlvLycgEgv/76a8DzBSB5eXlSV1cndXV1cvnyZdm/f7/YbDaZNm2adHZ2escWFBT0itwwDFmwYEHY1rsiIoNy2dDW1oakpCQYhuGzPSEhASaTybvO7JaUlNRrH93b2traMGLEiIDHu3r1KtavX4+JEyeivLzc75hVq1bh6NGj3sdZWVk+t/KioqKwZMkSLFmyBADw559/IicnBwcPHsT69evx9ttvB5xDt5aWFgAIOGeXywUAKCoqQlFRkd8xra2tQY9ls9kwceJE7+Pp06fDZrNh+fLlqKiowNq1a73PRUVF4fLlywAAi8WC1NRUDB06NPgJDaBBGW9cXBx++ukniIhPwPfv34fH40F8fLzP+Obm5l776N4WFxcX8FgtLS1YunQpYmJicOrUKVgsFr/jPv/8c2zatMn7eMiQIQH3O3LkSHz88cfYvHkzfvnll5Dj7V4/3r17t88x3edfXFyMpUuX+h0zbty4kI7X0/jx4wEATqfTZ3tERIRP6IPBoIx39uzZqKysxDfffONzj/TYsWPe55938eJFuFwu7xesrq4unDx5EmPGjAl4BfN4PMjNzcW9e/dw4cIFjBw5ss+xo0aN8vtF7OHDhzAMAzExMb2eu3HjBgDgjTfe6Ptke5g6dSpiY2Nht9uRn5/f69MH+DvMsWPHwul0YufOnSHvOxT19fUA/v6UG+wGZbwfffQRDh48iIKCAjQ1NeGdd97BlStXsHPnTixYsABz5szxGR8fH4/s7Gxs27bNe7fh5s2bQW+Xbd26FTU1Nfjwww8RHR2NH3/80e+47ltC/vz222+YN28e8vPzkZWVheTkZLS3t6OqqgqHDh3CzJkzMXXq1JDPPSYmBnv37sXq1asxZ84crFmzBomJiWhoaIDT6cSBAwcAAA6HAzk5OZg3bx4KCwuRkpKCBw8e4MaNG/j555/x9ddfBz2Wy+XynnNnZyfq6+tRVlaGYcOGYeXKlSHPOWzCttp+jr+b4G1tbbJu3TpJTk4Wk8kkqampUlxc7PNFQuSfm+3l5eUyZswYMZvNkpaWJsePHw963NTU1D6/dT//E0h7e7uUlZVJdna2pKSkSGRkpFitVpkwYYKUlZVJR0eHd2woX9i6ffvtt5KVlSVWq1Wio6MlIyNDdu/e7TPG6XRKXl6eJCQkiNlslqSkJMnOzg7pbknPczSbzTJ69GhZuXKlNDQ0+IwdbH9E6maI6P7fw4ZhYOPGjd4rEv3/GLR/YSMKhvGSWoPyC9u/oXzVQy+AV15Si/GSWoyX1GK8pFbIX9jG/veLgZzHK5dY1xXuKbw0lnZ3uKfwUl28VBzSOF55SS3GS2oxXlKL8ZJajJfUYrykFuMltRgvqcV4SS3GS2oxXlKL8ZJajJfUYrykFuMltRgvqcV4SS3GS2oxXlKL8ZJajJfUYrykFuMltRgvqcV4SS3GS2oxXlKL8ZJajJfUYrykFuMltRgvqcV4SS3GS2oxXlKL8ZJajJfUYrykFuMltRgvqcV4SS3GS2oxXlKL8ZJajJfUYrykFuMltRgvqcV4SS3GS2qZQh0Y73w2kPN45YZcuR3uKbw0npbWcE8hLHjlJbUYL6nFeEktxktqMV5Si/GSWoyX1GK8pBbjJbUYL6nFeEktxktqMV5Si/GSWoyX1GK8pBbjJbUYL6nFeEktxktqMV5Si/GSWoyX1GK8pBbjJbUYL6nFeEktxktqMV5Si/GSWoyX1GK8pBbjJbUYL6nFeEktxktqMV5Si/GSWoyX1GK8pBbjJbUYL6nFeEktxktqMV5Si/GSWoyX1GK8pBbjJbUYL6nFeEktU6gDf/jCMZDzeOX+48oM9xRemsbHtnBPISx45SW1GC+pxXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+pxXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+pxXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+pxXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+pxXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+pxXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+pZYiIhHsSRP3BKy+pxXhJLcZLajFeUovxklqMl9RivKQW4yW1GC+p9T+ty0XresM4RgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 200x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(2,2))\n",
    "plt.title('Top Z-Slice BP')\n",
    "plt.pcolormesh(sample_object_BP[:,:,-1].T)\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48019e0f",
   "metadata": {
    "papermill": {
     "duration": 0.010945,
     "end_time": "2024-11-19T22:42:38.879326",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.868381",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Part 2: Implementing The Forward/Backward Projections in the System Matrix Class"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3be6e522",
   "metadata": {
    "papermill": {
     "duration": 0.009537,
     "end_time": "2024-11-19T22:42:38.900053",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.890516",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We can now create a system matrix class that implements forward and back projection\n",
    "\n",
    "* We inherit from the `SystemMatrix` class of PyTomography\n",
    "* The method for forward projection should be called `forward`, while the method for back projection should be called `backward`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "65933c32",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.922184Z",
     "iopub.status.busy": "2024-11-19T22:42:38.921569Z",
     "iopub.status.idle": "2024-11-19T22:42:38.926977Z",
     "shell.execute_reply": "2024-11-19T22:42:38.926420Z"
    },
    "papermill": {
     "duration": 0.018037,
     "end_time": "2024-11-19T22:42:38.928262",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.910225",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSSystemMatrix(SystemMatrix):\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        projection_0degrees = object.sum(dim=0)\n",
    "        projection_90degrees = object.sum(dim=1)\n",
    "        projections = torch.stack([projection_0degrees, projection_90degrees], dim=0)\n",
    "        projections *= self.proj_meta.sensitivity_factor\n",
    "        return projections\n",
    "    def backward(self, projections, subset_idx = None):\n",
    "        object_BP_angle0 = (projections[0]*self.proj_meta.sensitivity_factor).unsqueeze(0).repeat(self.proj_meta.M,1,1)\n",
    "        object_BP_angle90 = (projections[1]*self.proj_meta.sensitivity_factor).unsqueeze(1).repeat(1,self.proj_meta.M,1)\n",
    "        object_BP = object_BP_angle0 + object_BP_angle90\n",
    "        return object_BP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2da28fc9",
   "metadata": {
    "papermill": {
     "duration": 0.010271,
     "end_time": "2024-11-19T22:42:38.948434",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.938163",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "This system matrix can now be used for forward/backward projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6902f897",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:38.970104Z",
     "iopub.status.busy": "2024-11-19T22:42:38.969751Z",
     "iopub.status.idle": "2024-11-19T22:42:38.973865Z",
     "shell.execute_reply": "2024-11-19T22:42:38.973296Z"
    },
    "papermill": {
     "duration": 0.016942,
     "end_time": "2024-11-19T22:42:38.975141",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.958199",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSSystemMatrix(object_meta=object_meta, proj_meta=proj_meta)\n",
    "FP = system_matrix.forward(sample_object)\n",
    "BP = system_matrix.forward(sample_object)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "003cc694",
   "metadata": {
    "papermill": {
     "duration": 0.00945,
     "end_time": "2024-11-19T22:42:38.994733",
     "exception": false,
     "start_time": "2024-11-19T22:42:38.985283",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In order to use the system matrix in reconstruction algorithms, one of the requirements is the computation of a normalization factor $H^T 1$. For this, we need to define the `compute_normalization_factor` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b60c3988",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.015722Z",
     "iopub.status.busy": "2024-11-19T22:42:39.015094Z",
     "iopub.status.idle": "2024-11-19T22:42:39.020598Z",
     "shell.execute_reply": "2024-11-19T22:42:39.020010Z"
    },
    "papermill": {
     "duration": 0.017123,
     "end_time": "2024-11-19T22:42:39.021884",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.004761",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSSystemMatrix(SystemMatrix):\n",
    "    def compute_normalization_factor(self):\n",
    "        # A clever implementation of this function will only compute the normalization factor once, and then store it for future use (e.g. using a boolean flag)\n",
    "        norm_projections = torch.ones((2,self.proj_meta.M, self.proj_meta.M))\n",
    "        return self.backward(norm_projections)\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        projection_0degrees = object.sum(dim=0)\n",
    "        projection_90degrees = object.sum(dim=1)\n",
    "        projections = torch.stack([projection_0degrees, projection_90degrees], dim=0)\n",
    "        projections *= self.proj_meta.sensitivity_factor\n",
    "        return projections\n",
    "    def backward(self, projections, subset_idx = None):\n",
    "        object_BP_angle0 = (projections[0]*self.proj_meta.sensitivity_factor).unsqueeze(0).repeat(self.proj_meta.M,1,1)\n",
    "        object_BP_angle90 = (projections[1]*self.proj_meta.sensitivity_factor).unsqueeze(1).repeat(1,self.proj_meta.M,1)\n",
    "        object_BP = object_BP_angle0 + object_BP_angle90\n",
    "        return object_BP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43a98f9c",
   "metadata": {
    "papermill": {
     "duration": 0.00932,
     "end_time": "2024-11-19T22:42:39.042131",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.032811",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Then we can compute the normalization factor as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0d12725e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.062539Z",
     "iopub.status.busy": "2024-11-19T22:42:39.061987Z",
     "iopub.status.idle": "2024-11-19T22:42:39.065919Z",
     "shell.execute_reply": "2024-11-19T22:42:39.065047Z"
    },
    "papermill": {
     "duration": 0.015581,
     "end_time": "2024-11-19T22:42:39.067143",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.051562",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSSystemMatrix(object_meta=object_meta, proj_meta=proj_meta)\n",
    "norm_factor = system_matrix.compute_normalization_factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "75a1f5ec",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.087365Z",
     "iopub.status.busy": "2024-11-19T22:42:39.086859Z",
     "iopub.status.idle": "2024-11-19T22:42:39.177656Z",
     "shell.execute_reply": "2024-11-19T22:42:39.177021Z"
    },
    "papermill": {
     "duration": 0.102685,
     "end_time": "2024-11-19T22:42:39.179269",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.076584",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANkAAADCCAYAAADEr2cuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAaBUlEQVR4nO3deVAUZ/oH8G8P4wgMSLw4REWWKMqR0qByKIYVF+ORxPWIGhe1cnjixmLZNfkZI7qUxmO3kuyKomtJrFXBCBIqugqKeERIVqNuYqnRqOsFosSI0YWB6ef3x8iEdnouM28w8HyqukrefqfnHcov79tv9/QrERGBMSaMprkbwFhLxyFjTDAOGWOCccgYE4xDxphgHDLGBOOQMSYYh4wxwThkjAnmcMgkSXJoKy0tFdhcYPr06Xbb0KNHD7vHqa+vR1ZWFgYMGIAOHTrA09MTQUFBeOmll7Bz505zvcuXL0OSJGRnZ5vLsrOzIUkSLl++7PoPaEPj53vvvfcs9jW26dixYz9rm5h9WkcrlpWVKX7+85//jAMHDqCkpERRHhYW5pqWWbFo0SLMmjVLdV92djaysrLw29/+1u5xkpOTkZ+fj/nz52PJkiVo27YtLl68iD179mDv3r02jzFq1CiUlZUhICDgsT/HT/Hee+9hxowZ6NChQ7O8P3MSPaZp06aRXq9/3Je7XFlZGel0OhoyZAjV19fbrHvx4kUCQO+++67qfqPRaP73pUuXCABt2rTJlc19LABo2LBhpNVqKTU1VbFv06ZNBID+/e9/u+S97t+/75LjMCKXnpN99913mDNnDgIDA6HT6fCrX/0KCxcuRF1dnaKeJElISUlBVlYWevXqhbZt2yIsLAw5OTmP9b6VlZUYN24cOnfujO3bt0Ortd1BV1dXA4DVnkijsf1rsTZc3LNnDxITE+Hj4wNPT0/06dMHy5cvV9Q5duwYXnzxRXTo0AHu7u7o168ftm/fbucT/ig0NBSvvfYa1qxZg//+97926xcWFiI2Nhaenp7w9vbGb37zG4tRSXp6OiRJwpdffonx48ejffv2CAkJAQD06NEDo0ePxqeffop+/frBw8MDffr0waeffmr+XfTp0wd6vR4DBw7k4aqax03noz3Z//73P3rmmWdIr9fT6tWrqaioiBYtWkRarZZGjhypeC0A6tatG4WFhdG2bduosLCQnn/+eQJAH3/8sVPtMBgMNHjwYNLpdFRWVubQa3744Qd66qmnyN/fn7KysujSpUtW66r1ZI29RtPX/eMf/yBJkighIYG2bt1K+/bto8zMTJozZ465TklJCel0OoqPj6fc3Fzas2cPTZ8+3eGeEgDNnTuXKioqyNPTk5KTky3a1LQn27JlCwGgpKQkKigooNzcXIqKiiKdTkeHDx8211u8eDEBoKCgIFqwYAEVFxdTQUEBEREFBQVR165dKSIigrZt20a7d++m6OhoatOmDb377rs0aNAgys/Pp507d1KvXr3Iz8+PHjx4YPeztCYuC9m6desIAG3fvl1Rb8WKFQSAioqKfnxTgDw8PKiystJc1tDQQL1796ann37aqXbMmTOHANC6deucet2uXbuoU6dOBIAAUMeOHWnChAlUWFioqOdIyO7du0ft2rWjwYMHkyzLVt+zd+/e1K9fP4vh7OjRoykgIEAxTFXTGDIiooULF5JGo6FTp04p2tQYMqPRSF26dKHIyEjFce/du0e+vr4UFxdnLmsMmdrwOSgoiDw8POjatWvmspMnTxIACggIUAwrCwoKCIDF77C1c9lwsaSkBHq9HuPHj1eUT58+HQCwf/9+RXliYiL8/PzMP7u5uWHixIm4cOECrl275tB7ZmdnIzMzE6+++ipmzpxpsV+WZTQ0NJg3o9Fo3jdy5EhcuXIFO3fuRFpaGsLDw1FQUIAXX3wRKSkpjn5sAMDRo0dRU1ODOXPmQJIk1ToXLlzA2bNnMWXKFABQtGvkyJGoqKjAuXPnHH7PP/3pT+jQoQMWLFiguv/cuXO4ceMGkpOTFcNfLy8vjBs3DuXl5Xjw4IHiNePGjVM9Vt++fREYGGj+uU+fPgCAhIQEeHp6WpQ7MoxtTVwWsurqavj7+1v8J/P19YVWqzWfBzXy9/e3OEZj2aN11Rw7dgyzZ89G//79kZmZqVrn1VdfRZs2bcxbYmKiYr+HhwfGjBmDVatW4eDBg7hw4QLCwsKwZs0anD592m4bGt26dQsA0LVrV6t1bt68CQBIS0tTtKlNmzaYM2cOAOD27dsOv2e7du3wzjvvYM+ePThw4IDFflvnnV26dIEsy7hz546i3No56qOzmDqdzmZ5bW2tg5+idXB4Ct+ejh074vPPPwcRKYJWVVWFhoYGdOrUSVG/srLS4hiNZR07drT5Xrdu3cLYsWPh5eWFvLw8tG3bVrVeenq6olfy9va2edzu3btjxowZmD9/Pk6fPo3w8HCb9Rt17twZAGz2wI2f/+2338bYsWNV64SGhjr0fo1mz56NDz74AAsWLMDs2bMV+xp/hxUVFRavu3HjBjQaDdq3b68ot9YLs5/GZSFLTEzE9u3bUVBQoLjGtHnzZvP+pvbv34+bN2+ah4xGoxG5ubkICQmx2SM0NDRgwoQJuHHjBoqKitC9e3erdXv06KF6YfrevXuQJAleXl4W+86cOQPA9NfeUXFxcfDx8cG6deswadIk1f+soaGh6NmzJ06dOoVly5Y5fGxbdDodMjIyMGXKFIs/YqGhoQgMDMTWrVuRlpZmbtP9+/eRl5dnnnFk4rksZFOnTsWaNWswbdo0XL58GZGRkThy5AiWLVuGkSNHYtiwYYr6nTp1wtChQ7Fo0SLo9XpkZmbi7Nmzdqfx//jHP+LgwYOYMmUKPD09UV5erlovJibG6jHOnTuH4cOHY9KkSXjuuecQEBCAO3fuYNeuXVi/fj0SEhIQFxfn8Gf38vLCX/7yF7z++usYNmwY3njjDfj5+eHChQs4deoU/v73vwMAsrKyMGLECAwfPhzTp09HYGAgvvvuO5w5cwZffvklPv74Y4ffs9HkyZOxevVq/Otf/1KUazQarFy5ElOmTMHo0aMxc+ZM1NXVYdWqVfj+++9V7xphgjzujInaxejq6mqaNWsWBQQEkFarpaCgIHr77beptrZWUQ8PZ8kyMzMpJCSE2rRpQ71796YtW7bYfd+goCDzjKCtzZY7d+5QRkYGDR06lAIDA0mn05Fer6e+fftSRkaGYgra0Sl8IqLdu3fTc889R3q9njw9PSksLIxWrFihqHPq1Cl6+eWXydfXl9q0aUP+/v40dOhQh2ZH0WR2samioiLz5370YnRBQQFFR0eTu7s76fV6SkxMpM8++0xRp3F28datWxbHDgoKolGjRjnUlsbf1apVq+x+ltZEIvr5n1YlSRLmzp1r/gvPWEvGd+EzJhiHjDHBXDbx4YxmGKEy1my4J2NMMA4ZY4JxyBgTrFnOyVjLV1tbC4PBYHW/TqeDu7v7z9ii5uOykE0un+GqQwn1+RfO3R/YHC68nNXcTbBL4/+N1X21tbUIDvJCZZXRah1/f39cunSpVQSNezLmcgaDAZVVRlw41g3tvC3PSGruyXi6/1UYDAYOGWM/hac3wdPb8nJNA1rXJRwOGRPGSASjyjVRtbKWjEPGhGmAjHor5a0Jh4wJU08y6lU6rXrikDHmEvLDTa28NeGQMWEMRDConH+plbVkHDImTAMk1MPyUQwNKmUtGYeMCSOTaVMrb0343kUmjAEaq5szli9fjgEDBsDb2xu+vr4YM2aM3WdUlpaWqq74c/bsWUW9vLw8hIWFmR8V33RFn0aZmZkIDg6Gu7s7oqKicPjwYafazyFjwtSTxurmjIMHD2Lu3LkoLy9HcXExGhoakJSUhPv379t97blz51BRUWHeevbsad5XVlaGiRMnIjk5GadOnUJycjJefvllfP755+Y6ubm5mD9/PhYuXIgTJ04gPj4eI0aMwJUrVxxuv8ue8cH3LrrOL/3exZqaGvj4+KDk627wUrmt6od7MoZGXMXdu3fRrl07p9/71q1b8PX1xcGDBzFkyBDVOqWlpfj1r3+NO3fu4KmnnlKtM3HiRNTU1Cie9PX888+jffv22LZtGwAgOjoazz77LNauXWuu06dPH4wZM8ZiMRFruCdjwjRY6cUaHvZkNTU1iu3R1X+suXv3LgDLJxir6devHwICApCYmGjxpOWysjIkJSUpyoYPH46jR48CMN2Defz4cYs6SUlJ5jqO4JAxYerJDfWkVdncAADdunWDj4+PeXOkZyAipKamYvDgwYiIiLBaLyAgAOvXr0deXh7y8/MRGhqKxMREHDp0yFynsrJSsR4DAPj5+ZmfZH379m0YjUabdRzBs4tMGCMkGFWm6xvLrl69qhguWnvcelMpKSn4z3/+gyNHjtisFxoaqnjseWxsLK5evYrVq1crhpiPPu2ZHnnMvKN1bOGejAlj6snUN8C0aEbTzV7I5s2bh8LCQhw4cMDmo9ytiYmJwfnz580/+/v7W/RIVVVV5p6rU6dOcHNzs1nHERwyJkw9aWFQ2erJuQEUESElJQX5+fkoKSlBcHDwY7XnxIkTipVrYmNjUVxcrKhTVFRkfkS7TqdDVFSURZ3i4mKnHuPOw0UmjAwNZJW/47KT3yebO3cutm7dik8++QTe3t7mnsXHxwceHh4ATKvlXL9+3bzAyfvvv48ePXogPDwcBoMB//znP5GXl4e8vDzzcd98800MGTIEK1aswEsvvYRPPvkE+/btUwxFU1NTkZycjP79+yM2Nhbr16/HlStXMGvWLIfbzyFjwhjIDdqHQ0NluXPHaZw+T0hIUJRv2rTJvMhkRUWF4tqVwWBAWloarl+/Dg8PD4SHh2PXrl0YOXKkuU5cXBxycnLwzjvvYNGiRQgJCUFubi6io6PNdSZOnIjq6mosXboUFRUViIiIwO7duxEUFORw+/k62ROopVwn2/BlFDy9LUP24J4Rbzx7/LGvk/3ScE/GhJEBGMlyFo6/6sKYi9STFlqVSQ61L3K2ZBwyJky9lXOyev4+GWOuYYQGRpXZRbWyloxDxoRpaHLhWVnOPRljLlFPbnDj4SKHjIkjkwayynfH1MpaMg4ZE6aeNFZ6stY1ic8hY8JYPyfjkDHmEkbSwKgyNFQra8k4ZEyYenKDhoeLHDImToOV2UUeLjLmIkaSVO9dVCtryThkTBij7IYG2bInM8rckzHmEvUkQVKZ5Kjnnowx1+CL0SYcMiZMPWms9GQcMsZcgnsyEw4ZE8aIH58W/Gh5a8IhY8I0yG6QVGYX1WYcWzIOGRNGJgmy2jM+eHaRMddosDLxoTaEbMk4ZEyYBlkDSVYJmUpZS8YhY8LwcNGEQ8aEMVq544PvXWTMRRpkDcDDxVZ2wYL9rBqHi2qbMx5nYfamPvvsM2i1WvTt21dRnpCQoLp4+6hRo8x10tPTLfb7+/s71X4OGROm8ZvRapszfsrC7Hfv3sXUqVORmJhosS8/P1+xaPvXX38NNzc3TJgwQVEvPDxcUe+rr75yqv08XGTCGK3MLhqdHC7u2bNH8fOmTZvg6+uL48ePW12YvdHMmTPxyiuvwM3NDQUFBYp9j645nZOTA09PT4uQabVap3uvprgnY8LYGy6KXph906ZN+Pbbb7F48WKHjrtx40ZMmjQJer1eUX7+/Hl06dIFwcHBmDRpEi5evOjQ8Rq5rCe79l5PVx1KqJDC8uZugl1PY2ZzN8Gui7+3X0eWNaq9lvywrFu3boryxYsXIz093eYxHV2Y/fz583jrrbdw+PBhaLX2/5t/8cUX+Prrr7Fx40ZFeXR0NDZv3oxevXrh5s2byMjIQFxcHE6fPo2OHTvaPS7Aw0UmkBESoPb4AcELsxuNRrzyyitYsmQJevXq5VBbN27ciIiICAwcOFBRPmLECPO/IyMjERsbi5CQEHz00UdITU116NgcMiYMkQRSCVljWeOC7I5qXJj90KFDNhdmv3fvHo4dO4YTJ04gJSUFACDLMogIWq0WRUVFGDp0qLn+gwcPkJOTg6VLl9ptg16vR2RkpGKBd3s4ZEwYoywBskpPplJmCxFh3rx52LlzJ0pLS+0uzN6uXTuLGcDMzEyUlJRgx44dFq/fvn076urq8Lvf/c5uW+rq6nDmzBnEx8c73H4OGRNGtjK7KDs5u+jswuwajcbifM3X1xfu7u6q53EbN27EmDFjVM+x0tLS8MILL6B79+6oqqpCRkYGampqMG3aNIfbzyFjwsgkQXLBvYuPszC7o7755hscOXIERUVFqvuvXbuGyZMn4/bt2+jcuTNiYmJQXl7u1MLsHDImjCwDksrQ0NknwpEDSy1lZ2fb3J+enq46c9mrVy+bx8/JybH73vZwyJgwrurJfuk4ZEwYe7OLrQWHjIkjSyC1mUQnZxd/6ThkTBjZyhS+zCFjzEVI/Y4P1bIWjEPGhCHZtKmVtyYcMiYMkfo5GU98MOYiPLtowiFj4vA5GQAOGRNJfriplbciHDImDvdkADhkTCCeXTThkDFhJFlSvUFYrawl45AxcejhplbeinDImDhWbqviexcZcxWeXQTAIWMi8ewiAA4ZE0iSTZtaeWvCIWPCSAAklUmO1tWPcciYSDxcBMAhYyLxxAcADhkTiM/JTDhkTBy+GA2AQ8YE4tuqTDhkTBgeLprwIoBMHLKxOUHUmtHZ2dmqa0bX1tYq6mVmZiI4OBju7u6IiorC4cOHnWo/h4yJI//YmzXdnJ1dFLVmNGBaAabpetAVFRVwd3c378/NzcX8+fOxcOFCnDhxAvHx8RgxYoRTz93n4SITxlXDRVFrRgOAJEk214P+61//itdeew2vv/46AOD999/H3r17sXbtWixfvtyh9nNPxprNk7Bm9A8//ICgoCB07doVo0ePxokTJ8z7DAYDjh8/jqSkJMVrkpKScPToUYfaCnDImEBqQ8WmvVu3bt3g4+Nj3hzpGZxdM3rLli1W14zu3bs3srOzUVhYiG3btsHd3R2DBg0yr6J5+/ZtGI1G+Pn5KV7n5+dnXiPNETxcZOIQ1M+/Hk58NPea0TExMYiJiTH/PGjQIDz77LP429/+hg8//NBcLknKSw5EZFFmC4eMCSORlRuEH5Y9KWtGN9JoNBgwYIC5J+vUqRPc3Nwseq2qqiqL3s0WHi4yYewNFx1FREhJSUF+fj5KSkocXjP65MmT5m3WrFkIDQ3FyZMnER0dbfV9Tp48iYCAAACATqdDVFQUiouLFfWKi4sRFxfncPu5J2PiuOgGYVFrRi9ZsgQxMTHo2bMnampq8OGHH+LkyZNYs2aNuU5qaiqSk5PRv39/xMbGYv369bhy5QpmzZrlcPs5ZEwYe8NFR4laM/r777/HjBkzUFlZCR8fH/Tr1w+HDh3CwIEDzXUmTpyI6upqLF26FBUVFYiIiMDu3budWjNaIkcW5HVA/JhVrjiMcO6FXzR3E+z69v0Y+5Wa2cXf/8HqvpqaGvj4+CD0zWVwa+tusd9YV4tzH/wf7t6969Q52S8V92RMGL530YRDxoRx1XDxl45DxsThb0YD4JAxgbgnM+GQMWE4ZCYcMiaOnduqWgsOGROGZxdNOGRMGB4umnDImDDck5m4LGQ34t1cdSihumCg/UrNLHqg48+veKLxFD4A7smYQDxcNOGQMWEkmSDJlolSK2vJOGRMGD4nM+GQMXH4Md0AOGRMIO7JTDhkTBgOmQmHjAnV2mYS1XDImDA8u2jCIWPCSEZAUnkemmT8+dvSnDhkTByeXQTAIWMC8XDRhEPGhOHbqkw4ZEwYnsI34ZAxYXi4aMIhY+LwxAcADhkTSDISJI1KT2ZsXSnjVV2YMBJZWdXlCVmYfcOGDYiPj0f79u3Rvn17DBs2DF98oXyMe3p6usXC7baWv1XDIWPiEFnfnCBqYfbS0lJMnjwZBw4cQFlZGbp3746kpCRcv35dUS88PFyxcPtXX33lVPt5uMiEedIXZt+yZYvi5w0bNmDHjh3Yv38/pk6dai7XarVO915NcU/GhJGMZHUDnoyF2Zt68OAB6uvrLY57/vx5dOnSBcHBwZg0aRIuXrzo0PEacciYMBKR1Q1o/oXZH/XWW28hMDAQw4YNM5dFR0dj8+bN2Lt3LzZs2IDKykrExcWhurraoWMCPFxkIslk2tTK0fwLsze1cuVKbNu2DaWlpXB3/3FNtREjRpj/HRkZidjYWISEhOCjjz5CamqqQ8fmkDFh7F2MflIWZl+9ejWWLVuGffv24ZlnnrHZBr1ej8jISPPi7Y7gkDFhXDXxQUSYN28edu7cidLSUocXZm8qMzMTJSUl2LFjh+L1q1atQkZGBvbu3Yv+/fvbbUtdXR3OnDmD+Ph4h9vPIWPi2BkuOkrUwuwrV67EokWLsHXrVvTo0cN8XC8vL3h5eQEA0tLS8MILL6B79+6oqqpCRkYGampqMG3aNIfbzxMfTBhJlq1uzli7di3u3r2LhIQEBAQEmLfc3FxzncdZmD0zMxMGgwHjx49XHHf16tXmOteuXcPkyZMRGhqKsWPHQqfToby83KmF2bknY+K4aOkkcuDidXZ2ts396enpSE9PV5RdvnzZ7nFzcnLs1rGHQ8aEkWSCpHICxnfhM+YqRiu34beyG4Q5ZEyYpheeHy1vTThkTBzZyhy+kxMfv3QcMiYODxcBcMiYQDxcNOGQMXGMVpbaNPJwkTHXIFn9/Is4ZIy5hrVvQfNwkTEXMRoBUnnwvdy6HobPIWPiGGX1oSFP4TPmIjxcBMAhYyLJVu4Q5nsXGXMRPicDwCFjIvFwEQCHjAlERiNIpScj7skYcxHZyuwiX4xmzEWs3YXPIWPMNchoBKmswq42hGzJOGRMHCP3ZACHjIlEVq6T8ewiY65hGi5aPnWQh4uMuUi9sRYEy0A1oL4ZWtN8JHLkoXaMOaG2thbBwcHmJ/Kq8ff3x6VLlxSLO7RUHDImRG1tLQwGg9X9Op2uVQQM4JAxJhw/C58xwThkjAnGIWNMMA4ZY4JxyBgTjEPGmGAcMsYE+3/Wbj9Z3JSJCQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 200x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(2,2))\n",
    "plt.title('Top Z-Slice Norm')\n",
    "plt.pcolormesh(norm_factor[:,:,-1].T)\n",
    "plt.colorbar()\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f464c7b7",
   "metadata": {
    "papermill": {
     "duration": 0.009722,
     "end_time": "2024-11-19T22:42:39.199623",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.189901",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now that the system matrix has been defined, we can reconstruct using the available reconstruction algorithms (up to this point, we haven't defined any subset functionality, so we can only use non-subset based algorithms) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "345794a5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.221804Z",
     "iopub.status.busy": "2024-11-19T22:42:39.221333Z",
     "iopub.status.idle": "2024-11-19T22:42:39.225220Z",
     "shell.execute_reply": "2024-11-19T22:42:39.224622Z"
    },
    "papermill": {
     "duration": 0.016722,
     "end_time": "2024-11-19T22:42:39.226446",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.209724",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample_object = torch.rand(object_meta.shape) # object has batch dimension\n",
    "sample_projections = system_matrix.forward(sample_object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dc1136c9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.247320Z",
     "iopub.status.busy": "2024-11-19T22:42:39.246801Z",
     "iopub.status.idle": "2024-11-19T22:42:39.250369Z",
     "shell.execute_reply": "2024-11-19T22:42:39.249859Z"
    },
    "papermill": {
     "duration": 0.015503,
     "end_time": "2024-11-19T22:42:39.251551",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.236048",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define system matrix\n",
    "system_matrix = EXSSystemMatrix(object_meta=object_meta, proj_meta=proj_meta)\n",
    "# Define likelihood that characterizes measured data (for SPECT/PET, this is PoissonLog, but here we'll use NegativeMSE)\n",
    "likelihood = NegativeMSELikelihood(system_matrix, projections=sample_projections, scaling_constant=0.01)\n",
    "# Define \n",
    "reconstruction_algorithm = MLEM(likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2385a26b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.272099Z",
     "iopub.status.busy": "2024-11-19T22:42:39.271619Z",
     "iopub.status.idle": "2024-11-19T22:42:39.285148Z",
     "shell.execute_reply": "2024-11-19T22:42:39.284502Z"
    },
    "papermill": {
     "duration": 0.025261,
     "end_time": "2024-11-19T22:42:39.286525",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.261264",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "recon = reconstruction_algorithm(n_iters=40)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb5e6b72",
   "metadata": {
    "papermill": {
     "duration": 0.010764,
     "end_time": "2024-11-19T22:42:39.307594",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.296830",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In this case our system matrix is underdetermined so we have no way of ensuring we get the right solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "30ff4c3c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.328919Z",
     "iopub.status.busy": "2024-11-19T22:42:39.328563Z",
     "iopub.status.idle": "2024-11-19T22:42:39.492258Z",
     "shell.execute_reply": "2024-11-19T22:42:39.491615Z"
    },
    "papermill": {
     "duration": 0.175661,
     "end_time": "2024-11-19T22:42:39.493643",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.317982",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fb3c05a3950>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWUAAAC8CAYAAACgw7/oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfXElEQVR4nO3df1AU5/0H8PcdhDtbv5wjNoCKSPwxYMivQqJAMU0TcUjiJJkQmPFbqIpJmGs1SnVaQlKN9RvHTKuHiaBOoDfOmIQmhsTMkOBlMgGptlMoOJ1Ca1JMRD1kMIlIEoG72+8f5Pa4u+eAPc+w7L1fM89k+Nze3rMzHz959tlnd3WSJEkgIiJV0E92B4iIyINFmYhIRViUiYhUhEWZiEhFWJSJiFSERZmISEVYlImIVIRFmYhIRViUiYhUhEWZiEhFWJQp5JqamrBq1SrMnj0bOp0O77zzzrjfaWxsRFpaGoxGI2655RYcOHDgxneUSIVYlCnkvv76a9xxxx145ZVXJrT92bNn8eCDDyI7OxttbW149tlnsXHjRhw9evQG95RIfXR8IBHdSDqdDnV1dXj00UcDbvOb3/wGx44dQ2dnpxwrKSnB6dOncerUqe+hl0TqwZEyTbpTp04hJyfHK7Zy5Uq0tLRgeHh4knpFNDkiJ7sDpA7Xrl3D0NCQ8DNJkqDT6bxiBoMBBoMhJL/d09OD2NhYr1hsbCwcDgf6+voQHx8fkt+hqW2sHAWAqKgoGI3G77FHN8aEi/Ktv917I/sRtIHF6htJJS+4ONldEGq41yKMX7t2DUmJ09HT6xR+Pn36dAwMDHjFtm3bhu3bt4esb75F3z2r5hsfS9r75SHrT0i9PWuye+DnyoqvJ7sLQp/mPyeMj5ejABAXF4ezZ89O+cLMkTJhaGgIPb1OnGmZi+j/8Z7R6r/qwuL08+ju7kZ0dLQcD9UoGRj5x9TT0+MV6+3tRWRkJGJiYkL2OzR1jZWjgCdPh4aGWJRJO6b/z0gbzfXdf6Ojo72KcihlZGTgvffe84odP34c6enpuOmmm27Ib9LUJMpRwJOnWsALfSQbllzCptTAwADa29vR3t4OYGTJW3t7O86dOwcAKCsrQ1FRkbx9SUkJPv/8c5SWlqKzsxM1NTWorq7Gli1bQnJcpB2BcjSYPFUrjpRJ5oALvjP0jiDGIC0tLbjvvvvkv0tLSwEAv/jFL2C1WmG32+UCDQBJSUmor6/H5s2bsX//fsyePRv79u3D448/HtRxkHaJctQd1woWZZI5JQlOn2Xrvn9PxE9/+lOMtfzdarX6xe6991784x//UPxbFF5EOeqOawWLMsmGIWEYkl+MSC1EOeqOawWLMsmGpZHmGyNSC1GOuuNawaJMMhd0cELnFyNSC1GOuuNawaJMsmFJh2FJ5xcjUgtRjrrjWsGiTLJhSY9hSe8Tm6TOEAmIcnQkPgmduUFYlEnmFJwaik4ViSaLKEfdca1gUSaZQ4rwG4U4NHRaSFOfKEdH4trJUxZlkg1JEbjJJ+GHNJTsNPWJcnQkrp08ZVEmmQs6uHzuvHdpaP0nTX2iHB2JaydPWZRJNjIKifCJTVJniAREOToSn4TO3CB8IBHJHIjAsE9zwP8fANFkEeVosHlaWVmJpKQkGI1GpKWl4cSJEwG3XbNmDXQ6nV+79dZb5W2sVqtwm2vXrinqF4syyZySXtiI1CJQjirN09raWmzatAnl5eVoa2tDdnY2cnNzvR6UNVpFRQXsdrvcuru7MXPmTDzxxBNe20VHR3ttZ7fbFT/fmdMXJBuWIjDsc2qopfWfNPWJcnQkrmw/e/bsQXFxMdavXw8AsFgsaGhoQFVVFXbt2uW3vclkgslkkv9+55138OWXX2Lt2rVe2+l0OsTFxSnrjA8Og0g2LEUKG5FaBMpRJXk6NDSE1tZWv5f15uTk4OTJkxPaR3V1NR544AEkJiZ6xQcGBpCYmIi5c+fi4YcfRltb24T75cZ/cSRzQg+nz/+nnRq6qk1TnyhHR+Ijedrf3+8VF73gt6+vD06nU/iyXt/XkonY7Xa8//77eO2117ziycnJsFqtuO2229Df34+KigpkZWXh9OnTWLRo0YSOD+BImUZxQC+fHrqbgylCKiLK0dF5mpCQIE81mEwm4VSEm+hlvRN5Ua/VasWMGTPw6KOPesWXLVuGn//857jjjjuQnZ2NP//5z1i8eDFefvllRcfIkTLJhqUIRPrNKXOkTOohytGR+EieTuQFv7NmzUJERITwZb2+o2dfkiShpqYGhYWFiIqKGnNbvV6Pu+++G5988smY2/l9T9HWpGlcfUFqN97qC/cLft1NVJSjoqKQlpYGm83mFbfZbMjMzBzz9xsbG/Hpp5+iuLh43L5KkoT29nbEx8crOEKOlGkUjpRJ7cYbKU9UaWkpCgsLkZ6ejoyMDBw6dAjnzp1DSUkJgJGX+164cAGHDx/2+l51dTWWLl2K1NRUv32+8MILWLZsGRYtWoT+/n7s27cP7e3t2L9/v6K+sSiTzCG4iu1gTSYVEeXoSFzZfgoKCnD58mXs2LEDdrsdqampqK+vl1dT+L7cFwCuXLmCo0ePoqKiQrjPr776Ck899RR6enpgMplw1113oampCffcc4+ivvHclGTuxyL6tmAouVsKAI4cOYI77rgDP/jBDxAfH4+1a9fi8uXLQf02aVegHA0mT81mMz777DMMDg6itbUVy5cvlz+zWq34+OOPvbY3mUz45ptv8OSTTwr3t3fvXnz++ecYHBxEb28vGhoakJGRobhfLMokG3mAuO+VbeUpovRuqebmZhQVFaG4uBj/+te/8Oabb+Lvf/+7vLCfyE2co8HlqVpp50joug27IoVNqdF3S6WkpMBisSAhIQFVVVXC7f/6179i/vz52LhxI5KSkvCTn/wETz/9NFpaWq73kEhjAuVoMHmqVizKJJOg++7RiJ4mfXda2N/f79UGBweF+wjmbqnMzEycP38e9fX1kCQJly5dwltvvYWHHnootAdIU54oR0fnqRawKJNs2BUhbMDEF+UHc7dUZmYmjhw5goKCAkRFRSEuLg4zZsxQvOietC9QjrrzVAu0M+an6zYsRUDvtyTOBWBii/JHU3K3VEdHBzZu3Ijf/e53WLlyJex2O7Zu3YqSkhJUV1cHcyikUaIcHYm7JqE3NwaLMslc0AvePOK9KH88wdwttWvXLmRlZWHr1q0AgNtvvx0//OEPkZ2djZ07dypefE/aJcpRd1wrtHMkdN2GXXphUyKYu6W++eYb6PXevxMRMTIaknjzCo0SKEeV5qmacaRMMofg1NARxGmh0rulVq1ahSeffBJVVVXy9MWmTZtwzz33YPbs2dd/YKQZohwdiXP6gjTIJQEun7cCu4IYqCq9W2rNmjW4evUqXnnlFfz617/GjBkz8LOf/Qy7d+++ruMh7RHlqDuuFSzKJBOPlIO7qm02m2E2m4WfWa1Wv9iGDRuwYcOGoH6LwkfgkTJXX5AGDbsioPNZWqSlpUY09Yly1B3XChZlkrmg85++0NCifJr6RDnqjmsFizLJnJIeDp9nCPB5yqQmohx1x7WCRZlkDsGpoUNDp4U09Yly1B3XChZlkrmfI+AbI1ILUY6641rBokwyh0sPnc8ifIeGFuXT1CfKUXdcK1iUScaiTGrHokxhRYL/aaCG1uSTBohy1B3XChZlkjlceoAjZVIxUY7KcY1gUSYZizKpHYsyhRVJ0kGSfJ+DrJ2r2jT1iXLUHdcKFmWSOSQ94LMIX7RQn2iyiHJUjmsEizLJnIIr204NnRbS1CfKUXdcK1iUScbpC1I7Tl+MYlh++Ub2I2gJW9W3GKY3O3GyuyB279gfuyQdnC7f5ylPnWTXvxkz2V0QmmEVv8V7Mhm/XDrZXRDLH/tjUY6641rBkTLJnIL5Oi096IWmPlGOynGN0M6R0HVzSTphI1KLQDkaTJ5WVlYiKSkJRqMRaWlpOHHiRMBt16xZA51O59duvfVWr+2OHj2KJUuWwGAwYMmSJairq1PcLxZlkrlcOmEjUotAOao0T2tra7Fp0yaUl5ejra0N2dnZyM3N9XpN2WgVFRWw2+1y6+7uxsyZM/HEE0/I25w6dQoFBQUoLCzE6dOnUVhYiPz8fPztb39T1DcWZZI5XXphI1KLQDmqNE/37NmD4uJirF+/HikpKbBYLEhISEBVVZVwe5PJhLi4OLm1tLTgyy+/xNq1a+VtLBYLVqxYgbKyMiQnJ6OsrAz3338/LBaLor7xXxzJJEnciNQiUI6687S/v9+rDQ4O+u1jaGgIra2tyMnJ8Yrn5OTg5MmJXZStrq7GAw88IL8MGBgZKfvuc+XKlRPepxuLMslGTgP1Pi246Qsl83UAMDg4iPLyciQmJsJgMGDBggWoqakJ6rdJu8Q56snThIQEmEwmue3atctvH319fXA6nYiNjfWKx8bGoqenZ9w+2O12vP/++1i/fr1XvKenJ+h9jsbVFyRzSTrofN/RF8QFFPd8XWVlJbKysnDw4EHk5uaio6MD8+bNE34nPz8fly5dQnV1NRYuXIje3l44HI6gjoO0S5Sj7jgAdHd3Izo6Wo4bDIaA+9LpfNfkS34xEavVihkzZuDRRx8N2T5HY1EmDwn+z0AMYvpi9HwdMDLX1tDQgKqqKuHI5YMPPkBjYyO6urowc+ZMAMD8+fOV/zBpnyhH4YlFR0d7FWWRWbNmISIiwm8E29vb6zfS9fsZSUJNTQ0KCwsRFRXl9VlcXFxQ+/TF6QuSSYIr2tJ3p4UTmasDgpuvO3bsGNLT0/HSSy9hzpw5WLx4MbZs2YJvv/02tAdIU54oR0fn6URERUUhLS0NNpvNK26z2ZCZmTnmdxsbG/Hpp5+iuLjY77OMjAy/fR4/fnzcffriSJlkkksPyecqtvvvhIQEr/i2bduwfft2v30EM1/X1dWF5uZmGI1G1NXVoa+vD2azGV988QXnlcmLKEfdcSVKS0tRWFiI9PR0ZGRk4NChQzh37hxKSkoAAGVlZbhw4QIOHz7s9b3q6mosXboUqampfvt85plnsHz5cuzevRuPPPII3n33XXz44Ydobm5W1DcWZZKJVlu4/1YyVwcom1tzuVzQ6XQ4cuQITCYTgJEpkLy8POzfvx/Tpk1TeCSkVYFWBCldJVRQUIDLly9jx44dsNvtSE1NRX19vbyawm63+61ZvnLlCo4ePYqKigrhPjMzM/HGG2/gueeew/PPP48FCxagtrYWS5cqu6WdRZlkkuA00P33RObqgODm6+Lj4zFnzhy5IANASkoKJEnC+fPnsWjRIqWHQholylF3XCmz2Qyz2Sz8zGq1+sVMJhO++eabMfeZl5eHvLw8xX0ZjXPKJJMknZz0clO4+iKY+bqsrCxcvHgRAwMDcuzMmTPQ6/WYO3eu8gMhzRLmaBB5qmYsyuQhBWgKlZaW4tVXX0VNTQ06OzuxefNmv/m6oqIiefvVq1cjJiYGa9euRUdHB5qamrB161asW7eOUxfkLVCOaugmJ05fkIekG2m+MYWUztdNnz4dNpsNGzZsQHp6OmJiYpCfn4+dO3de1+GQBoly1B3XCBZl8nDpRppvLAhK5+uSk5P9pjyI/Ihy1B3XCBZlko21+oJIDUK1+kLNWJTJI4QjZaIbgiNlCic610jzjRGphShH3XGtYFEmjxBd6CO6YXihj8KK67vmGyNSC1GOIkBsimJRJg/OKZPacU6ZwolOGmm+MSK1EOWoO64VLMrkEaLnKRPdMOM8T1kLWJRJpoNgpDwpPSESE+WoO64VLMrkwTllUjvOKVM44TplUjuuU6bwwjllUjvOKVM44UiZ1I4jZQovnFMmteOcMoUTrlMmteM6ZQovolNDDZ0WkgYEmL7QUp6yKJMHn31BasdnX1A44fQFqR2nLyi8cEkcqV0YLInj26xJppM8S47kFmSyV1ZWIikpCUajEWlpaThx4sSEvveXv/wFkZGRuPPOO4P7YdI0YY5eR56qEYsyeYTo1e21tbXYtGkTysvL0dbWhuzsbOTm5nq9wVrkypUrKCoqwv3336/8Ryk8BMpRFmXSIuEIJIgLKHv27EFxcTHWr1+PlJQUWCwWJCQkoKqqaszvPf3001i9ejUyMjKCPALSukA5qqWbR1iUSRaKZB8aGkJraytycnK84jk5OTh58mTA7/3pT3/Cf//7X2zbti2YrlOYCGVRVjrFNjg4iPLyciQmJsJgMGDBggWoqamRP7dardDpdH7t2rVrivrFC33kMcaFvv7+fq+wwWCAwWDw20VfXx+cTidiY2O94rGxsejp6RH+7CeffILf/va3OHHiBCIjmZI0hhBd6HNPsVVWViIrKwsHDx5Ebm4uOjo6MG/ePOF38vPzcenSJVRXV2PhwoXo7e2Fw+Hw2iY6Ohr/+c9/vGJGo1FR3/gvgGRjPfsiISHBK75t2zZs37498L503re9SpLkFwMAp9OJ1atX44UXXsDixYuD6jeFj1A9+2L0FBsAWCwWNDQ0oKqqCrt27fLb/oMPPkBjYyO6urowc+ZMAMD8+fP9+6HTIS4uTllnfHD6gmRjnRZ2d3fjypUrcisrKxPuY9asWYiIiPAbFff29vqNngHg6tWraGlpwa9+9StERkYiMjISO3bswOnTpxEZGYmPPvoo5MdJU9d40xf9/f1ebXBw0G8fwUyxHTt2DOnp6XjppZcwZ84cLF68GFu2bMG3337rtd3AwAASExMxd+5cPPzww2hra1N8jBwpk8cY0xfR0dGIjo4edxdRUVFIS0uDzWbDY489JsdtNhseeeQRv+2jo6Pxz3/+0ytWWVmJjz76CG+99RaSkpKUHgVp2TjTFxM5owtmiq2rqwvNzc0wGo2oq6tDX18fzGYzvvjiC3leOTk5GVarFbfddhv6+/tRUVGBrKwsnD59GosWLZrwIbIokyxUj+4sLS1FYWEh0tPTkZGRgUOHDuHcuXMoKSkBAJSVleHChQs4fPgw9Ho9UlNTvb5/8803w2g0+sWJxpu+6O7u9ho8iK57yN+Z4BQbALhcLuh0Ohw5cgQmkwnAyBRIXl4e9u/fj2nTpmHZsmVYtmyZ/J2srCz8+Mc/xssvv4x9+/ZN9BBZlMkjVLdZFxQU4PLly9ixYwfsdjtSU1NRX1+PxMREAIDdbh93zTKRyHi3WU/kjE7pFBsAxMfHY86cOXJBBoCUlBRIkoTz588LR8J6vR533303Pvnkk3GOyud7irYmbQvhonyz2YzPPvsMg4ODaG1txfLly+XPrFYrPv7444Df3b59O9rb24P7YdK2ENw8MnqKbTSbzYbMzEzhd7KysnDx4kUMDAzIsTNnzkCv12Pu3LnirkoS2tvbER8fP/HOgUWZRtH6onya+kK1Trm0tBSvvvoqampq0NnZic2bN/tNsRUVFcnbr169GjExMVi7di06OjrQ1NSErVu3Yt26dZg2bRoA4IUXXkBDQwO6urrQ3t6O4uJitLe3y/ucKE5fkMz9XAHfGJFaiHLUHVdC6RTb9OnTYbPZsGHDBqSnpyMmJgb5+fnYuXOnvM1XX32Fp556Cj09PTCZTLjrrrvQ1NSEe+65R1HfWJTJg0+JI7UL4VPizGYzzGaz8DOr1eoXS05O9pvyGG3v3r3Yu3ev8o74YFEmGV+cSmrHF6eOkju340b2I2jvvHj7ZHfBz8KY/052F4Iy1Yvy3148MNldELp11f9Odhf8HLt7z2R3IYAtY37KokzhhdMXpHZh8JB7FmWSTfWRMmkfR8oUVnQuCTqX5BcjUgtRjrrjWsGiTDK+OJXUji9OpbDC6QtSO05fUFhhUSa1Y1Gm8MLVF6R2XH1BYUUSXESRNJTtNPWJcvS7uFawKJOM0xekdpy+oLDCokxqx6JMYYVFmdSORZnCCm8eIbXjzSMUXrj6gtSOqy8onHCkTGrHkTKFFc4pk9qFw5wy39FHMvdzBXxbMCorK5GUlASj0Yi0tDScOHEi4LZvv/02VqxYgR/96EeIjo5GRkYGGhoagjwK0rJAOaqlZ1+wKJOHUxI3hWpra7Fp0yaUl5ejra0N2dnZyM3N9Xrn2WhNTU1YsWIF6uvr0draivvuuw+rVq1CW1vb9R4RaU2gHA0iT9WK0xck0wnultIFcafUnj17UFxcjPXr1wMALBYLGhoaUFVVhV27dvltb7FYvP5+8cUX8e677+K9997DXXfdpfj3SbtEOeqOawVHyiQLxWnh0NAQWltbkZOT4xXPycnByZMnJ7QPl8uFq1evYubMmcp+nDQvHKYvOFIm2VirL/r7+73iBoMBBoPBbx99fX1wOp2IjY31isfGxqKnp2dC/fjjH/+Ir7/+Gvn5+Uq6T2EgHFZfcKRMMp1TEjYASEhIgMlkkptoGsJrXzqd19+SJPnFRF5//XVs374dtbW1uPnmm4M/GNKkQDmq45wyadIYN490d3cjOjpaDotGyQAwa9YsRERE+I2Ke3t7/UbPvmpra1FcXIw333wTDzzwgNLeUzgIg5tHOFImmfvU0LcBQHR0tFcLVJSjoqKQlpYGm83mFbfZbMjMzAz426+//jrWrFmD1157DQ899FDoDoo0JVCOamn6giNlkumcEnQ+V0yCOS0sLS1FYWEh0tPTkZGRgUOHDuHcuXMoKSkBAJSVleHChQs4fPgwgJGCXFRUhIqKCixbtkweZU+bNg0mk+k6j4q0RJSj7rhWcKRMHpIkbgoVFBTAYrFgx44duPPOO9HU1IT6+nokJiYCAOx2u9ea5YMHD8LhcOCXv/wl4uPj5fbMM8+E7NBIIwLlaBB5quQGJwAYHBxEeXk5EhMTYTAYsGDBAtTU1Hhtc/ToUSxZsgQGgwFLlixBXV2d4n5xpEyyUD77wmw2w2w2Cz+zWq1ef3/88cdB/QaFn1CtvnDf4FRZWYmsrCwcPHgQubm56OjowLx584Tfyc/Px6VLl1BdXY2FCxeit7cXDodD/vzUqVMoKCjA73//ezz22GOoq6tDfn4+mpubsXTp0gn3jUWZPFyCO6M0NFdHGiDKUXdcAaU3OH3wwQdobGxEV1eXvH5+/vz5XttYLBasWLECZWVlAEam6RobG2GxWPD6669PuG+cviCZTpKEjUgtAuWoO0/7+/u92uDgoN8+grnB6dixY0hPT8dLL72EOXPmYPHixdiyZQu+/fZbeZtTp0757XPlypUTvmnKjSNl8nBJ/o/b4kiZ1ESUo+44RtbTj7Zt2zZs377dKxbMDU5dXV1obm6G0WhEXV0d+vr6YDab8cUXX8jzyj09Pdd105QbizLJdE4JOp8Fn1q6qk1TnyhH3XFg4uvpAWU3OLlcLuh0Ohw5ckReEbRnzx7k5eVh//79mDZtmuJ9BsKiTB6iq9icviA1CbTSQvJeTz+WYG5wio+Px5w5c7yWaKakpECSJJw/fx6LFi1CXFxcUDdN+eKcMnm4XOJGpBaBclRBngZzg1NWVhYuXryIgYEBOXbmzBno9XrMnTsXAJCRkeG3z+PHj49505QIizLJtP5MAZr6QvXsi9LSUrz66quoqalBZ2cnNm/e7HeDU1FRkbz96tWrERMTg7Vr16KjowNNTU3YunUr1q1bJ09dPPPMMzh+/Dh2796Nf//739i9ezc+/PBDbNq0SVHfOH1BHpy+ILUbZ/piogoKCnD58mXs2LEDdrsdqampY97gNH36dNhsNmzYsAHp6emIiYlBfn4+du7cKW+TmZmJN954A8899xyef/55LFiwALW1tYrWKAMsyjSa0wXAJYgRqYQoR+W4MkpucAKA5ORkv+kJX3l5ecjLy1Pcl9FYlMlDEszNSSzKpCKiHHXHNYJFmTxcgucicp0yqYkoR+W4NrAok4fLCcApiBGphChH5bg2sCiTh9PlfxrIJXGkJqIcBTSVpyzK5CFBsPpiUnpCJCbKUXdcI1iUycPpBCROX5CKiXIU0FSesiiTh0uw3EhDp4WkAaIclePawKJMHlx9QWrH1RcUTiSXE5LPqaHv30STSZSjgLbylEWZPJxOQOeT3BpKdtIAUY4CmspTFmXykASnhnz2BamJKEfluDbwKXEkk5xOYQuG0jcFNzY2Ii0tDUajEbfccgsOHDgQ1O+StgXK0WDzVI1YlMnD6Ro5PfRqyq9qu98UXF5ejra2NmRnZyM3N9frqVujnT17Fg8++CCys7PR1taGZ599Fhs3bsTRo0ev94hIa4Q5GlyeqhWLMskklyRsSo1+U3BKSgosFgsSEhJQVVUl3P7AgQOYN28eLBYLUlJSsH79eqxbtw5/+MMfrveQSGMC5WgweapWLMokC8VpYTBvCg70FuCWlhYMDw8rOwjStHCYvuCFPpI5pEG/5wo4MFIU+/v7veIGg0H4Uspg3hQc6C3ADocDfX19iI+PV3wspE2iHAU8eaoFEy7K/3f72zeyH0H7v9snuwdTX1RUFOLi4tDcUy/8fPr06RN6dftoSt/qK9peFB+LPu7MhLf9PnU+Ntk9ENk22R1QZLwcBYC4uDhERUV9j726MThSJhiNRpw9exZDQ0PCz0UFNdCr24N5U3CgtwBHRkYiJiZmoodBGjZejgIjhdtoNH6PvboxWJQJwEjShyKhR78p+LHHPENEm82GRx55RPidjIwMvPfee16x48ePIz09HTfddNN194m0IVQ5qna80Echp/RNwSUlJfj8889RWlqKzs5O1NTUoLq6Glu2bJmsQyCaNBwpU8gpfVNwUlIS6uvrsXnzZuzfvx+zZ8/Gvn378Pjjj0/WIRBNGp0kaej+RCKiKY7TF0REKsKiTESkIizKREQqwqJMRKQiLMpERCrCokxEpCIsykREKsKiTESkIizKREQqwqJMRKQiLMpERCrCokxEpCL/D5Jee+Zg7MDBAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 400x200 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4,2))\n",
    "plt.subplot(121)\n",
    "plt.pcolormesh(sample_object[:,:,0], vmin=0, vmax=1)\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "plt.subplot(122)\n",
    "plt.pcolormesh(recon[:,:,0])\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b6d7e80",
   "metadata": {
    "papermill": {
     "duration": 0.010405,
     "end_time": "2024-11-19T22:42:39.514559",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.504154",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Part 3: Incorporating Subsets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8d501ea",
   "metadata": {
    "papermill": {
     "duration": 0.010135,
     "end_time": "2024-11-19T22:42:39.534685",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.524550",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Up until now, we've ignored the `subset_idx` which specifies how projection data can be split into subsets in ordered subset based reconstruction algorithms. In such algorithms, the projection data is split up into disjoint subsets: in SPECT/PET, one typically partitions using the projection angle. We will do the same thing here, giving the option of `subset_idx=0` for the first projection angle and `subset_idx=1` for the second:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "880773c9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.556841Z",
     "iopub.status.busy": "2024-11-19T22:42:39.555995Z",
     "iopub.status.idle": "2024-11-19T22:42:39.562562Z",
     "shell.execute_reply": "2024-11-19T22:42:39.561969Z"
    },
    "papermill": {
     "duration": 0.018752,
     "end_time": "2024-11-19T22:42:39.563758",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.545006",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSSystemMatrix(SystemMatrix):\n",
    "    def compute_normalization_factor(self):\n",
    "        norm_projections = torch.ones((2,self.proj_meta.M, self.proj_meta.M))\n",
    "        return self.backward(norm_projections)\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        projection_0degrees = object.sum(dim=0)\n",
    "        projection_90degrees = object.sum(dim=0)\n",
    "        if subset_idx==0:\n",
    "            projections = projection_0degrees\n",
    "        elif subset_idx==1:\n",
    "            projections = projection_90degrees\n",
    "        else:\n",
    "            projections = torch.stack([projection_0degrees, projection_90degrees], dim=0)\n",
    "        projections *= self.proj_meta.sensitivity_factor\n",
    "        return projections\n",
    "    def backward(self, proj, subset_idx = None):\n",
    "        # Back projection expects projections in their subset\n",
    "        if subset_idx is not None:\n",
    "            object_BP = (proj[0]*self.proj_meta.sensitivity_factor).unsqueeze(subset_idx).repeat_interleave(self.proj_meta.M, subset_idx)\n",
    "        else:\n",
    "            object_BP_angle0 = (proj[0]*self.proj_meta.sensitivity_factor).unsqueeze(0).repeat(self.proj_meta.M,1,1)\n",
    "            object_BP_angle90 = (proj[1]*self.proj_meta.sensitivity_factor).unsqueeze(1).repeat(1,self.proj_meta.M,1)\n",
    "            object_BP = object_BP_angle0 + object_BP_angle90\n",
    "        return object_BP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0e7d5b2f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.585291Z",
     "iopub.status.busy": "2024-11-19T22:42:39.585062Z",
     "iopub.status.idle": "2024-11-19T22:42:39.589592Z",
     "shell.execute_reply": "2024-11-19T22:42:39.588866Z"
    },
    "papermill": {
     "duration": 0.017314,
     "end_time": "2024-11-19T22:42:39.590964",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.573650",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSSystemMatrix(object_meta=object_meta, proj_meta=proj_meta)\n",
    "FP_subset0 = system_matrix.forward(sample_object, subset_idx=0)\n",
    "FP_subset1 = system_matrix.forward(sample_object, subset_idx=1)\n",
    "BP_subset0 = system_matrix.backward(FP_subset0, subset_idx=0)\n",
    "BP_subset1 = system_matrix.backward(FP_subset0, subset_idx=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "847a197b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.613788Z",
     "iopub.status.busy": "2024-11-19T22:42:39.613142Z",
     "iopub.status.idle": "2024-11-19T22:42:39.692467Z",
     "shell.execute_reply": "2024-11-19T22:42:39.691797Z"
    },
    "papermill": {
     "duration": 0.092724,
     "end_time": "2024-11-19T22:42:39.693857",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.601133",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAADCCAYAAACYAHMjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAQAElEQVR4nO3de0zV9ePH8dfhIjcTh6CQTRhWCGo/dVpeMkzcVGrm/ZI3ZN51zszWzKlJdtHf+pouGDibaJn3UJtNy/uaZLIU08TG8DJTvKCmmDiQ9+8Pv5yfh+sRwXer52Njk895n8/lnDdPPudzjuowxhgBAKzwsL0DAPBvRoQBwCIiDAAWEWEAsIgIA4BFRBgALCLCAGAREQYAi4gwAFjkdoQdDodbX/v376/H3ZUSEhJq3IeIiIga11NcXKy0tDR16tRJQUFB8vf3V3h4uN544w1lZGQ4x509e1YOh0Pp6enOZenp6XI4HDp79mzdH2A1yh9nQECAoqOjtXDhQt25c8dlbPnHycfHR1FRUVqwYIGKiopq3NapU6c0evRoRUZGytfXV8HBwerQoYOmT5+uW7duPfK+79+/Xw6HQ5s3b3Zr/JUrV5SQkKDg4GD5+/urS5cu2rNnT7X3YY6mO5cxR+t3jl64cEEzZ85UbGysGjduXOHxfxRe7g7MzMx0+f6DDz7Qvn37tHfvXpflMTExtdoRd82bN0+TJ0+u9Lb09HSlpaVpwIABNa5n9OjR+uabbzRz5kwtXLhQPj4+ysvL086dO7Vr165q1/Haa68pMzNTYWFhtT6O2ho8eLDefvttSVJhYaEOHDigpKQkHT9+XFu2bHEZ6+fn53x+bty4oXXr1ikpKUk5OTnasGFDlds4evSounXrpujoaM2fP18RERG6du2asrOztX79es2ePVuNGjWqt2O8d++e4uLidPPmTS1btkxNmzZVcnKy+vTpo927dys2NrbS+zFH/x9ztH7naG5urtauXat27dopPj5e69atq/3KTC2NHTvWBAQE1PbudS4zM9M0aNDAvPLKK6a4uLjasXl5eUaSmT9/fqW3379/3/nnM2fOGElm1apVdbm7tSLJTJs2rcLy0aNHGw8PD3P37l3nsqqen+7duxtJ5sKFC1VuZ8yYMSYgIMDcunWr0ttLS0sfed/37dtnJJlNmzbVODY5OdlIMocOHXIuKy4uNjExMebFF190e5vM0Sfv3zJHH378jxw58liPf51eE75+/bqmTp2q5s2bq0GDBoqMjNTcuXN17949l3EOh0PTp09XWlqann/+efn4+CgmJkbr16+v1Xbz8/M1aNAghYSEaOPGjfLyqv4Ev6CgQJKqPEvw8Kj+Yanqpd7OnTsVFxenwMBA+fv7Kzo6Wh9//LHLmKysLPXr109BQUHy9fVV+/bttXHjxhqOsHqBgYFyOBzy9PSscWznzp0lSefOnatyTEFBgRo1aqSGDRtWervD4XD+OSIiQgkJCRXG9OjRQz169KiwvKioSLNmzVJoaKj8/PwUGxuro0ePuozJyMhQVFSUunTp4lzm5eWlUaNG6eeff9Yff/xR3SFWiznKHC3zOHO0psf/UdTZmoqKivTqq69qzZo1mjVrlnbs2KFRo0ZpyZIlGjhwYIXx27dv1/Lly5WUlKTNmzcrPDxcI0aMcPuaYZni4mINGTJE165d0+bNm9WsWbMa7xMdHa3GjRtr4cKFWrFiRZ1cN/viiy8UHx+v0tJSpaam6ttvv9WMGTN04cIF55h9+/apW7duunnzplJTU7Vt2za1a9dOw4YNc/t6kjFGJSUlKikp0c2bN7Vt2zatXr1aw4cPl7e3d433z83NlSSFhIRUOaZLly66dOmSRo4cqQMHDuju3btu7Zs73nvvPeXl5WnlypVauXKlLl68qB49eigvL8855sSJE3rhhRcq3Lds2cmTJ2u1beYoc9Qd7szROlWr82dT8aVEamqqkWQ2btzoMm7x4sVGkvn++++dyyQZPz8/k5+f71xWUlJiWrVqZZ599tlH2o+pU6caSSY1NfWR7rdjxw4THBxsJBlJpkmTJmbIkCFm+/btLuMqe6m3atUqI8mcOXPGGGPM7du3TaNGjczLL79c7cugVq1amfbt21d4Kfr666+bsLAwl5c4lSnb1/Jfffv2NYWFhS5jy56f4uJiU1xcbK5evWqWLVtmHA6H6dSpU7XbKSoqMv3793eu39PT07Rv397MnTvXXLlyxWVseHi4GTt2bIV1xMbGmtjYWOf3ZS/1OnTo4PIYnT171nh7e5vx48c7l3l7e5tJkyZVWOehQ4eMJPP1119Xu//lH4MyzFHm6MMeZ44+7G9zOWLv3r0KCAjQ4MGDXZaXvQwo/852XFycyxmBp6enhg0bptzcXJffzNVJT09XSkqKEhMTNWnSpAq3l5aWOn8jl5SU6P79+87b4uPjdf78eWVkZGj27Nlq3bq1tm7dqn79+mn69OnuHrYk6dChQ7p165amTp3q8jLoYbm5ucrJydHIkSMlyWW/4uPjdenSJZ0+fbrGbQ0dOlRHjhzRkSNHdPDgQS1fvlxZWVnq06dPhZfUd+7ckbe3t7y9vRUSEqKZM2eqb9++Lu+sV8bHx0cZGRn67bfftHTpUg0fPlxXr17Vhx9+qOjoaLf2sypvvvmmy2MUHh6url27at++fS7jqnoca7qtOsxR5qg73J2jdcXtT0fUpKCgQKGhoRWe4KZNm8rLy8t5jatMaGhohXWULSsoKNAzzzxT7faysrI0ZcoUdezYUSkpKZWOSUxM1OrVq53fx8bGunw8yc/PT/3791f//v0lSefPn1ffvn2VnJysKVOmqHXr1tXuQ5mrV69KUrX7fPnyZUnS7NmzNXv27ErHXLt2rcZthYSEqGPHjs7vu3fvrpCQEI0YMULp6ekuP+h+fn46ePCgpAeTNjw8/JHeMY6OjlZ0dLSkBy8xP/vsM82aNUvz5s2r9TXCqp737Oxs5/dNmjSpMF+kB9dzJSkoKKhW22aOMkfd4c4crUt1FuEmTZro8OHDMsa4TPIrV66opKREwcHBLuPz8/MrrKNsWZMmTard1tWrVzVw4EA1bNhQW7ZskY+PT6Xj3n//fZczhqeeeqra9bZo0UITJ07UzJkzdfLkSbcneNm1q+rOjsqOf86cOZVef5SkqKgot7ZXXtm10vKTxMPDw+WH4XE4HA699dZbSkpK0okTJ5zLfX19K5zdSA9+WMs/51LVz/vDz3nbtm3166+/VhhXtqxNmza1OgbmKHP0YY8zR+tSnV2OiIuLU2FhobZu3eqyfM2aNc7bH7Znzx7nb15Jun//vjZs2KCWLVtW+9u6pKREQ4YM0cWLF7Vhwwa1aNGiyrERERHq2LGj86tsAt2+fVuFhYWV3ufUqVOSpKeffrrqgy2na9euCgwMVGpqqkwV/1tUVFSUnnvuOWVnZ7vs08NfNf0AVuXYsWOSHpzR1YVLly5VuvzixYu6deuWy2MTERGh48ePu4z7/fffq3w5uG7dOpfH6Ny5czp06JDLu9QDBgxQTk6ODh8+7FxWUlKir776Si+99NIjPTcPY44yR8s87hytS3V2JjxmzBglJydr7NixOnv2rNq2basff/xRH330keLj49WrVy+X8cHBwerZs6fmzZungIAApaSkKCcnp8aPAL3zzjs6cOCARo4cKX9/f/3000+Vjiv7mEtlTp8+rd69e2v48OGKjY1VWFiYbty4oR07dmjFihXq0aOHunbt6vaxN2zYUJ9++qnGjx+vXr16acKECWrWrJlyc3OVnZ2tzz//XJKUlpamvn37qnfv3kpISFDz5s11/fp1nTp1Sr/88os2bdpU47YuX77sPOaioiIdO3ZMixYtUuPGjTVu3Di397k6EydO1M2bNzVo0CC1adNGnp6eysnJ0dKlS+Xh4aF3333XOXb06NEaNWqUpk6dqkGDBuncuXNasmRJle9sX7lyRQMGDNCECRP0559/asGCBfL19dWcOXOcYxITE5WcnKwhQ4bok08+UdOmTZWSkqLTp09r9+7dtT4u5ihztK7mqCTnp2TKPjWRlZXl/Mhc+fcdqlWrt/NM5R+0LigoMJMnTzZhYWHGy8vLhIeHmzlz5piioiKXcfrvB7pTUlJMy5Ytjbe3t2nVqpVZu3ZtjdsNDw+v8h3Yh7+qc+PGDbNo0SLTs2dP07x5c9OgQQMTEBBg2rVrZxYtWmT++usv51h33nku891335nY2FgTEBBg/P39TUxMjFm8eLHLmOzsbDN06FDTtGlT4+3tbUJDQ03Pnj3deue8/DF6e3ubyMhIM27cOJObm+sy9nH+osKuXbtMYmKiiYmJMYGBgcbLy8uEhYWZgQMHmszMTJexpaWlZsmSJSYyMtL4+vqajh07mr1791b5zvOXX35pZsyYYUJCQoyPj4/p3r27ycrKqrAP+fn5ZsyYMSYoKMj4+vqazp07mx9++OGRjoM5yhw1pv7maG2f2/Ic/13ZE+VwODRt2jTnb1/g74Y5iieFf0UNACwiwgBgkZXLEQCABzgTBgCLiDAAWESEAcAiIgwAFtXZ35grrzT/+fpaNf6hPEJ/f6Lba/npf57o9upS5DuZNQ/6m9p1sX7+IZwnoT7mKGfCAGAREQYAi4gwAFhEhAHAIiIMABYRYQCwiAgDgEVEGAAsIsIAYBERBgCLiDAAWESEAcAiIgwAFhFhALCICAOARUQYACwiwgBgEREGAIuIMABYRIQBwCIiDAAWEWEAsIgIA4BFRBgALCLCAGAREQYAi4gwAFhEhAHAIiIMABYRYQCwiAgDgEVEGAAsIsIAYBERBgCLiDAAWESEAcAiIgwAFhFhALCICAOARUQYACwiwgBgEREGAIuIMABYRIQBwCIiDAAWEWEAsIgIA4BFRBgALCLCAGAREQYAi7zqa8W9n/6f+lo1/qF+KLW9B8CTx5kwAFhEhAHAIiIMABYRYQCwiAgDgEVEGAAsIsIAYBERBgCLiDAAWESEAcAiIgwAFhFhALCICAOARUQYACwiwgBgEREGAIuIMABYRIQBwCIiDAAWEWEAsIgIA4BFRBgALCLCAGAREQYAi4gwAFhEhAHAIiIMABYRYQCwiAgDgEVEGAAsIsIAYBERBgCLiDAAWESEAcAiIgwAFhFhALCICAOARUQYACwiwgBgEREGAIuIMABYRIQBwCIiDAAWEWEAsIgIA4BFRBgALCLCAGAREQYAi4gwAFhEhAHAIiIMABZ51deK8/63S32tGgD+MTgTBgCLiDAAWESEAcAiIgwAFhFhALCICAOARUQYACwiwgBgEREGAIuIMABYRIQBwCIiDAAWEWEAsIgIA4BFRBgALCLCAGAREQYAi4gwAFhEhAHAIiIMABYRYQCwiAgDgEVEGAAsIsIAYBERBgCLiDAAWESEAcAiIgwAFhFhALCICAOARUQYACwiwgBgEREGAIuIMABYRIQBwCIiDAAWEWEAsIgIA4BFRBgALCLCAGAREQYAi4gwAFhEhAHAIiIMABYRYQCwiAgDgEVEGAAsIsIAYBERBgCLiDAAWOQwxhjbOwEA/1acCQOARUQYACwiwgBgEREGAIuIMABYRIQBwCIiDAAWEWEAsIgIA4BF/wfgH35dYY9RygAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4,2))\n",
    "plt.subplot(121)\n",
    "plt.title('Top Z-Slice BP Sub0')\n",
    "plt.pcolormesh(BP_subset0[:,:,-1].T)\n",
    "plt.axis('off')\n",
    "plt.subplot(122)\n",
    "plt.title('Top Z-Slice BP Sub1')\n",
    "plt.pcolormesh(BP_subset1[:,:,-1].T)\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "341ec6bb",
   "metadata": {
    "papermill": {
     "duration": 0.010928,
     "end_time": "2024-11-19T22:42:39.715543",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.704615",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "There are a couple other methods we need to implement for subset based image reconstruction algorithms:\n",
    "* `set_n_subsets`: Sets the number of \"subsets\" used in an iterative reconstruction algorithm. This requires partitioning the projection data into N distinct subsets. In general, the way in which this is done depends on the particular imaging modality (in SPECT/PET, one typically partitions using the projection angle).In our case, we have only set up our forward/back projectors to consider 2 subsets (we only have two possible angles). More advanced system matrices will need to consider implementation of arbitrary numbers of subsets \n",
    "* `get_projection_subset`: Returns the projection data corresponding to the ith subset \n",
    "* `get_weighting_subset`: This is used for scaling prior functions used in Bayesian reconstruction algorithms. It should be equal to the fraction of data elements in the particular subset. (For $N$ even subsets, this factor is $1/N$ for each subset, but sometimes the subsets aren't even)\n",
    "* We need to modify `compute_normalization_factor` so it can get the normalization factor of the $m$ th subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "c697007d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.738234Z",
     "iopub.status.busy": "2024-11-19T22:42:39.737994Z",
     "iopub.status.idle": "2024-11-19T22:42:39.745862Z",
     "shell.execute_reply": "2024-11-19T22:42:39.745215Z"
    },
    "papermill": {
     "duration": 0.020488,
     "end_time": "2024-11-19T22:42:39.747178",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.726690",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSSystemMatrix(SystemMatrix):\n",
    "    # ----\n",
    "    # NEW METHODS\n",
    "    # ----\n",
    "    def set_n_subsets(self, n_subsets):\n",
    "        self.n_subsets = n_subsets\n",
    "    def get_projection_subset(self, projections, subset_idx):\n",
    "        # Called when n_subsets>1 in internal pytomography code, in this case, assumes 2 subsets since thats the only possible number of subsets we have. In general, this should split data evenly (see SPECTSystemMatrix source code)\n",
    "        return projections[subset_idx].unsqueeze(0)\n",
    "    def compute_normalization_factor(self, subset_idx = None):\n",
    "        # This function generally looks the same for all system matrices\n",
    "        norm_projections = torch.ones((2,self.proj_meta.M, self.proj_meta.M))\n",
    "        if subset_idx is not None:\n",
    "            norm_projections = self.get_projection_subset(norm_projections, subset_idx)\n",
    "        return self.backward(norm_projections, subset_idx)\n",
    "    def get_weighting_subset(self, subset_idx):\n",
    "        if subset_idx is None:\n",
    "            return 1\n",
    "        elif self.n_subsets==2:\n",
    "            return 0.5 # equal weighting in this case, in general need to be careful with this\n",
    "    # ----\n",
    "    # SAME AS PREVIOUS\n",
    "    # ----\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        projection_0degrees = object.sum(dim=0)\n",
    "        projection_90degrees = object.sum(dim=0)\n",
    "        if subset_idx==0:\n",
    "            projections = projection_0degrees\n",
    "        elif subset_idx==1:\n",
    "            projections = projection_90degrees\n",
    "        else:\n",
    "            projections = torch.stack([projection_0degrees, projection_90degrees], dim=0)\n",
    "        projections *= self.proj_meta.sensitivity_factor\n",
    "        return projections\n",
    "    def backward(self, proj, subset_idx = None):\n",
    "        # Back projection expects projections in their subset\n",
    "        if subset_idx is not None:\n",
    "            object_BP = (proj[0]*self.proj_meta.sensitivity_factor).unsqueeze(subset_idx).repeat_interleave(self.proj_meta.M, subset_idx)\n",
    "        else:\n",
    "            object_BP_angle0 = (proj[0]*self.proj_meta.sensitivity_factor).unsqueeze(0).repeat(self.proj_meta.M,1,1)\n",
    "            object_BP_angle90 = (proj[1]*self.proj_meta.sensitivity_factor).unsqueeze(1).repeat(1,self.proj_meta.M,1)\n",
    "            object_BP = object_BP_angle0 + object_BP_angle90\n",
    "        return object_BP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "96b4bc58",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.772203Z",
     "iopub.status.busy": "2024-11-19T22:42:39.771521Z",
     "iopub.status.idle": "2024-11-19T22:42:39.774895Z",
     "shell.execute_reply": "2024-11-19T22:42:39.774272Z"
    },
    "papermill": {
     "duration": 0.019253,
     "end_time": "2024-11-19T22:42:39.776493",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.757240",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSSystemMatrix(object_meta=object_meta, proj_meta=proj_meta)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55228d8e",
   "metadata": {
    "papermill": {
     "duration": 0.010891,
     "end_time": "2024-11-19T22:42:39.798581",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.787690",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Our system matrix can now be used in reconstruction algorithms. Let's reconstruct our simulated data. Like the tutorials online, we now create a likelihood function and then a reconstruction algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "47d21aef",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.821845Z",
     "iopub.status.busy": "2024-11-19T22:42:39.821100Z",
     "iopub.status.idle": "2024-11-19T22:42:39.824782Z",
     "shell.execute_reply": "2024-11-19T22:42:39.824180Z"
    },
    "papermill": {
     "duration": 0.016297,
     "end_time": "2024-11-19T22:42:39.826047",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.809750",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample_object = torch.rand(object_meta.shape) # object has batch dimension\n",
    "sample_projections = system_matrix.forward(sample_object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9d4bfbea",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.848950Z",
     "iopub.status.busy": "2024-11-19T22:42:39.848181Z",
     "iopub.status.idle": "2024-11-19T22:42:39.851890Z",
     "shell.execute_reply": "2024-11-19T22:42:39.851347Z"
    },
    "papermill": {
     "duration": 0.016066,
     "end_time": "2024-11-19T22:42:39.853201",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.837135",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSSystemMatrix(object_meta=object_meta, proj_meta=proj_meta)\n",
    "likelihood = NegativeMSELikelihood(system_matrix, projections=sample_projections, scaling_constant=0.01)\n",
    "reconstruction_algorithm = OSEM(likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "fc7c44b2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.878293Z",
     "iopub.status.busy": "2024-11-19T22:42:39.876162Z",
     "iopub.status.idle": "2024-11-19T22:42:39.896956Z",
     "shell.execute_reply": "2024-11-19T22:42:39.896225Z"
    },
    "papermill": {
     "duration": 0.034393,
     "end_time": "2024-11-19T22:42:39.898351",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.863958",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "recon = reconstruction_algorithm(n_iters=40, n_subsets=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "a2b79907",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:39.925444Z",
     "iopub.status.busy": "2024-11-19T22:42:39.923635Z",
     "iopub.status.idle": "2024-11-19T22:42:40.090782Z",
     "shell.execute_reply": "2024-11-19T22:42:40.090000Z"
    },
    "papermill": {
     "duration": 0.1815,
     "end_time": "2024-11-19T22:42:40.092182",
     "exception": false,
     "start_time": "2024-11-19T22:42:39.910682",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7fb3c044fa90>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAADLCAYAAAC/DyLrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAMVElEQVR4nO3dX2xUdRqH8e8wbAeEzqxVy9r0zzYhoFhLIoVkDBrU2KSbJbJXXmxIE/EC0zYhvVmKFw0kbk1IjGzQBoNBb0iJwQJGaWwibUXSLGVpasguCYZsuy5I1O20He10W357sduG2go9p3P69myfT3IuZnKG8yovj+N0cog455wAACaWWQ8AAEsZEQYAQ0QYAAwRYQAwRIQBwBARBgBDRBgADBFhADBEhAHAEBEGAEOeItzc3Kzy8nLF43HF43Elk0mdPXs2qNkAz9hRhE3Ey70jPvroI0WjUa1du1aS9P777+vgwYO6fPmyHnvsscCGBOaKHUXYeIrwbPLy8nTw4EHt2rUrWzMBWcWOYjFb7veFExMT+uCDD5ROp5VMJn/2vEwmo0wmM/X49u3b+v777/XAAw8oEon4vTyWEOechoeHVVBQoGXL5v4JGjuKheJ3Rydf7ElfX59btWqVi0ajLpFIuI8//viu5zc2NjpJHBzzPgYGBthRjkV9zHVH7+T544ixsTH19/drcHBQJ0+e1NGjR9XZ2akNGzbMev5P32WkUikVFxfr73/5teKrw/fljN+te9x6BN+u/WmT9Qi+3P5xVP/8wx81ODioRCJxz/OztaPr3t2j6H2xrP1zLJRfHQzfn6tJtzbnWo/gy8TYqP527MCcd/ROnj+OyMnJmfqhR0VFhS5evKhDhw7pyJEjs54fi8UUi81c5PjqZYrnRr1e3tzyyC+sR/Bt2coV1iPMy1w/GsjWjkbvi4Uywsuj4ftzNSmaszR29E7z/k+mc27auwhgsWFHsZh5eie8b98+VVVVqaioSMPDw2ppaVFHR4fa2tqCmg/whB1F2HiK8DfffKOdO3fqxo0bSiQSKi8vV1tbm55//vmg5gM8YUcRNp4i/O677wY1B5AV7CjCJrw/RgWA/wNEGAAMEWEAMESEAcAQEQYAQ0QYAAwRYQAwRIQBwBARBgBDRBgADBFhADBEhAHAEBEGAENEGAAMEWEAMESEAcAQEQYAQ0QYAAwRYQAwRIQBwBARBgBDRBgADBFhADBEhAHAEBEGAENEGAAMEWEAMESEAcAQEQYAQ0QYAAwRYQAwRIQBwBARBgBDRBgADBFhADBEhAHAEBEGAENEGAAMEWEAMESEAcCQpwg3NTVp8+bNys3NVX5+vnbs2KGrV68GNRvgGTuKsPEU4c7OTtXU1Ki7u1vt7e0aHx9XZWWl0ul0UPMBnrCjCJvlXk5ua2ub9vjYsWPKz8/XpUuX9PTTT2d1MMAPdhRhM6/PhFOplCQpLy8vK8MA2caOYrHz9E74Ts451dfXa+vWrSorK/vZ8zKZjDKZzNTjoaEhv5cEPGFHEQa+I1xbW6u+vj6dP3/+ruc1NTVp//79M57/7V9/o+WrYn4vb+aHnQXWI/i2/pU/W4/gy7j7t/7h43Xz3dHb3QlFYit8XNmW67lgPYJva5Y9bj2CL+Pjo7ri87W+Po6oq6vTmTNndO7cORUWFt713IaGBqVSqaljYGDA16CAF+wowsLTO2HnnOrq6tTa2qqOjg6Vlpbe8zWxWEyxWPje8SKc2FGEjacI19TU6Pjx4zp9+rRyc3N18+ZNSVIikdDKlSsDGRDwgh1F2Hj6OKK5uVmpVErbtm3Tww8/PHWcOHEiqPkAT9hRhI3njyOAxYwdRdhw7wgAMESEAcAQEQYAQ0QYAAwRYQAwRIQBwBARBgBDRBgADBFhADBEhAHAEBEGAENEGAAMEWEAMESEAcAQEQYAQ0QYAAwRYQAwRIQBwBARBgBDRBgADBFhADBEhAHAEBEGAENEGAAMEWEAMESEAcAQEQYAQ0QYAAwRYQAwRIQBwBARBgBDRBgADBFhADBEhAHAEBEGAENEGAAMEWEAMESEAcAQEQYAQ0QYAAwRYQAw5DnCXV1d2r59uwoKChSJRHTq1KkAxgL8Y0cRJp4jnE6ntXHjRh0+fDiIeYB5Y0cRJsu9vqCqqkpVVVVBzAJkBTuKMPEcYa8ymYwymczU46GhoaAvCXjCjsJS4BFuamrS/v37Zzx/YzCu6NiKoC+fdb+0HmAeKvv+ZT2CL6Mj4zq3Jbhf/+d2tLvmqOK50eAuHJCNE69Yj+DbxNaU9Qi+TPzgpN/7e23g345oaGhQKpWaOgYGBoK+JOAJOwpLgb8TjsViisViQV8G8I0dhSW+JwwAhjy/Ex4ZGdG1a9emHl+/fl29vb3Ky8tTcXFxVocD/GBHESaeI9zT06Nnnnlm6nF9fb0kqbq6Wu+9917WBgP8YkcRJp4jvG3bNjnngpgFyAp2FGHCZ8IAYIgIA4AhIgwAhogwABgiwgBgiAgDgCEiDACGiDAAGCLCAGCICAOAISIMAIaIMAAYIsIAYIgIA4AhIgwAhogwABgiwgBgiAgDgCEiDACGiDAAGCLCAGCICAOAISIMAIaIMAAYIsIAYIgIA4AhIgwAhogwABgiwgBgiAgDgCEiDACGiDAAGCLCAGCICAOAISIMAIaIMAAYIsIAYIgIA4AhIgwAhogwABjyFeG3335bpaWlWrFihTZt2qTPP/8823MB88KOIiw8R/jEiRPas2ePXn31VV2+fFlPPfWUqqqq1N/fH8R8gGfsKMLEc4TfeOMN7dq1Sy+//LIeffRRvfnmmyoqKlJzc3MQ8wGesaMIk+VeTh4bG9OlS5e0d+/eac9XVlbqwoULs74mk8kok8lMPU6lUpKk2z9mZj1/sZsYi1qP4NvoyLj1CL5Mzu2cu+e52dzRoZHbfkc2NZEZtR7Bt4kfQtqF/809lx2dwXnw9ddfO0nuiy++mPb8a6+95tatWzfraxobG50kDo55H1999RU7yrGoj7ns6E95eic8KRKJTHvsnJvx3KSGhgbV19dPPR4cHFRJSYn6+/uVSCT8XN7M0NCQioqKNDAwoHg8bj2OJ2GePZVKqbi4WHl5eXN+DTsavt/nMM/uZ0cneYrwgw8+qGg0qps3b057/tatW1qzZs2sr4nFYorFYjOeTyQSofsXPSkejzO7gWXL7v0jDHb0v8L8+xzm2eeyozNe4+XknJwcbdq0Se3t7dOeb29v15NPPun54kC2saMIG88fR9TX12vnzp2qqKhQMpnUO++8o/7+fu3evTuI+QDP2FGEiecIv/jii/ruu+904MAB3bhxQ2VlZfrkk09UUlIyp9fHYjE1NjbO+r9/ix2z2/A6OzvK7AttPrNHnPPznQoAQDZw7wgAMESEAcAQEQYAQ0QYAAwtaITDenvBrq4ubd++XQUFBYpEIjp16pT1SHPS1NSkzZs3Kzc3V/n5+dqxY4euXr1qPdacNDc3q7y8fOqL+8lkUmfPng38uuzowmJHFzDCYb69YDqd1saNG3X48GHrUTzp7OxUTU2Nuru71d7ervHxcVVWViqdTluPdk+FhYV6/fXX1dPTo56eHj377LN64YUXdOXKlcCuyY4uPHZU8nQDn/nYsmWL271797TnHnnkEbd3796FGiErJLnW1lbrMXy5deuWk+Q6OzutR/Hl/vvvd0ePHg3s12dH7S3FHV2Qd8KTtxesrKyc9vzdbi+I7Ju8RaOfm4xYmpiYUEtLi9LptJLJZCDXYEcXh6W4o77uoubVt99+q4mJiRk3UFmzZs2MG60gGM451dfXa+vWrSorK7MeZ06+/PJLJZNJjY6OavXq1WptbdWGDRsCuRY7am+p7uiCRHiSl9sLIrtqa2vV19en8+fPW48yZ+vXr1dvb68GBwd18uRJVVdXq7OzM7AQS+yopaW6owsSYT+3F0T21NXV6cyZM+rq6lJhYaH1OHOWk5OjtWvXSpIqKip08eJFHTp0SEeOHMn6tdhRW0t5RxfkM2FuL2jDOafa2lp9+OGH+uyzz1RaWmo90rw456b9NUTZxI7aYEcX8OOIMN9ecGRkRNeuXZt6fP36dfX29iovL0/FxcWGk91dTU2Njh8/rtOnTys3N3fqXV4ikdDKlSuNp7u7ffv2qaqqSkVFRRoeHlZLS4s6OjrU1tYW2DXZ0YXHjmrhvqLmnHNvvfWWKykpcTk5Oe6JJ54IzddQzp07N+vfJ1VdXW092l3NNrMkd+zYMevR7umll16a2pWHHnrIPffcc+7TTz8N/Lrs6MJiR53jVpYAYIh7RwCAISIMAIaIMAAYIsIAYIgIA4AhIgwAhogwABgiwgBgiAgDgCEiDACGiDAAGCLCAGDoP7rQYhqe9GPjAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 400x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4,2))\n",
    "plt.subplot(121)\n",
    "plt.pcolormesh(sample_object[:,:,-1])\n",
    "plt.subplot(122)\n",
    "plt.pcolormesh(recon[:,:,-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33e97a20",
   "metadata": {
    "papermill": {
     "duration": 0.011359,
     "end_time": "2024-11-19T22:42:40.115325",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.103966",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Example 2: List Mode System Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb66c947",
   "metadata": {
    "papermill": {
     "duration": 0.010276,
     "end_time": "2024-11-19T22:42:40.137314",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.127038",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Modalities containing listmode data can also be implemented in PyTomography; their implementation may look slightly different than above, however. Regardless, if implemented properly, they are still usable with all available reconstruction algorithms\n",
    "\n",
    "For our simple imaging system, listmode data might be stored via a list of integers where each integer specifies a given pixel in projection space. Since we have 2 projections, each with $M \\times M$ pixels, the total number of detector elements is $2M^2$.\n",
    "\n",
    "We can simulate 100 random events as follows (in reality, a more sophisticated program, such as GATE, would be required to give a proper distribution according to a specific phantom/source distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "49218215",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.160154Z",
     "iopub.status.busy": "2024-11-19T22:42:40.159260Z",
     "iopub.status.idle": "2024-11-19T22:42:40.167831Z",
     "shell.execute_reply": "2024-11-19T22:42:40.167088Z"
    },
    "papermill": {
     "duration": 0.021004,
     "end_time": "2024-11-19T22:42:40.169176",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.148172",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([ 3,  3, 10, 14, 15, 14,  6, 14,  1, 16,  7, 14, 12,  2, 13, 13, 10,  6,\n",
       "        10,  5,  6,  1, 15,  9,  6,  2,  0, 17,  9,  2,  8,  9,  5, 10,  0,  1,\n",
       "         7, 16,  0, 11, 13,  5,  4,  9, 11, 17, 10,  9,  8, 10,  1, 13,  5, 17,\n",
       "         6,  2,  2, 10, 16, 16,  3, 10, 10, 10,  0,  9, 14, 10, 10,  4, 10, 14,\n",
       "        13, 13, 12,  1,  7, 13, 12,  1,  4, 16, 11,  6,  0, 17,  8, 15,  5, 12,\n",
       "         3, 16,  5,  8,  0,  9,  4,  3,  9,  6,  8, 17, 15,  5,  8,  5, 10, 12,\n",
       "         2,  1,  5,  7, 13,  8, 13, 17, 14, 12,  6,  0,  3, 15,  9,  7, 11,  0,\n",
       "         6,  4,  2, 11, 17,  9, 10,  6,  0,  1, 16, 16, 16,  2, 16,  7, 15, 12,\n",
       "        10, 17, 14, 14, 17,  1,  7,  9,  2, 10,  7, 12,  3, 14, 15, 14,  5,  7,\n",
       "        17,  8,  4, 17,  8,  1,  7, 16,  5, 16,  1, 17, 16, 14,  4, 11, 13,  6,\n",
       "        13,  1,  6,  3, 16,  0,  3, 11,  6,  8, 16,  6, 12,  7, 15, 10,  4, 13,\n",
       "        13,  2,  0,  7,  8, 10,  1, 10,  3, 13,  1,  9, 15, 16, 14,  9,  3,  0,\n",
       "         4,  4,  2, 10,  1, 13, 14, 17,  5, 12,  3, 17,  6,  3, 17,  9,  3,  7,\n",
       "        14,  5, 13, 13,  2, 13,  7,  4,  6, 13, 10,  9, 12,  9, 13, 10,  5, 13,\n",
       "        17,  0, 11,  4,  9,  3,  1, 17,  0, 13,  3, 11, 17, 10,  9, 15,  6, 13,\n",
       "        14, 13, 10,  5,  4,  9,  4, 10, 12,  9,  0,  3, 10, 13, 15,  3, 10,  3,\n",
       "         9, 17,  7,  6, 11,  2, 15, 12,  4,  6,  1, 10, 16, 16,  9, 17,  7, 11,\n",
       "         3, 11,  8,  3,  2, 14, 13,  7,  2,  7,  3, 15, 12,  2, 16, 16, 10,  2,\n",
       "         1,  4,  0, 10,  2, 16, 17,  6,  7, 10,  9,  1,  0,  7,  3,  0, 15,  9,\n",
       "         3, 10, 13,  6,  9,  1, 16, 10,  7,  6,  5, 12,  1,  0,  7,  1,  3,  4,\n",
       "        16, 12, 10,  5, 16,  3,  0, 11,  8, 16,  6,  7,  0, 11,  8, 12,  6, 10,\n",
       "        10,  7,  1, 14, 14, 13,  0,  3, 15,  4, 16,  6, 10,  9,  4, 13,  7, 16,\n",
       "         0, 14,  6,  2])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "detector_ids = torch.randint(low=0, high=2*M**2, size=(400,))\n",
    "detector_ids"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73f224b9",
   "metadata": {
    "papermill": {
     "duration": 0.010574,
     "end_time": "2024-11-19T22:42:40.190261",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.179687",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## MetaData"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4d87c3f",
   "metadata": {
    "papermill": {
     "duration": 0.010408,
     "end_time": "2024-11-19T22:42:40.211232",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.200824",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Each integer corresponds to a given detector element where a count has been measured.First, we may want to generate a lookup table for these integer values and corresponding detector coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "bd4f61ce",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.233944Z",
     "iopub.status.busy": "2024-11-19T22:42:40.233259Z",
     "iopub.status.idle": "2024-11-19T22:42:40.239252Z",
     "shell.execute_reply": "2024-11-19T22:42:40.238603Z"
    },
    "papermill": {
     "duration": 0.018948,
     "end_time": "2024-11-19T22:42:40.240639",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.221691",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[0, 0, 0],\n",
       "        [0, 0, 1],\n",
       "        [0, 0, 2],\n",
       "        [0, 1, 0],\n",
       "        [0, 1, 1],\n",
       "        [0, 1, 2],\n",
       "        [0, 2, 0],\n",
       "        [0, 2, 1],\n",
       "        [0, 2, 2],\n",
       "        [1, 0, 0],\n",
       "        [1, 0, 1],\n",
       "        [1, 0, 2],\n",
       "        [1, 1, 0],\n",
       "        [1, 1, 1],\n",
       "        [1, 1, 2],\n",
       "        [1, 2, 0],\n",
       "        [1, 2, 1],\n",
       "        [1, 2, 2]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scanner_LUT = torch.cartesian_prod(\n",
    "    torch.tensor([0,1]), # Angle\n",
    "    torch.arange(M), # row\n",
    "    torch.arange(M), # column\n",
    ")\n",
    "scanner_LUT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbc876df",
   "metadata": {
    "papermill": {
     "duration": 0.010483,
     "end_time": "2024-11-19T22:42:40.262475",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.251992",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "This lookup table gives the angle index, row index, and column index in projection space for each detector ID. For example, if we want the location corresponding to detector ID 11:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "6ae5406e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.285368Z",
     "iopub.status.busy": "2024-11-19T22:42:40.284980Z",
     "iopub.status.idle": "2024-11-19T22:42:40.289320Z",
     "shell.execute_reply": "2024-11-19T22:42:40.288791Z"
    },
    "papermill": {
     "duration": 0.017558,
     "end_time": "2024-11-19T22:42:40.290574",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.273016",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([0, 1, 2])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scanner_LUT[5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09dcf616",
   "metadata": {
    "papermill": {
     "duration": 0.011268,
     "end_time": "2024-11-19T22:42:40.312229",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.300961",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Similarily, if we wanted the location corresponding to the first 5 measured events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "5e4d9228",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.334899Z",
     "iopub.status.busy": "2024-11-19T22:42:40.334299Z",
     "iopub.status.idle": "2024-11-19T22:42:40.339176Z",
     "shell.execute_reply": "2024-11-19T22:42:40.338652Z"
    },
    "papermill": {
     "duration": 0.017621,
     "end_time": "2024-11-19T22:42:40.340466",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.322845",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[0, 1, 0],\n",
       "        [0, 1, 0],\n",
       "        [1, 0, 1],\n",
       "        [1, 1, 2],\n",
       "        [1, 2, 0]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scanner_LUT[detector_ids[0:5]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "240aa164",
   "metadata": {
    "papermill": {
     "duration": 0.010757,
     "end_time": "2024-11-19T22:42:40.363897",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.353140",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We can also get the sensitivity factor for each listmode event:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "dd743f9e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.390809Z",
     "iopub.status.busy": "2024-11-19T22:42:40.390350Z",
     "iopub.status.idle": "2024-11-19T22:42:40.399955Z",
     "shell.execute_reply": "2024-11-19T22:42:40.399405Z"
    },
    "papermill": {
     "duration": 0.026466,
     "end_time": "2024-11-19T22:42:40.401173",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.374707",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([1.0203, 1.0203, 1.2767, 1.2802, 1.2021, 1.2802, 1.2276, 1.2802, 1.2862,\n",
       "        1.2021, 1.2276, 1.2802, 1.2802, 1.2862, 1.2802, 1.2802, 1.2767, 1.2276,\n",
       "        1.2767, 1.0203, 1.2276, 1.2862, 1.2021, 1.2767, 1.2276, 1.2862, 1.2862,\n",
       "        1.2021, 1.2767, 1.2862, 1.2276, 1.2767, 1.0203, 1.2767, 1.2862, 1.2862,\n",
       "        1.2276, 1.2021, 1.2862, 1.2767, 1.2802, 1.0203, 1.0203, 1.2767, 1.2767,\n",
       "        1.2021, 1.2767, 1.2767, 1.2276, 1.2767, 1.2862, 1.2802, 1.0203, 1.2021,\n",
       "        1.2276, 1.2862, 1.2862, 1.2767, 1.2021, 1.2021, 1.0203, 1.2767, 1.2767,\n",
       "        1.2767, 1.2862, 1.2767, 1.2802, 1.2767, 1.2767, 1.0203, 1.2767, 1.2802,\n",
       "        1.2802, 1.2802, 1.2802, 1.2862, 1.2276, 1.2802, 1.2802, 1.2862, 1.0203,\n",
       "        1.2021, 1.2767, 1.2276, 1.2862, 1.2021, 1.2276, 1.2021, 1.0203, 1.2802,\n",
       "        1.0203, 1.2021, 1.0203, 1.2276, 1.2862, 1.2767, 1.0203, 1.0203, 1.2767,\n",
       "        1.2276, 1.2276, 1.2021, 1.2021, 1.0203, 1.2276, 1.0203, 1.2767, 1.2802,\n",
       "        1.2862, 1.2862, 1.0203, 1.2276, 1.2802, 1.2276, 1.2802, 1.2021, 1.2802,\n",
       "        1.2802, 1.2276, 1.2862, 1.0203, 1.2021, 1.2767, 1.2276, 1.2767, 1.2862,\n",
       "        1.2276, 1.0203, 1.2862, 1.2767, 1.2021, 1.2767, 1.2767, 1.2276, 1.2862,\n",
       "        1.2862, 1.2021, 1.2021, 1.2021, 1.2862, 1.2021, 1.2276, 1.2021, 1.2802,\n",
       "        1.2767, 1.2021, 1.2802, 1.2802, 1.2021, 1.2862, 1.2276, 1.2767, 1.2862,\n",
       "        1.2767, 1.2276, 1.2802, 1.0203, 1.2802, 1.2021, 1.2802, 1.0203, 1.2276,\n",
       "        1.2021, 1.2276, 1.0203, 1.2021, 1.2276, 1.2862, 1.2276, 1.2021, 1.0203,\n",
       "        1.2021, 1.2862, 1.2021, 1.2021, 1.2802, 1.0203, 1.2767, 1.2802, 1.2276,\n",
       "        1.2802, 1.2862, 1.2276, 1.0203, 1.2021, 1.2862, 1.0203, 1.2767, 1.2276,\n",
       "        1.2276, 1.2021, 1.2276, 1.2802, 1.2276, 1.2021, 1.2767, 1.0203, 1.2802,\n",
       "        1.2802, 1.2862, 1.2862, 1.2276, 1.2276, 1.2767, 1.2862, 1.2767, 1.0203,\n",
       "        1.2802, 1.2862, 1.2767, 1.2021, 1.2021, 1.2802, 1.2767, 1.0203, 1.2862,\n",
       "        1.0203, 1.0203, 1.2862, 1.2767, 1.2862, 1.2802, 1.2802, 1.2021, 1.0203,\n",
       "        1.2802, 1.0203, 1.2021, 1.2276, 1.0203, 1.2021, 1.2767, 1.0203, 1.2276,\n",
       "        1.2802, 1.0203, 1.2802, 1.2802, 1.2862, 1.2802, 1.2276, 1.0203, 1.2276,\n",
       "        1.2802, 1.2767, 1.2767, 1.2802, 1.2767, 1.2802, 1.2767, 1.0203, 1.2802,\n",
       "        1.2021, 1.2862, 1.2767, 1.0203, 1.2767, 1.0203, 1.2862, 1.2021, 1.2862,\n",
       "        1.2802, 1.0203, 1.2767, 1.2021, 1.2767, 1.2767, 1.2021, 1.2276, 1.2802,\n",
       "        1.2802, 1.2802, 1.2767, 1.0203, 1.0203, 1.2767, 1.0203, 1.2767, 1.2802,\n",
       "        1.2767, 1.2862, 1.0203, 1.2767, 1.2802, 1.2021, 1.0203, 1.2767, 1.0203,\n",
       "        1.2767, 1.2021, 1.2276, 1.2276, 1.2767, 1.2862, 1.2021, 1.2802, 1.0203,\n",
       "        1.2276, 1.2862, 1.2767, 1.2021, 1.2021, 1.2767, 1.2021, 1.2276, 1.2767,\n",
       "        1.0203, 1.2767, 1.2276, 1.0203, 1.2862, 1.2802, 1.2802, 1.2276, 1.2862,\n",
       "        1.2276, 1.0203, 1.2021, 1.2802, 1.2862, 1.2021, 1.2021, 1.2767, 1.2862,\n",
       "        1.2862, 1.0203, 1.2862, 1.2767, 1.2862, 1.2021, 1.2021, 1.2276, 1.2276,\n",
       "        1.2767, 1.2767, 1.2862, 1.2862, 1.2276, 1.0203, 1.2862, 1.2021, 1.2767,\n",
       "        1.0203, 1.2767, 1.2802, 1.2276, 1.2767, 1.2862, 1.2021, 1.2767, 1.2276,\n",
       "        1.2276, 1.0203, 1.2802, 1.2862, 1.2862, 1.2276, 1.2862, 1.0203, 1.0203,\n",
       "        1.2021, 1.2802, 1.2767, 1.0203, 1.2021, 1.0203, 1.2862, 1.2767, 1.2276,\n",
       "        1.2021, 1.2276, 1.2276, 1.2862, 1.2767, 1.2276, 1.2802, 1.2276, 1.2767,\n",
       "        1.2767, 1.2276, 1.2862, 1.2802, 1.2802, 1.2802, 1.2862, 1.0203, 1.2021,\n",
       "        1.0203, 1.2021, 1.2276, 1.2767, 1.2767, 1.0203, 1.2802, 1.2276, 1.2021,\n",
       "        1.2862, 1.2802, 1.2276, 1.2862])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sensitivity_factor[*scanner_LUT[detector_ids][:,:2].T]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5ccca9f",
   "metadata": {
    "papermill": {
     "duration": 0.010595,
     "end_time": "2024-11-19T22:42:40.423287",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.412692",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "With this in mind, we can define our new listmode projection metadata class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "412ad17d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.448302Z",
     "iopub.status.busy": "2024-11-19T22:42:40.447748Z",
     "iopub.status.idle": "2024-11-19T22:42:40.452609Z",
     "shell.execute_reply": "2024-11-19T22:42:40.451963Z"
    },
    "papermill": {
     "duration": 0.019801,
     "end_time": "2024-11-19T22:42:40.454048",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.434247",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSListmodeProjMeta(ProjMeta):\n",
    "    def __init__(self, shape, scanner_LUT, detector_ids, sensitivity_factor):\n",
    "        self.scanner_LUT = scanner_LUT\n",
    "        self.detector_ids = detector_ids\n",
    "        self.sensitivity_at_ids = sensitivity_factor[*scanner_LUT[:,1:].T]\n",
    "        self.shape = shape\n",
    "        if (sensitivity_factor.shape[0]!=M)*(sensitivity_factor.shape[1]!=M):\n",
    "            raise ValueError(\"sensitivity_factor should have side dimensions M\")\n",
    "proj_meta_listmode = EXSListmodeProjMeta((2,M,M), scanner_LUT, detector_ids, sensitivity_factor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ecf19f2",
   "metadata": {
    "papermill": {
     "duration": 0.011441,
     "end_time": "2024-11-19T22:42:40.476831",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.465390",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Lets now implement the forward/backward methods of the system matrix. We'll start with forward projection, which forward projects an object to each `id` form the listmode data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "49fe35e1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.500662Z",
     "iopub.status.busy": "2024-11-19T22:42:40.500213Z",
     "iopub.status.idle": "2024-11-19T22:42:40.505147Z",
     "shell.execute_reply": "2024-11-19T22:42:40.504537Z"
    },
    "papermill": {
     "duration": 0.018497,
     "end_time": "2024-11-19T22:42:40.506441",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.487944",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSListmodeSystemMatrix(SystemMatrix):\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        # There is probably a faster implementation, but I am trying to keep it simple for illustration purposes\n",
    "        projections = []\n",
    "        for i, detector_id in enumerate(self.proj_meta.detector_ids):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            sensitivity_factor_i = self.proj_meta.sensitivity_at_ids[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                projections.append(object[:,coord[1],coord[2]].sum() *  sensitivity_factor_i) # sum along x\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                projections.append(object[coord[1],:,coord[2]].sum() * sensitivity_factor_i)  # sum along y\n",
    "        return torch.tensor(projections)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "15285752",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.530694Z",
     "iopub.status.busy": "2024-11-19T22:42:40.530132Z",
     "iopub.status.idle": "2024-11-19T22:42:40.561530Z",
     "shell.execute_reply": "2024-11-19T22:42:40.561025Z"
    },
    "papermill": {
     "duration": 0.045043,
     "end_time": "2024-11-19T22:42:40.562732",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.517689",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([2.4181, 2.4181, 1.4410, 1.3367, 1.2952, 1.3367, 1.6440, 1.3367, 1.5810,\n",
       "        1.3257, 0.8895, 1.3367, 1.7677, 0.9422, 1.6551, 1.6551, 1.4410, 1.6440,\n",
       "        1.4410, 1.1332, 1.6440, 1.5810, 1.2952, 3.0902, 1.6440, 0.9422, 2.0810,\n",
       "        1.5294, 3.0902, 0.9422, 1.6234, 3.0902, 1.1332, 1.4410, 2.0810, 1.5810,\n",
       "        0.8895, 1.3257, 2.0810, 0.8263, 1.6551, 1.1332, 1.9618, 3.0902, 0.8263,\n",
       "        1.5294, 1.4410, 3.0902, 1.6234, 1.4410, 1.5810, 1.6551, 1.1332, 1.5294,\n",
       "        1.6440, 0.9422, 0.9422, 1.4410, 1.3257, 1.3257, 2.4181, 1.4410, 1.4410,\n",
       "        1.4410, 2.0810, 3.0902, 1.3367, 1.4410, 1.4410, 1.9618, 1.4410, 1.3367,\n",
       "        1.6551, 1.6551, 1.7677, 1.5810, 0.8895, 1.6551, 1.7677, 1.5810, 1.9618,\n",
       "        1.3257, 0.8263, 1.6440, 2.0810, 1.5294, 1.6234, 1.2952, 1.1332, 1.7677,\n",
       "        2.4181, 1.3257, 1.1332, 1.6234, 2.0810, 3.0902, 1.9618, 2.4181, 3.0902,\n",
       "        1.6440, 1.6234, 1.5294, 1.2952, 1.1332, 1.6234, 1.1332, 1.4410, 1.7677,\n",
       "        0.9422, 1.5810, 1.1332, 0.8895, 1.6551, 1.6234, 1.6551, 1.5294, 1.3367,\n",
       "        1.7677, 1.6440, 2.0810, 2.4181, 1.2952, 3.0902, 0.8895, 0.8263, 2.0810,\n",
       "        1.6440, 1.9618, 0.9422, 0.8263, 1.5294, 3.0902, 1.4410, 1.6440, 2.0810,\n",
       "        1.5810, 1.3257, 1.3257, 1.3257, 0.9422, 1.3257, 0.8895, 1.2952, 1.7677,\n",
       "        1.4410, 1.5294, 1.3367, 1.3367, 1.5294, 1.5810, 0.8895, 3.0902, 0.9422,\n",
       "        1.4410, 0.8895, 1.7677, 2.4181, 1.3367, 1.2952, 1.3367, 1.1332, 0.8895,\n",
       "        1.5294, 1.6234, 1.9618, 1.5294, 1.6234, 1.5810, 0.8895, 1.3257, 1.1332,\n",
       "        1.3257, 1.5810, 1.5294, 1.3257, 1.3367, 1.9618, 0.8263, 1.6551, 1.6440,\n",
       "        1.6551, 1.5810, 1.6440, 2.4181, 1.3257, 2.0810, 2.4181, 0.8263, 1.6440,\n",
       "        1.6234, 1.3257, 1.6440, 1.7677, 0.8895, 1.2952, 1.4410, 1.9618, 1.6551,\n",
       "        1.6551, 0.9422, 2.0810, 0.8895, 1.6234, 1.4410, 1.5810, 1.4410, 2.4181,\n",
       "        1.6551, 1.5810, 3.0902, 1.2952, 1.3257, 1.3367, 3.0902, 2.4181, 2.0810,\n",
       "        1.9618, 1.9618, 0.9422, 1.4410, 1.5810, 1.6551, 1.3367, 1.5294, 1.1332,\n",
       "        1.7677, 2.4181, 1.5294, 1.6440, 2.4181, 1.5294, 3.0902, 2.4181, 0.8895,\n",
       "        1.3367, 1.1332, 1.6551, 1.6551, 0.9422, 1.6551, 0.8895, 1.9618, 1.6440,\n",
       "        1.6551, 1.4410, 3.0902, 1.7677, 3.0902, 1.6551, 1.4410, 1.1332, 1.6551,\n",
       "        1.5294, 2.0810, 0.8263, 1.9618, 3.0902, 2.4181, 1.5810, 1.5294, 2.0810,\n",
       "        1.6551, 2.4181, 0.8263, 1.5294, 1.4410, 3.0902, 1.2952, 1.6440, 1.6551,\n",
       "        1.3367, 1.6551, 1.4410, 1.1332, 1.9618, 3.0902, 1.9618, 1.4410, 1.7677,\n",
       "        3.0902, 2.0810, 2.4181, 1.4410, 1.6551, 1.2952, 2.4181, 1.4410, 2.4181,\n",
       "        3.0902, 1.5294, 0.8895, 1.6440, 0.8263, 0.9422, 1.2952, 1.7677, 1.9618,\n",
       "        1.6440, 1.5810, 1.4410, 1.3257, 1.3257, 3.0902, 1.5294, 0.8895, 0.8263,\n",
       "        2.4181, 0.8263, 1.6234, 2.4181, 0.9422, 1.3367, 1.6551, 0.8895, 0.9422,\n",
       "        0.8895, 2.4181, 1.2952, 1.7677, 0.9422, 1.3257, 1.3257, 1.4410, 0.9422,\n",
       "        1.5810, 1.9618, 2.0810, 1.4410, 0.9422, 1.3257, 1.5294, 1.6440, 0.8895,\n",
       "        1.4410, 3.0902, 1.5810, 2.0810, 0.8895, 2.4181, 2.0810, 1.2952, 3.0902,\n",
       "        2.4181, 1.4410, 1.6551, 1.6440, 3.0902, 1.5810, 1.3257, 1.4410, 0.8895,\n",
       "        1.6440, 1.1332, 1.7677, 1.5810, 2.0810, 0.8895, 1.5810, 2.4181, 1.9618,\n",
       "        1.3257, 1.7677, 1.4410, 1.1332, 1.3257, 2.4181, 2.0810, 0.8263, 1.6234,\n",
       "        1.3257, 1.6440, 0.8895, 2.0810, 0.8263, 1.6234, 1.7677, 1.6440, 1.4410,\n",
       "        1.4410, 0.8895, 1.5810, 1.3367, 1.3367, 1.6551, 2.0810, 2.4181, 1.2952,\n",
       "        1.9618, 1.3257, 1.6440, 1.4410, 3.0902, 1.9618, 1.6551, 0.8895, 1.3257,\n",
       "        2.0810, 1.3367, 1.6440, 0.9422])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system_matrix = EXSListmodeSystemMatrix(object_meta=object_meta, proj_meta=proj_meta_listmode)\n",
    "sample_object = torch.rand(object_meta.shape) # object has batch dimension\n",
    "sample_projections = system_matrix.forward(sample_object)\n",
    "sample_projections"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aa08b40",
   "metadata": {
    "papermill": {
     "duration": 0.012837,
     "end_time": "2024-11-19T22:42:40.589418",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.576581",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Note that our \"projections\" are now in listmode format: namely, these correspond to forward projected values at all the \"detector_ids\" which formed our listmode dataset. The backward method is analogous:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "43e643fa",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.612970Z",
     "iopub.status.busy": "2024-11-19T22:42:40.612596Z",
     "iopub.status.idle": "2024-11-19T22:42:40.619728Z",
     "shell.execute_reply": "2024-11-19T22:42:40.619074Z"
    },
    "papermill": {
     "duration": 0.020883,
     "end_time": "2024-11-19T22:42:40.621604",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.600721",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSListmodeSystemMatrix(SystemMatrix):\n",
    "    # ----\n",
    "    # SAME AS ABOVE\n",
    "    # ----\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        # There is probably a faster implementation, but I am trying to keep it simple for illustration purposes\n",
    "        projections = []\n",
    "        for i, detector_id in enumerate(self.proj_meta.detector_ids):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            sensitivity_factor_i = self.proj_meta.sensitivity_at_ids[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                projections.append(object[:,coord[1],coord[2]].sum() *  sensitivity_factor_i) # sum along x\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                projections.append(object[coord[1],:,coord[2]].sum() * sensitivity_factor_i)  # sum along y\n",
    "        return torch.tensor(projections)\n",
    "    # ---\n",
    "    # NEW CODE\n",
    "    # ---\n",
    "    def backward(self, projections, subset_idx = None):\n",
    "        # There is probably a faster implementation, but I am trying to keep it simple for illustration purposes\n",
    "        object = torch.zeros(object_meta.shape)\n",
    "        projections *= self.proj_meta.sensitivity_at_ids[self.proj_meta.detector_ids]\n",
    "        for i, detector_id in enumerate(detector_ids):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                object[:,coord[1],coord[2]] += projections[i]\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                object[coord[1],:,coord[2]] += projections[i]\n",
    "        return object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "0cf8e99a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.646181Z",
     "iopub.status.busy": "2024-11-19T22:42:40.645617Z",
     "iopub.status.idle": "2024-11-19T22:42:40.690809Z",
     "shell.execute_reply": "2024-11-19T22:42:40.690237Z"
    },
    "papermill": {
     "duration": 0.058901,
     "end_time": "2024-11-19T22:42:40.692265",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.633364",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSListmodeSystemMatrix(object_meta=object_meta, proj_meta=proj_meta_listmode)\n",
    "sample_object = torch.rand(object_meta.shape) # object has batch dimension\n",
    "FP = system_matrix.forward(sample_object)\n",
    "BP = system_matrix.backward(FP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e69772e",
   "metadata": {
    "papermill": {
     "duration": 0.012116,
     "end_time": "2024-11-19T22:42:40.719315",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.707199",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The computation of the normalization factor for listmode reconstruction is slightly more nuanced since it requires defining a \"projection of all 1s\". In this case, the projection of all 1s corresponds to a value of 1 in each location of the scanner LUT.\n",
    "\n",
    "* This is distinct from the regular forward/back projection which projects to all the measured detector IDs. For example, if we measure 1000 events, then forward projection produces an array of size 1000, and back projection takes in projections of size 1000 to produce an object. Computation of the normalization factor, however, takes in detector IDs corresponding to each detector, in the case $M=3$, it would take in 18 detector IDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "b9aeacd0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.743404Z",
     "iopub.status.busy": "2024-11-19T22:42:40.743064Z",
     "iopub.status.idle": "2024-11-19T22:42:40.750563Z",
     "shell.execute_reply": "2024-11-19T22:42:40.749948Z"
    },
    "papermill": {
     "duration": 0.02127,
     "end_time": "2024-11-19T22:42:40.751872",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.730602",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSListmodeSystemMatrix(SystemMatrix):\n",
    "    # ----\n",
    "    # SAME AS ABOVE\n",
    "    # ----\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        # There is probably a faster implementation, but I am trying to keep it simple for illustration purposes\n",
    "        projections = []\n",
    "        for i, detector_id in enumerate(self.proj_meta.detector_ids):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            sensitivity_factor_i = self.proj_meta.sensitivity_at_ids[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                projections.append(object[:,coord[1],coord[2]].sum() *  sensitivity_factor_i) # sum along x\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                projections.append(object[coord[1],:,coord[2]].sum() * sensitivity_factor_i)  # sum along y\n",
    "        return torch.tensor(projections)\n",
    "    def backward(self, projections, subset_idx = None):\n",
    "        # There is probably a faster implementation, but I am trying to keep it simple for illustration purposes\n",
    "        object = torch.zeros(object_meta.shape)\n",
    "        projections *= self.proj_meta.sensitivity_at_ids[self.proj_meta.detector_ids]\n",
    "        for i, detector_id in enumerate(detector_ids):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                object[:,coord[1],coord[2]] += projections[i]\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                object[coord[1],:,coord[2]] += projections[i]\n",
    "        return object\n",
    "    # ----\n",
    "    # NEW CODE\n",
    "    # ----\n",
    "    def compute_normalization_factor(self):\n",
    "        norm_BP = torch.zeros(object_meta.shape)\n",
    "        # Now we loop through unique detector ids instead\n",
    "        unique_detector_ids = torch.arange(scanner_LUT.shape[0])\n",
    "        for i, detector_id in enumerate(unique_detector_ids):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            sensitivity_factor_i = self.proj_meta.sensitivity_at_ids[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                norm_BP[:,coord[1],coord[2]] += sensitivity_factor_i\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                norm_BP[coord[1],:,coord[2]] += sensitivity_factor_i\n",
    "        return norm_BP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "92cb175c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.775441Z",
     "iopub.status.busy": "2024-11-19T22:42:40.775079Z",
     "iopub.status.idle": "2024-11-19T22:42:40.779920Z",
     "shell.execute_reply": "2024-11-19T22:42:40.779346Z"
    },
    "papermill": {
     "duration": 0.01808,
     "end_time": "2024-11-19T22:42:40.781203",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.763123",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[1.2862, 1.0203, 1.2276],\n",
       "        [1.2767, 1.2802, 1.2021],\n",
       "        [1.2675, 1.1575, 1.2548]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sensitivity_factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "8fa09e9a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.805665Z",
     "iopub.status.busy": "2024-11-19T22:42:40.805249Z",
     "iopub.status.idle": "2024-11-19T22:42:40.809736Z",
     "shell.execute_reply": "2024-11-19T22:42:40.809206Z"
    },
    "papermill": {
     "duration": 0.018185,
     "end_time": "2024-11-19T22:42:40.811126",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.792941",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSListmodeSystemMatrix(object_meta=object_meta, proj_meta=proj_meta_listmode)\n",
    "norm_BP = system_matrix.compute_normalization_factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ed422116",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:40.835237Z",
     "iopub.status.busy": "2024-11-19T22:42:40.834855Z",
     "iopub.status.idle": "2024-11-19T22:42:40.958759Z",
     "shell.execute_reply": "2024-11-19T22:42:40.957875Z"
    },
    "papermill": {
     "duration": 0.138322,
     "end_time": "2024-11-19T22:42:40.960605",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.822283",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fb3c0637910>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANEAAADLCAYAAADwd3YSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAPvElEQVR4nO3df2hbVf8H8PdNapOiTaBK1pXWURiTuTJ1P9COOTZ1lSKD/aUw2AQrOEiL0r82faCdoEGGsoGuTJTOgZtluq0bzmJhtrXzKV9XWur3D4ePDBpK60+W1tLlx73n+aNrHrKcJPfmNL25y/sF54+c3pt7aPPp53PO/RFNCCFARHlz2T0AIqdjEBEpYhARKWIQESliEBEpYhARKWIQESliEBEpYhARKWIQESmyFERdXV3YuHEjfD4ffD4fGhsb8c033xRqbESOoFm5du7y5ctwu91Yu3YtAOCzzz7D0aNHMTY2hg0bNhRskETFzFIQyVRVVeHo0aNoaWlZrjEROUpZvjvquo5z585hfn4ejY2NGbeLRqOIRqPJ14Zh4O+//8aDDz4ITdPyPTw5hBACc3NzqKmpgcuVPnu4ffs2YrFYxv3Ly8vh9XoLOUR1wqKJiQlx//33C7fbLfx+v/j666+zbt/R0SEAsJV4C4fDaZ+NhYUFUR1wZ92vurpaLCwsWP2YrijL5VwsFsPk5CRu3bqFr776Cp988gkGBwfx6KOPSre/OxNFIhE8/PDD2PTCv+C+r3j/wzz6+v/bPYScro4V/zzUuH0bU2++g1u3bsHv96f8bHZ2Fn6/H/+5XgdfZXqWmp0zsHZLGJFIBD6fb6WGbJnlcq68vDy5sLBlyxb8+OOPOH78OE6ePCnd3uPxwOPxpPW77/OirIiDqPyBcruHkJOronh/f3fLVro/UKnhgcr0nxtwRrmf95xoiRAiJdMQWRUXOuKSgiguDBtGY52lIHrzzTfR3NyMuro6zM3N4YsvvsDAwAD6+voKNT4qAQkYiGfodwJLQfTbb79h//79mJ6eht/vx8aNG9HX14fdu3cXanxUAgwIGEjPRLK+YmQpiD799NNCjYNKWFyIDOXcPRhERIUQEwIxScDI+ooRg4hsZ9xpsn4nYBCR7RJCQ1ykL2cnJH3FiEFEttOhQZecE5L1FSMGEdkuLlyIi/QrFuLOmBIxiMh+MbgRk9zaFmMmIjJHCA2GZP4jOCciMicm3LhPUs7FGERE5sThQhxuSb8zMIjIdrpwQZdkIp0nW4nMScAtzUQJG8aSDwYR2S4uyhAXknKOcyIic3ShQZcEjKyvGDGIyHbMRESK4nAhJguie/F+IqJCMOCCIbliQdZXjBhEZLu4cKNMWs4xExGZEhNlcIv0j2LMGTHEICL7GRmunZP1FSMGEdkuIcoQl2SiBDMRkTlx4YabcyKi/OmQ38Wqr/xQ8sIgItvFjTK4jfSPYtxgJiIyJSHc0isWEvfiY4SJCoGrc0SK4sINl3RhgZmIyBQGEZEilnNEihIZMhEXFohMShguuAxJEBnOOFPEICLbGdCkXy1ZMl83SaQqbrihSTJRXNJXjBhEZDsdLiRkj8ziTXlE5nB1jkhRIkM5l2A5R2ROQrigSco5WYlXjBhEZDunl3POCHW6pyUMV8ZmRSgUwtatW1FZWYlAIIC9e/fixo0bpve/du0aysrK8Pjjj1s6LoOIbKcLDQnhSmtWn4A6ODiIYDCIkZER9Pf3I5FIoKmpCfPz8zn3jUQiOHDgAJ599lnL42c5R7ZbrnKur68v5XV3dzcCgQBGR0exY8eOrPu+9tpr2LdvH9xuNy5evGjpuJYykWq6JJLJVc7Nzs6mtGg0aup9I5EIAKCqqirrdt3d3fj111/R0dGR1/gtBZFKuiTKZOn7iWQNAOrq6uD3+5MtFArlfE8hBNrb27F9+3Y0NDRk3O6XX37BoUOH8Pnnn6OsLL/CzNJeKumSKJNc5Vw4HIbP50v2ezyenO/Z2tqKiYkJDA8PZ9xG13Xs27cPR44cwbp16/IY+SKlOZHZdEmUjW64oElW4vQ7fT6fLyWIcmlra8OlS5cwNDSE2trajNvNzc3h+vXrGBsbQ2trKwDAMAwIIVBWVoZvv/0WzzzzTM7j5R1EZtNlNBpNqWFnZ2fzPSTdowzDlQyYu/utEEKgra0NFy5cwMDAAOrr67Nu7/P58NNPP6X0nThxAlevXsWXX36Zc/8leQeRmXQJLC5GHDlyJK3/96c0uLzFezJtuPbfdg8hp4N2D8CE2D8xfJZjGwFA9pxGqw/MCgaDOHPmDHp7e1FZWYmZmRkAgN/vR0VFBQDg8OHDmJqawunTp+FyudISQCAQgNfrzZoY7pbXeaKldPndd99lTZdLg45EIskWDofzOSTdw3ItLJjV1dWFSCSCnTt3YvXq1cnW09OT3GZ6ehqTk5PLOn5LmchqugQWJ4FmJoJUunRDAwzJE1AlfdkIE48dPnXqVNafd3Z2orOz09JxLQWRmXRJZJUQGoRkdU7WV4ws5Usz6ZLIKv3OwoKsOYHlco5ouRkGoElKN8MZD/vhtXNkP6eXcwwisp0hNGgOvp+IQUT2MzQI2UqcxdU5uzCIyHZCZDjZ6pApOIOIbCcMF4RkJU7WV4wYRGQ7YSw2Wb8TMIjIdlydI1IkhHxhgUFEZJbQFpus3wEYRGQ/Afl9D1ydIzIpw1XcPE9EZBJX54hUcU5EpEYzFpus3wkYRGQ/zomIFHF1jkgRMxGRGqfPiZxxmSxREWMmIttpQpM+Y0F2t2sxYhCR/Yw7TdbvAAwisp0mFpus3wkYRGQ/ZiIiNZqRYU7EJW4ik3iylUiN088TMYjIfhmCiHMiIrNYzhGp4RI3kSpmIiI1msiwsMAgIjKHq3NEqljOEalhJiJSxCAiUsVyjkgNMxGRIgYRkSqWc0RqeNkPkSKnl3OWH5k1NDSEPXv2oKamBpqm4eLFiwUYFpUUkaU5gOUgmp+fx2OPPYYPP/ywEOOhErR07Vxac0gQWS7nmpub0dzcXIixUIlyejlX8DlRNBpFNBpNvp6dnS30IclpuDqXXSgUwpEjR9L66w/9H8q0+wp9+LytxUG7h5BT9bXi/5Ql4rdzbrNcmSgUCuH8+fP4+eefUVFRgW3btuG9997DI488knGf8+fPo6urC+Pj44hGo9iwYQM6Ozvx/PPPmz5uwZ/FffjwYUQikWQLh8OFPiQ5jHQ+lOm5C1kMDg4iGAxiZGQE/f39SCQSaGpqwvz8fMZ9hoaGsHv3bly5cgWjo6PYtWsX9uzZg7GxMdPHLXgm8ng88Hg8hT4MOdkylXN9fX0pr7u7uxEIBDA6OoodO3ZI9zl27FjK63fffRe9vb24fPkynnjiCVPH5Xkisp1mCGhGesQs9d09jzb7jzkSiQAAqqqqTI/FMAzMzc1Z2sdyOffPP/9gfHwc4+PjAICbN29ifHwck5OTVt+KCEDucq6urg5+vz/ZQqFQzvcUQqC9vR3bt29HQ0OD6bG8//77mJ+fx4svvmh6H8uZ6Pr169i1a1fydXt7OwDg5ZdfxqlTp6y+HVHOy37C4TB8Pl+y30wWam1txcTEBIaHh02P4+zZs+js7ERvby8CgYDp/SwH0c6dOyFE8a8KkXPkWp3z+XwpQZRLW1sbLl26hKGhIdTW1prap6enBy0tLTh37hyee+4508cCOCeiYpDhaT9WFxaEEGhra8OFCxcwMDCA+vp6U/udPXsWr7zyCs6ePYsXXnjB2kHBIKJiIMRik/VbEAwGcebMGfT29qKyshIzMzMAAL/fj4qKCgCLp1ympqZw+vRpAIsBdODAARw/fhxPPfVUcp+Kigr4/X5Tx+V3tpLtlus8UVdXFyKRCHbu3InVq1cnW09PT3Kb6enplEWwkydPIpFIIBgMpuzz+uuvmz4uMxHZTtMBTfLvXNOtvY+Zufrdi18DAwPWDiLBICLb8aY8IkW5TrYWOwYR2Y63QhCpWqbVObswiMh2zEREijRdQHNJ5kQ6MxGRObyzlUiNJjKsznFORGQO50REijQhpFmHmYjIJE0X0CSXJ3BhgcgsQyw2Wb8DMIjIdrx2jkgRyzkiVSzniNRwdY5IlSEAWenGTERkjmYIaJIzq7yfiMgs3gpBpEbTBTTJ1aZcnSMyy8hw8ZzhjIvnGERkP5ZzRGpYzhGp0g0AktJNZzlHZA7LOSJFwpAvIghmIiJzdB0QkmcGGxafI2wTBhHZj+UckSLdkJduPE9EZJIhIF2d47VzRCYZGZa4mYmITGIQEakRug4hWZ0TXJ0jMklkuD2cq3NEJum6/LslZeeOihCDiGwndB1CEkSyEq8Y5fXt4SdOnEB9fT28Xi82b96M77//frnHRaVk6WSrrDmA5SDq6enBG2+8gbfeegtjY2N4+umn0dzcnPK15kSW6MZiSZfWnLE6ZzmIPvjgA7S0tODVV1/F+vXrcezYMdTV1aGrq6sQ46MSIHQ9Y3MCS3OiWCyG0dFRHDp0KKW/qakJP/zwg3SfaDSKaDSafB2JRAAACcSL+kucjNu37R5CTol4Ef8C79Dji79HkaU0ixsxCMmHIYF4wca1rIQFU1NTAoC4du1aSv8777wj1q1bJ92no6Nj6XvQ2Eq4hcPhtM/GwsKCqK6uzrpfdXW1WFhYsPIxXXF5rc5pmpbyWgiR1rfk8OHDaG9vT76+desW1qxZg8nJSfj9/nwOTwBmZ2dRV1eHcDgMn89n93AyEkJgbm4ONTU1aT/zer24efMmYrFYxv3Ly8vh9XoLOURlloLooYcegtvtxszMTEr/77//jlWrVkn38Xg88Hg8af1+v7+o//hO4fP5iv73mO2fpdfrLfogycXSwkJ5eTk2b96M/v7+lP7+/n5s27ZtWQdG5BSWy7n29nbs378fW7ZsQWNjIz7++GNMTk7i4MGDhRgfUdGzHEQvvfQS/vrrL7z99tuYnp5GQ0MDrly5gjVr1pja3+PxoKOjQ1rikXn8PRYPTQiHnBYmKlJ5XfZDRP/DICJSxCAiUsQgIlK0okHEWyjUhUIhbN26FZWVlQgEAti7dy9u3Lhh97BK2ooFEW+hWB6Dg4MIBoMYGRlBf38/EokEmpqaMD8/b/fQStaKLXE/+eST2LRpU8otE+vXr8fevXsRCoVWYgj3pD/++AOBQACDg4PYsWOH3cMpSSuSiZZuoWhqakrpz3YLBZmzdGtJVVWVzSMpXSsSRH/++Sd0XU+7SHXVqlVpF7OSeUIItLe3Y/v27WhoaLB7OCVrRR9UYuUWCsqttbUVExMTGB4etnsoJW1FgiifWygou7a2Nly6dAlDQ0Oora21ezglbUXKOd5CsXyEEGhtbcX58+dx9epV1NfX2z2kkrdi5RxvoVgewWAQZ86cQW9vLyorK5PZ3e/3o6KiwubRlaiVvBf9o48+EmvWrBHl5eVi06ZNYnBwcCUPf09AhmcRdHd32z20ksVbIYgU8do5IkUMIiJFDCIiRQwiIkUMIiJFDCIiRQwiIkUMIiJFDCIiRQwiIkUMIiJFDCIiRf8Fzq45g1SjDOIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 200x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(2,2))\n",
    "plt.pcolormesh(norm_BP[:,:,1])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce5513ff",
   "metadata": {
    "papermill": {
     "duration": 0.011441,
     "end_time": "2024-11-19T22:42:40.985941",
     "exception": false,
     "start_time": "2024-11-19T22:42:40.974500",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now we can use the listmode system matrix with the rest of the library.\n",
    "\n",
    "**For list mode system matrices, we no longer need to provide `projections` to the likelihood function**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "90bce2c6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:41.014743Z",
     "iopub.status.busy": "2024-11-19T22:42:41.014359Z",
     "iopub.status.idle": "2024-11-19T22:42:41.018207Z",
     "shell.execute_reply": "2024-11-19T22:42:41.017578Z"
    },
    "papermill": {
     "duration": 0.017812,
     "end_time": "2024-11-19T22:42:41.019481",
     "exception": false,
     "start_time": "2024-11-19T22:42:41.001669",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSListmodeSystemMatrix(object_meta=object_meta, proj_meta=proj_meta_listmode)\n",
    "likelihood = NegativeMSELikelihood(system_matrix, scaling_constant=0.01)\n",
    "reconstruction_algorithm = OSEM(likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "9893b652",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:41.043160Z",
     "iopub.status.busy": "2024-11-19T22:42:41.042823Z",
     "iopub.status.idle": "2024-11-19T22:42:42.681511Z",
     "shell.execute_reply": "2024-11-19T22:42:42.680612Z"
    },
    "papermill": {
     "duration": 1.652942,
     "end_time": "2024-11-19T22:42:42.683447",
     "exception": false,
     "start_time": "2024-11-19T22:42:41.030505",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "recon = reconstruction_algorithm(n_iters=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "83b34a34",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:42.708361Z",
     "iopub.status.busy": "2024-11-19T22:42:42.707928Z",
     "iopub.status.idle": "2024-11-19T22:42:42.793314Z",
     "shell.execute_reply": "2024-11-19T22:42:42.792639Z"
    },
    "papermill": {
     "duration": 0.098963,
     "end_time": "2024-11-19T22:42:42.794892",
     "exception": false,
     "start_time": "2024-11-19T22:42:42.695929",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fb3c0397590>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAL8AAAC2CAYAAACI5VzOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAMzUlEQVR4nO3dX2hUVx4H8O/8u5lkbAJNYlJqjNZCTf/YNROwqQQpVIsVWXGhoQW1oFQhtqR9aLVCG/rSUvpPWLUWAqUF02xbQWGFVhZqI9oHQ9J9SB8WGk1rY2NMNzExk8nMvfswMZs5dya998yYOyfn+wkH65lh7il8PTlz5s75+SzLskCkIb/XAyDyCsNP2mL4SVsMP2mL4SdtMfykLYaftMXwk7aCXg+A1BSLxRCPxzM+ZhgGwuHwAo/IPcfhn/T57uQ48mLsv16PwJmqf3k9Age2Z//gPxaLYdm99+LGyEjGx6urq9Hf31/w/wA485Nr8XgcN0ZG8M9330WkuDjtsYnJSWx59VXE43GGnxavSDyOJX7hbWOWpVAhYvhJ3sAAYBjpfQw/aeHKFSAUSu+bnvZmLBIYfpI3OgoEhQglEt6MRQLDT/IGBgBxzW+a3oxFAsNP8n77zesR5IThJ2l+2G8RUOmWAYafpDH8pK3ATBP7VMHwkzTfTBP7VMHwkzTO/KQthp+0ptIyR8TwkzTO/KQtY6aJfapg+Eka9/lJW5z5SVtBAKEMfapQaaxUYPiGl7QVgn3mF/9eyBh+ksbwk7Z8YUA80cZnAYh5MhzXGH6SVwX73qYJ4IoHY5Gg0rYsFZqlAKqFtlTupY4ePYqVK1ciHA4jGo2iq6sr63PPnz+P9evXo7y8HMXFxVi9ejU+/PBD19fkzE/ylsCeIInvr3d2dqK1tRVHjx7F+vXrcfz4cWzevBl9fX1Yvny57fmRSAT79+/HmjVrEIlEcP78eezduxeRSAQvvPCC4+v6nBak43GF+aP6cYVjY2MoKyvDaDNQKnyqNRYHyjqB0dFRlJaWOrrUunXrUF9fj2PHjs321dXVYdu2bXj77bedDXf7dkQiEXz++eeOng9w2UO5qIJ92VOVemhsbCytTU1NZXyJeDyO7u5ubNq0Ka1/06ZNuHDhgqNh9PT04MKFC9iwYYOr4TP8JC+C1NJnboukHqqpqUFZWdlsyzaDDw8PI5lMoqqqKq2/qqoK165dm/fyy5YtQ1FRERoaGtDS0oI9e/a4Gj7X/CSvCkCx0DeZ+uOXX35JW/YUFRXN+1I+YVltWZatT9TV1YXx8XH88MMPOHDgAO6//348++yzTkfP8FMOwjNtrpm3CqWlpY7W/BUVFQgEArZZfmhoyPbbQLRy5UoAwCOPPILff/8dbW1trsLPZQ/JCyA1fc5tLm/uMQwD0WgUZ8+eTes/e/YsHn/8ccevY1lW1vcV2XDmJ3mZZn6J0wpfeeUV7NixAw0NDWhsbMQnn3yCgYEB7Nu3DwBw8OBBXL16FZ999hkA4MiRI1i+fDlWr14NILXv/9577+HFF190dV2Gn+TlKfzNzc24ceMG3nrrLQwODuLhhx/GmTNnUFtbCwAYHBzEwMDA/y9hmjh48CD6+/sRDAaxatUqvPPOO9i7d6+r63Kf3wOLZp//JFAaER6bAMq2u9vn9wpnfpKXaeZPejEQOQw/ySuCPfzqHM/P8FMObu/wiH2KUGioVGiSoVQT+1TB8JM0hp+0ZQVSTexTBcNP0pLBDDO/QolSaKhUaEwj1cQ+VTD8JM30p5rYpwqGn6SZoQwzP9/wkg7iAWAqYO9TBcNP0hK+VBP7VMHwk7RYAAgF7H2qYPhJ2lSG8IvLoELG8JM0bZY9xZ13chj50V7m9Qic2brpz5/jtVoHz5kKAEHO/KSjmB8I+O19qmD4SVrCn2pinyoYfpI2iKzH9iiB4Sdpv8L+RS5FjuYHwPBTDsYAiCfluDs5x1sMP0m7CnsZomkvBiKJ4SdpvyIvx/N7huEnaWOwn06o0MklDD/JuwpA/EDX0QloBYLhJ2ljSdineoWmfoaf5Jmwn80pcVanVxh+khefaWKfIhh+kjcN+96mQnudDD/JM2Ff4yu07FHoNiQqONNZmgQ3RahPnjyJjRs3orKyEqWlpWhsbMQ333zj+poMP8nLU/hvF6E+dOgQenp60NTUhM2bN6cVpJjr+++/x8aNG3HmzBl0d3fjiSeewNatW9HT0+Pquo6LU+Afhf8Vnb8/4/UInNk67vUI/lztkj8vToF/A7hLePAmgDULX4T6oYceQnNzM9544w1Hzwc481MuErDP+jP3NyxkEWrTNHHz5k3cfffdrobP8JO8eJaGhSlCfdv777+PiYkJPPOMu1/93O0hefN8yLUQRagBoKOjA21tbTh16hSWLl3qZNSzGH6SN8+HXAtRhLqzsxO7d+/Gl19+iSeffNLFwFO47CF586z5nZItQt3R0YHnn38eJ06cwJYtW1wOPIUzP8nL041tbotQd3R0YOfOnTh8+DAee+yx2d8axcXFqV0ohxh+kpen2xvcFqE+fvw4EokEWlpa0NLSMtu/a9cufPrpp46vy31+Dyyaff5TAIQi1JgA8FcWoabFjje2kbZuv+EV+xTB8JM8fpOLtMVlD2krDnuC+E0u0oIF++0NCh3fwPCTvDjsB/dw5ictmEWAKXz+Y1pQ5cROhp/k+UoAn3B7mM8Ew0+LX6ACCAjrnkASwB+eDMcthp/kBSqAgBChQALAfzwZjlsMP8nzRQCfECGfOh/xMvwkz18OBAyhT53tHoaf5PkrAb/w9US/Gm92AYafcuEvyRB+dQrxMvwkz3cX4BNK0vnUKUnH8FMOQrBX5VLntk6Gn3Lgh/3+BnXORGD4KQeZZn5udZIWGH7SVgD2ZQ93e0gLQahchprhpxwYM20uLntIB2Yw1cQ+RagzUio8ZhFgCh9ymYtwn7971aU7OY682H+yweshODLxN69H4ICT7+KaYSCpQfiJbKxgqol9ilBnpFR4kkX2mT/JN7ykg2SGZQ/DT1rgsoe0lQwBScPepwh1bsGjwjNPNUa33FRgHxwcxHPPPYcHHngAfr8fra2tUtdk+Ene7eMK5zaJ4wrdVmCfmppCZWUlDh06hEcffVR6+Aw/ycvTzP/BBx9g9+7d2LNnD+rq6vDRRx+hpqYmrSL7XCtWrMDhw4exc+dOVzW4RAw/yRMrMc45snwhK7DLYvhJXjJLw8JWYJfF3R6SN09ZooWqwJ4Lhp/kTcO+dpj5x7AQFdhzxWUPyYtlaS7IVmDPB878JG8E9u+ySOz2uK3ADgC9vb0AgPHxcVy/fh29vb0wDAMPPvig4+sy/CTvD+TlW4xuK7ADwNq1a2f/u7u7GydOnEBtbS0uX77s+LqOK7B3d3c7flGvRK/wfv58icwTi9kK7E29QPCu9AcTN4Guv7ACOy1yf8B+WIM632Vh+CkHI7BvmYjVGQsYw0/yJsHwk6bUKL2VFcNPOZiCfbuHxSlICwnYD6ni1xhJCzHYI8TiFKSFKdgjxGUPaYHLHtLWFOyfcnHmJy3EYN/o55qfNOBHEj7hfgYLSWU+52L4SVoFkvAL4TeRxJBH43GL4SdpFUgiIIQ/yfCTDkoABIWDetTZ62H4KQcVSCIkzPzTCt3TzPCTtEqYMIS3t3Fl3u4y/JSDkgzhDzL8pIMKmAgLYY8x/KSDQIbdHvHvhYzhJ2nWzI/YpwqGn6QlZ37EPlUw/CTNhGkLu8k1P+mAyx7SFpc9pC1z5kfsUwXDT9IYftLWKEZgIL3oRJzf5CIdjGAYIeGM8mnZWqQeYPhJ2iRuISGcST4tc0a5Rxh+kjaCYQSFCCUUuqOfZYlI2giGcQPX09oIhqVey00FdgA4d+4cotEowuEw7rvvPnz88ceur8nwk7QYbmESE2kthluuX8dtBfb+/n48/fTTaGpqQk9PD15//XW89NJL+Prrr11dl5VZPLBYKrOEUQIfhBKisBDDLVeVWdatW4f6+vq0iut1dXXYtm1bxvq9r732Gk6fPo2ffvpptm/fvn348ccfcfHiRUfXBDjzUw7mm/nvZAX2ixcv2p7/1FNP4dKlS5iedv6G2/Eb3mg06vhFPRNV476SiBrDzMowDFRXV2etkL5kyRLU1NSk9b355ptoa2uzPVemAvu1a9cyPj+RSGB4eBj33HOPo/8P7vaQa+FwGP39/YjHM+/pZ6qenu8K7Jmen6l/Pgw/SQmHwwiHwzm/jkwF9ky/dYaGhhAMBlFeXu742lzzk6dkKrA3Njbanv/tt9+ioaEBoZBYKWYeFpHHvvjiCysUClnt7e1WX1+f1draakUiEevy5cuWZVnWgQMHrB07dsw+/+eff7ZKSkqsl19+2err67Pa29utUChkffXVV66uy/BTQThy5IhVW1trGYZh1dfXW+fOnZt9bNeuXdaGDRvSnv/dd99Za9eutQzDsFasWGEdO3bM9TUd7/MTLTZc85O2GH7SFsNP2mL4SVsMP2mL4SdtMfykLYaftMXwk7YYftIWw0/a+h+/OIAI9x76bgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 200x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(2,2))\n",
    "plt.pcolormesh(recon[:,:,1], vmin=0, cmap='nipy_spectral')\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc2a8e5b",
   "metadata": {
    "papermill": {
     "duration": 0.011324,
     "end_time": "2024-11-19T22:42:42.819291",
     "exception": false,
     "start_time": "2024-11-19T22:42:42.807967",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Part 2: Incorporating Subsets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ac55a2f",
   "metadata": {
    "papermill": {
     "duration": 0.011585,
     "end_time": "2024-11-19T22:42:42.843916",
     "exception": false,
     "start_time": "2024-11-19T22:42:42.832331",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Incorporating subsets is slightly more intuitive with listmode data. In this case, we split up all the detected events into $N$ approximately equal subsets using the `torch.tensor_split` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "7ddc7653",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:42.870503Z",
     "iopub.status.busy": "2024-11-19T22:42:42.870070Z",
     "iopub.status.idle": "2024-11-19T22:42:42.880572Z",
     "shell.execute_reply": "2024-11-19T22:42:42.879949Z"
    },
    "papermill": {
     "duration": 0.025889,
     "end_time": "2024-11-19T22:42:42.881877",
     "exception": false,
     "start_time": "2024-11-19T22:42:42.855988",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class EXSListmodeSystemMatrix(SystemMatrix):\n",
    "    # ----\n",
    "    # NEW CODE\n",
    "    # ----\n",
    "    def set_n_subsets(self, n_subsets):\n",
    "        self.n_subsets = n_subsets\n",
    "        idx = torch.arange(proj_meta_listmode.detector_ids.shape[0])\n",
    "        self.subset_indices_array = torch.tensor_split(idx, self.n_subsets)\n",
    "    def get_projection_subset(self, projections, subset_idx):\n",
    "        # Needs to consider cases where projection is simply a 1 element tensor in the numerator, but also cases of scatter where it is a longer tensor\n",
    "        if (projections.shape[0]>1)*(subset_idx is not None):\n",
    "            subset_indices = self.subset_indices_array[subset_idx]\n",
    "            proj_subset = projections[subset_indices]\n",
    "        else:\n",
    "            proj_subset = projections\n",
    "        return proj_subset\n",
    "    def get_weighting_subset(self, subset_idx):\n",
    "        if subset_idx is None:\n",
    "            return 1\n",
    "        else:\n",
    "            # Fraction of events in the subset\n",
    "            return self.subset_indices_array[subset_idx].shape[0] / proj_meta_listmode.detector_ids.shape[0]\n",
    "    def forward(self, object, subset_idx = None):\n",
    "        # There is probably a faster implementation, but I am trying to keep it simple for illustration purposes\n",
    "        projections = []\n",
    "        if subset_idx is None:\n",
    "            detector_ids_subset = self.proj_meta.detector_ids\n",
    "        else:\n",
    "            detector_ids_subset = self.proj_meta.detector_ids[self.subset_indices_array[subset_idx]]\n",
    "        for i, detector_id in enumerate(detector_ids_subset):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            sensitivity_factor_i = self.proj_meta.sensitivity_at_ids[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                projections.append(object[:,coord[1],coord[2]].sum() *  sensitivity_factor_i) # sum along x\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                projections.append(object[coord[1],:,coord[2]].sum() * sensitivity_factor_i)  # sum along y\n",
    "        return torch.tensor(projections)\n",
    "    def backward(self, projections, subset_idx = None):\n",
    "        # There is probably a faster implementation, but I am trying to keep it simple for illustration purposes\n",
    "        object = torch.zeros(object_meta.shape)\n",
    "        if subset_idx is None:\n",
    "            detector_ids_subset = self.proj_meta.detector_ids\n",
    "        else:\n",
    "            detector_ids_subset = self.proj_meta.detector_ids[self.subset_indices_array[subset_idx]]\n",
    "        projections *= self.proj_meta.sensitivity_at_ids[detector_ids_subset]\n",
    "        for i, detector_id in enumerate(detector_ids_subset):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                object[:,coord[1],coord[2]] += projections[i]\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                object[coord[1],:,coord[2]] += projections[i]\n",
    "        return object\n",
    "    def compute_normalization_factor(self, subset_idx=None):\n",
    "        fraction_considered = self.get_weighting_subset(subset_idx)\n",
    "        norm_BP = torch.zeros(object_meta.shape)\n",
    "        # Now we loop through unique detector ids instead\n",
    "        unique_detector_ids = torch.arange(scanner_LUT.shape[0])\n",
    "        for i, detector_id in enumerate(unique_detector_ids):\n",
    "            coord = self.proj_meta.scanner_LUT[detector_id]\n",
    "            sensitivity_factor_i = self.proj_meta.sensitivity_at_ids[detector_id]\n",
    "            if coord[0]==0: # If angle 0:\n",
    "                norm_BP[:,coord[1],coord[2]] += sensitivity_factor_i\n",
    "            elif coord[0]==1: # If angle 90:\n",
    "                norm_BP[coord[1],:,coord[2]] += sensitivity_factor_i\n",
    "        return norm_BP * fraction_considered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "a2debf5a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:42.906196Z",
     "iopub.status.busy": "2024-11-19T22:42:42.905627Z",
     "iopub.status.idle": "2024-11-19T22:42:42.909300Z",
     "shell.execute_reply": "2024-11-19T22:42:42.908754Z"
    },
    "papermill": {
     "duration": 0.017123,
     "end_time": "2024-11-19T22:42:42.910590",
     "exception": false,
     "start_time": "2024-11-19T22:42:42.893467",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "system_matrix = EXSListmodeSystemMatrix(object_meta=object_meta, proj_meta=proj_meta_listmode)\n",
    "likelihood = NegativeMSELikelihood(system_matrix, scaling_constant=0.01)\n",
    "reconstruction_algorithm = OSEM(likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "66a7cea9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:42.935140Z",
     "iopub.status.busy": "2024-11-19T22:42:42.934681Z",
     "iopub.status.idle": "2024-11-19T22:42:45.515625Z",
     "shell.execute_reply": "2024-11-19T22:42:45.514558Z"
    },
    "papermill": {
     "duration": 2.595581,
     "end_time": "2024-11-19T22:42:45.518115",
     "exception": false,
     "start_time": "2024-11-19T22:42:42.922534",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "recon_2subsets = reconstruction_algorithm(n_iters=20, n_subsets=2)\n",
    "reconstruction_algorithm = OSEM(likelihood)\n",
    "recon_1subset = reconstruction_algorithm(n_iters=40)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1371e63c",
   "metadata": {
    "papermill": {
     "duration": 0.011251,
     "end_time": "2024-11-19T22:42:45.542647",
     "exception": false,
     "start_time": "2024-11-19T22:42:45.531396",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now lets plot the reconstructions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "8e190a1d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-19T22:42:45.569088Z",
     "iopub.status.busy": "2024-11-19T22:42:45.568444Z",
     "iopub.status.idle": "2024-11-19T22:42:45.730212Z",
     "shell.execute_reply": "2024-11-19T22:42:45.729582Z"
    },
    "papermill": {
     "duration": 0.17736,
     "end_time": "2024-11-19T22:42:45.731917",
     "exception": false,
     "start_time": "2024-11-19T22:42:45.554557",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fb3c02a7150>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAC2CAYAAACVvjV8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUt0lEQVR4nO3db2gc5b4H8O9mk822W3cvNnUrmsaoYKNVb7OBmpRQBLsSxXuLh2NQSCu02EBVYl8cGwNa+iZy8F/hNLGFglRozT0qKJyAhgPWlFQOLqn3RX1xwdTt6UmMSXuytjbZZHfui2RnsjPPJLvbZJ95dr6fMpzT3042z9IvP5+ZeWbWo2maBiIiWnVlsgdAROQWbLhEREXChktEVCRsuERERcKGS0RUJGy4RERFwoZLRFQkbLhEREVSLnsA5AzT09NIJpPC13w+H/x+f5FHRJRtqYwCauQ054Z70+NZzXEULPFv2SOwCv9d9ghsPCu+qXB6ehp333UXJq9eFb6+ceNGjIyMOD7MCYdmdPaa7BFYrf+H7BHYiBaWUUCNnHKGS0gmk5i8ehV/+/OfEVizJuu1Gzdv4uk//QnJZNLRQabStlRGAXVyyoZLukAyiXVlptP6SxzCERWbMKOAMjllwyVDPA74fNk1RYJMLiHKKKBMTtlwyfDzz0BFRXZtdlbOWIhERBkFlMkpGy4ZpqaAclMk5ubkjIVIRJRRQJmcsuGSIR4HzOfH0mk5YyESEWUUUCanbLhk+Ne/ZI+AaGmKZ5QNl3RlsN56yFsRyUlEGYVNzYnYcEnHhktOx4ZLJcO7sJlrRE4hyihsak7Ehks6z8JmrhE5hSijsKk5ERsu6TjDJafjDJdKBhsuOR0bLpUUVQ7NyL1UzigbLuk4wyWn4wyXSoZvYTPXiJxClFHY1JyIDZd0XIdLTsd1uFQyOMMlp+MMl0pGOQDzg+8YEHISUUYzdRWoMk4qAl40I6fjRTMqGRWwzh5EswkiWUQZhU3NidhwSceGS07Hhkslw+MHzN807tEATEsZDpGFKKOAOjllwyVDGNb1NWkAP0sYC5GIKKOAMjlVZfkaFcMdADaatjsKe6uenh7U1tbC7/cjEolgcHDQdt9z585h+/btWL9+PdasWYPNmzfj/fffL+wXU2kTZbTAnMrIKGe4ZFgHayIK+G6+vr4+dHR0oKenB9u3b8fx48fR0tKCixcvYtOmTZb9A4EAXn75ZTzyyCMIBAI4d+4c9u/fj0AggJdeeqmgj0IlSpRRIO+cysqoR9M0LZcdb4pOnDhA4t+yR2AV/rvsEdh4VvxPnUgkEAqFMNUKBE0ryBNJINQHTE1NIRgM5vRrtm3bhvr6evT29uq1uro67Nq1C93d3bkN9dlnEQgE8PHHH+e0PwAkHJrR2WuyR2C1/h+yR2Ajmn9GgfxzKiujPKVAhjCsh2rh+ZcSiUTWNjMzI3yLZDKJWCyGaDSaVY9GoxgaGsppGMPDwxgaGsKOHTsK/SRUqkQZzTOnMjPKhkuGAOYP2RZvgfmXqqurEQqF9M1uFjAxMYFUKoVwOJxVD4fDGBsbW/LX33333aisrERDQwMOHDiAffv23eonolIjymieOZWZUZ7DJUMYwBpT7eb8/1y+fDnrUK2ysnLJt/KYDu81TbPUzAYHB3H9+nV89913OHToEO6//348//zzuY6e3ECUUaCgnMrIKBsuGfwL22ILp9SCwWBO58aqqqrg9XotM4Xx8XHLjMKstrYWAPDwww/jl19+weHDh9lwKZsoo0BeOZWZUZ5SIIMX8/8JXrzleZO6z+dDJBLBwMBAVn1gYABNTU05v4+mabbnicnFRBnNM6cyM8oZLhlEs4d0/m9z8OBBtLW1oaGhAY2NjThx4gTi8Tja29sBAJ2dnbhy5QpOnToFADh27Bg2bdqEzZs3A5hf8/jOO+/glVdeuYUPQyXJboabZ05lZZQNlwwr1HBbW1sxOTmJI0eOYHR0FFu2bEF/fz9qamoAAKOjo4jH48avSKfR2dmJkZERlJeX47777sPbb7+N/fv3F/5ZqDStUMOVlVGuw10Fyq7D/RwIBkyv3QBCz+a3DlcWrsPNnbLrcAUZBdTJKWe4ZBDNHlIyBkJkw26Gq0hO2XDJUAlrmAu4tZdo1YgyCiiTUzZcMmSu+JprRE4hyihsag6kyDCpGFIV85u5RuQUooxm6ipgwyUdGy45HRsulQzNO7+Za0ROIcpopq4CNlzSpcoFM1wmhBxElNFMXQWKDJOKIe2b38w1IqcQZTRTVwEbLunSZfObuUbkFKKMZuoqYMMlXbpCMMNV5GIEuYMoo5m6CthwSZf0AjNea43IKUQZzdRVwIZLujnP/GauETmFKKOZugrYcEk37QUqvNYakVOIMpqpq4ANl3QzgjCLDt+IZBFlNFNXARsu6XhKgZzONacU1vSt5jAKdzIkewRWz0SX30eGmmVen/EC5QrPcINfyB6B2F/+Q/YIrP64TfYIxJb+RjFxRjN1FXCGS7rpMsBbZq0ROYUoo5m6CthwSTdXNr+Za0ROIcpopq4CNlzSjQJYY6rdlDEQIhuijALq5JQNl3T/hPVh+tMyBkJkQ5RRQJ2csuGSLgFgxlQz/51IJlFGYVNzIjZc0l0BYL4lfVbGQIhsiDIKqJNTNlzS/RPWQCjy3XzkEqKMAurklA2XdAkA5uWMinz7NLmEKKOAOjllwyXdFQDmG3Y0GQMhsiHKKKBOTtlwSZdIwTpVUGXqQK4gzChsag7EhkuG9MJmrhE5hSijsKk5EBsuGZILm7lG5BSijMKm5kBsuGSYhXV9jSrrbcgdRBmFTc2B2HDJkIb1XJgih2rkEqKMZuoKUOSRD1QUszZbAXp6elBbWwu/349IJILBwUHbfT///HPs3LkTGzZsQDAYRGNjI7766qvCfjGVNruMFpBTGRllwyXDCgW5r68PHR0d6OrqwvDwMJqbm9HS0oJ4PC7c/9tvv8XOnTvR39+PWCyGxx9/HM888wyGh4cL/yxUmlao4crKqEfTtNyWsP2PMx+p/pfnZI/A6pnrskcgVrNO/E+dSCQQCoWA/wVwm+nF3wA8AkxNTSEYDOb0e7Zt24b6+nr09vbqtbq6OuzatQvd3d05vcdDDz2E1tZWvPnmmzntDwD40qEZ/S/ZI7D645TsEYiFQwVkFMg7p7IyyhkuGeZgnTUs3DOZSCSytpkZ8eNCkskkYrEYotHsr72IRqMYGhrKaRjpdBq//fYbbr/99kI/CZUqUUbzzKnMjLLhkiFpswGorq5GKBTSN7tZwMTEBFKpFMLh7C9LCYfDGBsby2kY7777Lm7cuIHnnnPg4QvJZZfRPHIqM6NcpUCGJW58uHz5ctahWmVl5ZJv5fFkH95rmmapiZw5cwaHDx/GF198gTvuuCOXUZObLHPjQz45lZFRNlwyLHHjQzAYzOncWFVVFbxer2WmMD4+bplRmPX19WHv3r3461//iieeeCKPgZNrLHPjQy45lZlRnlIgwxLncHPl8/kQiUQwMDCQVR8YGEBTU5Ptz505cwYvvvgiTp8+jaeffjrPgZNrLHMONxcyM8oZLhlW6OE1Bw8eRFtbGxoaGtDY2IgTJ04gHo+jvb0dANDZ2YkrV67g1KlTAOaDvHv3bhw9ehSPPfaYPvNYs2bN/JVpoowVeniNrIyy4ZJhhW7tbW1txeTkJI4cOYLR0VFs2bIF/f39qKmpAQCMjo5mrXc8fvw45ubmcODAARw4cECv79mzBx999FH+A6DStUK39srKKNfhrgJl1+F+ASBgevEGgP/Obx2uNFyHmzNl1+GKMgook1POcMnAh9eQ0/HhNVQyMhckzDUipxBlNFNXABsuGfiND+R0/MYHKhk8pUBOx1MKVDKSsCZCkSfpk0uIMpqpK4ANlwwarLdNqvJ1qOQOooxm6gpgwyVDEoBXUCNyClFGM3UFsOGSIV0JpE1rWdMaAPGjGImKTpRRQJmcsuGSwbMW8Jger+FJQ4Ugk0uIMgook1M2XDJ4qwCv6XjNmwJwTcpwiCxEGQWUySkbLhm8VYDXFAnvHID/kzIcIgtRRgFlcsqGSwZPAPCYIuFR5BYecgdRRgFlcsqGS4ay9YDXZ6opcvmX3EGUUUCZnLLhkqFsA1Bm+kqSMudfiCAXEWUUUCanbLhkKFsraLiiRY9EkogyCiiTUzZcMnhuAzx+U21azliIREQZBZTJKRsuLVKxsC2myGOYyCVEGQVUySkbLi1SBut9k/yeUXISUUYzdedjw6VFRLMHNZbbkFvYzXDVyCkbLi3ChktOx4ZLJcML6+GaGld/yS1EGYVNzXnYcGmRclhnD4o8Sp9cQpRRQJWcsuHSIr6FbTE1DtXILUQZBVTJKRsuGdLl85u5RuQUooxm6gpQY5RUHOlKIG1aVJ5WY30juYQoo4AyOc254cbu+341x1Gwlz9vkD0Eixt/kD0CG8t971PaD6TUbbixu5jRXCVKKaOAMjnlDJcMWvn8Zq4ROYUoo5m6AtQYJRVHqtI6e0ipcTGCXEKUUUCZnLLhkiElOFxTJMjkEqKMAsrklA2XDDylQE7HUwpUMlIVQMpnrRE5hSijmboC1HjEDhVH0mYrQE9PD2pra+H3+xGJRDA4OGi77+joKF544QU88MADKCsrQ0dHR2G/lEqfXUYLyKmMjLLhkkEDkDZtyy3TEejr60NHRwe6urowPDyM5uZmtLS0IB6PC/efmZnBhg0b0NXVhUcffbTw8VPpE2W0gJzKyigbLhlWaObw3nvvYe/evdi3bx/q6urwwQcfoLq6Gr29vcL977nnHhw9ehS7d+9GKBQqfPxU+lZohisro2y4ZJi12QAkEomsbWZG/KV9yWQSsVgM0Wg0qx6NRjE0NLSKgydXsMtoHjmVmVE2XDKkbDYA1dXVCIVC+tbd3S18i4mJCaRSKYTD4ax6OBzG2NjYKg6eXMEuo3nkVGZGuUqBDHOwPuVuYXnj5cuXEQwG9XJlpeCbUxfxeDxZf9c0zVIjypsoo5k68supjIyy4ZJhFtZjnoVwB4PBrCDbqaqqgtfrtcwUxsfHLTMKoryJMpqpI7ecyswoTymQYdpmy4PP50MkEsHAwEBWfWBgAE1NTSszTnIvu4zmkVOZGeUMlwxXYX22cwGrFA4ePIi2tjY0NDSgsbERJ06cQDweR3t7OwCgs7MTV65cwalTp/SfuXDhAgDg+vXr+PXXX3HhwgX4fD48+OCDhX0WKk2ijAJ551RWRtlwyXANK/INO62trZicnMSRI0cwOjqKLVu2oL+/HzU1NQDmF5Gb1ztu3bpV//+xWAynT59GTU0NLl26lP8AqHSJMgrknVNZGfVompbTkuFYLJbzmxZT5GfnPWvUqc/DDdj8UycSifm1hc0XgPLbsl+c+w0Y/E9MTU3ldA5XJmY0d059Hm6wkIwCyuSUM1wyXIP1y0/VeK4zuYUoo4AyOWXDJcNVWC+jpmUMhMiGKKOAMjllwyXDTbDhkrOJMgook1M2XDJckz0AomUonlE2XFpkBtZLwOJnJhDJIcpopu58bLi0yBz0eySzakROIcoobGrOw4ZLi0zDGok8bzUjWlWijGbqzseGS4vMwBoJNQ7VyC1EGc3UnY8NlxbhKQVyOp5SoJIxA+uqcjVmDuQWooxm6s7HhkuLTMO6yFGNc2PkFqKMZurOx4ZLujKk4DHdI6khpcqacnIBUUYBdXLKhku6KqRQZgpzGimMSxoPkZkoo4A6OWXDJV0VUvCawpxSJMjkDqKMAurklA2XdGsBlCP78XhqXPsltxBlFFAnp2y4pKtCChWm2cOsKs+9I1cQZRRQJ6dsuKTbgDR8pksPSSUuRZBbiDIKqJNTNlzSrRWEuVyRIJM7iDIKqJNTNlzSVSENvym404oEmdxBlFFAnZyy4ZLOK7gCLLoiTCSLKKOZugrYcEmnLfwx14icQpTRTF0FbLikSy38MdeInEKU0UxdBWy4pEsjbQluWpFzY+QOooxm6ipgwyUdTymQ0/GUApUMnlIgp+MpBSoZ6YU/5hqRU4gymqmrgA2XdGy45HRsuFQypnAVPlRm1ZKKPEmf3EGUUUCdnLLhku4qJlABX1ZtFklJoyGyEmUUUCenbLiku4nfMYfZrNqs6e9EMokyCqiTUzZc0l3FBMpNkZhT5kmj5AaijALq5FT0bWzkUlcxgUn8mrVdxURB79XT04Pa2lr4/X5EIhEMDg4uuf/Zs2cRiUTg9/tx77334sMPPyzo91JpE2W00JzKyCgbLumm8Ttu4kbWNo3f836fvr4+dHR0oKurC8PDw2hubkZLSwvi8bhw/5GRETz11FNobm7G8PAw3njjDbz66qv47LPPbvUjUYkRZbSQnMrKqEfTtJxu0YjFYnm9cbFEfm6QPQSLG3+QPQKxgM0/dSKRQCgUgh9r4YEn6zUNGqbxO6amphAMBnP6Pdu2bUN9fT16e3v1Wl1dHXbt2oXu7m7L/q+//jq+/PJL/Pjjj3qtvb0dP/zwA86fP5/T7wSY0XwkHJrRYAEZBfLPqayMcoZLuqVmuIlEImubmREvw0kmk4jFYohGo1n1aDSKoaEh4c+cP3/esv+TTz6J77//HrOzalwMoeJYboabS05lZjTni2aRSCTnNy2qiPPuoQ44b0hL8vl82LhxI8bGxoSvr1u3DtXV1Vm1t956C4cPH7bsOzExgVQqhXA4nFUPh8O27z82Nibcf25uDhMTE7jzzjtz+hzMaO6CzhvSkpbLKJB7TmVmlKsUCH6/HyMjI0gmxWsZNU2Dx5N9GFdZaV18vph5f9F7LLe/qE7utFxGgfxzKiOjbLgEYD7Qfr//lt+nqqoKXq/XMlMYHx+3zBAyRDOX8fFxlJeXY/369bc8JioNpZBRnsOlFeXz+RCJRDAwMJBVHxgYQFNTk/BnGhsbLft//fXXaGhoQEVFxaqNldxJakY1ohX2ySefaBUVFdrJkye1ixcvah0dHVogENAuXbqkaZqmHTp0SGtra9P3/+mnn7S1a9dqr732mnbx4kXt5MmTWkVFhfbpp5/K+ghU4mRllA2XVsWxY8e0mpoazefzafX19drZs2f11/bs2aPt2LEja/9vvvlG27p1q+bz+bR77rlH6+3tLfKIyW1kZDTndbhERHRreA6XiKhI2HCJiIqEDZeIqEjYcImIioQNl4ioSNhwiYiKhA2XiKhI2HCJiIqEDZeIqEjYcImIioQNl4ioSP4fOH7famCt9SwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x200 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "max_value = max([recon_1subset[:,:,1].max().item(), recon_1subset[:,:,1].max().item()])\n",
    "plt.figure(figsize=(4,2))\n",
    "plt.subplot(121)\n",
    "plt.pcolormesh(recon_1subset[:,:,1], vmin=0, vmax=max_value, cmap='nipy_spectral')\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "plt.subplot(122)\n",
    "plt.pcolormesh(recon_2subsets[:,:,1], vmin=0, vmax=max_value, cmap='nipy_spectral')\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b08003",
   "metadata": {
    "papermill": {
     "duration": 0.011503,
     "end_time": "2024-11-19T22:42:45.756843",
     "exception": false,
     "start_time": "2024-11-19T22:42:45.745340",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "For a better imaging system (we're not collecting nearly enough angles here), these two images should look very similar."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pytomo_install_test",
   "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.11.10"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 13.046144,
   "end_time": "2024-11-19T22:42:46.889182",
   "environment_variables": {},
   "exception": null,
   "input_path": "t_examplesystemmatrix.ipynb",
   "output_path": "t_examplesystemmatrix.ipynb",
   "parameters": {},
   "start_time": "2024-11-19T22:42:33.843038",
   "version": "2.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}