Gamma between 3D ref and 2D evaluation

Hello and thank you for your high-quality work!

I would like to simulate a patient specific qa (PSQA) as it is done in radiotherapy.
To do so, I would like to compare a 3D reference dose (from a dicom RT DOSE) to a 2D evaluation dose Normaly, the 2D distribution goes from a 2D array detector. For me, it’s just a chosen plane from a 3D distribution.

I want to keep my 3 dimensions evaluation distribution because I would like the gamma to be calculated around the chosen slice and not only 2D/2D.

Is there an easy way to do it?

Thanks in advance!



I tried to find a solution but the result is not what I expected.

My first idea was: As I can compare two distributions of different sizes, maybe I can try to reduce the size of the evaluation distribution until it reaches the thickness of a single slice. Thus, the dimension of my new “hybrid” 3D distribution is (1,y,x) and the axis_evaluation is ([chosen_slice(mm)],[y],[ x ]).

My code can be found there. It creates a new dose_evaluation and a new dose_axis. (7.7 KB)

Everything seems more or less coherent :

My new dose_evaluation looks good

All the verification assertions pass


But the gamma calculation is extremely slow and some errors appear.

I hope that all these details will help to understand my initial request. I hope I didn’t miss something at the beginning.
Now I plan to try to use more than one slice.

Thanks in advance


@Matthew-Jennings, might you be in a position to give @julien103 a hand?

Hi @julien103,

Apologies for the delay! A couple of thoughts:

  • Since the last official release, we’ve introduced a performance enhancement for Gamma. You can access this now by installing the 0.40.0.dev2 development release of pymedphys as follows:
pip install pymedphys==0.40.0.dev2
  • If you don’t need to obtain an accurate gamma value for values greater than unity (as implied by your code, which sets max_gamma=1.1), you could also try cropping your reference dataset in a similar manner to the cropping you’ve already performed to obtain your evaluation (2D) dataset. As long as you include all voxels in the reference dataset that are within, say, 2× your DTA criterion (i.e., 2 mm for the 1 mm DTA set in the file you provided) of all voxels in the evaluation dataset, you should be good I think. I say 2× instead of 1.1× since I’ve seen some funky stuff happen near such thresholds. Make sure you consider your voxel size in all of this, of course.

I am in no hurry at all and am grateful for your response! Thank you very much.

I tried what you suggested (with the new version of pymedphys). First I added a method that creates a segmented version of the reference distribution. The list containing the axis values is consistent as you can see in the exemple below (sagittal axis chosen)

array([-155.8484, -152.8484, -149.8484])


Here, evaluation and reference distribution are the same. Voxel size is 3 mm x 3 mm x 3 mm. As the dta is set to 1 mm, I chose to keep one slice on each side of the main slice.

The function pymedphys.gamma returns an Assertion error:

  File "X:\DEV\Gamma\", line 100, in calculate_hybrid_2D_gamma
    self.gamma = pymedphys.gamma(self.axis_reference, self.dose_reference,
  File "C:\Users\M11142\Anaconda3\envs\gamma\lib\site-packages\pymedphys\_gamma\implementation\", line 166, in gamma_shell
    current_gamma = gamma_loop(options)
  File "C:\Users\M11142\Anaconda3\envs\gamma\lib\site-packages\pymedphys\_gamma\implementation\", line 353, in gamma_loop
    min_relative_dose_difference = calculate_min_dose_difference(
  File "C:\Users\M11142\Anaconda3\envs\gamma\lib\site-packages\pymedphys\_gamma\implementation\", line 490, in calculate_min_dose_difference
    evaluation_dose = interpolate_evaluation_dose_at_distance(
  File "C:\Users\M11142\Anaconda3\envs\gamma\lib\site-packages\pymedphys\_gamma\implementation\", line 547, in interpolate_evaluation_dose_at_distance
    coords_evaluation_grid = interpolation.splines.CGrid(*options.axes_evaluation)
  File "C:\Users\M11142\Anaconda3\envs\gamma\lib\site-packages\interpolation\splines\", line 27, in CGrid
    assert a.shape[0] >= 2

I tried with 'max_gamma': 2, same error with 'max_gamma': 1.1.

Maybe you can understand where this error comes from :upside_down_face:

my new code there : (12.1 KB)

Thanks in advance!


Hi @julien103,

Try switching the reference and evaluation distributions. The reference distribution is the one whose voxels are used to calculate the pass rate (its voxel coordinates form the centres of each gamma ellipsoid). The evaluation distribution is the one that is searched over. Switching them seems to apply to your situation better. Sorry, my last post may have contributed to this error.

Searching requires interpolation, and your code is getting an error because the underlying interpolator needs the distribution to have at least two voxels along each axis (to interpolate between).