Session 1

Introduction to X-ray attenuation and its implementation in

Authors: Franck Vidal and Jean-Michel Létang

(version 1.0, 22 Sep 2022)

Aims of this session

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

Acknowledgements

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

What is gVXR and why has it been developed?

Medical Training using VR

Timeline

• 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;

Trend

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

Why gVXR?

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

is not a replacement for Monte Carlo methods, it's a complement when a trade-off between speed and accuracy must met.

In a nutshell, gVXR is (1/2)

• 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

In a nutshell, gVXR is (2/2)

• 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)

Simulation supports:

• 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).

More recently

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

Fact sheet

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

Examples of applications that use gVirtualXRay

• Teaching particle physics to undergraduates;
• Virtual Reality Simulation;
• Virtual Testing Lab
• Design new clinical imaging modality to reduce dose exposure
• Teaching radiography in medicine

Synchrotron μ-tomography

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

Use of fast realistic simulations on GPU to extract CAD models from microtomographic data in the presence of strong CT artefacts

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

Registration of the TiAlV matrix

Evolution of the core and fibre radii during the optimisation

Evolution of beam spectrum parameters during the optimisation.

Evolution of the core and fibre radii and phase contrast parameters during the optimisation

Principles of rojection X-ray imaging

X-ray production: X-ray tube

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.

X-ray production: Synchrotron

ESRF, Grenoble, France	Diamond Light Source, Didcot, UK
Source: Photograph by German Wikipedian Christian Hendrich, October 2004	Copyright of Diamond Light Source Ltd

How do X-rays interact with matter?

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

1. Directly transmitted photons (no interaction);
2. Absorbed photons;
3. Scattered photons;
4. 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).

Beer-Lambert law to compute the attenuation of X-rays by matter

• 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_{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.

Linear attenuation coefficient

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

%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'}));

Attenuation of 500 photons of 100 keV by 1cm of copper

In the previous cell,

1. Select copper as material
2. Click on Run Interact
3. Using the mouse cursor, find what the total linear attenuation coefficient of copper is for 100 keV,
4. Replace ??? in the cell below with the corresponding value
5. 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

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

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")

Attenuation of 500 photons of 100 keV by ??cm of water

• 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:

$$N_{\mathrm{H_2O}}= 500 \times \exp(-\mu_{H_2O} \times x) = 8$$

• What is the value of $x$?

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}}$$

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")

Check if we made a mistake

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

What does it have to do with gVirtualXray?

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.

Path Length: Naive Approach

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

Path Length: L-Buffer

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-Buffer Implementation

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

Multipass Rendering Pipeline

pixel = E x N_out

pixel = E x N_in(E) e^{(-Σ_i μ_i L_p(i))}

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

• For each object of the scene:
  1. Compute L_p(i);
  2. Update results of Σ μ_i L_p(i).
• For the final image only:
  1. Compute N_out;
  2. (Optional when direct display only is needed).

Quantitative Validation

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

  

Mass attenuation coefficients from elements

Mass attenuation coefficients from mixtures

Going back to previous slide

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

Going back to previous slide

• 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

More testing

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

  

Test case

• 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).

Gate vs. gVirtualXRay

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.

Gate vs. gVirtualXRay

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%

Real digital radiograph (DR) taken with a clinically utilised X-ray equipment vs. gVirtualXRay

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

Normalised cross-correlation (NCC) = 98.92%

Acknowledgements

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