Using Gamma for np.arrays

Hi Simon,

I am trying to use the pymedphys.gamma function on a numpy array. I have done this by following the steps used to get gamma from DICOM images. However, I have run into an issue when I do this. I am trying to generate gamma analysis maps from 3D numpy arrays, however, the array that is being pulled from the reference and evaluation is outputted in the wrong axis and in the wrong orientation. I was wondering how I would go about re-orienting the array so that the data can be outputted correctly. I have tried to fix the array myself but I seem to be getting nowhere. I have attached an image of how the data is currently being outputted, it looks like a dose profile on its side at the moment.

Please note, I also noticed that the axes_reference and axes_evaluation are ordered so that the axes are pulled in the order z,y,x whereas I ordered my arrays x,y,z but when I tried to change this in the code, the gamma analysis stopped printing as well?



Attached below is my code

x = (np.linspace(-1.9875,1.9875,x_voxel))
y = (np.linspace(-1.9875,1.9875,y_voxel))
z = (np.linspace(0.0125,29.9875,z_voxel))

coords = (x,y,z)

axes_reference = coords
axes_evaluation = coords

dose_reference = reduced_dose
dose_evaluation = sens_combined

(z_ref, y_ref, x_ref) = axes_reference
(z_eval, y_eval, x_eval) = axes_evaluation

gamma_options = {
    'dose_percent_threshold': 3,
    'distance_mm_threshold': 2,
    'lower_percent_dose_cutoff': 2.8,
    'interp_fraction': 10,  # Should be 10 or more for more accurate results
    'max_gamma': 3,
    'random_subset': None,
    'local_gamma': True,
    'ram_available': 2**29  # 1/2 GB
gamma = pymedphys.gamma(
    axes_reference, dose_reference, 
    axes_evaluation, dose_evaluation, 

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.ylabel("Probability Density")

# plt.savefig('gamma_hist.png', dpi=300)

########### GAMMA AT VARIOUS DEPTHS ##############
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, 0)) > 
slice_start = np.max([
slice_end = np.min([

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)
    fig.colorbar(c00, ax=ax[0,0], label='Dose (Gy)')
    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)
    fig.colorbar(c01, ax=ax[0,1], label='Dose (Gy)')
    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].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'))
        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].set_xlabel('x (mm)')
    ax[1,1].set_ylabel('z (mm)')

Hi Nevin,

Would you be able to trim down your example until it is a minimal amount of code that reproduces your issue? Also, if you could make it self contained, as in, have the code example generate the initial conditions so that it can be run without any external data or external code.


TestFile for (4.4 KB)

Hi Simon,

I have attached the python file, I haven’t been able to trim it down but I have been able to figure out where I think the problem is occurring, I have note this down on the python file itself.