{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Speeding up gamma calculations\n",
    "\n",
    "The PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function allows the user to compute gamma ($\\gamma$) index distributions between a reference and an evaluation dose distribution in 1D, 2D and 3D, and obtain gamma pass rates for these distributions.\n",
    "\n",
    "In this tutorial, you will compute the gamma pass rate using the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function. You will tweak its parameters to speed up the computation of the gamma index distribution.\n",
    "\n",
    "You will test the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) on two 3D (and 2D) dose distributions, and you will focus on two of the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) parameters, <code>max_gamma</code> and <code>random_subset</code>, which are particularly useful when performing pass rate calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we import the required modules and download the DICOM RTDOSE data that you will use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pydicom\n",
    "import pymedphys\n",
    "import numpy as np\n",
    "\n",
    "reference_filepath = pymedphys.data_path(\"original_dose_beam_4.dcm\") # reference (e.g., calculated) dose\n",
    "evaluation_filepath = pymedphys.data_path(\"logfile_dose_beam_4.dcm\") # evaluation (e.g., measured) dose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both the reference and evaluation datasets are 3D dose distributions.\n",
    "\n",
    "Then, we read the DICOM files using pydicom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reference = pydicom.dcmread(reference_filepath)\n",
    "evaluation = pydicom.dcmread(evaluation_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to extract the x, y, and z axes, as well as the dose grids from the reference and evaluation RTDOSE datasets in a form readable by the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function.\n",
    "\n",
    "We can obtain them using the <code>zyx_and_dose_from_dataset()</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axes_reference, dose_reference = pymedphys.dicom.zyx_and_dose_from_dataset(reference)\n",
    "axes_evaluation, dose_evaluation = pymedphys.dicom.zyx_and_dose_from_dataset(evaluation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<code>axes_reference</code> and <code>axes_evaluation</code> are tuples of numpy arrays, whose values define the grid coordinates in the DICOM frame of reference.\n",
    "\n",
    "Let's check the shape of the reference image axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_ref, y_ref, x_ref = axes_reference\n",
    "\n",
    "print(z_ref.shape, y_ref.shape, x_ref.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, in this case the 3D reference dose image is composed of 152 z-slices, 158 y-slices, and 254 x-slices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<code>dose_reference</code> and <code>dose_evaluation</code> are matrices containing the actual dose values at each grid point.\n",
    "\n",
    "Their shape will match that of the axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(dose_reference)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let' have a look at the dose distributions above and below the z-slice where they disagree the most"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "dose_difference = dose_evaluation - dose_reference\n",
    "\n",
    "max_diff = np.max(np.abs(dose_difference))\n",
    "\n",
    "# get the z-slice with the maximum dose difference\n",
    "z_max_diff, _, _ = np.unravel_index(np.argmax(np.abs(dose_difference), axis=None), dose_difference.shape)\n",
    "\n",
    "# we consider 10 z-slices above and below the maximum dose difference\n",
    "z_start = z_max_diff - 10\n",
    "z_end = z_max_diff + 10\n",
    "\n",
    "fig, axes = plt.subplots(figsize=(15,10), nrows=4, ncols=5, sharex=True, sharey=True)\n",
    "ax = axes.ravel().tolist()\n",
    "ax[0].invert_yaxis() # invert just once as y axes are shared\n",
    "\n",
    "for i, dose_diff_slice in enumerate(dose_difference[z_start:z_end]):\n",
    "\n",
    "    im = ax[i].contourf(x_ref, y_ref, dose_diff_slice, 100, cmap=plt.get_cmap(\"seismic\"), vmin=-max_diff, vmax=max_diff)\n",
    "    ax[i].set_title(f\"Slice Z_{z_start + i}\")\n",
    "    \n",
    "    if i >= 15:\n",
    "        ax[i].set_xlabel(\"x (mm)\")\n",
    "    \n",
    "    if i % 5 == 0:\n",
    "        ax[i].set_ylabel(\"y (mm)\")\n",
    "\n",
    "fig.tight_layout()\n",
    "cbar = fig.colorbar(im, ax=ax, label=\"[Dose Eval] - [Dose Ref] (Gy)\", aspect=40)\n",
    "cbar.formatter.set_powerlimits((0, 0))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the largest dose differences are concentrated in small regions of the slices, which will end up in a fast gamma evaluation, as the computation will be immediate for those points where the dose difference is zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we move to the actual calculation.\n",
    "\n",
    "First of all, we define the acceptance criteria as 1% dose difference and 1mm distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dose_percent_threshold = 1\n",
    "distance_mm_threshold = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the 1% of <code>dose_percent_threshold</code> refers to the 1% of the maximum value of the reference dose, i.e."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Maximum reference dose: {np.max(dose_reference): .2f} Gy\")\n",
    "print(f\"Dose threshold: {np.max(dose_reference) * .01:.4f} Gy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to see in detail the messages traced by the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html), we import the logging module and set the logging level to DEBUG:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.DEBUG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we compute the gamma using the default settings of the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html), and measure its computation time with <code>%%time</code>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "gamma = pymedphys.gamma(\n",
    "    axes_reference, dose_reference, \n",
    "    axes_evaluation, dose_evaluation, \n",
    "    dose_percent_threshold, distance_mm_threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's breakdown the log trace of the gamma function:\n",
    "\n",
    "1. _Global normalisation set to 1.946_: the global normalization is set by default to the maximum value of the reference dose\n",
    "2. _Global dose threshold set to 0.019 (1.00% of normalisation)_: the dose threshold value is (<code>dose_percent_threshold</code> * the global normalization / 100)\n",
    "3. _Distance threshold set to [1]_: the <code>distance_mm_threshold</code>\n",
    "4. _Lower dose cutoff set to 0.389 (20.0% of normalisation)_: the gamma is computed only for those points where the dose value is bigger than this cutoff (20% of the maximum reference dose)\n",
    "5. _Current distance: 0.50 mm | Number of reference points remaining: 128_:\n",
    "    - _Current distance_: the current distance between reference grid point and evaluated grid points. The computation stops once a distance is reached for which the value of the gamma cannot improve further\n",
    "    - _Number of reference points remaining_: the number of points left to be evaluated in the reference dose grid\n",
    "6. _Points tested per reference point: 348 | RAM split count: 1_:\n",
    "    - _Points tested per reference point_: the number of points considered in the evaluation dose grid for the current reference point\n",
    "    - _RAM split count_: in case the number of points to check for the gamma value doesn't fit into the designated RAM memory, the computation is split and performed in stages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output stored in the <code>gamma</code> variable is a matrix containing the gamma index value for each point of the reference dose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the distribution of the gamma values to get a sense of the function output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# remove NaN grid points that were not evaluated as the dose value was below the dose threshold\n",
    "valid_gamma = gamma[~np.isnan(gamma)]\n",
    "\n",
    "out = plt.hist(valid_gamma, bins=20, density=True)\n",
    "plt.xlabel(\"gamma index\")\n",
    "_ = plt.ylabel(\"probability density\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the gamma index distribution has been computed, one can consider the gamma pass rate, defined as the fraction of points having $\\gamma\\le 1$, to assess the agreement between the reference and evaluated dose distributions.\n",
    "\n",
    "Suppose that, for our purposes, we consider a gamma pass rate above 95% in order to have \"passed\" this comparison test.\n",
    "\n",
    "Then, we check if the test passes (probably, you already guessed the answer by looking at the histogram above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Pass rate(\\u03B3<=1): {len(valid_gamma[valid_gamma <= 1]) / len(valid_gamma) * 100}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, all reference points have $\\gamma\\le 1$, thus the pass rate is 100%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the computation of the gamma was quite fast because most of the grid points have zero dose difference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = dose_difference[dose_difference == 0]\n",
    "b = dose_difference[dose_difference != 0]\n",
    "\n",
    "print(f\"Fraction of grid points where [Dose Eval] = [Dose Ref]: {np.sum(dose_difference == 0) / len(dose_difference.flat) * 100:.1f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice, calculations of 3D gamma index distributions can be very time consuming. The problem is further exacerbated in applications where multiple comparisons are required.\n",
    "\n",
    "For many applications, it is sufficient to simply calculate an accurate pass rate.\n",
    "\n",
    "The PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function accepts optional parameters that can speed up the calculation in such cases:\n",
    "- <code>max_gamma</code>: stop the calculation when a certain value of gamma is reached around the reference point\n",
    "- <code>random_subset</code>: use only a random subset of the reference dose grid (range: 0-x*y*z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how these parameters affect the gamma calculation time.\n",
    "\n",
    "For simplicity, we restrict ourself on the z-slice with the highest disagreement.\n",
    "\n",
    "We also add a small shift to the evaluation dose distribution to obtain a worse gamma agreement which will make the computation slower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add 1% of the evaluation dose standard deviation where the evaluation dose in non zero\n",
    "dose_evaluation_Z = np.where(\n",
    "    dose_evaluation[z_max_diff] != 0,\n",
    "    dose_evaluation[z_max_diff] + .1 * np.std(dose_evaluation[z_max_diff]),\n",
    "    dose_evaluation[z_max_diff]\n",
    ")\n",
    "dose_reference_Z = np.array(dose_reference[z_max_diff])\n",
    "\n",
    "# keep only the y, x axes\n",
    "axes_reference_subset = axes_reference[1:]\n",
    "axes_evaluation_subset = axes_evaluation[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# increase the logging level to silence the traces of the PyMedPhys gamma\n",
    "logger.setLevel(logging.ERROR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, let's say we are now interested in evaluating all those points above 1% of the maximum value of the reference dose (instead of the default 20%): we pass the <code>lower_percent_dose_cutoff</code> to the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "gamma = pymedphys.gamma(\n",
    "    axes_reference_subset, dose_reference_Z, \n",
    "    axes_evaluation_subset, dose_evaluation_Z, \n",
    "    dose_percent_threshold, distance_mm_threshold,\n",
    "    lower_percent_dose_cutoff=1 # 1% lower threshold\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_gamma = gamma[~np.isnan(gamma)]\n",
    "\n",
    "out = plt.hist(valid_gamma, bins=20, density=True)\n",
    "plt.xlabel(\"gamma index\")\n",
    "plt.ylabel(\"probability density\")\n",
    "title = plt.title(f\"Gamma passing rate: {np.sum(valid_gamma <= 1) / len(valid_gamma) * 100:.1f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to evaluate the gamma pass rate, we are uninterested in accurate gamma index values for $\\gamma > 1$.\n",
    "\n",
    "We can suppress calculations of gamma values above 1 by passing a suitable value (e.g., 1.1) to the <code>max_gamma</code> parameter of the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "gamma = pymedphys.gamma(\n",
    "    axes_reference_subset, dose_reference_Z, \n",
    "    axes_evaluation_subset, dose_evaluation_Z, \n",
    "    dose_percent_threshold, distance_mm_threshold,\n",
    "    lower_percent_dose_cutoff=1, # 1% lower threshold\n",
    "    max_gamma=1.1 # stop when gamma > 1.1\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_gamma = gamma[~np.isnan(gamma)]\n",
    "\n",
    "out = plt.hist(valid_gamma, bins=20, density=True)\n",
    "plt.xlabel(\"gamma index\")\n",
    "plt.ylabel(\"probability density\")\n",
    "title = plt.title(f\"Gamma passing rate: {np.sum(valid_gamma <= 1) / len(valid_gamma) * 100:.1f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the gamma passing rate value (shown in the plot title) did not change. Of course now all the points with $\\gamma > 1.1$ are collapsed to the 1.1 bin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, random sampling of a subset of the dose distributions for gamma index calculations may also produce a sufficiently accurate pass rate. This can also introduce significant decreases in calculation time.\n",
    "\n",
    "We can compute the gamma index distribution from a random subset of grid points by providing a suitable fraction to the <code>random_subset</code> parameter of the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "gamma = pymedphys.gamma(\n",
    "    axes_reference_subset, dose_reference_Z, \n",
    "    axes_evaluation_subset, dose_evaluation_Z, \n",
    "    dose_percent_threshold, distance_mm_threshold,\n",
    "    lower_percent_dose_cutoff=1,\n",
    "    random_subset=int(len(dose_reference_Z.flat) // 10) # sample only 1/10 of the grid points\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\gamma$ index distribution should look something like the following plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_gamma = gamma[~np.isnan(gamma)]\n",
    "\n",
    "out = plt.hist(valid_gamma, bins=20, density=True)\n",
    "plt.xlabel(\"gamma index\")\n",
    "plt.ylabel(\"probability density\")\n",
    "title = plt.title(f\"Gamma passing rate: {np.sum(valid_gamma <= 1) / len(valid_gamma) * 100:.1f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the computation time decreased and, despite the the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html) function sampled one 10th of the points, the pass rate deviates by less than 1% from the true value. Therefore, you can use <code>random_subset</code> to get a sense of what could be the pass rate of the dose distributions you are assessing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To conclude, you have performed several computations of the gamma index and the pass rate for 3D and 2D dose distributions using the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html).\n",
    "\n",
    "Moreover, you have obtained fast pass rate computations by tweaking the <code>max_gamma</code> and <code>random_subset</code> parameters of the PyMedPhys [gamma](https://docs.pymedphys.com/users/ref/lib/gamma.html), independently.\n",
    "\n",
    "Of course, both <code>max_gamma</code> and <code>random_subset</code> can be used at the same time. You are now ready to test them by yourself!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.10.4"
  },
  "nbsphinx": {
   "timeout": 600
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}