Gamma from DICOM¶
PyMedPhys has multiple ways to calculate Gamma. There are also a range of interfaces that can be used. Presented here is a simplified interface which receives as its input two DICOM file paths for the purpose of directly calculating Gamma from a pair of RT DICOM dose files.
import numpy as np
import matplotlib.pyplot as plt
!pip install pydicom pymedphys
import pydicom
import pymedphys
Getting the demo DICOM files¶
Let’s download some demo files for the purpose of demonstrating gamma_dicom
usage.
reference_filepath = pymedphys.data_path("original_dose_beam_4.dcm")
reference_filepath
evaluation_filepath = pymedphys.data_path("logfile_dose_beam_4.dcm")
evaluation_filepath
reference = pydicom.read_file(str(reference_filepath), force=True)
evaluation = pydicom.read_file(str(evaluation_filepath), force=True)
axes_reference, dose_reference = pymedphys.dicom.zyx_and_dose_from_dataset(reference)
axes_evaluation, dose_evaluation = pymedphys.dicom.zyx_and_dose_from_dataset(evaluation)
(z_ref, y_ref, x_ref) = axes_reference
(z_eval, y_eval, x_eval) = axes_evaluation
Importantly the shape of the coordinates are in the same order as the dose axis order
np.shape(z_ref)
np.shape(y_ref)
np.shape(x_ref)
np.shape(dose_reference)
Calculate and display Gamma¶
gamma_options = {
'dose_percent_threshold': 1,
'distance_mm_threshold': 1,
'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
}
gamma = pymedphys.gamma(
axes_reference, dose_reference,
axes_evaluation, dose_evaluation,
**gamma_options)
valid_gamma = gamma[~np.isnan(gamma)]
num_bins = (
gamma_options['interp_fraction'] * gamma_options['max_gamma'])
bins = np.linspace(0, gamma_options['max_gamma'], num_bins + 1)
plt.hist(valid_gamma, bins, density=True)
plt.xlim([0, gamma_options['max_gamma']])
pass_ratio = np.sum(valid_gamma <= 1) / len(valid_gamma)
plt.title(f"Local Gamma ({gamma_options['dose_percent_threshold']}%/{gamma_options['distance_mm_threshold']}mm) | Percent Pass: {pass_ratio*100:.2f} %")
# plt.savefig('gamma_hist.png', dpi=300)
Plotting the Dose and the Gamma¶
max_ref_dose = np.max(dose_reference)
lower_dose_cutoff = gamma_options['lower_percent_dose_cutoff'] / 100 * max_ref_dose
relevant_slice = (
np.max(dose_reference, axis=(1, 2)) >
lower_dose_cutoff)
slice_start = np.max([
np.where(relevant_slice)[0][0],
0])
slice_end = np.min([
np.where(relevant_slice)[0][-1],
len(z_ref)])
z_vals = z_ref[slice(slice_start, slice_end, 5)]
eval_slices = [
dose_evaluation[np.where(z_i == z_eval)[0][0], :, :]
for z_i in z_vals
]
ref_slices = [
dose_reference[np.where(z_i == z_ref)[0][0], :, :]
for z_i in z_vals
]
gamma_slices = [
gamma[np.where(z_i == z_ref)[0][0], :, :]
for z_i in z_vals
]
diffs = [
eval_slice - ref_slice
for eval_slice, ref_slice
in zip(eval_slices, ref_slices)
]
max_diff = np.max(np.abs(diffs))
for i, (eval_slice, ref_slice, diff, gamma_slice) in enumerate(zip(eval_slices, ref_slices, diffs, gamma_slices)):
fig, ax = plt.subplots(figsize=(13,10), nrows=2, ncols=2)
c00 = ax[0,0].contourf(
x_eval, y_eval, eval_slice, 100,
vmin=0, vmax=max_ref_dose)
ax[0,0].set_title("Evaluation")
fig.colorbar(c00, ax=ax[0,0], label='Dose (Gy)')
ax[0,0].invert_yaxis()
ax[0,0].set_xlabel('x (mm)')
ax[0,0].set_ylabel('z (mm)')
c01 = ax[0,1].contourf(
x_ref, y_ref, ref_slice, 100,
vmin=0, vmax=max_ref_dose)
ax[0,1].set_title("Reference")
fig.colorbar(c01, ax=ax[0,1], label='Dose (Gy)')
ax[0,1].invert_yaxis()
ax[0,1].set_xlabel('x (mm)')
ax[0,1].set_ylabel('z (mm)')
c10 = ax[1,0].contourf(
x_ref, y_ref, diff, 100,
vmin=-max_diff, vmax=max_diff, cmap=plt.get_cmap('seismic'))
ax[1,0].set_title("Dose difference")
fig.colorbar(c10, ax=ax[1,0], label='[Dose Eval] - [Dose Ref] (Gy)')
ax[1,0].invert_yaxis()
ax[1,0].set_xlabel('x (mm)')
ax[1,0].set_ylabel('z (mm)')
c11 = ax[1,1].contourf(
x_ref, y_ref, gamma_slice, 100,
vmin=0, vmax=2, cmap=plt.get_cmap('coolwarm'))
ax[1,1].set_title(
f"Local Gamma ({gamma_options['dose_percent_threshold']} % / {gamma_options['distance_mm_threshold']} mm)")
fig.colorbar(c11, ax=ax[1,1], label='Gamma Value')
ax[1,1].invert_yaxis()
ax[1,1].set_xlabel('x (mm)')
ax[1,1].set_ylabel('z (mm)')
plt.show()
print("\n")