diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index ecc95fb5..026b73ed 100644 --- a/docs/source/changelog.rst +++ b/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 ^^^^ diff --git a/pylinac/ct.py b/pylinac/ct.py index 55646fdc..7b31c9f7 100644 --- a/pylinac/ct.py +++ b/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]) diff --git a/pyproject.toml b/pyproject.toml index 23becfe8..0aa06b8f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -98,6 +98,9 @@ filterwarnings = [ "ignore::UserWarning:pydicom", "ignore:FigureCanvasAgg.*:UserWarning" ] +markers = [ + 'catphan604', +] [tool.coverage.run] source = [ diff --git a/tests_basic/test_cbct.py b/tests_basic/test_cbct.py index 10003a90..7683cab9 100644 --- a/tests_basic/test_cbct.py +++ b/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