Merged in bugfix/RAM-3754_reverse_order_gamma_profile (pull request #474

)

RAM-3754 Handle reversed-order x-values for geometric gamma.

Approved-by: Randy Taylor

docs/source/changelog.rst

@@ -44,6 +44,10 @@ Gamma
 * :bdg-warning:`Fixed` Evaluating gamma on a profile that had at least one x-spacing that was double the distance-to-agreement
   parameter would result in an error. This is because for computational speed, only nearby simplexes were included in the
   evaluation. However, if the x-spacing was too large, the nearest simplex would end up being the same point as the evaluation.
+* :bdg-warning:`Fixed` The geometric gamma function would fail if the x-values if the reference coordinates or evaluation coordinates
+  were not in ascending order. Coordinates can now be in any order.
+* :bdg-warning:`Fixed` The geometric gamma function now validates that the x-coordinates are monotonically increasing or decreasing and
+  will error out if not.
 
 Core
 ^^^^

docs/source/topics/gamma.rst

@@ -161,7 +161,9 @@ of the simplex support. In the 1D case, we can simply take the minimum distance
 Logic
 ~~~~~
 
-#. The reference and evaluation distributions are both normalized by the maximum reference value * the dose to agreement parameter.
+#. The reference and evaluation distributions are both normalized in the y-axis by the maximum reference value * the dose to agreement parameter.
+   Note the axes in Fig 1 of the paper above.
+#. The reference and evaluation distributions are both normalized in the x-axis by the distance to agreement parameter.
 #. For each evaluation point :math:`p`:
 
    #. The vertices of each simplex are found. In the 1D scenario these will always be line segments defined by 2 points. All the points that are within and just beyond

pylinac/core/array_utils.py

@@ -416,3 +416,21 @@ def _rt_image_position(array: np.ndarray, dpmm: float) -> list[float, float]:
     x_position = -(width_mm / 2) + (pixel_size_mm / 2)
     y_position = -(height_mm / 2) + (pixel_size_mm / 2)
     return [x_position, y_position]
+
+
+@validate(array=(array_not_empty, single_dimension))
+def is_monotonically_increasing(array: np.ndarray) -> bool:
+    """Check if an array is monotonically increasing."""
+    return np.all(np.diff(array) > 0)
+
+
+@validate(array=(array_not_empty, single_dimension))
+def is_monotonically_decreasing(array: np.ndarray) -> bool:
+    """Check if an array is monotonically decreasing."""
+    return np.all(np.diff(array) < 0)
+
+
+@validate(array=(array_not_empty, single_dimension))
+def is_monotonic(array: np.ndarray) -> bool:
+    """Check if an array is monotonically increasing or decreasing."""
+    return is_monotonically_increasing(array) or is_monotonically_decreasing(array)

pylinac/core/gamma.py

@@ -7,6 +7,10 @@
 from skimage.draw import disk
 
 from . import validators
+from .array_utils import (
+    is_monotonic,
+    is_monotonically_decreasing,
+)
 
 
 def _construct_matrices(
@@ -150,12 +154,20 @@ def gamma_geometric(
         raise ValueError("Distance to agreement must be greater than 0")
     if reference_coordinates is None:
         reference_coordinates = np.arange(len(reference), dtype=float)
+    if not is_monotonic(reference_coordinates):
+        raise ValueError(
+            "Reference x-values must be monotonically increasing or decreasing"
+        )
     if len(reference) != len(reference_coordinates):
         raise ValueError(
             f"Reference and reference_x_values must be the same length. Got reference: {len(reference)} and reference_x_values: {len(reference_coordinates)}"
         )
     if evaluation_coordinates is None:
         evaluation_coordinates = np.arange(len(evaluation), dtype=float)
+    if not is_monotonic(evaluation_coordinates):
+        raise ValueError(
+            "Evaluation x-values must be monotonically increasing or decreasing"
+        )
     if len(evaluation) != len(evaluation_coordinates):
         raise ValueError(
             f"Evaluation and evaluation_x_values must be the same length. Got evaluation: {len(evaluation)} and evaluation_x_values: {len(evaluation_coordinates)}"
@@ -199,9 +211,15 @@ def gamma_geometric(
         # not change the gamma calculation.
         # This sub-sampling of the vertices is all for computational efficiency.
         left_diffs = np.abs(normalized_reference_x - (eval_x - distance_to_agreement))
+        right_diffs = np.abs(normalized_reference_x - (eval_x + distance_to_agreement))
+        if is_monotonically_decreasing(normalized_reference_x):
+            # it can be the case that the x-values go from high to low vs low to high.
+            # If the profiles are decreasing, we need to swap the left and right differences
+            # We cannot sort/flip the array itself because then the gamma order
+            # would be reversed from the profile order.
+            left_diffs, right_diffs = right_diffs, left_diffs
         # we need to ensure we don't go out of bounds if evaluating at the edge, hence the max/min
         left_idx = max(np.argmin(left_diffs) - 1, 0)
-        right_diffs = np.abs(normalized_reference_x - (eval_x + distance_to_agreement))
         right_idx = min(np.argmin(right_diffs) + 1, len(normalized_reference) - 1)
         # the vertices are the (x, y) pairs of the reference profile
         vertices_x = normalized_reference_x[left_idx : right_idx + 1].tolist()

tests_basic/core/test_array_utils.py

@@ -16,6 +16,9 @@
     geometric_center_value,
     ground,
     invert,
+    is_monotonic,
+    is_monotonically_decreasing,
+    is_monotonically_increasing,
     normalize,
     stretch,
 )
@@ -330,3 +333,35 @@ def test_dtypes(self, input_dtype, output_dtype):
             ds_reload = pydicom.dcmread(f)
         self.assertTrue(np.allclose(ds_reload.pixel_array, arr))
         self.assertEqual(ds_reload.pixel_array.dtype, output_dtype)
+
+
+class TestMonotonic(TestCase):
+    @parameterized.expand(
+        [
+            [np.array([1, 2, 3, 4, 5]), True],
+            [np.array([5, 4, 3, 2, 1]), False],
+            [np.array([1, 1, 2, 3, 4]), False],
+        ]
+    )
+    def test_monotonic_increasing(self, arr, expected):
+        self.assertEqual(is_monotonically_increasing(arr), expected)
+
+    @parameterized.expand(
+        [
+            [np.array([1, 2, 3, 4, 5]), False],
+            [np.array([5, 4, 3, 2, 1]), True],
+            [np.array([4, 3, 3, 2, 1]), False],
+        ]
+    )
+    def test_monotonic_decreasing(self, arr, expected):
+        self.assertEqual(is_monotonically_decreasing(arr), expected)
+
+    @parameterized.expand(
+        [
+            [np.array([1, 2, 3, 4, 5]), True],
+            [np.array([5, 4, 3, 2, 1]), True],
+            [np.array([4, 3, 3, 2, 1]), False],
+        ]
+    )
+    def test_is_monotonic(self, arr, expected):
+        self.assertEqual(is_monotonic(arr), expected)

tests_basic/core/test_gamma.py

@@ -521,6 +521,68 @@ def test_at_left_edge(self, ref, eval, expected):
         )
         self.assertTrue(np.array_equal(g, expected, equal_nan=True))
 
+    @parameterized.expand([(np.arange(4, -5, -1),), (np.arange(-4, 5, 1),)])
+    def test_reversed_x_values_have_same_result(self, x):
+        # test when the x values are reversed
+        ref = np.asarray([0, 1, 3, 4, 5, 4, 4, 1, 0])
+        eval = np.asarray([0, 1, 3, 4, 5, 4, 3, 1, 0])
+        g = gamma_geometric(
+            reference=ref,
+            reference_coordinates=x,
+            evaluation=eval,
+            evaluation_coordinates=x,
+            dose_to_agreement=1,
+            distance_to_agreement=1,
+        )
+        self.assertTrue(
+            np.allclose(
+                g,
+                [np.nan, 0, 0, 0, 0, 0, 0.3332, 0, np.nan],
+                equal_nan=True,
+                atol=0.001,
+            )
+        )
+
+    @parameterized.expand([(np.arange(4, -5, -1),), (np.arange(-4, 5, 1),)])
+    def test_reversed_x_values_for_one_profile(self, x):
+        # test that when one profile is reversed, the gamma is still the same as if they were both in the same order
+        ref = np.asarray([0, 1, 3, 4, 5, 4, 4, 1, 0])
+        eval = np.asarray([0, 1, 3, 4, 5, 4, 3, 1, 0])
+        g = gamma_geometric(
+            reference=ref,
+            reference_coordinates=x,
+            evaluation=eval,
+            evaluation_coordinates=np.flip(x),
+            dose_to_agreement=1,
+            distance_to_agreement=1,
+        )
+        self.assertTrue(
+            np.allclose(
+                g,
+                [np.nan, 0, 0.3332, 0, 0, 0, 0, 0, np.nan],
+                equal_nan=True,
+                atol=0.001,
+            )
+        )
+
+    @parameterized.expand(
+        [
+            ("bad eval", [0, 1, 2, 3, 4], [1, 3, 4, 3, 6]),
+            ("bad ref", [1, 3, 4, 3, 6], [0, 1, 2, 3, 4]),
+            ("bad both", [1, 3, 4, 3, 6], [1, 3, 4, 3, 6]),
+        ]
+    )
+    def test_non_monotonic_fails(self, _, ref_x, eval_x):
+        ref = eval = np.ones(5)
+        with self.assertRaises(ValueError) as e:
+            gamma_geometric(
+                reference=ref,
+                reference_coordinates=np.array(ref_x),
+                evaluation=eval,
+                evaluation_coordinates=np.array(eval_x),
+            )
+        self.assertIn("monotonically", str(e.exception))
+
 
 class TestGamma1D(TestCase):
     def test_resolution_below_1mm(self):