{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Effects of Noise on Gamma\n",
    "\n",
    "Here is an example of the effects noise can have on gamma. \n",
    "Simply by adding noise to the evaluation distribution the Gamma **goes from a passing rate of about `50%` to about `99%`**.\n",
    "For a Monte Carlo planning system this dose distribution includes statistical noise potentially pushing the Gamma pass rate artificially higher than should that noise be absent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install pymedphys\n",
    "\n",
    "import pymedphys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_options = {\n",
    "    'dose_percent_threshold': 3,\n",
    "    'distance_mm_threshold': 3,\n",
    "    'lower_percent_dose_cutoff': 20,\n",
    "    'interp_fraction': 10,  # Should be 10 or more for more accurate results\n",
    "    'max_gamma': 2,\n",
    "    'random_subset': None,\n",
    "    'local_gamma': True,\n",
    "    'ram_available': 2**29,  # 1/2 GB\n",
    "    'quiet': True\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = 0.5\n",
    "scale_factor = 1.035\n",
    "noise = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin = -28\n",
    "xmax = 28\n",
    "ymin = -25\n",
    "ymax = 25\n",
    "\n",
    "extent = [xmin-grid/2, xmax+grid/2, ymin-grid/2, ymax+grid/2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining an example dose reference\n",
    "\n",
    "Here we are defining an idealised square field using an exponential function. We are also creating a `coords` variable which we will be passing to the `pymedphys.gamma` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(xmin, xmax + grid, grid)\n",
    "y = np.arange(ymin, ymax + grid, grid)\n",
    "\n",
    "coords = (y, x)\n",
    "\n",
    "xx, yy = np.meshgrid(x, y)\n",
    "dose_ref = np.exp(-((xx/15)**20 + (yy/15)**20))\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Reference dose')\n",
    "\n",
    "plt.imshow(\n",
    "    dose_ref, clim=(0, 1.04), extent=extent)\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of importance the length of the first coordinate set in `coords` matches the first dimension of `dose_ref`, and the length of the second coordinate set in `coords` matches the second dimension of `dose_ref`. This is required because `pymedphys.gamma` internally uses [`scipy.interpolate.RegularGridInterpolator`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html) and the gamma implimentation has been chosen to align with this scipy coordinate convention uses here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(coords[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(coords[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(dose_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dimensions_of_dose_ref = np.shape(dose_ref)\n",
    "assert dimensions_of_dose_ref[0] == len(coords[0])\n",
    "assert dimensions_of_dose_ref[1] == len(coords[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining an example evaluation dose\n",
    "\n",
    "Let's scale our reference dose by a scaling factor. We have chosen above to scale by 3.5% purposefully so that the dose will go larger than the dose percent evaluation criterion of 3% that was chosen above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dose_eval = dose_ref * scale_factor\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Evaluation dose')\n",
    "\n",
    "plt.imshow(\n",
    "    dose_eval, clim=(0, 1.04), extent=extent)\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Seeing a dose difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dose_diff = dose_eval - dose_ref\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Dose Difference')\n",
    "\n",
    "plt.imshow(\n",
    "    dose_diff, \n",
    "    clim=(-0.1, 0.1), extent=extent,\n",
    "    cmap='seismic'\n",
    ")\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation Gamma\n",
    "\n",
    "Now lets evaluate gamma for the distributions defined above. As expected, due to the dose being scaled larger than the dose threshold, quite a few points fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_no_noise = pymedphys.gamma(\n",
    "    coords, dose_ref,\n",
    "    coords, dose_eval,\n",
    "    **gamma_options)\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Gamma Distribution')\n",
    "\n",
    "plt.imshow(\n",
    "    gamma_no_noise, clim=(0, 2), extent=extent,\n",
    "    cmap='coolwarm')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "valid_gamma_no_noise = gamma_no_noise[~np.isnan(gamma_no_noise)]\n",
    "no_noise_passing = 100 * np.round(np.sum(valid_gamma_no_noise <= 1) / len(valid_gamma_no_noise), 4)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(f'Gamma Histogram | Passing rate = {round(no_noise_passing, 2)}%')\n",
    "plt.xlabel('Gamma')\n",
    "plt.ylabel('Number of pixels')\n",
    "\n",
    "plt.hist(valid_gamma_no_noise, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigating the effect of noise\n",
    "\n",
    "Let's now slightly adjust our evaluation distribution by adding some random normal noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "dose_eval_noise = (\n",
    "    dose_eval + \n",
    "    np.random.normal(loc=0, scale=noise, size=np.shape(dose_eval))\n",
    ")\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Evaluation dose with Noise')\n",
    "\n",
    "plt.imshow(\n",
    "    dose_eval_noise, clim=(0, 1.04), extent=extent)\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what our dise difference looks like with this noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dose_diff_with_noise = dose_eval_noise - dose_ref\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Dose Difference with Noise')\n",
    "\n",
    "plt.imshow(\n",
    "    dose_diff_with_noise, \n",
    "    clim=(-0.1, 0.1), extent=extent,\n",
    "    cmap='seismic'\n",
    ")\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But what about our gamma?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_with_noise = pymedphys.gamma(\n",
    "    coords, dose_ref,\n",
    "    coords, dose_eval_noise,\n",
    "    **gamma_options)\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Gamma Distribution')\n",
    "\n",
    "plt.imshow(\n",
    "    gamma_with_noise, clim=(0, 2), extent=extent,\n",
    "    cmap='coolwarm')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()\n",
    "valid_gamma_with_noise = gamma_with_noise[~np.isnan(gamma_with_noise)]\n",
    "with_noise_passing = 100 * np.round(np.sum(valid_gamma_with_noise <= 1) / len(valid_gamma_with_noise), 4)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(f'Gamma Histogram | Passing rate = {round(with_noise_passing, 2)}%')\n",
    "plt.xlabel('Gamma')\n",
    "plt.ylabel('Number of pixels')\n",
    "\n",
    "plt.hist(valid_gamma_with_noise, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "So, `pymedphys` does provide a gamma implementation, but should it be a defacto standard of validation? I highly recommend giving the following paper a read to highlight some other tools we can use to augment our usage of Gamma.\n",
    "\n",
    "> [Evaluating IMRT and VMAT dose accuracy: Practical examples of failure to detect systematic errors when applying a commonly used metric and action levels](http://download.xuebalib.com/xuebalib.com.42814.pdf)\n",
    "\n",
    "A key take home message from that paper relevant to this `pymedphys.gamma` utility is the following:\n",
    "\n",
    "> Using 3%G/3 mm gamma passing rates as a basis of gauging “sufficient accuracy” as per TG-119 has knock-on effects\n",
    "outside of commissioning. One such effect is that it is often\n",
    "used as the de facto standard to “validate” a new or modified\n",
    "product. For example, even very recently 3%G/3 mm passing rates have been quoted to prove the accuracy of TPS dose\n",
    "algorithms, delivery systems, or dosimetry devices.\n",
    ">\n",
    "> In this work, we show a variety of clear examples of the\n",
    "metric’s insensitivity, supporting the argument that it is inappropriate as a sole basis for commissioning. The same conclusion can be made for the yet further upstream task of product\n",
    "validation by medical device manufacturers and their clinical\n",
    "partners. "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}