Merged in bugfix/RAM-3200_catphan_604_slice_thickness (pull request #327

)

RAM-3200 Add refinement for 604 analyses

Approved-by: Randy Taylor

docs/source/changelog.rst

@@ -24,6 +24,50 @@ Nuclear
 * A new module has been created. This module is a Python implementation of the NMQC toolkit for SPECT.
   It contains 9 tests that are very similar to the ImageJ toolkit. See :ref:`nuclear` for more.
 
+CT
+^^
+
+* Analysis of Catphan 604 datasets often did not find the HU module center correctly. This had to do with some of the
+  HU plugs being longer than the rest of the features in the 604 model. This was not causing issues and was left as-is
+  for quite some time. However, several RadMachine customers had noticed the slice thickness may be different because of this.
+  The algorithm has been adjusted to find the center
+  of the HU plugs more accurately by performing a second pass over the center slices using the relative angle between the wire ramps. This only affects the Catphan 604.
+  Users may notice a small change in HU values since the slice may now be different by 1-3 slices. Users may also notice
+  a change in the slice thickness value. All test dataset results either stayed the same or were closer to the nominal value.
+  Contrast values may also change slightly. Each of the modules are now almost always centered on the top bright marker
+  above the module.
+
+  .. figure:: images/604_old.png
+     :width: 600
+     :align: center
+
+     The old algorithm. Note the wire ramp is on the left side of the ROI for the top position. This indicates we are
+     not at the center of the HU module. Also note the side view line is barely off-center to the left for the HU module.
+
+  .. figure:: images/604_new.png
+     :width: 600
+     :align: center
+
+     The new algorithm. Note the wire ramp is now in the center of the ROI for the top position. This indicates we are
+     at the center of the HU module.
+
+* Due to the above change, a new method is available to override if desired: ``refine_origin_slice()``. This method
+  will perform the second pass over the center slices to find the HU module center. This method is available for all
+  Catphan analyses and will be empty for all phantoms besides the 604 for the time being.
+
+  If the old behavior is desired, the ``refine_origin_slice()`` method can be overridden to simply pass the
+  initial slice number:
+
+  .. code-block:: python
+
+    from pylinac import CatPhan604
+
+
+    class MyCatPhan604(CatPhan604):
+        def refine_origin_slice(self, initial_slice_num):
+            return initial_slice_num
+
+
 Profiles
 ^^^^^^^^
 

pylinac/ct.py

@@ -1770,6 +1770,9 @@ def localize(self) -> None:
         self._phantom_center_func = self.find_phantom_axis()
         self.origin_slice = self.find_origin_slice()
         self.catphan_roll = self.find_phantom_roll()
+        self.origin_slice = self.refine_origin_slice(
+            initial_slice_num=self.origin_slice
+        )
         # now that we have the origin slice, ensure we have scanned all linked modules
         if not self._ensure_physical_scan_extent():
             raise ValueError(
@@ -1905,9 +1908,13 @@ def find_origin_slice(self) -> int:
         hu_slices = hu_slices[((c + ln / 2) >= hu_slices) & (hu_slices >= (c - ln / 2))]
         center_hu_slice = int(round(float(np.median(hu_slices))))
         if self._is_within_image_extent(center_hu_slice):
-            # print(center_hu_slice)
             return center_hu_slice
 
+    def refine_origin_slice(self, initial_slice_num: int) -> int:
+        """Apply a refinement to the origin slice. This was added to handle
+        the catphan 604 at least due to variations in the length of the HU plugs."""
+        return initial_slice_num
+
     def _is_right_area(self, region: RegionProperties):
         thresh = np.pi * ((self.air_bubble_radius_mm / self.mm_per_pixel) ** 2)
         return thresh * 2 > region.filled_area > thresh / 2
@@ -2370,7 +2377,7 @@ class CatPhan604(CatPhanBase):
     modules = {
         CTP404CP604: {"offset": 0},
         CTP486: {"offset": -80},
-        CTP528CP604: {"offset": 42},
+        CTP528CP604: {"offset": 40},
         CTP515: {"offset": -40},
     }
 
@@ -2382,6 +2389,82 @@ def run_demo(show: bool = True):
         print(cbct.results())
         cbct.plot_analyzed_image(show)
 
+    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.
+
+        Starting with the initial slice, go +/- 5 slices to find the slice with the least angle
+        between the left and right wires.
+
+        This suffers from a weakness that the roll is not yet determined.
+        This will thus return the slice that has the least ABSOLUTE
+        roll. If the phantom has an inherent roll, this will not be detected and may be off by a slice or so.
+        Given the angle of the wire, the error would be small and likely only 1-2 slices max.
+        """
+        angles = []
+        # make a CTP module for the purpose of easily extracting the thickness ROIs
+        ctp404, offset = self._get_module(CTP404CP604, raise_empty=True)
+        # we don't want to set up our other ROIs (they sometimes fail) so we temporarily override the method
+        original_setup = ctp404._setup_rois
+        ctp404._setup_rois = lambda x: x
+        ctp = ctp404(
+            self,
+            offset=offset,
+            clear_borders=self.clear_borders,
+            hu_tolerance=0,
+            scaling_tolerance=0,
+            thickness_tolerance=0,
+        )
+        # we have to reset the method after we're done for future calls
+        ctp404._setup_rois = original_setup
+        for slice_num in range(initial_slice_num - 5, initial_slice_num + 5):
+            slice = Slice(self, slice_num, clear_borders=self.clear_borders)
+            # make a slice and add the wire ROIs to it.
+            troi = {}
+            for name, setting in ctp.thickness_roi_settings.items():
+                troi[name] = ThicknessROI(
+                    slice.image,
+                    setting["width_pixels"],
+                    setting["height_pixels"],
+                    setting["angle_corrected"],
+                    setting["distance_pixels"],
+                    slice.phan_center,
+                )
+            # now find the angle between the left and right and top and bottom wires via the long profile
+            left_wire = troi["Left"].long_profile.center_idx
+            right_wire = troi["Right"].long_profile.center_idx
+            h_angle = abs(left_wire - right_wire)
+            top_wire = troi["Top"].long_profile.center_idx
+            bottom_wire = troi["Bottom"].long_profile.center_idx
+            v_angle = abs(top_wire - bottom_wire)
+            angle = (h_angle + v_angle) / 2
+
+            angles.append(
+                {
+                    "slice": slice_num,
+                    "angle": angle,
+                    "left width": troi["Left"].long_profile.field_width_px,
+                    "right width": troi["Right"].long_profile.field_width_px,
+                }
+            )
+
+        # 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
+        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
+        for angle_set in angles.copy():
+            if (
+                angle_set["left width"] < median_width * 0.7
+                or angle_set["right width"] < median_width * 0.7
+            ):
+                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])
+        refined_slice = angles[m_slice_num[0]]["slice"]
+        return refined_slice
+
 
 class CatPhan600(CatPhanBase):
     """A class for loading and analyzing CT DICOM files of a CatPhan 600.

tests_basic/test_cbct.py

@@ -554,7 +554,7 @@ def test_vial_roi(self):
 class CatPhan604Test(CatPhanMixin, TestCase):
     catphan = CatPhan604
     file_name = "CBCTCatPhan604.zip"
-    origin_slice = 45
+    origin_slice = 47
     hu_values = {
         "Poly": -47,
         "Acrylic": 105,
@@ -1444,7 +1444,7 @@ class CatPhan604Sen(CatPhan604Mixin, TestCase):
 
 class CatPhan604Som(CatPhan604Mixin, TestCase):
     file_name = "SiemensSomCatPhan604.zip"
-    origin_slice = 175
+    origin_slice = 179
     hu_values = {
         "Poly": -34,
         "Acrylic": 120,
@@ -1471,9 +1471,9 @@ class CatPhan604wJig(CatPhan604Mixin, TestCase):
         "Acrylic": 130,
         "Delrin": 371,
         "Air": -999,
-        "Teflon": 967,
+        "Teflon": 975,
         "PMP": -179,
-        "LDPE": -85,
+        "LDPE": -91,
         "50% Bone": 717,
         "20% Bone": 246,
     }
@@ -1482,6 +1482,28 @@ class CatPhan604wJig(CatPhan604Mixin, TestCase):
     lowcon_visible = 1
 
 
+class CatPhan604wJig2(CatPhan604Mixin, TestCase):
+    file_name = "Catphan604wJig2.zip"
+    expected_roll = -0.61
+    origin_slice = 49
+    slice_thickness = 1.95
+    avg_line_length = 49.9
+    hu_values = {
+        "Poly": -46,
+        "Acrylic": 120,
+        "Delrin": 361,
+        "Air": -999,
+        "Teflon": 965,
+        "PMP": -192,
+        "LDPE": -107,
+        "50% Bone": 691,
+        "20% Bone": 233,
+    }
+    unif_values = {"Center": 10, "Left": -2, "Right": 9, "Top": 7, "Bottom": 1}
+    mtf_values = {50: 0.28}
+    lowcon_visible = 2
+
+
 class CatPhan503SliceOverlap(CatPhan503Mixin, TestCase):
     """This dataset is 6mm slices with a 2mm overlap. The slice thickness
     default algorithm gives incorrect values unless the straddle is set explicitly"""