01-Introduction-to-Xray-attenuation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Session 1\n",
    "\n",
    "## Introduction to X-ray attenuation and its implementation in ![gVXR](img/gvxr_logo.png)\n",
    "\n",
    "## Authors: Franck Vidal and Jean-Michel L&eacute;tang\n",
    "\n",
    "(version 1.0, 22 Sep 2022)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Aims of this session\n",
    "\n",
    "- Explain what gVirtualXray (gVXR) is and why it has been developed;\n",
    "- Introduce projection X-ray imaging and how X-rays are produced;\n",
    "- Understand how X-rays interact with matter;\n",
    "- Become familiar with the Beer-Lambert law to compute the attenuation of X-rays by matter.\n",
    "- Describe how the Beer-Lambert law is implemented in [gVirtualXray](https://gvirtualxray.sourceforge.io/);\n",
    "- Compare images simulated using [gVirtualXray](https://gvirtualxray.sourceforge.io/) with ground truth images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Acknowledgements\n",
    "\n",
    "This notebook re-uses material from Jupyter Book by Jean-Michel L&eacute;tang: [*X-ray imaging: Physics, Instrumentation & Applications*](https://www.creatis.insa-lyon.fr/~letang/x-ray-book/physics.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## What is gVXR and why has it been developed?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Medical Training using VR\n",
    "\n",
    "### Timeline\n",
    "\n",
    "- MSc on *Constructing a GUI using 3D reconstruction for a radiographer's training tool* (2002)\n",
    "- MRes on *Modelling the response of x-ray detectors and removing artefacts in 3D tomography* (2003)\n",
    "- PhD on *Simulation of image guided needle puncture: contribution to real-time ultrasound and fluoroscopic rendering* (start 2003)\n",
    "- Postdoc *Development and validation of a virtual reality simulator for training in interventional radiological visceral needle puncture procedures*;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Trend\n",
    "\n",
    "- Create large amount of realistic simulated X-ray images;\n",
    "- Development of VR apps to train Radiographers or Interventional Radiologists;\n",
    "    - Take X-ray radiographs;\n",
    "    - Stick needles and catheters in the human body whilst looking at images\n",
    "    - Such as fluoroscopy, i.e. \n",
    "    - Real-time X-ray images on a TV screen.\n",
    "- Not interested in scattering;\n",
    "- Noise is not an issue;\n",
    "- Ray-tracing is viable in this context."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Why gVXR?\n",
    "\n",
    "- Simulation of X-Ray attenuation extensively studied in physics;\n",
    "- Different physically-based simulation code available;\n",
    "- Physically-based simulation usually performed using Monte Carlo methods on CPU (often used in dosimetry for radiotherapy);\n",
    "    - Very accurate; but\n",
    "    - Computing an image requires a very very very long time (e.g. days or weeks);\n",
    "- Ray-tracing is an alternative, but\n",
    "    - Still relatively slow on CPU;\n",
    "    - Does not easily take into account scattering;\n",
    "        - Scattering does not necessarily matter in X-ray transmission imaging;\n",
    "    - Does not include Poisson counting noise;\n",
    "        - Poisson noise can be added as a post-process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "![gVXR](img/gvxr_logo.png) is not  a replacement for Monte Carlo methods, it's a complement when a trade-off between speed and accuracy must met."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## In a nutshell, gVXR is (1/2)\n",
    "\n",
    "- an application programming interface (API);\n",
    "- a C++ X-Ray simulation library;\n",
    "- Open-source (with a flexible license):\n",
    "    - Enable reproducible research;\n",
    "    - Can be used in open-source projects;\n",
    "    - Can be used in closed source commercial applications too.\n",
    "- Realtime:\n",
    "    - GPU implementation.\n",
    "- Portable:\n",
    "    - Run on GNU/Linux, Mac OS and Microsoft Windows\n",
    "- Validated:\n",
    "    - against Monte Carlo simulation when scattering can be ignored"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## In a nutshell, gVXR is (2/2)\n",
    "\n",
    "- Flexible\n",
    "    - Available for most popular programming languages (Python3, R, Ruby, Tcl, C#, Java, and GNU Octave);\n",
    "    - Its Python package [\"gVXR\"](https://pypi.org/project/gVXR/) is listed on the [Python Package Index ![Pypi's logo](https://pypi.org/static/images/logo-small.95de8436.svg)](https://pypi.org/)\n",
    "    - Uses polygon meshes to model 3-D geometries from popular file formats (eg. STL, PLY, 3DS, OBJ, DXF, X3D, DAE)\n",
    "\n",
    "![Example of wireframe model](img/wireframe_model2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Simulation supports:\n",
    "\n",
    "- Various source shapes\n",
    "    - Point source;\n",
    "    - Cube source;\n",
    "    - Parallel beam.\n",
    "- Incident beams:        \n",
    "    - Monochromatic;\n",
    "    - Polychromatic.\n",
    "- Geometries:\n",
    "    - Surface meshes (triangles);\n",
    "    - Volume meshes (tetrahedrons). **NEW**\n",
    "- Material properties:\n",
    "    - Chemical elements (e.g. W for tungsten);\n",
    "    - Compounds (e.g. H2O for water); **NEW**\n",
    "    - Mixture (e.g. Ti90/Al6/V4); **NEW**\n",
    "    - Hounsfield units (for medical applications)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## More recently\n",
    "\n",
    "- Improved installation from source code;\n",
    "- Supports [xraylib](https://github.com/tschoonj/xraylib) from [![ESRF](img/ESRF-logo.jpg)](https://www.esrf.eu/) for Photon Cross Sections;\n",
    "- Installation from the Python Package Index using `pip install gvxr`;\n",
    "- Containerisation using ![Docker](img/horizontal-logo-monochromatic-white.png);\n",
    "- Converter from Abaqus files 2 STL files: **EXPERIMENTAL**\n",
    "    - The simulation works with volumetric meshes (tetrahedrons),\n",
    "    - With surface meshes (triangles)\n",
    "- Always working on improvements (when time allows)\n",
    "    - If you want a new feature, just drop me an email or open a [ticket on SourceForge](https://sourceforge.net/p/gvirtualxray/tickets/).\n",
    "- Realistic noise model for synchrotron CT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Fact sheet\n",
    "\n",
    "- SVN repository hosted by [![Sourceforge](img/sourceforge-logo.png)](https://sourceforge.net/)\n",
    "    - [https://sourceforge.net/projects/gvirtualxray/](https://sourceforge.net/projects/gvirtualxray/)\n",
    "    - Registered on 2013-12-01;\n",
    "    - 2298 commits on 2022-08-10 (v2.0.1).\n",
    "- Implemented in ![C++](img/cpp-logo.png) using ![OpenGL](img/OpenGL-logo.png)\n",
    "\n",
    "| Language                            | files      |    blank  |      comment    |       code |\n",
    "|-------------------------------------|------------|-----------|-----------------|------------|\n",
    "|  C/C++ Header                       |     88     |   5585    | 12720           |  38687     |\n",
    "|  C++                                |    117     |  11724    | 14024           |  33256     |\n",
    "|  HTML  |  50  |  592  |  474  |  13100   | \n",
    "|  XML  |  5  |  0  |  0  |  12121   |  \n",
    "| ...  | ... | ... | ... | ... |\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "| Language                            | files      |    blank  |      comment    |       code |\n",
    "|-------------------------------------|------------|-----------|-----------------|------------|\n",
    "| ...  | ... | ... | ... | ... |\n",
    "|  TeX  |                                   17      |      986        |   747        |   5441 |\n",
    "|  CMake  |  69  |  1410  |  1448  |  4490   | \n",
    "|  SWIG  |  2  |  736  |  516  |  4010   | \n",
    "|  Python  |  31  |  1383  |  1232  |  3541   | \n",
    "|  Markdown  |  22  |  767  |  0  |  2155   | \n",
    "|  GLSL  |  56  |  1029  |  2961  |  2060   | \n",
    "|  CSS  |  3  |  307  |  80  |  1481   | \n",
    "|  JavaScript |                             24    |        105         |   210          | 1107 |\n",
    "|  XSD  |  2  |  101  |  10  |  1086   | \n",
    "|  CSV  |  7  |  1  |  0  |  379   | \n",
    "|  make  |  5  |  36  |  18  |  234   | \n",
    "|  Windows Module Definition |               1 |             0     |         0        |    183 |\n",
    "|  Text  |  7  |  12  |  0  |  166   | \n",
    "|  Bourne Shell   |                         17   |          51       |      67     |       128 |\n",
    "| ...  | ... | ... | ... | ... |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "| Language                            | files      |    blank  |      comment    |       code |\n",
    "|-------------------------------------|------------|-----------|-----------------|------------|\n",
    "| ...  | ... | ... | ... | ... |\n",
    "|  Java  |  1  |  15  |  24  |  51   | \n",
    "|  TOML  |  1  |  4  |  2  |  44   | \n",
    "|  Ruby  |  1  |  17  |  26  |  41   | \n",
    "|  C#  |    1  |           12         |    23        |     39 |\n",
    "|  MATLAB  |  1  |  16  |  25  |  39   | \n",
    "|  R  |  1  |  14  |  26  |  35   | \n",
    "|  Perl  |  1  |  15  |  21  |  34   | \n",
    "|  Tcl/Tk         |                          1    |         15        |     24      |       33 |\n",
    "|  DOS Batch                 |               1        |      4        |      0     |        20 |\n",
    "|  XSLT  |  1  |  0  |  5  |  10   | \n",
    "|  **SUM:**  |                                 **533**    |      **24937**       |   **34683**      |   **123971** |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### ![OpenGL is](img/OpenGL-logo.png)\n",
    "\n",
    "- cross-language, \n",
    "- cross-platform application programming interface (API)\n",
    "- for rendering 2D and 3D vector graphics. \n",
    "\n",
    "The API is typically used to interact with a graphics processing unit (GPU), to achieve hardware-accelerated rendering.\n",
    "\n",
    "In other words, takes polygon meshes as an input to produce 2D images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "source": [
    "```bash\n",
    "cloc --report-file=/home/fpvidal/PROGRAMMING/gvirtualxray-code/gvxr-code.txt ~/Downloads/gvirtualxray-code-r2298-trunk.zip\n",
    "\n",
    "sed \"s/\\([A-Z]*[a-z]*\\)[   ]*\\([0-9]*\\)[   ]*\\([0-9]*\\)[   ]*\\([0-9]*\\)[   ]*\\([0-9]*\\)/|  \\1  |  \\2  |  \\3  |  \\4  |  \\5   | /\" < /home/fpvidal/PROGRAMMING/gvirtualxray-code/gvxr-code.txt\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Examples of applications that use gVirtualXRay\n",
    "\n",
    "- [Teaching particle physics to undergraduates](https://doi.org/10.1088/1361-6404/ab5011);\n",
    "- Virtual Reality Simulation;\n",
    "- [Virtual Testing Lab](https://doi.org/10.5281/zenodo.1452506)\n",
    "- [Design new clinical imaging modality to reduce dose exposure](https://doi.org/10.1109/ACCESS.2019.2895925)\n",
    "- [Teaching radiography in medicine](https://doi.org/10.2312/cgvc.20191267)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/sXB-9fG2AbU\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Synchrotron &mu;-tomography\n",
    "\n",
    "| View from Mount Jalla on the ESRF and ILL in Grenoble | Scanned object  |  CT slice |\n",
    "|-------------------------------------------------------|---|---|\n",
    "| ![Source: Photograph by German Wikipedian Christian Hendrich, October 2004](img/Esrf_grenoble.jpg) |  ![](img/sample.png)   |   ![](img/CT_ref_annotated.png)   |\n",
    "\n",
    "- Studied how/why artefacts in some CT volumes occurred.\n",
    "- \"Investigation of artefact sources in synchrotron microtomography via virtual X-ray imaging\". Vidal, F.P.; L&eacute;tang, J.M.; Peix, G.; Cloetens, P. In: *Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms*, Vol. 234, No. 3, 01.06.2005, p. 333-348. [DOI: 10.1016/j.nimb.2005.02.003](https://doi.org/10.1016/j.nimb.2005.02.003)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Use of fast realistic simulations on GPU to extract CAD models from microtomographic data in the presence of strong CT artefacts\n",
    "\n",
    "Vidal, Franck; Mitchell, Iwan; Letang, J.M. In: *Precision Engineering*, Vol. 74, 01.03.2022, p. 110-125. [DOI: 10.1016/j.precisioneng.2021.10.014](https://doi.org/10.1016/j.precisioneng.2021.10.014)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Registration of the TiAlV matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<video width=\"320\" height=\"240\" controls>\n",
    "  <source src=\"videos/1-s2.0-S0141635921002609-mmc7.mp4\" type=\"video/mp4\">\n",
    "</video>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Evolution of the core and fibre radii during the optimisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<video width=\"320\" height=\"240\" controls>\n",
    "  <source src=\"videos/1-s2.0-S0141635921002609-mmc6.mp4\" type=\"video/mp4\">\n",
    "</video>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Evolution of beam spectrum parameters during the optimisation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<video width=\"320\" height=\"240\" controls>\n",
    "  <source src=\"videos/1-s2.0-S0141635921002609-mmc2.mp4\" type=\"video/mp4\">\n",
    "</video>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Evolution of the core and fibre radii and phase contrast parameters during the optimisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<video width=\"320\" height=\"240\" controls>\n",
    "  <source src=\"videos/1-s2.0-S0141635921002609-mmc4.mp4\" type=\"video/mp4\">\n",
    "</video>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Principles of rojection X-ray imaging\n",
    "\n",
    "![](img/dragon2.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## X-ray production: X-ray tube\n",
    "\n",
    "| X-ray tube principle: William Coolidge in 1913 | Old X-ray tube (glass wall chamber) |\n",
    "|-----------------------------------------------------------------------|-----------|\n",
    "| ![X-ray tube principle: William Coolidge in 1913](img/x-ray_tube.jpg) | ![](img/GE-D17.jpg) |\n",
    "\n",
    "- Electrons in a heated filament (with a current of a few mA typically)\n",
    "- are extracted using a very high voltage (several tens of kV typically).\n",
    "- Their kinetic energy is close to 0 in the filament (they can be considered at rest).\n",
    "- Therefore their kinetic energy reach keV just before impinging onto the anode target.\n",
    "- From electron-matter interactions in the target follow X-ray mostly produced by Bremsstrahlung but also to a lesser extent Fluorescence. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "# X-ray production: Synchrotron\n",
    "\n",
    "| ESRF, Grenoble, France | Diamond Light Source, Didcot, UK |\n",
    "|------------------------|----------------------------------|\n",
    "| ![ESRF](img/Esrf_grenoble.jpg) |  ![](img/diamond.jpg) |\n",
    "| Source: Photograph by German Wikipedian Christian Hendrich, October 2004 | Copyright of Diamond Light Source Ltd |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## How do X-rays interact with matter?\n",
    "\n",
    "\n",
    "- X-photons cross matter;\n",
    "- During their path into any material, they can interact with matter.\n",
    "\n",
    "![*Illustration of X-ray photon/matter interaction*](img/interactions.png)\n",
    "\n",
    "1. Directly transmitted photons (no interaction);\n",
    "2. Absorbed photons;\n",
    "3. Scattered photons;\n",
    "4. Absorbed scattered photons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "**For most X-rays imaging modalities, only directly transmitted photons are essential:**\n",
    "\n",
    "- Scattered photons decrease the image quality;\n",
    "- Absorbed photons do not reach the detector;\n",
    "- Scattered photons may be ignored (but not necessarily).\n",
    "\n",
    "![*Illustration of X-ray photon/matter interaction*](img/interactions.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Beer-Lambert law to compute the attenuation of X-rays by matter\n",
    "\n",
    "- also known as (Attenuation Law)\n",
    "\n",
    "![*Illustration of the Beer-Lambert law*](img/attenuation.png)\n",
    "\n",
    "$$N_{out}(E) = N_{in}(E) \\times \\mathrm{e}^{-\\sum_i \\mu_i(E, \\rho, Z) L_p(i)}$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "$$N_{out}(E) = N_{in}(E) \\times \\mathrm{e}^{-\\sum_i \\mu_i(E, \\rho, Z) L_p(i)}$$  \n",
    "\n",
    "- $N_{in}(E)$ the number of incident photons at energy $E$;\n",
    "- $N_{out}(E)$ the number of transmitted photons of energy $E$;\n",
    "- $\\mu_i$ the linear attenuation coefficient (in cm<sup>-1</sup>) of the $i$<sup>th</sup> object. It depends on:\n",
    "    - $E$ the energy of incident photons;\n",
    "    - $\\rho$ the material density of the object;\n",
    "    - $Z$ the atomic number of the object material.\n",
    "    - $L_p(i)$ the path length of the ray in the $i$<sup>th</sup> object.\n",
    "- $E_{out} = N_{out}(E) \\times E$\n",
    "    - $E_{out}$ the energy received by the pixel, i.e. as recorded in the X-ray image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Linear attenuation coefficient\n",
    "\n",
    "The linear [attenuation coefficient](https://www.creatis.insa-lyon.fr/~letang/x-ray-book/physics.html#attenuation-coefficient) $\\mu$ appears in the Beer-Lambert [Attenuation law](https://www.creatis.insa-lyon.fr/~letang/x-ray-book/physics.html#attenuation-law) which gives the number of directly transmitted photons $N_{out}$ (i.e. without interaction) in terms of the number of incident photons $N_{in}$:\n",
    "$$\n",
    "N_{\\mathrm{dt}}(E)=N_0(E)\\exp(-\\mu(E)X)\n",
    "$$\n",
    "where $X$ is the thickness of the traversed material. This expression is only valid for photons which have the same energy $E$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "import utilities\n",
    "import ipywidgets as widgets\n",
    "widgets.interact_manual(utilities.mu,material=widgets.Dropdown(options=[('Polyethylene','H2C'),('Water','H2O'),('Aluminium','Al'),('Copper','Cu'),('Yttrium','Y'),('Tin','Sn'),('Lead','Pb')],value=\"H2C\",layout={'width': 'max-content'},description='Material:',style={'description_width': 'initial'}));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Attenuation of 500 photons of 100 keV by 1cm of copper\n",
    "\n",
    "In the previous cell, \n",
    "1. Select copper as material\n",
    "2. Click on `Run Interact`\n",
    "3. Using the mouse cursor, find what the total linear attenuation coefficient of copper is for 100 keV,\n",
    "4. Replace `???` in the cell below with the corresponding value\n",
    "5. In the cell below, compute $$N_{\\mathrm{dt}}(E)=N_0(E)\\exp(-\\mu(E)X)$$\n",
    "    - with $E$ = 100 keV,\n",
    "    - $N_0(E)$ = 500\n",
    "    - $X$ = 1 and\n",
    "    - $\\mu(E)$ the value you wrote down"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "mu_Cu = ??? # Linear attenuation of copper at 100 keV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "N_0 = 500 # Number of incident photons\n",
    "X_Cu = 1 # Thickness in cm\n",
    "\n",
    "N_dt = N_0 * math.exp(-mu_Cu * X_Cu)\n",
    "print(\"Out of\", N_0, \"photons, only\", round(N_dt), \"are transmitted\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Attenuation of 500 photons of 100 keV by ??cm of water\n",
    "\n",
    "- We know that $$N_{\\mathrm{Cu}}=500 \\times \\exp(-\\mu_{Cu} \\times 1) = 8$$\n",
    "- We want to illustrate than water does not attenuate X-rays of 100 keV as much as copper.\n",
    "- Repeat some of the previous steps to find the what the total linear attenuation coefficient of water is for 100 keV.\n",
    "- The Beer-Lambert law in this case is:\n",
    "$$N_{\\mathrm{H_2O}}= 500 \\times \\exp(-\\mu_{H_2O} \\times x) = 8$$\n",
    "- What is the value of $x$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "X_H2O = None # Thickness in cm, it is unknown\n",
    "mu_H2O = ??? # Linear attenuation of water at 100 keV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "$$500 \\times \\exp(-\\mu_{H_2O} \\times x_{H_2O}) = N_{\\mathrm{H_2O}} = 8$$\n",
    "$$\\exp(-\\mu_{H_2O} \\times x_{H_2O}) = \\frac{N_{\\mathrm{H_2O}}}{500}$$\n",
    "$$-\\mu_{H_2O} \\times x_{H_2O} = \\ln\\left(\\frac{N_{\\mathrm{H_2O}}}{500}\\right)$$\n",
    "$$\\mu_{H_2O} \\times x_{H_2O} = -\\ln\\left(\\frac{N_{\\mathrm{H_2O}}}{500}\\right)$$\n",
    "$$x_{H_2O} = \\frac{-\\ln\\left(\\frac{N_{\\mathrm{H_2O}}}{500}\\right)}{\\mu_{H_2O}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_H2O = (-math.log(N_dt / 500)) / mu_H2O\n",
    "print(round(X_H2O), \"cm of water is needed to stop as many photons of 100 keV as\", X_Cu,\"cm of copper\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Check if we made a mistake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_dt = N_0 * math.exp(-mu_H2O * X_H2O)\n",
    "print(\"Out of\", N_0, \"photons, only\", round(N_dt), \"are transmitted\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## What does it have to do with gVirtualXray?\n",
    "\n",
    "**Given $E$, $L_p$ and $\\mu$, we can compute X-ray images.**\n",
    "\n",
    "- Computational bottlenecks are:\n",
    "    - computing $L_p$ from polygon meshes\n",
    "    - the exponential,\n",
    "    - the accumulation $(\\sum)$\n",
    "- gVirtualXray provides real-time mechanisms to compute them on GPU."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Path Length: Naive Approach\n",
    "\n",
    "![*Is finding intersections in the right order important?*](img/intersections.png)\n",
    "\n",
    "1. Detect every intersection between a ray and the objects;\n",
    "2. Sort intersection (can be handled by GPUs using depth-peeling, a multi-pass rendering technique for semi-transparent polygonal objects without sorting polygons);\n",
    "3. Compute path length."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Path Length: L-Buffer\n",
    "\n",
    "![*Finding intersections in any order does not matter*](img/l-buffer.png)\n",
    "\n",
    "**Intersection sorting is actually not needed!**\n",
    "\n",
    "- By convention normals are outward;\n",
    "- A ray penetrates into an object when the dot product between the view vector (V) and the normal (N<sub>i</sub>) at the intersection point is positive;\n",
    "- It leaves an object when the dot product is negative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## L-Buffer Implementation\n",
    "\n",
    "L<sub>p</sub>=&Sigma;<sub>i</sub> - sng(V &middot; N<sub>i</sub>) x d<sub>i</sub>\n",
    "\n",
    "\n",
    "- i refers to i<sup>th</sup> intersection in an arbitrary order;\n",
    "- d<sub>i</sub> distance from X-ray source to intersection point;\n",
    "- sgn(V &middot; N<sub>i</sub>) stands for the sign of the dot product between V and N<sub>i</sub>;\n",
    "- In a shader program [^1] , compute:\n",
    "    - sgn(V &middot; N<sub>i</sub>);\n",
    "    - d<sub>i</sub> the distance between the X-ray source and the intersection;\n",
    "    - Assign -sng(V &middot; N<sub>i</sub>) x d<sub>i</sub> as the fragment value.\n",
    "- For each pixel, compute L<sub>p</sub> thanks to high-dynamic range and OpenGL blending function (pixel values may not be between 0 and 1)[^2].\n",
    "\n",
    "[^1]: In OpenGL a shader program is a small program that is usually run on a GPU in parallel for every pixel of an image.\n",
    "\n",
    "[^2]: See [DOI: 10.2312/LocalChapterEvents/TPCG/TPCG09/025-032](http://dx.doi.org/10.2312/LocalChapterEvents/TPCG/TPCG09/025-032) for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Multipass Rendering Pipeline\n",
    "\n",
    "pixel = E x N<sub>out</sub>\n",
    "\n",
    "pixel =  <span style=\"color:green\">E x N<sub>in</sub>(E) e</span><sup>(<span style=\"color:red\">-&Sigma;<sub>i</sub> &mu;<sub>i</sub> </span> <span style=\"color:blue\">  L<sub>p</sub>(i)</span>)</sup>\n",
    "\n",
    "The Beer-Lambert law is broken down in several rendering passes:\n",
    "\n",
    "<!-- - Needs 3 FBOs with high-dynamic range capability for off-line rendering:\n",
    " -->\n",
    "\n",
    "- For each object of the scene:\n",
    "  1. Compute L<sub>p</sub>(i);\n",
    "  2. Update results of &Sigma; &mu;<sub>i</sub> L<sub>p</sub>(i).\n",
    "- For the final image only:\n",
    "  1. Compute N<sub>out</sub>;\n",
    "  2. (Optional when direct display only is needed)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "![*OpenGL pipeline to compute the Beer-Lambert law (monochromatic case).*](img/pipeline1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Quantitative Validation\n",
    "\n",
    "- Simulating an image relies on a Beer-Lambert law implementation;\n",
    "- Solving the Beer-Lambert law relies on Linear Attenuation Coefficients; (&mu;)\n",
    "<!-- - &mu; is not known for given incident energies; -->\n",
    "- &mu; is computed using Mass Attenuation Coefficients (&mu;/&rho;) and material density (&rho;).\n",
    "- **Are the values used in gVirtualXRay accurate?**\n",
    "    - Compare values computed in gVirtualXRay with those from the literature.\n",
    "- **Are the Beer-Lambert law implementations accurate?**\n",
    "    - Compare values computed in gVirtualXRay with theoretical ones.\n",
    "- **Are the simulated images accurate?**\n",
    "    - Compare images computed using gVirtualXRay with those using a state-of-the-art Monte Carlo software, e.g.\n",
    "\n",
    "    [![Geant4](./img/g4logo-web.png)](https://geant4.web.cern.ch/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "source": [
    "## Mass Attenuation Coefficients (medical)\n",
    "\n",
    "- Any tissue can be described by its Hounsfiled Unit (HU):\n",
    "  - $HU(material) = 1000 \\times \\frac{\\mu(material) - \\mu_{H2O}}{\\mu_{H2O}}$\n",
    "- Given a HU value for any simulated object;\n",
    "  - If $\\mu_{H2O}$ is known for any energy,\n",
    "  - then $\\mu(material)$ for any HU and for any energy can be computed:\n",
    "      - $\\mu(material) = \\frac{HU(material) \\times \\mu_{H2O}}{1000} + \\mu_{H2O}$\n",
    "- Mass attenuation coefficients ($\\frac{\\mu}{\\rho}$) for various human tissues can be found in the literature;\n",
    "- The density ($\\rho$) for various human tissues can be found in the literature.\n",
    "- **Linear attenuation coefficients can therefore be computed for various human tissues** and\n",
    "- **used to solve the Beer-Lambert law**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "source": [
    "## Mass Attenuation Coefficients: Tissue, Soft (ICRU-44)\n",
    "\n",
    "- ![*Image from [https://physics.nist.gov/PhysRefData/XrayMassCoef/ComTab/tissue.html](https://physics.nist.gov/PhysRefData/XrayMassCoef/ComTab/tissue.html)*](./img/tissue.png)\n",
    "\n",
    "- &mu;(water) from literature is provided at given energies only: Interpolation needed for missing energies\n",
    "\n",
    "  - ![*Image computed with gVirtualXRay using linear interpolation*](./img/test_soft_tissue_1.png)\n",
    "\n",
    "  - ![*Image computed with gVirtualXRay using interpolation in the log scale*](./img/test_soft_tissue_2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "source": [
    "## Mass Attenuation Coefficients: Bone, Cortical (ICRU-44)\n",
    "\n",
    "- ![*Image from [https://physics.nist.gov/PhysRefData/XrayMassCoef/ComTab/bone.html](https://physics.nist.gov/PhysRefData/XrayMassCoef/ComTab/bone.html)*](img/bone.png)\n",
    "\n",
    "- &mu;(water) from literature is provided at given energies only: Interpolation needed for missing energies\n",
    "\n",
    "  - ![*Image computed with gVirtualXRay using linear interpolation*](./img/test_bone_1.png)\n",
    "\n",
    "  - ![*Image computed with gVirtualXRay using interpolation in the log scale*](./img/test_bone_2.png)\n",
    "\n",
    "**Not a good match as the peaks are not visible in &mu;(water)**\n",
    "\n",
    "- ![*Image computed with gVirtualXRay using NIST's XCOM database*](./img/test_bone_3.png)\n",
    "\n",
    "**This is now a good match**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Mass attenuation coefficients from elements\n",
    "\n",
    "![](img/mass_attenuation_coefficients_from_elements.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Mass attenuation coefficients from mixtures\n",
    "\n",
    "![](img/mass_attenuation_coefficients_from_ICRU44.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Going back to previous slide\n",
    "\n",
    "- Simulating an image relies on a Beer-Lambert law implementation;\n",
    "- Solving the Beer-Lambert law relies on Linear Attenuation Coefficients; (&mu;)\n",
    "- &mu; is not known for given incident energies;\n",
    "- &mu; is computed using Mass Attenuation Coefficients (&mu;/&rho;) and material density (&rho;).\n",
    "- **Are the values used in gVirtualXRay accurate?**\n",
    "    - Compare values computed in gVirtualXRay with those from the literature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Going back to previous slide\n",
    "\n",
    "- Simulating an image relies on a Beer-Lambert law implementation;\n",
    "- Solving the Beer-Lambert law relies on Linear Attenuation Coefficients; (&mu;)\n",
    "- &mu; is not known for given incident energies;\n",
    "- &mu; is computed using Mass Attenuation Coefficients (&mu;/&rho;) and material density (&rho;).\n",
    "- **Are the values used in gVirtualXRay accurate?**\n",
    "    - Compare values computed in gVirtualXRay with those from the literature.\n",
    "    - **YES**\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## More testing\n",
    "\n",
    "- **Are the Beer-Lambert law implementations accurate?**\n",
    "    - Compare values computed in gVirtualXRay with theoretical ones.\n",
    "- **Are the simulated images accurate?**\n",
    "    - Compare images computed using gVirtualXRay with those using a state-of-the-art Monte Carlo software, e.g.\n",
    "\n",
    "    [![Geant4](img/g4logo-web.png)](https://geant4.web.cern.ch/)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Test case\n",
    "\n",
    "![*Simulated object*](./img/validation_3D_scene.png)\n",
    "\n",
    "- Cube: edge length of 3 cm, made of soft tissue (HU = 52).\n",
    "- Cylinder: height of 3 m, diameter of 2 cm, made of bone (HU = 1330).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Gate vs. gVirtualXRay\n",
    "\n",
    "We simulate a test case twice:\n",
    "\n",
    "- Using a Monte Carlo method for particle physics implemented in [GATE](http://www.opengatecollaboration.org/);\n",
    "- Using our GPU implementation.\n",
    "\n",
    "*GATE is an opensource software developed by an international collaboration. Its focus is on Monte Carlo simulation in medical imaging and radiotherapy. GATE makes use of the Geant4 libraries. Geant 4 is CERN's Monte Carlo simulation platform dedicated to particle physics in nuclear research. CERN is the European Organization for Nuclear Research.*\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Gate vs. gVirtualXRay\n",
    "\n",
    "| Simulation parameters | Image computed with GATE (2 weeks of computations on supercomputer) | Image computed with gVirtualXRay (less than 1 sec. of computations on GPU) |\n",
    "|---|---|---|\n",
    "| ![*Simulation parameters*](./img/POINT_SOURCE/scene.png) | ![*Image computed with GATE (2 weeks of computations on supercomputer)*](./img/POINT_SOURCE/gate_norm.png) | ![*Image computed with gVirtualXRay (less than 1 sec. of computations on GPU)*](./img/POINT_SOURCE/gpu_norm.png) |\n",
    "\n",
    "Normalised cross-correlation (NCC) = 99.747&percnt;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "![*Profiles*](img/POINT_SOURCE/test_profile.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Real digital radiograph (DR) taken with a clinically utilised X-ray equipment vs. gVirtualXRay\n",
    "\n",
    "| Lungman phantom | 3D model of the phantom, source and detector|\n",
    "|-----------------|----------| \n",
    "| ![](img/lungman.png) | ![](img/lungman_model.png) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "![](./img/lungman-projection-harder.png)\n",
    "\n",
    "![](img/lungman-profiles-projection-postprocessing.png)\n",
    "\n",
    "Normalised cross-correlation (NCC) = 98.92&percnt;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Acknowledgements\n",
    "\n",
    "This research and/or training day would not have been possible without Jean Michel L&eacute;tang, Jean Yves Buffi&egrave;re, Nicolas Freud, Nigel John, Tianci Wen, Iwan Mitchell, Jammie Pointon, Llion Evans, Ben Thorpes, David Fairbrother, ..."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "585e9fc9b818a92a95ff682530f87f6de2c4282b8239c594bde282f49476a8c7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}