- Explain what gVirtualXray (gVXR) is and why it has been developed;
- Introduce projection X-ray imaging and how X-rays are produced;
- Understand how X-rays interact with matter;
- Become familiar with the Beer-Lambert law to compute the attenuation of X-rays by matter.
- Describe how the Beer-Lambert law is implemented in gVirtualXray;
- Compare images simulated using gVirtualXray with ground truth images.

This notebook re-uses material from Jupyter Book by Jean-Michel Létang: *X-ray imaging: Physics, Instrumentation & Applications*.

- MSc on
*Constructing a GUI using 3D reconstruction for a radiographer's training tool*(2002) - MRes on
*Modelling the response of x-ray detectors and removing artefacts in 3D tomography*(2003) - PhD on
*Simulation of image guided needle puncture: contribution to real-time ultrasound and fluoroscopic rendering*(start 2003) - Postdoc
*Development and validation of a virtual reality simulator for training in interventional radiological visceral needle puncture procedures*;

- Create large amount of realistic simulated X-ray images;
- Development of VR apps to train Radiographers or Interventional Radiologists;
- Take X-ray radiographs;
- Stick needles and catheters in the human body whilst looking at images
- Such as fluoroscopy, i.e.
- Real-time X-ray images on a TV screen.

- Not interested in scattering;
- Noise is not an issue;
- Ray-tracing is viable in this context.

- Simulation of X-Ray attenuation extensively studied in physics;
- Different physically-based simulation code available;
- Physically-based simulation usually performed using Monte Carlo methods on CPU (often used in dosimetry for radiotherapy);
- Very accurate; but
- Computing an image requires a very very very long time (e.g. days or weeks);

- Ray-tracing is an alternative, but
- Still relatively slow on CPU;
- Does not easily take into account scattering;
- Scattering does not necessarily matter in X-ray transmission imaging;

- Does not include Poisson counting noise;
- Poisson noise can be added as a post-process.

- an application programming interface (API);
- a C++ X-Ray simulation library;
- Open-source (with a flexible license):
- Enable reproducible research;
- Can be used in open-source projects;
- Can be used in closed source commercial applications too.

- Realtime:
- GPU implementation.

- Portable:
- Run on GNU/Linux, Mac OS and Microsoft Windows

- Validated:
- against Monte Carlo simulation when scattering can be ignored

- Flexible
- Available for most popular programming languages (Python3, R, Ruby, Tcl, C#, Java, and GNU Octave);
- Its Python package "gVXR" is listed on the Python Package Index
- Uses polygon meshes to model 3-D geometries from popular file formats (eg. STL, PLY, 3DS, OBJ, DXF, X3D, DAE)

- Various source shapes
- Point source;
- Cube source;
- Parallel beam.

- Incident beams:
- Monochromatic;
- Polychromatic.

- Geometries:
- Surface meshes (triangles);
- Volume meshes (tetrahedrons).
**NEW**

- Material properties:
- Chemical elements (e.g. W for tungsten);
- Compounds (e.g. H2O for water);
**NEW** - Mixture (e.g. Ti90/Al6/V4);
**NEW** - Hounsfield units (for medical applications).

- Improved installation from source code;
- Supports xraylib from for Photon Cross Sections;
- Installation from the Python Package Index using
`pip install gvxr`

; - Containerisation using ;
- Converter from Abaqus files 2 STL files:
**EXPERIMENTAL**- The simulation works with volumetric meshes (tetrahedrons),
- With surface meshes (triangles)

- Always working on improvements (when time allows)
- If you want a new feature, just drop me an email or open a ticket on SourceForge.

- Realistic noise model for synchrotron CT.

- SVN repository hosted by
- https://sourceforge.net/projects/gvirtualxray/
- Registered on 2013-12-01;
- 2298 commits on 2022-08-10 (v2.0.1).

- Implemented in using

Language | files | blank | comment | code |
---|---|---|---|---|

C/C++ Header | 88 | 5585 | 12720 | 38687 |

C++ | 117 | 11724 | 14024 | 33256 |

HTML | 50 | 592 | 474 | 13100 |

XML | 5 | 0 | 0 | 12121 |

... | ... | ... | ... | ... |

Language | files | blank | comment | code |
---|---|---|---|---|

... | ... | ... | ... | ... |

TeX | 17 | 986 | 747 | 5441 |

CMake | 69 | 1410 | 1448 | 4490 |

SWIG | 2 | 736 | 516 | 4010 |

Python | 31 | 1383 | 1232 | 3541 |

Markdown | 22 | 767 | 0 | 2155 |

GLSL | 56 | 1029 | 2961 | 2060 |

CSS | 3 | 307 | 80 | 1481 |

JavaScript | 24 | 105 | 210 | 1107 |

XSD | 2 | 101 | 10 | 1086 |

CSV | 7 | 1 | 0 | 379 |

make | 5 | 36 | 18 | 234 |

Windows Module Definition | 1 | 0 | 0 | 183 |

Text | 7 | 12 | 0 | 166 |

Bourne Shell | 17 | 51 | 67 | 128 |

... | ... | ... | ... | ... |

Language | files | blank | comment | code |
---|---|---|---|---|

... | ... | ... | ... | ... |

Java | 1 | 15 | 24 | 51 |

TOML | 1 | 4 | 2 | 44 |

Ruby | 1 | 17 | 26 | 41 |

C# | 1 | 12 | 23 | 39 |

MATLAB | 1 | 16 | 25 | 39 |

R | 1 | 14 | 26 | 35 |

Perl | 1 | 15 | 21 | 34 |

Tcl/Tk | 1 | 15 | 24 | 33 |

DOS Batch | 1 | 4 | 0 | 20 |

XSLT | 1 | 0 | 5 | 10 |

SUM: |
533 |
24937 |
34683 |
123971 |

- cross-language,
- cross-platform application programming interface (API)
- for rendering 2D and 3D vector graphics.

The API is typically used to interact with a graphics processing unit (GPU), to achieve hardware-accelerated rendering.

In other words, takes polygon meshes as an input to produce 2D images.

View from Mount Jalla on the ESRF and ILL in Grenoble | Scanned object | CT slice |
---|---|---|

- Studied how/why artefacts in some CT volumes occurred.
- "Investigation of artefact sources in synchrotron microtomography via virtual X-ray imaging". Vidal, F.P.; Lé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

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

X-ray tube principle: William Coolidge in 1913 | Old X-ray tube (glass wall chamber) |
---|---|

- Electrons in a heated filament (with a current of a few mA typically)
- are extracted using a very high voltage (several tens of kV typically).
- Their kinetic energy is close to 0 in the filament (they can be considered at rest).
- Therefore their kinetic energy reach keV just before impinging onto the anode target.
- From electron-matter interactions in the target follow X-ray mostly produced by Bremsstrahlung but also to a lesser extent Fluorescence.

ESRF, Grenoble, France | Diamond Light Source, Didcot, UK |
---|---|

Source: Photograph by German Wikipedian Christian Hendrich, October 2004 | Copyright of Diamond Light Source Ltd |

- X-photons cross matter;
- During their path into any material, they can interact with matter.

- Directly transmitted photons (no interaction);
- Absorbed photons;
- Scattered photons;
- Absorbed scattered photons.

**For most X-rays imaging modalities, only directly transmitted photons are essential:**

- Scattered photons decrease the image quality;
- Absorbed photons do not reach the detector;
- Scattered photons may be ignored (but not necessarily).

- also known as (Attenuation Law)

$$N_{out}(E) = N_{in}(E) \times \mathrm{e}^{-\sum_i \mu_i(E, \rho, Z) L_p(i)}$$

- $N_{in}(E)$ the number of incident photons at energy $E$;
- $N_{out}(E)$ the number of transmitted photons of energy $E$;
- $\mu_i$ the linear attenuation coefficient (in cm
^{-1}) of the $i$^{th}object. It depends on:- $E$ the energy of incident photons;
- $\rho$ the material density of the object;
- $Z$ the atomic number of the object material.
- $L_p(i)$ the path length of the ray in the $i$
^{th}object.

- $E_{out} = N_{out}(E) \times E$
- $E_{out}$ the energy received by the pixel, i.e. as recorded in the X-ray image.

The linear attenuation coefficient $\mu$ appears in the Beer-Lambert 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_{\mathrm{dt}}(E)=N_0(E)\exp(-\mu(E)X) $$ where $X$ is the thickness of the traversed material. This expression is only valid for photons which have the same energy $E$.

In [ ]:

```
%matplotlib widget
import utilities
import ipywidgets as widgets
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'}));
```

In the previous cell,

- Select copper as material
- Click on
`Run Interact`

- Using the mouse cursor, find what the total linear attenuation coefficient of copper is for 100 keV,
- Replace
`???`

in the cell below with the corresponding value - In the cell below, compute $$N_{\mathrm{dt}}(E)=N_0(E)\exp(-\mu(E)X)$$
- with $E$ = 100 keV,
- $N_0(E)$ = 500
- $X$ = 1 and
- $\mu(E)$ the value you wrote down

In [ ]:

```
mu_Cu = ??? # Linear attenuation of copper at 100 keV
```

In [ ]:

```
N_0 = 500 # Number of incident photons
X_Cu = 1 # Thickness in cm
N_dt = N_0 * math.exp(-mu_Cu * X_Cu)
print("Out of", N_0, "photons, only", round(N_dt), "are transmitted")
```

- We know that $$N_{\mathrm{Cu}}=500 \times \exp(-\mu_{Cu} \times 1) = 8$$
- We want to illustrate than water does not attenuate X-rays of 100 keV as much as copper.
- Repeat some of the previous steps to find the what the total linear attenuation coefficient of water is for 100 keV.
- The Beer-Lambert law in this case is:

- What is the value of $x$?

In [ ]:

```
X_H2O = None # Thickness in cm, it is unknown
mu_H2O = ??? # Linear attenuation of water at 100 keV
```

$$500 \times \exp(-\mu_{H_2O} \times x_{H_2O}) = N_{\mathrm{H_2O}} = 8$$$$\exp(-\mu_{H_2O} \times x_{H_2O}) = \frac{N_{\mathrm{H_2O}}}{500}$$$$-\mu_{H_2O} \times x_{H_2O} = \ln\left(\frac{N_{\mathrm{H_2O}}}{500}\right)$$$$\mu_{H_2O} \times x_{H_2O} = -\ln\left(\frac{N_{\mathrm{H_2O}}}{500}\right)$$$$x_{H_2O} = \frac{-\ln\left(\frac{N_{\mathrm{H_2O}}}{500}\right)}{\mu_{H_2O}}$$

In [ ]:

```
X_H2O = (-math.log(N_dt / 500)) / mu_H2O
print(round(X_H2O), "cm of water is needed to stop as many photons of 100 keV as", X_Cu,"cm of copper")
```

In [ ]:

```
N_dt = N_0 * math.exp(-mu_H2O * X_H2O)
print("Out of", N_0, "photons, only", round(N_dt), "are transmitted")
```

**Given $E$, $L_p$ and $\mu$, we can compute X-ray images.**

- Computational bottlenecks are:
- computing $L_p$ from polygon meshes
- the exponential,
- the accumulation $(\sum)$

- gVirtualXray provides real-time mechanisms to compute them on GPU.

- Detect every intersection between a ray and the objects;
- Sort intersection (can be handled by GPUs using depth-peeling, a multi-pass rendering technique for semi-transparent polygonal objects without sorting polygons);
- Compute path length.

**Intersection sorting is actually not needed!**

- By convention normals are outward;
- A ray penetrates into an object when the dot product between the view vector (V) and the normal (N
_{i}) at the intersection point is positive; - It leaves an object when the dot product is negative.

L_{p}=Σ_{i} - sng(V · N_{i}) x d_{i}

- i refers to i
^{th}intersection in an arbitrary order; - d
_{i}distance from X-ray source to intersection point; - sgn(V · N
_{i}) stands for the sign of the dot product between V and N_{i}; - In a shader program [^1] , compute:
- sgn(V · N
_{i}); - d
_{i}the distance between the X-ray source and the intersection; - Assign -sng(V · N
_{i}) x d_{i}as the fragment value.

- sgn(V · N
- For each pixel, compute L
_{p}thanks to high-dynamic range and OpenGL blending function (pixel values may not be between 0 and 1)[^2].

[^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.

[^2]: See DOI: 10.2312/LocalChapterEvents/TPCG/TPCG09/025-032 for more details.

pixel = E x N_{out}

pixel = E x N_{in}(E) e^{(-Σi μi Lp(i))}

The Beer-Lambert law is broken down in several rendering passes:

- For each object of the scene:
- Compute L
_{p}(i); - Update results of Σ μ
_{i}L_{p}(i).

- Compute L
- For the final image only:
- Compute N
_{out}; - (Optional when direct display only is needed).

- Compute N

- Simulating an image relies on a Beer-Lambert law implementation;
- Solving the Beer-Lambert law relies on Linear Attenuation Coefficients; (μ)

μ is computed using Mass Attenuation Coefficients (μ/ρ) and material density (ρ).

**Are the values used in gVirtualXRay accurate?**- Compare values computed in gVirtualXRay with those from the literature.

**Are the Beer-Lambert law implementations accurate?**- Compare values computed in gVirtualXRay with theoretical ones.

**Are the simulated images accurate?**- Compare images computed using gVirtualXRay with those using a state-of-the-art Monte Carlo software, e.g.

- Simulating an image relies on a Beer-Lambert law implementation;
- Solving the Beer-Lambert law relies on Linear Attenuation Coefficients; (μ)
- μ is not known for given incident energies;
- μ is computed using Mass Attenuation Coefficients (μ/ρ) and material density (ρ).
**Are the values used in gVirtualXRay accurate?**- Compare values computed in gVirtualXRay with those from the literature.

- Simulating an image relies on a Beer-Lambert law implementation;
- Solving the Beer-Lambert law relies on Linear Attenuation Coefficients; (μ)
- μ is not known for given incident energies;
- μ is computed using Mass Attenuation Coefficients (μ/ρ) and material density (ρ).
**Are the values used in gVirtualXRay accurate?**- Compare values computed in gVirtualXRay with those from the literature.
**YES**

- Cube: edge length of 3 cm, made of soft tissue (HU = 52).
- Cylinder: height of 3 m, diameter of 2 cm, made of bone (HU = 1330).

We simulate a test case twice:

- Using a Monte Carlo method for particle physics implemented in GATE;
- Using our GPU implementation.

*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.*

Simulation parameters | Image computed with GATE (2 weeks of computations on supercomputer) | Image computed with gVirtualXRay (less than 1 sec. of computations on GPU) |
---|---|---|

Normalised cross-correlation (NCC) = 99.747%

Lungman phantom | 3D model of the phantom, source and detector |
---|---|

Normalised cross-correlation (NCC) = 98.92%

This research and/or training day would not have been possible without Jean Michel Létang, Jean Yves Buffière, Nicolas Freud, Nigel John, Tianci Wen, Iwan Mitchell, Jammie Pointon, Llion Evans, Ben Thorpes, David Fairbrother, ...