An Example showing 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.

[1]: import numpy as np import matplotlib.pyplot as plt import pymedphys

[2]: 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 }

[3]: grid = 0.5 scale_factor = 1.035 noise = 0.01

[4]: 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. [5]: 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 (); 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. [6]: len ( coords [ 0 ]) [6]: 101 [7]: len ( coords [ 1 ]) [7]: 113 [8]: np . shape ( dose_ref ) [8]: (101, 113) [9]: 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. [10]: 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¶ [11]: 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. [12]: 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 = 100 * np . round ( np . sum ( valid_gamma_no_noise <= 1 ) / len ( valid_gamma_no_noise ), 4 ) 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 );

Investigating the effect of noise¶ Let’s now slightly adjust our evaluation distribution by adding some random normal noise. [13]: 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 dise difference looks like with this noise [14]: 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? [15]: 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 = 100 * np . round ( np . sum ( valid_gamma_with_noise <= 1 ) / len ( valid_gamma_with_noise ), 4 ) 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 );