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
Requirement already satisfied: pydicom in /opt/buildhome/python3.7/lib/python3.7/site-packages (2.1.2)
Requirement already satisfied: pymedphys in /opt/buildhome/python3.7/lib/python3.7/site-packages (0.36.0.dev1)
Requirement already satisfied: typing-extensions in /opt/buildhome/python3.7/lib/python3.7/site-packages (from pymedphys) (3.7.4.3)
WARNING: You are using pip version 20.3.1; however, version 20.3.3 is available.
You should consider upgrading via the '/opt/buildhome/python3.7/bin/python3.7 -m pip install --upgrade pip' command.
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
original_dose_beam_4.dcm?download=1: 0.00B [00:00, ?B/s]
original_dose_beam_4.dcm?download=1: 0%| | 8.19k/24.4M [00:01<58:04, 7.00kB/s]
original_dose_beam_4.dcm?download=1: 0%| | 16.4k/24.4M [00:01<27:17, 14.9kB/s]
original_dose_beam_4.dcm?download=1: 0%| | 41.0k/24.4M [00:01<09:07, 44.5kB/s]
original_dose_beam_4.dcm?download=1: 0%| | 90.1k/24.4M [00:01<03:37, 112kB/s]
original_dose_beam_4.dcm?download=1: 1%| | 197k/24.4M [00:01<01:28, 274kB/s]
original_dose_beam_4.dcm?download=1: 2%|▏ | 393k/24.4M [00:01<00:40, 592kB/s]
original_dose_beam_4.dcm?download=1: 3%|▎ | 811k/24.4M [00:01<00:17, 1.32MB/s]
original_dose_beam_4.dcm?download=1: 7%|▋ | 1.63M/24.4M [00:01<00:08, 2.79MB/s]
original_dose_beam_4.dcm?download=1: 13%|█▎ | 3.24M/24.4M [00:02<00:03, 5.77MB/s]
original_dose_beam_4.dcm?download=1: 19%|█▉ | 4.75M/24.4M [00:02<00:02, 7.66MB/s]
original_dose_beam_4.dcm?download=1: 31%|███▏ | 7.67M/24.4M [00:02<00:01, 12.5MB/s]
original_dose_beam_4.dcm?download=1: 36%|███▋ | 8.90M/24.4M [00:02<00:01, 11.7MB/s]
original_dose_beam_4.dcm?download=1: 48%|████▊ | 11.7M/24.4M [00:02<00:00, 15.1MB/s]
original_dose_beam_4.dcm?download=1: 56%|█████▌ | 13.7M/24.4M [00:02<00:00, 15.8MB/s]
original_dose_beam_4.dcm?download=1: 60%|█████▉ | 14.6M/24.4M [00:02<00:00, 13.3MB/s]
original_dose_beam_4.dcm?download=1: 71%|███████ | 17.3M/24.4M [00:02<00:00, 16.2MB/s]
original_dose_beam_4.dcm?download=1: 79%|███████▉ | 19.3M/24.4M [00:03<00:00, 12.7MB/s]
original_dose_beam_4.dcm?download=1: 92%|█████████▏| 22.4M/24.4M [00:03<00:00, 16.2MB/s]
original_dose_beam_4.dcm?download=1: 95%|█████████▌| 23.3M/24.4M [00:03<00:00, 13.7MB/s]
original_dose_beam_4.dcm?download=1: 24.4MB [00:03, 7.14MB/s]
PosixPath('/opt/buildhome/.pymedphys/data/original_dose_beam_4.dcm')
evaluation_filepath = pymedphys.data_path("logfile_dose_beam_4.dcm")
evaluation_filepath
logfile_dose_beam_4.dcm?download=1: 0.00B [00:00, ?B/s]
logfile_dose_beam_4.dcm?download=1: 0%| | 8.19k/24.4M [00:01<57:14, 7.10kB/s]
logfile_dose_beam_4.dcm?download=1: 0%| | 16.4k/24.4M [00:01<26:57, 15.1kB/s]
logfile_dose_beam_4.dcm?download=1: 0%| | 41.0k/24.4M [00:01<09:01, 45.0kB/s]
logfile_dose_beam_4.dcm?download=1: 0%| | 81.9k/24.4M [00:01<04:03, 99.8kB/s]
logfile_dose_beam_4.dcm?download=1: 1%| | 188k/24.4M [00:01<01:31, 265kB/s]
logfile_dose_beam_4.dcm?download=1: 2%|▏ | 393k/24.4M [00:01<00:39, 601kB/s]
logfile_dose_beam_4.dcm?download=1: 3%|▎ | 803k/24.4M [00:01<00:18, 1.31MB/s]
logfile_dose_beam_4.dcm?download=1: 7%|▋ | 1.63M/24.4M [00:01<00:08, 2.80MB/s]
logfile_dose_beam_4.dcm?download=1: 13%|█▎ | 3.26M/24.4M [00:02<00:03, 5.80MB/s]
logfile_dose_beam_4.dcm?download=1: 26%|██▌ | 6.38M/24.4M [00:02<00:01, 11.6MB/s]
logfile_dose_beam_4.dcm?download=1: 39%|███▊ | 9.41M/24.4M [00:02<00:00, 15.7MB/s]
logfile_dose_beam_4.dcm?download=1: 50%|█████ | 12.3M/24.4M [00:02<00:00, 18.0MB/s]
logfile_dose_beam_4.dcm?download=1: 63%|██████▎ | 15.3M/24.4M [00:02<00:00, 20.0MB/s]
logfile_dose_beam_4.dcm?download=1: 75%|███████▍ | 18.2M/24.4M [00:02<00:00, 21.5MB/s]
logfile_dose_beam_4.dcm?download=1: 87%|████████▋ | 21.1M/24.4M [00:02<00:00, 22.3MB/s]
logfile_dose_beam_4.dcm?download=1: 99%|█████████▉| 24.2M/24.4M [00:02<00:00, 23.5MB/s]
logfile_dose_beam_4.dcm?download=1: 24.4MB [00:02, 8.35MB/s]
PosixPath('/opt/buildhome/.pymedphys/data/logfile_dose_beam_4.dcm')
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)
(152,)
np.shape(y_ref)
(158,)
np.shape(x_ref)
(254,)
np.shape(dose_reference)
(152, 158, 254)
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)
Calcing using local normalisation point for gamma
Global normalisation set to 1.9462822416083623
Global dose threshold set to [0.01946282] ([1]% of normalisation)
Distance threshold set to [1]
Lower dose cutoff set to 0.3892564483216725 (20% of normalisation)
Current distance: 0.00 mm | Number of reference points remaining: 106910 | Points tested per reference point: 1 | RAM split count: 1
Current distance: 0.10 mm | Number of reference points remaining: 106910 | Points tested per reference point: 23 | RAM split count: 1
Current distance: 0.20 mm | Number of reference points remaining: 74257 | Points tested per reference point: 67 | RAM split count: 2
Current distance: 0.30 mm | Number of reference points remaining: 16220 | Points tested per reference point: 135 | RAM split count: 1
Current distance: 0.40 mm | Number of reference points remaining: 2732 | Points tested per reference point: 227 | RAM split count: 1
Current distance: 0.50 mm | Number of reference points remaining: 760 | Points tested per reference point: 348 | RAM split count: 1
Current distance: 0.60 mm | Number of reference points remaining: 257 | Points tested per reference point: 485 | RAM split count: 1
Complete!
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("Local Gamma (0.5%/0.5mm) | Percent Pass: {0:.2f} %".format(pass_ratio*100))
# plt.savefig('gamma_hist.png', dpi=300)
Text(0.5, 1.0, 'Local Gamma (0.5%/0.5mm) | Percent Pass: 100.00 %')
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")
