Merged in bugfix/RAM-3829_starshot_stabilization (pull request #431)

RAM-3829 stabilize starshot algorithm

Approved-by: Randy Taylor
Approved-by: Hasan Ammar

docs/source/changelog.rst

@@ -20,6 +20,30 @@ Picket Fence
   to be **relative to the CAX**. Previously, the positions were relative to the left/top of the image. This attribute
   was requested by French sites and this change was also requested by them.
 
+Starshot
+^^^^^^^^
+
+TL;DR: The algorithm is more robust, but manually passing bad starting points will now error.
+
+* :bdg-warning:`Fixed` Sometimes images analyzed with a large radius and/or with recursive analysis would "find" spokes
+  where there were clearly no spokes. This was due to a normalization issue of the sampled ring. The normalization
+  now takes into account the mean value of the central third of the image. This should prevent "ghost" spokes.
+* :bdg-warning:`Fixed` The image is now "grounded" (lowest value set to 0) before analysis. This should prevent
+  issues with images that have a high background value and relatively low signal values.
+* :bdg-primary:`Refactor` Detected spokes that are not within 10mm of the initial starting point will raise an error.
+  This is to prevent poor analyses, most often for gantry starshots, where some spokes are not detected.
+  In combination with ``recursive``, this should provide a more robust analysis.
+
+  .. image:: images/bad_line_match.png
+
+* :bdg-danger:`Change` Due to the internal changes required to fix the above issues, providing a ``start_point`` to the analysis method was usually robust enough to iteratively
+  correct itself. This is no longer the case. Passing a ``start_point`` will now be the only and final point that the ring is
+  centered around. I.e. if you pass a bad starting point you'll get an error most likely. The intent of the
+  start point was to provide a workaround if the automatic centering failed. This is still true but assumes that
+  the user is providing a reasonable starting point (<5mm from the radiation center). A future feature will provide
+  the distance from a start point to the radiation isocenter, allowing the workflow of passing the start point which
+  represents the mechanical isocenter.
+
 Gamma
 ^^^^^
 

docs/source/starshot_docs.rst

@@ -159,30 +159,39 @@ Algorithm
 * The image must have at least 6 spokes (3 angles).
 * The center of the "star" must be in the central 1/3 of the image.
 * The radiation spokes must extend to both sides of the center. I.e. the spokes must not end at the center of the circle.
+* Markings like pin pricks or marker dots can cause wild pixel values and analysis errors and should be cropped out if possible.
 
 **Pre-Analysis**
 
-
-* **Check for image noise** -- The image is checked for unreasonable noise by comparing the min and max
-  to the 1/99th percentile pixel values respectively. If there is a large difference then there is likely an artifact
-  and a median filter is applied until the min/max and 1/99th percentiles are similar.
+* **Ground** -- The image is "grounded" by setting the lowest pixel value to 0. This is done to avoid issues with
+  high-background pixel values.
 * **Check image inversion** -- The image is checked for proper inversion using histogram analysis.
-* **Set algorithm starting point** -- Unless the user has manually set the pixel location of the start point,
-  it is automatically found by summing the image along each axis and finding the
-  center of the full-width, 80%-max of each sum. The maximum value point is also located. Of the two points, the
-  one closest to the center of the image is chosen as the starting point.
+* **Set algorithm starting point** -- If the user provided a ``start_point`` to the ``analyze`` method,
+  this value is used. Otherwise, the start point is automatically found by examining the central 1/3 of the image and taking the max of that image window along each axis and finding the
+  center of the full-width, 80%-max of each profile.
 
 **Analysis**
 
 * **Extract circle profile** -- A circular profile is extracted from the image centered around the starting point
   and at the radius given.
-* **Find spokes** -- The circle profile is analyzed for peaks. Optionally, the profile is reanalyzed to find the center
-  of the FWHM. An even number of spokes must be found (1 for each side; e.g. 3 collimator angles should produce 6
-  spokes, one for each side of the CAX).
+* **Find spokes** -- The circle profile is analyzed for peaks (either simple peaks or FWHM peaks depending on the ``fwhm`` parameter).
+  An even number of spokes (1 for each side; e.g. 3 collimator angles should produce 6
+  spokes, one for each side of the CAX) and at least 6 spokes (3 angles) must be found or an error is raised. If ``recursive`` is set to true, this process
+  is repeated for the radii range 0.1-0.95 and from 0.05-0.95 ``min_peak_height``. If no combination of these
+  parameters produces a valid set of spokes, an error is raised.
 * **Match peaks** -- Peaks are matched to their counterparts opposite the CAX to compose a line using a simple peak number offset.
+  After matching, the distances from the lines to the start point is measured. If the distance is greater than 10mm of
+  any line, an error is raised and follows the same logic regarding ``recursive`` as the previous step. This is to
+  avoid issues specifically with gantry starshots. This can also be corrected ahead of time by using a lower ``min_peak_height``.
+
+  .. figure:: images/bad_line_match.png
+
+      The above image has an even number of spokes, but has not caught the low-dose end of the gantry spokes on the right.
+
+
 * **Find wobble** -- Starting at the initial starting point, a Nelder-Mead gradient method is utilized to find the
-  point of minimum distance to all lines. If recursive is set to True and a "reasonable" wobble (<2mm) is not found
-  using the passes settings, the peak height and radius are iterated until a reasonable wobble is found.
+  point of minimum distance to all lines. The wobble is assumed to be found if the wobble diameter is <2mm and is <10 mm from the
+  starting point. If not, the same logic as the previous steps is followed regarding the ``recursive`` parameter.
 
 **Post-Analysis**
 

pylinac/core/profile.py

@@ -2080,7 +2080,7 @@ def find_fwxm_peaks(
         _, peak_props = find_peaks(
             self.values,
             threshold=threshold,
-            min_width=min_distance,
+            peak_separation=min_distance,
             max_number=max_number,
             search_region=search_region,
             peak_sort=peak_sort,

pylinac/starshot.py

@@ -182,7 +182,7 @@ def from_zip(cls, zip_file: str, **kwargs):
             else:
                 return cls(image_files[0], **kwargs)
 
-    def _get_reasonable_start_point(self) -> Point:
+    def _get_reasonable_start_point(self) -> (Point, float):
         """Set the algorithm starting point automatically.
 
         Notes
@@ -198,8 +198,11 @@ def _get_reasonable_start_point(self) -> Point:
         right_third = int(left_third * 2)
         central_array = self.image.array[top_third:bottom_third, left_third:right_third]
 
-        x_sum = np.sum(central_array, 0)
-        y_sum = np.sum(central_array, 1)
+        # use max because there is an edge case where multiple spokes
+        # that are nearly vertical or horizontal can cause the sum to be higher than the center
+        # allowing for a dip that can affect the initial center point
+        x_sum = np.max(central_array, 0)
+        y_sum = np.max(central_array, 1)
 
         # Calculate Full-Width, 80% Maximum center
         fwxm_x_point = (
@@ -209,7 +212,7 @@ def _get_reasonable_start_point(self) -> Point:
             round(FWXMProfile(values=y_sum, fwxm_height=80).center_idx) + top_third
         )
         center_point = Point(fwxm_x_point, fwxm_y_point)
-        return center_point
+        return center_point, np.percentile(central_array, 90)
 
     @argue.bounds(radius=(0.2, 0.95), min_peak_height=(0.05, 0.95))
     def analyze(
@@ -267,18 +270,20 @@ def analyze(
         """
         self.tolerance = tolerance
         self.image.check_inversion_by_histogram(percentiles=[4, 50, 96])
+        self.image.ground()
         if invert:
             self.image.invert()
 
+        auto_point, local_max = self._get_reasonable_start_point()
         if start_point is None:
-            start_point = self._get_reasonable_start_point()
+            start_point = auto_point
 
         self._get_reasonable_wobble(
-            start_point, fwhm, min_peak_height, radius, recursive
+            start_point, fwhm, min_peak_height, radius, recursive, local_max
         )
 
     def _get_reasonable_wobble(
-        self, start_point, fwhm, min_peak_height, radius, recursive
+        self, start_point, fwhm, min_peak_height, radius, recursive, local_max
     ):
         """Determine a wobble that is "reasonable". If recursive is false, the first iteration will be passed,
         otherwise the parameters will be tweaked to search for a reasonable wobble."""
@@ -288,14 +293,20 @@ def _get_reasonable_wobble(
         radius_gen = get_radius()
         while wobble_unreasonable:
             try:
+                min_height = max(min_peak_height * local_max, 1.01)
                 self.circle_profile = StarProfile(
-                    self.image, focus_point, radius, min_peak_height, fwhm
+                    self.image, focus_point, radius, min_height, fwhm
                 )
+                # must have at least 3 spokes and we must detect an even number (e.g. gantry starshot spokes can fade off)
                 if (len(self.circle_profile.peaks) < 6) or (
                     len(self.circle_profile.peaks) % 2 != 0
                 ):
                     raise ValueError
-                self.lines = LineManager(self.circle_profile.peaks)
+                self.lines = LineManager(
+                    self.circle_profile.peaks,
+                    focus_point=focus_point,
+                    dpmm=self.image.dpmm,
+                )
                 self._find_wobble_minimize()
             except ValueError:
                 if not recursive:
@@ -318,23 +329,19 @@ def _get_reasonable_wobble(
                             )
 
             else:  # if no errors are raised
-                # set the focus point to the wobble minimum
-                # focus_point = self.wobble.center
-                # finally:
                 # stop after first iteration if not recursive
                 if not recursive:
                     wobble_unreasonable = False
                 # otherwise, check if the wobble is reasonable
                 else:
-                    # if so, stop
-                    if self.wobble.diameter_mm < 2:
-                        focus_near_center = (
-                            self.wobble.center.distance_to(focus_point) < 5
-                        )
-                        if focus_near_center:
-                            wobble_unreasonable = False
-                        else:
-                            focus_point = self.wobble.center
+                    # see if the wobble is reasonable and near the center
+                    # if so, we can assume we found the wobble
+                    focus_near_center = (
+                        self.wobble.center.distance_to(focus_point)
+                        < 10 * self.image.dpmm
+                    )
+                    if self.wobble.diameter_mm < 2 and focus_near_center:
+                        wobble_unreasonable = False
                     # otherwise, iterate through peak height
                     else:
                         try:
@@ -669,14 +676,16 @@ def diameter_mm(self) -> float:
 class LineManager:
     """Manages the radiation lines found."""
 
-    def __init__(self, points: list[Point]):
+    def __init__(self, points: list[Point], focus_point: Point, dpmm: float):
         """
         Parameters
         ----------
         points :
             The peak points found by the StarProfile
         """
         self.lines = []
+        self.focus_point = focus_point
+        self.dpmm = dpmm
         self.construct_rad_lines(points)
 
     def __getitem__(self, item):
@@ -704,6 +713,14 @@ def construct_rad_lines(self, points: list[Point]):
         geometry.Line : returning object
         """
         self.match_points(points)
+        # check that the lines pass near the focus point. If not, we can
+        # still be in a situation where we're missing spoke-halves like a gantry starshot.
+        for line in self.lines:
+            if line.distance_to(self.focus_point) > 10 * self.dpmm:
+                raise ValueError(
+                    "The radiation lines are not near the center of the image. "
+                    "This could be due to missing spoke halves, such as in a gantry starshot."
+                )
 
     def match_points(self, points: list[Point]):
         """Match the peaks found to the same radiation lines.

tests_basic/test_starshot.py

@@ -101,6 +101,7 @@ class StarMixin(CloudFileMixin):
     min_peak_height = 0.25
     radius = 0.85
     test_all_radii = True
+    radii_range = np.linspace(0.9, 0.25, 8)
     fwxm = True
     wobble_tolerance = 0.2
     kwargs = {"sid": 1000}
@@ -167,24 +168,25 @@ def test_all_radii_give_same_wobble(self):
         """Test that the wobble stays roughly the same for all radii."""
         if self.test_all_radii:
             star = self.construct_star()
-            radii = []
-            for radius in np.linspace(0.9, 0.25, 8):
+            diameters = []
+            for radius in self.radii_range:
                 star.analyze(
                     radius=float(radius),
                     min_peak_height=self.min_peak_height,
                     recursive=self.recursive,
                     fwhm=self.fwxm,
                 )
+                diameters.append(star.wobble.diameter_mm)
+            if self.verbose:
+                print(
+                    f"Diameter mean: {np.mean(diameters):2.2f}, range: {np.max(diameters) - np.min(diameters):2.2f}"
+                )
+            for diameter in diameters:
                 self.assertAlmostEqual(
-                    star.wobble.diameter_mm,
+                    diameter,
                     self.wobble_diameter_mm,
                     delta=self.wobble_tolerance,
                 )
-                radii.append(star.wobble.diameter_mm)
-            if self.verbose:
-                print(
-                    f"Radii mean: {np.mean(radii):2.2f}, range: {np.max(radii) - np.min(radii):2.2f}"
-                )
 
 
 class Demo(StarMixin, TestCase):
@@ -263,7 +265,8 @@ def test_fails_with_tight_tol(self):
         self.assertFalse(star.passed)
 
     def test_bad_inputs_still_recovers(self):
-        self.star.analyze(radius=0.3, min_peak_height=0.1)
+        star = Starshot.from_demo_image()
+        star.analyze(radius=0.3, min_peak_height=0.1)
         self.test_wobble_center()
         self.test_wobble_diameter()
 
@@ -275,13 +278,6 @@ def test_image_inverted(self):
         top_left_corner_val_after = star.image.array[0, 0]
         self.assertNotEqual(top_left_corner_val_before, top_left_corner_val_after)
 
-    def test_bad_start_point_recovers(self):
-        """Test that even at a distance start point, the search algorithm recovers."""
-        self.star.analyze(start_point=(1000, 1000))
-        self.test_passed()
-        self.test_wobble_center()
-        self.test_wobble_diameter()
-
     def test_publish_pdf(self):
         with tempfile.TemporaryFile() as t:
             self.star.publish_pdf(t, notes="stuff", metadata={"Unit": "TB1"})
@@ -355,10 +351,14 @@ class Starshot5(StarMixin, TestCase):
 
 
 class Starshot6(StarMixin, TestCase):
+    # for the radii comparison, the wobble at 0.25 is very high due to a bad spoke
+    # detection. Setting FWHM to false will fix this. We thus clip the lower radius
+    # to 0.3 instead of 0.25.
     file_name = "Starshot#6.tif"
     wobble_center = Point(528, 607)
     wobble_diameter_mm = 0.3
     num_rad_lines = 7
+    radii_range = np.linspace(0.9, 0.3, 8)
 
 
 class Starshot7(StarMixin, TestCase):
@@ -472,7 +472,7 @@ class Starshot21(StarMixin, TestCase):
 class Starshot22(StarMixin, TestCase):
     file_name = "Starshot#22.tiff"
     wobble_center = Point(1305, 1513)
-    wobble_diameter_mm = 0.9
+    wobble_diameter_mm = 0.95
     num_rad_lines = 9
     # outside 0.93mm
 
@@ -524,3 +524,12 @@ class ChicagoSet(StarMixin, TestCase):
     def get_filename(cls):
         """Return the canonical path to the file."""
         return get_folder_from_cloud_test_repo([*cls.dir_path, cls.file_name])
+
+
+class MarkerDots(StarMixin, TestCase):
+    file_name = "marker_dots.tif"
+    wobble_center = Point(566, 559)
+    wobble_diameter_mm = 1.7
+    wobble_tolerance = 0.25
+    num_rad_lines = 3
+    passes = False