diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 8233f8ef..9036b93b 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -13,6 +13,8 @@ Metrics without knowing where they might be. This is relatively efficient if there are multiple BBs in the image compared with using the :class:`~pylinac.core.metrics.DiskLocator` class multiple times, even when the BB locations are known. +* The metric :class:`~pylinac.core.metrics.GlobalSizedFieldLocator`` is also available. This metric + will find a number of open fields within an image. See :ref:`global_sized_field_locator` for more. Planar Imaging ^^^^^^^^^^^^^^ diff --git a/docs/source/images/global_sized_field_locator.png b/docs/source/images/global_sized_field_locator.png new file mode 100644 index 00000000..bf04db6a Binary files /dev/null and b/docs/source/images/global_sized_field_locator.png differ diff --git a/docs/source/topics/image_metrics.rst b/docs/source/topics/image_metrics.rst index b3b7a9b8..7e35108a 100644 --- a/docs/source/topics/image_metrics.rst +++ b/docs/source/topics/image_metrics.rst @@ -221,6 +221,8 @@ This will look for the disk/BB 30mm right and 30mm down from the center of the i Global Disk Locator ^^^^^^^^^^^^^^^^^^^ +.. versionadded:: 3.17 + The :class:`~pylinac.core.metrics.GlobalDiskLocator` metric is similar to the :class:`~pylinac.core.metrics.DiskLocator` metric except that it searches the entire image for disks/BB, not just a small window. This is useful for finding the BB in images where the BB is not in the expected location or unknown. This is also efficient for finding BBs in images, @@ -249,7 +251,65 @@ This will result in an image like so: :width: 600 :align: center +.. _global_sized_field_locator: + +Global Sized Field Locator +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. versionadded:: 3.17 + +The :class:`~pylinac.core.metrics.GlobalSizedFieldLocator` metric is similar to the :class:`~pylinac.core.metrics.GlobalDiskLocator` metric +except that it searches the entire image for fields of a given size. This is useful for finding one or more fields in images +where the field is not in the expected location or unknown. This is also efficient when multiple fields are present in the image. + +The locator will find the weighted center of the field(s) and return the location(s) as a :class:`~pylinac.core.geometry.Point` objects. +The boundary of the detected field(s) will be plotted on the image in addition to the center. + +The locator will use pixels by default, but also has a ``from_physical`` class method to use physical units. + +An example plot of finding multiple fields can be seen below: + +.. image:: ../images/global_sized_field_locator.png + :width: 600 + :align: center + +For example: + +.. code-block:: python + :caption: Search for at least 2 fields of size 30x30 pixels with a tolerance of 4 pixels & plot + img = DicomImage("my_image.dcm") + img.compute( + metrics=GlobalSizedFieldLocator( + field_width_px=30, field_height_px=30, field_tolerance_px=4, max_number=2 + ) + ) + img.plot() # this will plot the image with the fields overlaid + +Using physical units +#################### + +To perform a similar field location using mm instead of pixels: + +.. code-block:: python + :caption: Search for at least 2 fields of size 30x30mm with a tolerance of 4mm + + img = DicomImage("my_image.dcm") + img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=30, field_height_mm=30, field_tolerance_mm=4, max_number=2 + ) + ) + +Usage tips +########## + +* Whenever possible, set the ``max_number`` parameter. This can **greatly** speed up the computation for several reasons. + First, it will stop searching once the number of fields is found. Second, the thresholding algorithm will have a much + better initial guess and also a better step size. This is because the approximate area of the field is known relative + to the total image size. +* The ``field_tolerance_`` parameter can be relatively tight if the ``max_number`` parameter is set. Without a + ``max_number`` parameter, you may have to increase the field tolerance to find all fields. Writing Custom Plugins ---------------------- @@ -315,3 +375,7 @@ API .. autoclass:: pylinac.core.metrics.GlobalDiskLocator :inherited-members: :members: + +.. autoclass:: pylinac.core.metrics.GlobalSizedFieldLocator + :inherited-members: + :members: diff --git a/pylinac/core/metrics.py b/pylinac/core/metrics.py index 68c6db11..95db96e6 100644 --- a/pylinac/core/metrics.py +++ b/pylinac/core/metrics.py @@ -76,9 +76,11 @@ def is_right_circumference(region: RegionProperties, *args, **kwargs) -> bool: def is_right_square_perimeter(region: RegionProperties, *args, **kwargs) -> bool: """Test the regionprop's perimeter attr to see if it matches - that of an equivalent square""" + that of an equivalent square. In reality, edges aren't perfectly straight, so + the real perimeter is always going to be higher than the theoretical perimeter. + We thus add a larger tolerance (20%) to the upper perimeter""" actual_perimeter = region.perimeter / kwargs["dpmm"] - upper_perimeter = 2 * ( + upper_perimeter = 1.20 * 2 * ( kwargs["field_width_mm"] + kwargs["field_tolerance_mm"] ) + 2 * (kwargs["field_height_mm"] + kwargs["field_tolerance_mm"]) lower_perimeter = 2 * ( @@ -496,37 +498,158 @@ def plot(self, axis: plt.Axes) -> None: class GlobalSizedFieldLocator(MetricBase): fields: list[Point] + boundaries: list[np.ndarray] + is_from_physical: bool = False def __init__( self, - field_with_mm: float, - field_height_mm: float, - field_tolerance_mm: float, + field_width_px: float, + field_height_px: float, + field_tolerance_px: float, min_number: int = 1, max_number: int | None = None, - min_separation_mm: float = 5, name: str = "Field Finder", detection_conditions: list[callable] = ( is_right_square_perimeter, is_right_area_square, ), + default_threshold_step_size: float = 2, ): - self.field_width_mm = field_with_mm - self.field_height_mm = field_height_mm - self.field_tolerance_mm = field_tolerance_mm + """Finds fields globally within an image. + + Parameters + ---------- + field_width_px : float + The width of the field in px. + field_height_px : float + The height of the field in px. + field_tolerance_px : float + The tolerance of the field size in px. + min_number : int + The minimum number of fields to find. If not found, an error is raised. + max_number : int, None + The maximum number of fields to find. If None, no maximum is set. + name : str + The name of the metric. + detection_conditions : list[callable] + A list of functions that take a regionprops object and return a boolean. + default_threshold_step_size : float + The default step size for the threshold iteration. This is based on the max number of fields and the field size. + """ + self.field_width_mm = field_width_px + self.field_height_mm = field_height_px + self.field_tolerance_mm = field_tolerance_px self.min_number = min_number self.max_number = max_number or 1e6 self.name = name self.detection_conditions = detection_conditions - self.min_separation_mm = min_separation_mm + self.default_threshold_step_size = default_threshold_step_size + + @classmethod + def from_physical( + cls, + field_width_mm: float, + field_height_mm: float, + field_tolerance_mm: float, + min_number: int = 1, + max_number: int | None = None, + name: str = "Field Finder", + detection_conditions: list[callable] = ( + is_right_square_perimeter, + is_right_area_square, + ), + default_threshold_step_size: float = 2, + ): + """Construct an instance using physical dimensions. + + Parameters + ---------- + field_width_mm : float + The width of the field in mm. + field_height_mm : float + The height of the field in mm. + field_tolerance_mm : float + The tolerance of the field size in mm. + min_number : int + The minimum number of fields to find. If not found, an error is raised. + max_number : int, None + The maximum number of fields to find. If None, no maximum is set. + name : str + The name of the metric. + detection_conditions : list[callable] + A list of functions that take a regionprops object and return a boolean. + default_threshold_step_size : float + The default step size for the threshold iteration. This is based on the max number of fields and the field size. + """ + instance = cls( + field_width_px=field_width_mm, + field_height_px=field_height_mm, + field_tolerance_px=field_tolerance_mm, + min_number=min_number, + max_number=max_number, + name=name, + detection_conditions=detection_conditions, + default_threshold_step_size=default_threshold_step_size, + ) + instance.is_from_physical = True + return instance + + @property + def threshold_step_size(self) -> float: + """Set the step size for the threshold. This is based on the max number of fields and the field size.""" + if not self.max_number: + return self.default_threshold_step_size + else: + # usually the threshold is actually very small + # since the field is very small compared to the + # image size. In this case, we want to increase + # the threshold much slower than the default. + # In combination with the threshold_start, + # this is actually quite sensitive and quick. + # In effect, we are shifting the threshold to whatever + # 10% of the expected total field area is or 2, whichever is smaller. + # For larger fields, this can be quite large, thus the 2 max. + calculated_step_size = ( + self.max_number + * (self.field_width_mm * self.field_height_mm) + * (self.image.dpmm**2) + / self.image.size + * 10 + ) + return min((calculated_step_size, self.default_threshold_step_size)) + + @property + def threshold_start(self) -> float: + """The starting percentile for the threshold. This is based on the max number of fields and the field size.""" + if not self.max_number: + return 5 + else: + # start at a higher threshold if we have a max number + # by using the expected total area of the fields / image size + # this offset from 100 and adds a 1.5 safety margin + # E.g. for a 10x10 field, this might result in a starting threshold of 99.6 + return ( + 100 + - 100 + * 1.5 + * self.max_number + * (self.field_width_mm * self.field_height_mm) + * (self.image.dpmm**2) + / self.image.size + ) def calculate(self) -> list[Point]: - """Find up to N BBs/disks in the image. This will look for BBs at every percentile range. - Multiple BBs may be found at different threshold levels.""" + """Find up to N fields in the image. This will look for fields at every percentile range. + Multiple fields may be found at different threshold levels.""" + if not self.is_from_physical: + self.field_width_mm /= self.image.dpmm + self.field_height_mm /= self.image.dpmm + self.field_tolerance_mm /= self.image.dpmm fields = [] + boundaries = [] sample = self.image.array # search for multiple BBs by iteratively raising the high-pass threshold value. - threshold_percentile = 5 + threshold_percentile = self.threshold_start while threshold_percentile < 100 and len(fields) < self.max_number: try: binary_array = sample > np.percentile(sample, threshold_percentile) @@ -558,32 +681,44 @@ def calculate(self) -> list[Point]: Point(region.centroid[1], region.centroid[0]) for region in fields_regions ] + # find the boundaries of the fields + # this is solely for plotting purposes + # these will be bool arrays + # we pad the boundaries to offset the ROI to the right + # position on the image. boundaries = [ - find_boundaries( - region.image, - connectivity=region.image.ndim, - mode="inner", - background=0, + np.pad( + find_boundaries( + # padding is needed as boundary edges aren't detected otherwise + np.pad( + region.image, + pad_width=1, + mode="constant", + constant_values=0, + ), + connectivity=region.image.ndim, + mode="inner", + background=0, + ), + ((region.bbox[0] - 1, 0), (region.bbox[1] - 1, 0)), + mode="constant", + constant_values=0, ) for region in fields_regions ] - - # print(boundaries) - # print(marked) - # the separation is the the minimum value + field size + # the separation is the minimum value + field size fields = deduplicate_points( original_points=fields, new_points=points, - min_separation_px=( - self.min_separation_mm - + min((self.field_height_mm, self.field_height_mm)) + min_separation_px=min( + (self.field_height_mm, self.field_width_mm) ) * self.image.dpmm, ) except (IndexError, ValueError): pass finally: - threshold_percentile += 2 + threshold_percentile += self.threshold_step_size if len(fields) < self.min_number: # didn't find the number we needed raise ValueError( @@ -593,33 +728,25 @@ def calculate(self) -> list[Point]: self.boundaries = boundaries return fields - def plot(self, axis: plt.Axes) -> None: - """Plot the BB centers""" + def plot( + self, + axis: plt.Axes, + show_boundaries: bool = True, + color: str = "red", + markersize: float = 3, + alpha: float = 0.25, + ) -> None: + """Plot the BB centers and boundary of detection.""" for point in self.fields: - axis.plot(point.x, point.y, "g+") - for boundary in self.boundaries: - axis.imshow(boundary) - - -# TODO -# class GlobalFieldLocator(MetricBase): -# def __init__( -# self, -# low_threshold_percentile: float = 5, -# high_threshold_percentile: float = 99.9, -# name: str = "Field Finder", -# ): -# self.low_threshold_percentile = low_threshold_percentile -# self.high_threshold_percentile = high_threshold_percentile -# self.name = name -# -# def calculate(self, image: BaseImage) -> Point: -# min, max = np.percentile( -# image.array, [self.low_threshold_percentile, self.high_threshold_percentile] -# ) -# threshold_img = image.as_binary((max - min) / 2 + min) -# filled_img = ndimage.binary_fill_holes(threshold_img) -# coords = ndimage.center_of_mass(filled_img) -# return Point(x=coords[-1], y=coords[0]) -# -# + axis.plot(point.x, point.y, color=color, marker="+", alpha=alpha) + if show_boundaries: + for boundary in self.boundaries: + boundary_y, boundary_x = np.nonzero(boundary) + axis.scatter( + boundary_x, + boundary_y, + c=color, + marker="s", + alpha=alpha, + s=markersize, + ) diff --git a/tests_basic/core/test_profile_metrics.py b/tests_basic/core/test_profile_metrics.py index c05baf96..6e8f60b0 100644 --- a/tests_basic/core/test_profile_metrics.py +++ b/tests_basic/core/test_profile_metrics.py @@ -10,7 +10,12 @@ GaussianFilterLayer, Layer, ) -from pylinac.core.metrics import DiskLocator, GlobalDiskLocator, deduplicate_points +from pylinac.core.metrics import ( + DiskLocator, + GlobalDiskLocator, + GlobalSizedFieldLocator, + deduplicate_points, +) from tests_basic.utils import get_file_from_cloud_test_repo @@ -34,6 +39,43 @@ def create_bb_image( return as1000.as_dicom() +def create_open_field_image( + field_size=(50, 50), + offset=(0, 0), + field: type[Layer] = FilteredFieldLayer, +) -> Dataset: + as1000 = AS1000Image( + sid=1000 + ) # this will set the pixel size and shape automatically + as1000.add_layer( + field(field_size_mm=field_size, cax_offset_mm=offset) + ) # create a 50x50mm square field + as1000.add_layer( + GaussianFilterLayer(sigma_mm=2) + ) # add an image-wide gaussian to simulate penumbra/scatter + as1000.add_layer(RandomNoiseLayer()) + return as1000.as_dicom() + + +def create_multi_open_field( + field_sizes=((50, 50), (50, 50)), + offsets=((0, 0), (40, 40)), + field: type[Layer] = FilteredFieldLayer, +) -> Dataset: + as1000 = AS1000Image( + sid=1000 + ) # this will set the pixel size and shape automatically + for field_size, offset in zip(field_sizes, offsets): + as1000.add_layer( + field(field_size_mm=field_size, cax_offset_mm=offset) + ) # create a 50x50mm square field + as1000.add_layer( + GaussianFilterLayer(sigma_mm=2) + ) # add an image-wide gaussian to simulate penumbra/scatter + as1000.add_layer(RandomNoiseLayer()) + return as1000.as_dicom() + + class TestGlobalDiskLocator(TestCase): @classmethod def setUpClass(cls) -> None: @@ -350,3 +392,96 @@ def test_shifted_bb(self): ) ] ) + + +class TestGlobalSizedFieldLocator(TestCase): + def test_perfect_image(self): + ds = create_open_field_image(field_size=(60, 60)) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=60, + field_height_mm=60, + field_tolerance_mm=2, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + + def test_small_image(self): + ds = create_open_field_image(field_size=(10, 10)) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=10, + field_height_mm=10, + field_tolerance_mm=2, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + + def test_large_image(self): + ds = create_open_field_image(field_size=(250, 250)) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=250, + field_height_mm=250, + field_tolerance_mm=5, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + + def test_multiple_medium_fields(self): + ds = create_multi_open_field( + field_sizes=((30, 30), (30, 30)), offsets=((0, 0), (40, 40)) + ) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=30, + field_height_mm=30, + field_tolerance_mm=5, + max_number=2, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + self.assertAlmostEqual(fields[1].x, 613.6, delta=1) + self.assertAlmostEqual(fields[1].y, 485.6, delta=1) + + def test_lots_of_fields(self): + ds = create_multi_open_field( + field_sizes=((30, 30), (30, 30), (30, 30), (30, 30)), + offsets=((0, 0), (40, 40), (-40, 0), (-60, -60)), + ) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=30, + field_height_mm=30, + field_tolerance_mm=5, + max_number=4, + ) + ) + self.assertEqual(len(fields), 4) + + def test_not_from_physical(self): + ds = create_open_field_image(field_size=(10, 10)) + img = DicomImage.from_dataset(ds) + dpmm = img.dpmm + fields = img.compute( + metrics=GlobalSizedFieldLocator( + field_width_px=10 * dpmm, + field_height_px=10 * dpmm, + field_tolerance_px=2 * dpmm, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1)