Gamma analysis fails simple unit test

[TimosBrain]

The gamma function fails for a very simple unit test in 2D. Please see the test code:

import numpy as np
import pymedphys
import matplotlib.pyplot as plt
refdata = np.full([11, 11, 11], 1., dtype=np.float32)
refx = np.arange(-5., 5.1, 1, dtype=np.float32)
refy = np.arange(-5., 5.1, 1, dtype=np.float32)
refz = np.arange(-5., 5.1, 1, dtype=np.float32)
# test data with twice the resolution and one voxel with dose
evaldata = np.zeros([25, 25, 25], dtype=np.float32)
evaldata[12, 12, 12] = 1.
assert evaldata[12, 12, 12] == 1.
evalx = np.arange(-6, 6.1, 0.5, dtype=np.float32)
evaly = np.arange(-6, 6.1, 0.5, dtype=np.float32)
evalz = np.arange(-6, 6.1, 0.5, dtype=np.float32)
nx=3
ny=3
refplane = np.full([11,11], 1., dtype=np.float32)
evalplane = np.zeros([25,25], dtype=np.float32)
evalplane[12, 12] = 1
axes_ref = (refz, refy, refx)
axes_eval = (evalz, evaly, evalx)
fig2d, ax2d = plt.subplots(3, 3)
im = []
for i in range(nx):
    for j in range(ny):
        nint = 1 + i*nx +j
        gamma_options = {
                    "dose_percent_threshold": 3,
                    "distance_mm_threshold": np.sqrt(2),
                    "lower_percent_dose_cutoff": 0,
                    "interp_fraction": nint,
                    "local_gamma": False,
                    "max_gamma": 5}

        gamma_map = pymedphys.gamma(axes_reference=(refy, refx),
                                    dose_reference= refplane,
                                    axes_evaluation = (evaly, evalx),
                                    dose_evaluation= evalplane,
                                    **gamma_options)
        print(gamma_map)

        im.append(ax2d[i][j].imshow(gamma_map))
        ax2d[i][j].set_title("interpolation="+str(nint))
plt.show()

It compares a homogeneous dose distribution of 1Gy (reference) with grid size 1mm to an evaluation dose distribution with twice the resolution which is 0 everywhere except for the central point (0,0), where it is 1Gy. The distance criterion is sqrt(2)mm. So I would expect that gamma=0 for the central point (which it is) and monotonically increasing with increasing distance up to the maximum gamma set as parameter. The direct neighbors along the axis and diagonally should pass. In any case, the gamma map should be isotropic around the central point. For increasing interpolation, the gamma map should converge against the correct value. In the attached unit test, I try interpolation factors between 1 and 9 and the resulting gamma maps look all totally different and not a single one of them fulfills my expectations. It seems like something is going totally wrong here...

[SimonBiggs]

Would you be able to copy some exemplar figures into this issue

[SimonBiggs]

I also edited your issue to include the file inline.

[SimonBiggs]

Can you also provide a result where you set the interp fraction to something like 20?

[SimonBiggs]

As well as provide the geometric result that you expect, and plot the difference?

[TimosBrain]

I slightly modified the above code to show higher interpolation values. By simply looping over the evaluation points in the environment of the reference point grids, I get to following, expected gamma map: with 0.0 in the center, 1/sqrt(2) for the straight neighbors and 1.0 for the neighbors on the diagonals

[SimonBiggs]

May you include colour bars on your figures as well as set vmin and vmax so that the same colour scale is used across all of the figures?

[TimosBrain]

For interactive figures, that even show you the exact value of each pixel when hovering with the mouse over it, please simply execute the above minimal code snippet or modify it to your needs. This is even better than color scales. By adapting the color range to each plot individually, it is even more obvious from the plots, that the gamma map does not fulfill the expected symmetries or that the gamma often decreases when the distance to agreement increases. As long as those issues exist, it makes little sense to compare single pixel values among different plots with different interpolation factors.

[SimonBiggs]

Kk, might you be in a position to dig into the code to determine what the issue might be?

On Mon, 20 Nov 2023, 10:08 pm TimosBrain, @.***> wrote:

For interactive figures, that even show you the exact value of each pixel when hovering with the mouse over it, please simply execute the above minimal code snippet or modify it to your needs. This is even better than color scales. By adapting the color range to each plot individually, it is even more obvious from the plots, that the gamma map does not fulfill the expected symmetries or that the gamma often decreases when the distance to agreement increases. As long as those issues exist, it makes little sense to compare single pixel values among different plots with different interpolation factors.

— Reply to this email directly, view it on GitHub https://github.com/pymedphys/pymedphys/issues/1813#issuecomment-1818841646, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABSBK65L4JMARBVPKNVC2DLYFM22FAVCNFSM6AAAAAA67SKTFWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMJYHA2DCNRUGY . You are receiving this because you commented.Message ID: @.***>

[SimonBiggs]

Is the issue simply that because you have chosen eval to be 0 everywhere and 1 at the centre that when the interpolation is searching for the closest gamma value within the eval plane that it is effectively able to interpolate to be anywhere between 0 or 1 and find the closest point to the dose reference within that range?

And then the gamma value just becomes impacted by the distance away from the centre point with some interpolation artifacts?

[SimonBiggs]

Hi @TimosBrain,

Does the above clarify this issue? Would it be okay to close this issue report? It appears that in fact the gamma implementation is working as expected when interpolation is taken into account.

[TimosBrain]

Hi Simon, unfortunately, I don't have time to further look into the code. But I don't see any clear evidence, that the gamma implementation really works except for interpolation artifacts that would be negligible in a real case. Could you please provide some?

[SimonBiggs]

Here are the tests that replicate the testing suite from the following paper: http://dx.doi.org/10.1016/j.radonc.2015.11.034

https://github.com/pymedphys/pymedphys/blob/f6acdf9bd2e8a32e372966879284fbd71c612358/lib/pymedphys/tests/gamma/test_agnew_mcgarry.py

Data downloadable from here: https://zenodo.org/record/3870436/files/gamma_test_data.zip?download=1

An indication for why interpolation is helpful in some clinical scenarios is provided in the following paper: http://dx.doi.org/10.1118/1.2721657

[SimonBiggs]

Given the gamma algorithm has been implemented as intended and the original concern has been addressed I will close this issue unless there is a desire for it to be reopened.