Gamma from DICOM

This example reads two DICOM RT Dose files with pydicom, extracts their coordinate arrays and dose grids with pymedphys.dicom.zyx_and_dose_from_dataset, and passes those arrays to pymedphys.gamma. The gamma function does not take file paths directly.

To run this notebook yourself, for example on Colab, install PyMedPhys first. The next cell does that and is skipped when the documentation is built.

For an illustrated explanation of patient directions, array reordering and gamma overlays, see DICOM coordinates and gamma, illustrated. The explanation is written for physicists and can be read without following the Python code.

%pip install "pymedphys[user]"

import numpy as np
import matplotlib.pyplot as plt

import pydicom
import pymedphys

Getting the demo DICOM files

Let’s download two demo RT Dose files for this example.

reference_filepath = pymedphys.data_path("original_dose_beam_4.dcm")
reference_filepath

PosixPath('/home/docs/.pymedphys/data/original_dose_beam_4.dcm')

evaluation_filepath = pymedphys.data_path("logfile_dose_beam_4.dcm")
evaluation_filepath

PosixPath('/home/docs/.pymedphys/data/logfile_dose_beam_4.dcm')

reference = pydicom.dcmread(str(reference_filepath), force=True)
evaluation = pydicom.dcmread(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

The axes and dose are returned as a matching pair in patient (z, y, x) order. For example, dose_reference[k, j, i] is at (z_ref[k], y_ref[j], x_ref[i]).

The DICOM conversion makes each axis ascending and reverses or transposes the dose array to match. This preserves every dose value at its physical position; it does not interpolate or resample the dose. For non-HFS orientations or reversed slice order, the returned array can differ in order or shape from reference.pixel_array. Keep the returned axes and dose together.

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)

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 = pymedphys.gamma_pass_rate(gamma) / 100

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)

Text(0.5, 1.0, 'Dose cut: 20% | Local gamma (1%/1mm) | Pass Rate(γ<=1): 100.00% \n ref pts: 6100064 | valid γ pts: 106910')

Plotting the Dose and the Gamma

Gamma is returned on the reference grid: gamma[k, j, i] belongs to (z_ref[k], y_ref[j], x_ref[i]), just like dose_reference[k, j, i]. Plot gamma[k] against x_ref and y_ref, at the plane z_ref[k]. The ascending axes do not prevent physical alignment with the original dose; the corresponding dose values have been reordered with them.

Gamma only normalises the evaluation grid internally. If reference axes are supplied directly in descending order, gamma keeps that reference order in its output. The ascending reference axes here come from the DICOM conversion.

Use the converted reference dose for an overlay with gamma. Do not overlay gamma directly on raw reference.pixel_array by index unless their order is known to agree. To restore raw storage order, apply the inverse reversals and dimension permutation. Changing the plot’s viewing direction, such as with invert_yaxis(), is a display choice and does not change voxel positions.

The cells below display slices of an already calculated 3D gamma result. They do not perform a separate planar gamma calculation.

These demo datasets have coincident axes, allowing direct dose subtraction and selection of evaluation slices by matching z coordinates. For other inputs, use the evaluation dose’s own axes, choose slices by physical position rather than a shared index, and interpolate onto a common grid before subtracting dose arrays when sample positions differ. Gamma itself does not require coincident grids. The gamma panels always use the reference axes.

For the derivation and affected historical DICOM gamma workflows, see the illustrated guide; for the validation evidence and known limitations, see the DICOM coordinate validation record.

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")