From fd6043bf8c7cd1fc6b80134e6f224e28b937525f Mon Sep 17 00:00:00 2001 From: James Kerns Date: Thu, 7 Nov 2024 10:45:13 -0600 Subject: [PATCH] handle reversed-order x-values for geometric gamma. --- docs/source/changelog.rst | 4 ++ docs/source/topics/gamma.rst | 4 +- pylinac/core/array_utils.py | 18 ++++++++ pylinac/core/gamma.py | 20 ++++++++- tests_basic/core/test_array_utils.py | 35 ++++++++++++++++ tests_basic/core/test_gamma.py | 62 ++++++++++++++++++++++++++++ 6 files changed, 141 insertions(+), 2 deletions(-) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 9e1c11a5..d76fc2ca 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -43,6 +43,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 ^^^^ diff --git a/docs/source/topics/gamma.rst b/docs/source/topics/gamma.rst index 93d9af29..3145cb20 100644 --- a/docs/source/topics/gamma.rst +++ b/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 diff --git a/pylinac/core/array_utils.py b/pylinac/core/array_utils.py index 282acde2..07cffaf2 100644 --- a/pylinac/core/array_utils.py +++ b/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) diff --git a/pylinac/core/gamma.py b/pylinac/core/gamma.py index db8db850..c1e6511b 100644 --- a/pylinac/core/gamma.py +++ b/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() diff --git a/tests_basic/core/test_array_utils.py b/tests_basic/core/test_array_utils.py index feb5bceb..da32f8f7 100644 --- a/tests_basic/core/test_array_utils.py +++ b/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) diff --git a/tests_basic/core/test_gamma.py b/tests_basic/core/test_gamma.py index 5e036d2c..3b3c3887 100644 --- a/tests_basic/core/test_gamma.py +++ b/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):