This page was generated from docs/user/examples/gamma_from_dicom.ipynb. Interactive online version: Binder badge

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.

[1]:
import os
import warnings
from glob import glob
import zipfile

import numpy as np
import matplotlib.pyplot as plt

import pydicom

from pymedphys.gamma import gamma_dicom
from pymedphys.dicom import zyx_and_dose_from_dataset

Getting the demo DICOM files

Over at https://app.pymedphys.com there are some demo files which can be used. Let’s grab those here for the purpose of demonstrating gamma_dicom usage.

[2]:
pymedphys_app_demo_files_zip = '../../../app/src/data/demo-files.zip'

with zipfile.ZipFile(pymedphys_app_demo_files_zip) as myzip:
    with myzip.open('original_dose_beam_4.dcm', 'r') as myfile:
        reference = pydicom.dcmread(myfile, force=True)

    with myzip.open('logfile_dose_beam_4.dcm', 'r') as myfile:
        evaluation = pydicom.dcmread(myfile, force=True)

Calculate and display Gamma

[3]:
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 = gamma_dicom(
    reference, evaluation, **gamma_options)
Calcing using local normalisation point for gamma
Global normalisation set to 1.9462822416083623 Gy
Global dose threshold set to [0.01946282] Gy ([1]%)
Distance threshold set to [1] mm
Lower dose cutoff set to 0.3892564483216725 Gy (20%)

Current distance: 0.60 mm | Number of reference points remaining: 257 | Points tested per reference point: 485 | RAM split count: 1
Complete!
[4]:
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)
[4]:
Text(0.5, 1.0, 'Local Gamma (0.5%/0.5mm) | Percent Pass: 100.00 %')
../../_images/user_examples_gamma_from_dicom_6_1.png

Plotting the Dose and the Gamma

[5]:
(z_ref, y_ref, x_ref), dose_reference = zyx_and_dose_from_dataset(reference)
(z_eval, y_eval, x_eval), dose_evaluation = zyx_and_dose_from_dataset(evaluation)

dose_reference = dose_reference * 100
dose_evaluation= dose_evaluation * 100
[6]:
lower_dose_cutoff = gamma_options['lower_percent_dose_cutoff'] / 100 * np.max(dose_reference)

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)])
[7]:
max_ref_dose = np.max(dose_reference)

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, cmap=plt.get_cmap('viridis'))
    ax[0,0].set_title("Evaluation")
    fig.colorbar(c00, ax=ax[0,0], label='Dose (cGy)')
    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, cmap=plt.get_cmap('viridis'))
    ax[0,1].set_title("Reference")
    fig.colorbar(c01, ax=ax[0,1], label='Dose (cGy)')
    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] (cGy)')
    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("Local Gamma (0.5%/0.5mm)")
    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.savefig('{}.png'.format(i), dpi=300)
    plt.show()
    print("\n")
../../_images/user_examples_gamma_from_dicom_10_0.png


../../_images/user_examples_gamma_from_dicom_10_2.png


../../_images/user_examples_gamma_from_dicom_10_4.png


../../_images/user_examples_gamma_from_dicom_10_6.png


../../_images/user_examples_gamma_from_dicom_10_8.png


../../_images/user_examples_gamma_from_dicom_10_10.png


../../_images/user_examples_gamma_from_dicom_10_12.png