The Effects of Noise on Gamma

Here is an example of the effects noise can have on gamma. Simply by adding noise to the evaluation distribution the Gamma goes from a passing rate of about 50% to about 99%. 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.

To run this notebook yourself, for example on Colab, install PyMedPhys first. The next cell does that and is skipped when the documentation is built.

%pip install "pymedphys[user]"

import numpy as np
import matplotlib.pyplot as plt

import pymedphys

gamma_options = {
    'dose_percent_threshold': 3,
    'distance_mm_threshold': 3,
    'lower_percent_dose_cutoff': 20,
    'interp_fraction': 10,  # Should be 10 or more for more accurate results
    'max_gamma': 2,
    'random_subset': None,
    'local_gamma': True,
    'ram_available': 2**29,  # 1/2 GB
    'quiet': True
}

grid = 0.5
scale_factor = 1.035
noise = 0.01

xmin = -28
xmax = 28
ymin = -25
ymax = 25

extent = [xmin-grid/2, xmax+grid/2, ymin-grid/2, ymax+grid/2]

Defining an example dose reference

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.

x = np.arange(xmin, xmax + grid, grid)
y = np.arange(ymin, ymax + grid, grid)

coords = (y, x)

xx, yy = np.meshgrid(x, y)
dose_ref = np.exp(-((xx/15)**20 + (yy/15)**20))

plt.figure()
plt.title('Reference dose')

plt.imshow(
    dose_ref, clim=(0, 1.04), extent=extent)
plt.colorbar();

The first coordinate array in coords must match the first axis of dose_ref, and the second must match the second axis. PyMedPhys uses the same coordinate-array convention as SciPy’s RegularGridInterpolator, but the default gamma calculation now uses PyMedPhys’ own interpolation implementation. The optional interp_algo="scipy" selects SciPy instead.

len(coords[0])

101

len(coords[1])

113

np.shape(dose_ref)

(101, 113)

dimensions_of_dose_ref = np.shape(dose_ref)
assert dimensions_of_dose_ref[0] == len(coords[0])
assert dimensions_of_dose_ref[1] == len(coords[1])

Defining an example evaluation dose

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.

dose_eval = dose_ref * scale_factor

plt.figure()
plt.title('Evaluation dose')

plt.imshow(
    dose_eval, clim=(0, 1.04), extent=extent)
plt.colorbar();

Seeing a dose difference

dose_diff = dose_eval - dose_ref

plt.figure()
plt.title('Dose Difference')

plt.imshow(
    dose_diff, 
    clim=(-0.1, 0.1), extent=extent,
    cmap='seismic'
)
plt.colorbar();

Evaluation Gamma

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.

gamma_no_noise = pymedphys.gamma(
    coords, dose_ref,
    coords, dose_eval,
    **gamma_options)

plt.figure()
plt.title('Gamma Distribution')

plt.imshow(
    gamma_no_noise, clim=(0, 2), extent=extent,
    cmap='coolwarm')
plt.colorbar()

plt.show()
valid_gamma_no_noise = gamma_no_noise[~np.isnan(gamma_no_noise)]
no_noise_passing = np.round(pymedphys.gamma_pass_rate(gamma_no_noise), 2)

plt.figure()
plt.title(f'Gamma Histogram | Passing rate = {round(no_noise_passing, 2)}%')
plt.xlabel('Gamma')
plt.ylabel('Number of pixels')

plt.hist(valid_gamma_no_noise, 20);

/tmp/ipykernel_1401/1841687323.py:1: DeprecationWarning: Parameter `quiet` will be deprecated in the future
  gamma_no_noise = pymedphys.gamma(

Investigating the effect of noise

Let’s now slightly adjust our evaluation distribution by adding some random normal noise.

np.random.seed(0)
dose_eval_noise = (
    dose_eval + 
    np.random.normal(loc=0, scale=noise, size=np.shape(dose_eval))
)

plt.figure()
plt.title('Evaluation dose with Noise')

plt.imshow(
    dose_eval_noise, clim=(0, 1.04), extent=extent)
plt.colorbar();

Let’s see what our dose difference looks like with this noise.

dose_diff_with_noise = dose_eval_noise - dose_ref

plt.figure()
plt.title('Dose Difference with Noise')

plt.imshow(
    dose_diff_with_noise, 
    clim=(-0.1, 0.1), extent=extent,
    cmap='seismic'
)
plt.colorbar();

But what about our gamma?

gamma_with_noise = pymedphys.gamma(
    coords, dose_ref,
    coords, dose_eval_noise,
    **gamma_options)

plt.figure()
plt.title('Gamma Distribution')

plt.imshow(
    gamma_with_noise, clim=(0, 2), extent=extent,
    cmap='coolwarm')
plt.colorbar()

plt.show()
valid_gamma_with_noise = gamma_with_noise[~np.isnan(gamma_with_noise)]
with_noise_passing = np.round(pymedphys.gamma_pass_rate(gamma_with_noise), 2)

plt.figure()
plt.title(f'Gamma Histogram | Passing rate = {round(with_noise_passing, 2)}%')
plt.xlabel('Gamma')
plt.ylabel('Number of pixels')

plt.hist(valid_gamma_with_noise, 20);

/tmp/ipykernel_1401/2350559230.py:1: DeprecationWarning: Parameter `quiet` will be deprecated in the future
  gamma_with_noise = pymedphys.gamma(

Conclusion

This example illustrates how noise can increase a gamma pass rate. It does not establish that a high pass rate alone validates a dose calculation.

For examples of systematic errors missed by commonly used gamma criteria, see Nelms et al. (2013), Evaluating IMRT and VMAT dose accuracy. The authors compare gamma results with additional diagnostic methods and show why those methods can reveal errors that the pass rate misses.