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.



Hi Simon,

Thank you for the sscce example, that’s actually really helpful in isolating my problem. I tried to cut down my code based on the short, self-contained, correct example. Unfortunately, every time I tried to cut down any of the code in the example, python kept throwing up error messages.

I then realised there was a simpler way to explain my issue, I am using the exact same code as used in the gamma from DICOM example, Gamma from DICOM — PyMedPhys.

I have only made changes to the first 12 lines of code in the example to get the gamma function to work for numpy arrays instead of DICOM images. i.e.

x = (np.linspace(-1.9875,1.9875,32))
y = (np.linspace(-1.9875,1.9875,80))
z = (np.linspace(0.0125,29.9875,240))

coords = (x,y,z)

axes_reference = coords
axes_evaluation = coords

reduced_dose = np.load('reduced_dose.npy')
sens_combined = np.load('sens_combined.npy')

dose_reference = reduced_dose
dose_evaluation = sens_combined

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

I also changed the plane defined by the axis from (1,2) to (1,0) in the line:

relevant_slice = ( 
    np.max(dose_reference, axis=(1, 0)) >  

Would this make it easier to analyse?



Yup, that is certainly much simpler. And something I can quickly parse in a short period of time. Much more likely to get a prompt response :slightly_smiling_face:

Would you be able to copy in some of the error messages that you believe would be most relevant? It helps others who are “Google Searching” for that error message, and it also helps me diagnose the issue.

I suspect that the issue is a “coords order” issue, but do please still copy in the error messages.

So that others following along can know what @nko41 was referring to, I asked Nevin to try and create a “Short, Self Contained, Correct Example”. Below is the link:

Hi Simon,

This is the error message that comes up when I change the value of coords to zyx
coords = (z,y,x)

ValueError: Length of items in axes_reference ((240, 80, 32)) does not match the shape of dose_reference ((32, 80, 240))

This is the error message that comes up when I change the order of variables in axes_reference and axes_evaluation while keeping coords =(x,y,z)

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

IndexError: index 35 is out of bounds for axis 0 with size 32

I tried to remove the gamma pass rate plot by the deleting the code for the pass rate plot, however this threw the error

NameError: name ‘gamma’ is not defined

Although I think this error only occurred because I removed the gamma variable in the gamma pass rate section. I doubt that its actually contributing to the problem, I agree that the issue must be a coords order issue as well.



What happens if you leave the coordinate order as z, y, x but run:

reduced_dose = np.swapaxis(reduced_dose, 0, 2)
sens_combined = np.swapaxis(sens_combined, 0, 2)

Can you copy in the full traceback for the following error:

Also, when copying it in write the following around it so that it formats nicely:

<your traceback goes here>

When I leave the coordinate order as

coords = (z,y,x)

and use the np.swapaxes (swapaxis doesn’t exist), I get the following error

Traceback (most recent call last):

  File "C:\Users\Nevin Koshy\AppData\Local\Packages\CanonicalGroupLimited.Ubuntu20.04onWindows_79rhkp1fndgsc\LocalState\rootfs\home\nevinkoshy\topas\masters_practice\TestFile for", line 125, in <module>
    c00 = ax[0,0].contourf(

  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\", line 1361, in inner
    return func(ax, *map(sanitize_sequence, args), **kwargs)

  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\axes\", line 6434, in contourf
    contours = mcontour.QuadContourSet(self, *args, **kwargs)

  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\", line 777, in __init__
    kwargs = self._process_args(*args, **kwargs)

  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\", line 1366, in _process_args
    x, y, z = self._contour_args(args, kwargs)

  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\", line 1424, in _contour_args
    x, y, z = self._check_xyz(args[:3], kwargs)

  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\", line 1465, in _check_xyz
    raise TypeError(f"Length of x ({nx}) must match number of "

TypeError: Length of x (240) must match number of columns in z (32)```

Also for the index error.

Traceback (most recent call last):

  File "C:\Users\Nevin Koshy\AppData\Local\Packages\CanonicalGroupLimited.Ubuntu20.04onWindows_79rhkp1fndgsc\LocalState\rootfs\home\nevinkoshy\topas\masters_practice\TestFile for", line 93, in <module>
    eval_slices = [

  File "C:\Users\Nevin Koshy\AppData\Local\Packages\CanonicalGroupLimited.Ubuntu20.04onWindows_79rhkp1fndgsc\LocalState\rootfs\home\nevinkoshy\topas\masters_practice\TestFile for", line 94, in <listcomp>
    dose_evaluation[np.where(z_i == z_eval)[0][0], :, :]   #where I think the issue is occuring

IndexError: index 35 is out of bounds for axis 0 with size 32

Hi Nevin,

It looks like order (z, y, x) is what you want. Both of those errors you’ve provided actually have nothing to do with the pymedphys gamma calculation, instead they are to do with indexing numpy arrays and utilisation of matplotlib.

You’ll need to investigate your usage of those libraries and whether or not the code used is relevant to the dataset you have provided.

Cheers :slight_smile:,

Hi Simon,

Ah you’re right, my axes reference and axes evaluation had the wrong order, I had ordered it as x,y,z and should have ordered them z,y,x. When I do this, the gamma functionality works perfectly.


1 Like

Awesome to hear :slight_smile:. If you do end up with a nice simple working self contained 3D example we can put your example up on the docs if you like to help others?

Hi Simon,

I do have a working example now and I would be happy to put it up on the docs page. I have attached the test file below. It does require the external files for reduced dose and sens_combined that I have sent earlier, I can’t attach the files as they exceed the 4 mb size limit. TestFile for (4.2 KB)

Hi @nko41,

Can you first verify that there is no identifying information within the DICOM files, and then would you be able to upload your files to the PyMedPhys data repo:

Using the following button:

You can then make your code download those files at the very beginning, therefore making the script self contained.

Lastly, to be able to include it within the docs, for these sorts of walk-through tutorials we use Jupyter Notebook. Might you be able to write up your script as a notebook making sure to include explanatory text with the intent to make it clear what is going on at each step. The target audience would be a PyMedPhys new comer.

Also, feel free to add a link to your own personal website or social media/GitHub profile somewhere in there so readers of the document can recognise who put the work into writing it :slight_smile:.

Cheers :slightly_smiling_face:,

Hi Simon,

Apologies for the delay, have been busy with some university work, I can definitely work on a Jupyter notebook tutorial. I’ll start working on it now. I have tried to upload the files onto github but it says that uploads are disabled?



No need to apologise at all :slight_smile:.

Hmm, I’m not sure why it’s saying it’s disabled. @cedricbelanger just did it two days ago… @cedricbelanger do you have any tips for @nko41 to get it working?