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?

Thanks,

Nevin

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, 
    **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.figure(2)
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.xlabel("Gamma")
plt.ylabel("Probability Density")
plt.show()
# 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)) > 
    lower_dose_cutoff)
slice_start = np.max([
        np.where(relevant_slice)[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)
    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('z (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('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].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(
        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].invert_yaxis()
    ax[1,1].set_xlabel('x (mm)')
    ax[1,1].set_ylabel('z (mm)')
    plt.show()
    print("\n")
              
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.
Cheers,

Simon
              
TestFile for Simon.py (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.
Thanks,
Nevin
              
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)) >  
    lower_dose_cutoff) 
Would this make it easier to analyse?
Thanks,
Nevin
 nko41:
Unfortunately, every time I tried to cut down any of the code in the example, python kept throwing up error messages.
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.
              
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.
Thanks,
Nevin
 nko41:
ValueError: Length of items in axes_reference ((240, 80, 32)) does not match the shape of dose_reference ((32, 80, 240))
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:
 nko41:
IndexError: index 35 is out of bounds for axis 0 with size 32
Also, when copying it in write the following around it so that it formats nicely:
```python
<your traceback goes here>
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 Simon.py", line 125, in <module>
    c00 = ax[0,0].contourf(
  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\__init__.py", line 1361, in inner
    return func(ax, *map(sanitize_sequence, args), **kwargs)
  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\axes\_axes.py", line 6434, in contourf
    contours = mcontour.QuadContourSet(self, *args, **kwargs)
  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\contour.py", line 777, in __init__
    kwargs = self._process_args(*args, **kwargs)
  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\contour.py", line 1366, in _process_args
    x, y, z = self._contour_args(args, kwargs)
  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\contour.py", line 1424, in _contour_args
    x, y, z = self._check_xyz(args[:3], kwargs)
  File "C:\Users\Nevin Koshy\anaconda3\lib\site-packages\matplotlib\contour.py", 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)```
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 Simon.py", 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 Simon.py", 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 ,

Simon
              
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.
Cheers!
              
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 Simon.py (4.2 KB)
image910×288 28.6 KB
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 .
Cheers ,

Simon
              
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?
Thanks,
Nevin