Merged in bugfix/RAM-3947_catphan_604_wire_detection (pull request #442)

Bugfix/RAM-3947 catphan 604 wire detection

Approved-by: Randy Taylor

docs/source/changelog.rst

@@ -46,6 +46,8 @@ CT
 
 * :bdg-success:`Feature` Documentation was added for how to handle the Siemens DirectDensity reconstruction algorithm which causes
   drastically different HU values. See :ref:`cbct_siemens_directdensity`.
+* :bdg-warning:`Fixed` A handful of CatPhan 604 datasets were failing to properly find the
+  "origin slice". This was caused by incorrect wire ramp detection.
 
 Core
 ^^^^

pylinac/ct.py

@@ -2718,6 +2718,15 @@ def refine_origin_slice(self, initial_slice_num: int) -> int:
         """The HU plugs are longer than the 'wire section'. This applies a refinement to find the
         slice that has the least angle between the centers of the left and right wires.
 
+        Under normal conditions, we would simply apply an offset to the
+        initial slice. The rods extend 4-5mm past the wire section.
+        Unfortunately, we also sometimes need to account for users with the
+        RM R1-4 jig which will cause localization issues due to the base of the
+        jig touching the phantom.
+        This solution is robust to the jig being present, but can suffer
+        from images where the angle of the wire, due to noise,
+        appears small but doesn't actually represent the wire ramp.
+
         Starting with the initial slice, go +/- 5 slices to find the slice with the least angle
         between the left and right wires.
 
@@ -2770,20 +2779,48 @@ def refine_origin_slice(self, initial_slice_num: int) -> int:
                     "angle": angle,
                     "left width": troi["Left"].long_profile.field_width_px,
                     "right width": troi["Right"].long_profile.field_width_px,
+                    "left profile": troi["Left"].long_profile.values,
+                    "right profile": troi["Right"].long_profile.values,
                 }
             )
 
         # some slices might not include the wire
         # we need to drop those; we do so by dropping pairs that have a field width well below the median
+        # or by fields who don't appear to have the wire in them (max value is near median)
         median_width_l = np.median([angle["left width"] for angle in angles])
         median_width_r = np.median([angle["right width"] for angle in angles])
         median_width = (median_width_l + median_width_r) / 2
+        # get median and max pixel values of all the profiles
+        median_left_pixel_val = np.median(
+            np.concatenate([a["left profile"] for a in angles])
+        )
+        median_right_pixel_val = np.median(
+            np.concatenate([a["right profile"] for a in angles])
+        )
+        median_pixel_val = (median_left_pixel_val + median_right_pixel_val) / 2
+        max_left_pixel_val = np.max(np.concatenate([a["left profile"] for a in angles]))
+        max_right_pixel_val = np.max(
+            np.concatenate([a["right profile"] for a in angles])
+        )
+        max_pixel_val = (max_left_pixel_val + max_right_pixel_val) / 2
+
         for angle_set in angles.copy():
+            # field width is well below the median; probably not in the slice; drop it
             if (
                 angle_set["left width"] < median_width * 0.7
                 or angle_set["right width"] < median_width * 0.7
             ):
                 angles.remove(angle_set)
+                continue
+            # if the max pixel value of the angle set is closer to the overall median than the max
+            # it means the wire isn't in the slice; drop it
+            max_pixel = max(
+                angle_set["left profile"].max(), angle_set["right profile"].max()
+            )
+            delta_median = abs(median_pixel_val - max_pixel)
+            delta_max = abs(max_pixel_val - max_pixel)
+            if delta_median < delta_max:
+                angles.remove(angle_set)
 
         # now find the slice with the least angle, accounting for the phantom roll
         m_slice_num = np.argsort([a["angle"] - self.catphan_roll for a in angles])

pyproject.toml

@@ -98,6 +98,9 @@ filterwarnings = [
     "ignore::UserWarning:pydicom",
     "ignore:FigureCanvasAgg.*:UserWarning"
     ]
+markers = [
+    'catphan604',
+]
 
 [tool.coverage.run]
 source = [

tests_basic/test_cbct.py

@@ -7,6 +7,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+from pytest import mark
 
 from pylinac import CatPhan503, CatPhan504, CatPhan600, CatPhan604
 from pylinac.core.geometry import Point
@@ -616,6 +617,7 @@ def test_vial_roi(self):
         self.assertIn("Vial", self.cbct.ctp404.rois)
 
 
+@mark.catphan604
 class CatPhan604Test(CatPhanMixin, TestCase):
     catphan = CatPhan604
     file_name = "CBCTCatPhan604.zip"
@@ -1488,6 +1490,7 @@ class Siemens5MM(CatPhan504Mixin, TestCase):
     lowcon_visible = 0
 
 
+@mark.catphan604
 @skip("Investigate later")
 class CatPhan604Sen(CatPhan604Mixin, TestCase):
     expected_roll = -0.5
@@ -1510,6 +1513,7 @@ class CatPhan604Sen(CatPhan604Mixin, TestCase):
     lowcon_visible = 2
 
 
+@mark.catphan604
 class CatPhan604Som(CatPhan604Mixin, TestCase):
     file_name = "SiemensSomCatPhan604.zip"
     origin_slice = 179
@@ -1530,6 +1534,7 @@ class CatPhan604Som(CatPhan604Mixin, TestCase):
     slice_thickness = 1
 
 
+@mark.catphan604
 class CatPhan604wJig(CatPhan604Mixin, TestCase):
     file_name = "Catphan604-with-jig.zip"
     origin_slice = 43
@@ -1551,6 +1556,7 @@ class CatPhan604wJig(CatPhan604Mixin, TestCase):
     avg_noise_power = 0.287
 
 
+@mark.catphan604
 class CatPhan604wJig2(CatPhan604Mixin, TestCase):
     file_name = "Catphan604wJig2.zip"
     expected_roll = -0.61
@@ -1622,6 +1628,7 @@ class CatPhan503SliceOverlap(CatPhan503Mixin, TestCase):
     mtf_values = {50: 0.40}
 
 
+@mark.catphan604
 class CatPhan604NegativeSliceOverlap(CatPhan604Mixin, TestCase):
     """Has a negative value for slice overlap. Has to do with stack order, but is irrelevant for our needs:
     https://dicom.innolitics.com/ciods/nm-image/nm-reconstruction/00180088"""
@@ -1647,6 +1654,60 @@ class CatPhan604NegativeSliceOverlap(CatPhan604Mixin, TestCase):
     lowcon_visible = 6
 
 
+@mark.catphan604
+class CatPhan604WireLocalization(CatPhan604Mixin, TestCase):
+    """Caused detection issues due to the wire not being present in some slices
+    causing noise to be detected as the wire and messing up the localization refinement
+    """
+
+    file_name = "CatPhan 604 - wire localization.zip"
+    expected_roll = 0.4
+    origin_slice = 45
+    hu_values = {
+        "Poly": -40,
+        "Acrylic": 123,
+        "Delrin": 353,
+        "Air": -1000,
+        "Teflon": 953,
+        "PMP": -189,
+        "LDPE": -101,
+        "50% Bone": 679,
+        "20% Bone": 235,
+    }
+    thickness_slice_straddle = 0
+    slice_thickness = 2.7
+    unif_values = {"Center": 0, "Left": 5, "Right": 11, "Top": 6, "Bottom": 11}
+    mtf_values = {50: 0.30}
+    lowcon_visible = 2
+
+
+@mark.catphan604
+class CatPhan604CoCr(CatPhan604Mixin, TestCase):
+    """Caused detection issues due to the wire not being present in some slices
+    causing noise to be detected as the wire and messing up the localization refinement
+    """
+
+    file_name = "CatPhan604CT - CoCr.zip"
+    expected_roll = -0.4
+    origin_slice = 62
+    hu_values = {
+        "Poly": -46,
+        "Acrylic": 129,
+        "Delrin": 398,
+        "Air": -1023,
+        "Teflon": 1075,
+        "PMP": -216,
+        "LDPE": -112,
+        "50% Bone": 740,
+        "20% Bone": 226,
+    }
+    thickness_slice_straddle = 0
+    slice_thickness = 2.35
+    unif_values = {"Center": 14, "Left": 14, "Right": 14, "Top": 14, "Bottom": 14}
+    mtf_values = {50: 0.30}
+    lowcon_visible = 2
+
+
 class CatPhan504NearEdge(CatPhan504Mixin, TestCase):
     file_name = "phantom_edge.zip"
     expected_roll = 1.4