Merged in bugfix/RAM-3754_geometric_gamma_dta (pull request #458)

RAM-3754 fix x-scaling of gamma distances.

Approved-by: Randy Taylor

docs/source/changelog.rst

@@ -16,13 +16,13 @@ v 3.28.0
 General
 ^^^^^^^
 
-* The minimum version of Python supported is now 3.9 as the end-of-life of Python 3.8 is October 2024.
-* Along with this, minimum versions of numpy, scipy, and matplotlib have been bumped to versions that support Python 3.9.
+* :bdg-danger:`Change` The minimum version of Python supported is now 3.9 as the end-of-life of Python 3.8 is October 2024.
+* :bdg-danger:`Change` Along with this, minimum versions of numpy, scipy, and matplotlib have been bumped to versions that support Python 3.9.
 
 Plan Generator
 ^^^^^^^^^^^^^^
 
-* The patient name and patient ID can now be passed to the :class:`~pylinac.plan_generator.dicom.PlanGenerator` and :class:`~pylinac.plan_generator.dicom.HalcyonPlanGenerator` classes. This is useful for
+* :bdg-success:`Feature` The patient name and patient ID can now be passed to the :class:`~pylinac.plan_generator.dicom.PlanGenerator` and :class:`~pylinac.plan_generator.dicom.HalcyonPlanGenerator` classes. This is useful for
   setting the patient name and ID in the generated DICOM files.
 
 Starshot
@@ -33,11 +33,22 @@ Starshot
 CT
 ^^
 
-* The Tomo and CatPhan localization algorithm has changed slightly to handle very high HU values such as
+* :bdg-warning:`Fixed` The Tomo and CatPhan localization algorithm has changed slightly to handle very high HU values such as
   bbs and metal rods. In some instances, the HU values were so high the localization algorithm did not
   detect the phantom as it appeared at a background level. Values are now capped to -1000 and +1000 HU.
   Note this is only relevant for localization. HU values are not changed or capped in the analysis itself.
 
+Gamma
+^^^^^
+
+* :bdg-warning:`Fixed` Geometric gamma was incorrectly calculating the minimum threshold to evaluate over due to scaling of the dataset by the
+  dose to agreement. E.g. if the dose to agreement was 2mm, the minimum threshold for evaluating gamma
+  would be 20% of the dose. This would cause the evaluation points calculated over to be fewer than desired.
+* :bdg-warning:`Fixed` Geometric gamma was incorrectly calculating the gamma value for the distance to agreement due to a lack of scaling
+  by the distance to agreement. This would cause gamma values to be much higher than expected and thus gamma
+  pass rates to be lower than expected.
+* :bdg-warning:`Fixed` Dose to agreement and distance to agreement must be greater than 0.
+
 Core
 ^^^^
 

pylinac/core/gamma.py

@@ -144,6 +144,10 @@ def gamma_geometric(
         raise ValueError(
             f"Reference and evaluation arrays must be 1D. Got reference: {reference.ndim} and evaluation: {evaluation.ndim}"
         )
+    if distance_to_agreement <= 0:
+        raise ValueError("Dose to agreement must be greater than 0")
+    if dose_to_agreement <= 0:
+        raise ValueError("Distance to agreement must be greater than 0")
     if reference_coordinates is None:
         reference_coordinates = np.arange(len(reference), dtype=float)
     if len(reference) != len(reference_coordinates):
@@ -165,18 +169,22 @@ def gamma_geometric(
         raise ValueError(
             "The evaluation x-values must be within the range of the reference x-values"
         )
-    threshold = float(dose_threshold)
-    # convert dose to normalized distance of dose to agreement. I.e. D/delta(D)
+    # normalize the dose threshold by the DTA
+    threshold = float(dose_threshold) / float(dose_to_agreement)
+    # convert dose to normalized distance of dose to agreement. I.e. D/delta(D) in Figure 1.
     normalized_reference = (
         reference.astype(float) * 100 / (reference.max() * dose_to_agreement)
     )
     normalized_evaluation = (
         evaluation.astype(float) * 100 / (reference.max() * dose_to_agreement)
     )
+    # normalize the x-values; i.e. X/delta(d) in Figure 1.
+    normalized_reference_x = reference_coordinates / distance_to_agreement
+    normalized_evaluation_x = evaluation_coordinates / distance_to_agreement
 
     gamma = np.full(len(evaluation), fill_value)
     for idx, (eval_x, eval_point) in enumerate(
-        zip(evaluation_coordinates, normalized_evaluation)
+        zip(normalized_evaluation_x, normalized_evaluation)
     ):
         # skip if below dose threshold
         if eval_point < threshold:
@@ -185,13 +193,13 @@ def gamma_geometric(
         # we need to grab the vertices just beyond the edge of the DTA
         # so we evaluate within the entire DTA range
         left_idx = np.argmin(
-            np.abs(reference_coordinates - (eval_x - distance_to_agreement))
+            np.abs(normalized_reference_x - (eval_x - distance_to_agreement))
         )
         right_idx = np.argmin(
-            np.abs(reference_coordinates - (eval_x + distance_to_agreement))
+            np.abs(normalized_reference_x - (eval_x + distance_to_agreement))
         )
         # the vertices are the (x, y) pairs of the reference profile
-        vertices_x = reference_coordinates[left_idx : right_idx + 1].tolist()
+        vertices_x = normalized_reference_x[left_idx : right_idx + 1].tolist()
         vertices_y = normalized_reference[left_idx : right_idx + 1].tolist()
         vertices = list(np.array((x, y)) for x, y in zip(vertices_x, vertices_y))
         # iterate in pairs of x, y

tests_basic/core/test_gamma.py

@@ -2,6 +2,7 @@
 from unittest import TestCase, skip
 
 import numpy as np
+from parameterized import parameterized
 
 from pylinac.core.gamma import _compute_distance, gamma_1d, gamma_2d, gamma_geometric
 from pylinac.core.image import DicomImage
@@ -440,6 +441,53 @@ def test_non_1d_array(self):
         with self.assertRaises(ValueError):
             gamma_geometric(reference=ref, evaluation=eval)
 
+    @parameterized.expand(
+        [
+            (1, 50),
+            (2, 50),
+            (3, 100),
+            (4, 100),
+            (5, 100),
+        ]
+    )
+    def test_distance_scaling(self, dist: float, pass_rate: float):
+        # see RAM-3754
+        ref = np.asarray([0, 0.5, 1, 1, 1, 1])
+        # eval is 3mm off. We want to test that gamma is <100 <3mm, then goes to 100 when we reach >=3mm
+        eval = np.asarray([0, 0, 0, 0, 0.5, 1])
+        # for dist in range(1, 5):
+        gamma = gamma_geometric(
+            reference=ref,
+            evaluation=eval,
+            dose_to_agreement=0.1,  # large dose to ensure it's not the limiting factor
+            distance_to_agreement=dist,
+        )
+        valid_elements = gamma[~np.isnan(gamma)]
+        measured_pass_rate = 100 * (np.sum(valid_elements < 1) / len(valid_elements))
+        self.assertAlmostEqual(measured_pass_rate, pass_rate, delta=1)
+
+    def test_positive_distance_to_agreement(self):
+        ref = np.ones(5)
+        eval = np.ones(5)
+        with self.assertRaises(ValueError):
+            gamma_geometric(
+                reference=ref,
+                evaluation=eval,
+                dose_to_agreement=1,
+                distance_to_agreement=-1,
+            )
+
+    def test_postive_dose_to_agreement(self):
+        ref = np.ones(5)
+        eval = np.ones(5)
+        with self.assertRaises(ValueError):
+            gamma_geometric(
+                reference=ref,
+                evaluation=eval,
+                dose_to_agreement=-1,
+                distance_to_agreement=1,
+            )
+
 
 class TestGamma1D(TestCase):
     def test_resolution_below_1mm(self):
@@ -667,4 +715,4 @@ def test_different_epids(self):
             reference_profile=p1200_prof, dose_to_agreement=1, gamma_cap_value=2
         )
         # gamma is very low; just pixel noise from the image generator
-        self.assertAlmostEqual(np.nanmean(gamma), 0.948, delta=0.01)
+        self.assertAlmostEqual(np.nanmean(gamma), 0.938, delta=0.01)