Gamma from 1D ASCII data

PyMedPhys can also be used to calculate gamma index for one dimensional dose data. Presented here is an example which reads in two .csv files (reference and evaluation dose distributions) as input for the purpose of directly calculating the gamma index of the reference dose distribution. In this example, the reference dose distribution is from film measurement and the evaluation dose distribution is from Monte Carlo simulation.

import numpy as np
import matplotlib.pyplot as plt
import csv
import pymedphys
import urllib.request

Getting the demo .csv files

Let’s download the reference and evaluation dose data in .csv files for the purpose of demonstrating gamma analysis for 1D dose data. Users, please note that, there are not explict links of these two .csv files.

url = 'https://github.com/pymedphys/data/raw/200a79623e5c8b226e039e45a606512129a9df97/dose_film_1D_z0mm.csv'
save_filepath = 'dose_film_1D_z0mm.csv'

urllib.request.urlretrieve(url, save_filepath)

('dose_film_1D_z0mm.csv', <http.client.HTTPMessage at 0x7ff2a56aaf40>)

url = 'https://github.com/pymedphys/data/raw/200a79623e5c8b226e039e45a606512129a9df97/dose_MC_1D_z0mm.csv'
save_filepath = 'dose_MC_1D_z0mm.csv'

urllib.request.urlretrieve(url, save_filepath)

('dose_MC_1D_z0mm.csv', <http.client.HTTPMessage at 0x7ff2a56aa6a0>)

Calculate gamma

reference = np.genfromtxt('dose_film_1D_z0mm.csv',delimiter=',',skip_header=1)
evaluation = np.genfromtxt('dose_MC_1D_z0mm.csv',delimiter=',',skip_header=1)

axis_reference = reference[:, 0] #1st column is x in mm
dose_reference = reference[:, 1] #2nd column is dose in Gy/MU

axis_evaluation = evaluation[:, 0]
dose_evaluation = evaluation[:, 1]
gamma_options = {
    'dose_percent_threshold': 1,
    'distance_mm_threshold': 1,
    'lower_percent_dose_cutoff': 10,
    'interp_fraction': 10,  # Should be 10 or more for more accurate results
    'max_gamma': 2,
    'random_subset': None,
    'local_gamma': False, # False indicates global gamma is calculated
    'ram_available': 2 ** 29  # 1/2 GB
}
#for global dose normalization, the maximum reference dose is used
#but in TG218, it said usually the prescribed dose or the maximum dose in a plan (evaluation) is used
gamma = pymedphys.gamma(
    axis_reference, dose_reference,
    axis_evaluation, dose_evaluation,
    **gamma_options)

Plot the histogram of gamma index

valid_gamma = gamma[~np.isnan(gamma)]
print('# of reference points with a valid gamma {0}'.format( len(valid_gamma)) )
num_bins = (
    gamma_options['interp_fraction'] * gamma_options['max_gamma'])
bins = np.linspace(0, gamma_options['max_gamma'], num_bins + 1)

if gamma_options['local_gamma']:
    gamma_norm_condition = 'Local gamma'
else:
    gamma_norm_condition = 'Global gamma'

pass_ratio = np.sum(valid_gamma <= 1) / len(valid_gamma)
    
figure_1 = plt.figure(1, figsize=(8, 6), dpi=300, facecolor='w', edgecolor='k')
figure_1.suptitle(f"Dose cut: {gamma_options['lower_percent_dose_cutoff']}% | {gamma_norm_condition} ({gamma_options['dose_percent_threshold']}%/{gamma_options['distance_mm_threshold']}mm) | Pass Rate(\u03B3<=1): {pass_ratio*100:.2f}% \n ref pts: {len(dose_reference)} | valid \u03B3 pts: {len(valid_gamma)}")

plt.hist(valid_gamma, bins, density=True) #y value is probability density in each bin
#plt.hist(valid_gamma, bins, density=False) #y value is counts in each bin
plt.xlim([0, gamma_options['max_gamma']])
plt.xlabel('gamma index of reference point') #FG
plt.ylabel('probability density')
#plt.ylabel('counts')
plt.vlines(x=[1], ymin=0, ymax=1, colors='purple', ls='--', lw=2, label='target')
plt.savefig('1D_{0}_histogram.png'.format(gamma_norm_condition), dpi=300) #plt.savefig() must be before plt.show()

# of reference points with a valid gamma 401

Display gamma together with dose points

#show gamma index together with dose points
max_ref_dose = np.max(dose_reference) #max reference dose
max_eva_dose = np.max(dose_evaluation) #max evaluation dose
lower_dose_cutoff = gamma_options['lower_percent_dose_cutoff']/100*max_ref_dose

figure_2 = plt.figure(2, figsize=(8, 6), dpi=120, facecolor='w', edgecolor='k')
figure_2.suptitle('Reference and evaluation dose curves, and {0} index of reference dose points'.format(gamma_norm_condition),fontsize=12)

ax_1 = figure_2.add_subplot(111)
ax_1.tick_params(direction='in')
ax_1.tick_params(axis='x', bottom='on', top='on')  
ax_1.tick_params(labeltop='on')
ax_1.minorticks_on()
ax_1.set_xlabel('x(mm)')
ax_1.set_ylabel('dose(Gy/MU)')
ax_1.set_ylim([0, max(max_ref_dose,max_eva_dose) * 1.1])


ax_2 = ax_1.twinx()
ax_2.minorticks_on()
ax_2.set_ylabel('gamma index')
ax_2.set_ylim([0, gamma_options['max_gamma'] * 2.0 ])

curve_0 = ax_1.plot(axis_reference, dose_reference,'k-',label='reference dose')
curve_1 = ax_1.plot(axis_evaluation, dose_evaluation,'bo', mfc='none', markersize=5, label='evaluation dose')
curve_2 = ax_2.plot(axis_reference, gamma, 'r*', label=f"gamma ({gamma_options['dose_percent_threshold']}%/{gamma_options['distance_mm_threshold']}mm)")
curves = curve_0 + curve_1 + curve_2

labels = [l.get_label() for l in curves] 
ax_1.legend(curves, labels, loc='best') 
ax_1.grid(True)

# save figure first and show it
figureName = '1D dose_reference_evaluation_{0} index.png'.format(gamma_norm_condition)
plt.savefig(figureName) #plt.savefig() must be before plt.show()
plt.show()