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
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.8/site-packages/IPython/utils/_process_posix.py:148, in ProcessHandler.system(self, cmd)
147 else:
--> 148 child = pexpect.spawn(self.sh, args=['-c', cmd]) # Vanilla Pexpect
149 flush = sys.stdout.flush
File ~/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.8/site-packages/pexpect/pty_spawn.py:205, in spawn.__init__(self, command, args, timeout, maxread, searchwindowsize, logfile, cwd, env, ignore_sighup, echo, preexec_fn, encoding, codec_errors, dimensions, use_poll)
204 else:
--> 205 self._spawn(command, args, preexec_fn, dimensions)
206 self.use_poll = use_poll
File ~/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.8/site-packages/pexpect/pty_spawn.py:303, in spawn._spawn(self, command, args, preexec_fn, dimensions)
300 self.args = [a if isinstance(a, bytes) else a.encode(self.encoding)
301 for a in self.args]
--> 303 self.ptyproc = self._spawnpty(self.args, env=self.env,
304 cwd=self.cwd, **kwargs)
306 self.pid = self.ptyproc.pid
File ~/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.8/site-packages/pexpect/pty_spawn.py:315, in spawn._spawnpty(self, args, **kwargs)
314 '''Spawn a pty and return an instance of PtyProcess.'''
--> 315 return ptyprocess.PtyProcess.spawn(args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.8/site-packages/ptyprocess/ptyprocess.py:315, in PtyProcess.spawn(cls, argv, cwd, env, echo, preexec_fn, dimensions, pass_fds)
314 os.close(exec_err_pipe_write)
--> 315 exec_err_data = os.read(exec_err_pipe_read, 4096)
316 os.close(exec_err_pipe_read)
KeyboardInterrupt:
During handling of the above exception, another exception occurred:
UnboundLocalError Traceback (most recent call last)
Cell In[2], line 1
----> 1 get_ipython().system('pip install pydicom pymedphys')
3 import pydicom
4 import pymedphys
File ~/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.8/site-packages/ipykernel/zmqshell.py:657, in ZMQInteractiveShell.system_piped(self, cmd)
655 self.user_ns["_exit_code"] = system(cmd)
656 else:
--> 657 self.user_ns["_exit_code"] = system(self.var_expand(cmd, depth=1))
File ~/checkouts/readthedocs.org/user_builds/pymedphys/envs/latest/lib/python3.8/site-packages/IPython/utils/_process_posix.py:164, in ProcessHandler.system(self, cmd)
159 out_size = len(child.before)
160 except KeyboardInterrupt:
161 # We need to send ^C to the process. The ascii code for '^C' is 3
162 # (the character is known as ETX for 'End of Text', see
163 # curses.ascii.ETX).
--> 164 child.sendline(chr(3))
165 # Read and print any more output the program might produce on its
166 # way out.
167 try:
UnboundLocalError: local variable 'child' referenced before assignment
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)
#if density is True, y value is probability density; otherwise, it is count in a bin
plt.xlim([0, gamma_options['max_gamma']])
plt.xlabel('gamma index')
plt.ylabel('probability density')
pass_ratio = np.sum(valid_gamma <= 1) / len(valid_gamma)
if gamma_options['local_gamma']:
gamma_norm_condition = 'Local gamma'
else:
gamma_norm_condition = 'Global gamma'
plt.title(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(z_ref)*len(y_ref)*len(x_ref)} | valid \u03B3 pts: {len(valid_gamma)}")
# 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)
fig.suptitle('Slice Z_{0}'.format(slice_start+i*5), fontsize=12)
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('y (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('y (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")
cbar = fig.colorbar(c10, ax=ax[1,0], label='[Dose Eval] - [Dose Ref] (Gy)')
cbar.formatter.set_powerlimits((0, 0)) #use scientific notation
ax[1,0].invert_yaxis()
ax[1,0].set_xlabel('x (mm)')
ax[1,0].set_ylabel('y (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"{gamma_norm_condition} ({gamma_options['dose_percent_threshold']} % / {gamma_options['distance_mm_threshold']} mm)")
fig.colorbar(c11, ax=ax[1,1], label='gamma index')
ax[1,1].invert_yaxis()
ax[1,1].set_xlabel('x (mm)')
ax[1,1].set_ylabel('y (mm)')
plt.show()
print("\n")