diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index aa70893d..01c193d3 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -47,6 +47,14 @@ Winston-Lutz Metrics ^^^^^^^ +* There is a new ``metrics`` module in pylinac. Existing metrics have been moved into this module. + + E.g. instead of ``from pylinac.core.metrics import SizedDiskLocator`` you would now do ``from pylinac.metrics.image import SizedDiskLocator``. + Image-based metrics are now under ``pylinac.metrics.image``. Profile-based metrics are now under ``pylinac.metrics.profile``. + Individual feature detection functions are now under ``pylinac.metrics.features``. + + For backward compatibility (even though metrics are relatively new feature), the old import locations will still work + but will raise a deprecation warning. * The documentation for metrics has been updated considerably. See :ref:`image-metrics`. * The detection algorithm for disk/field metrics has been written out; see :ref:`image_metric_algorithm`. * The ``DiskLocator`` class was renamed to ``SizedDiskLocator``. diff --git a/docs/source/topics/image_metrics.rst b/docs/source/topics/image_metrics.rst index 4af8ea56..0e290dc7 100644 --- a/docs/source/topics/image_metrics.rst +++ b/docs/source/topics/image_metrics.rst @@ -1,7 +1,7 @@ .. _image-metrics: -Custom Image Metrics -==================== +Images & 2D Metrics +=================== .. versionadded:: 3.16 @@ -21,23 +21,23 @@ Use Cases Tool Legend ----------- -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ -| Use Case | Constraint | Class | -+=========================================================+=====================================================+====================================================================+ -| Find the location of a BB in the image | The BB size and location is known approximately | :class:`~pylinac.core.metrics.SizedDiskLocator` | -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ -| Find the ROI properties of a BB in the image | The BB size and location is known approximately | :class:`~pylinac.core.metrics.SizedDiskRegion` | -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ -| Find the location of *N* BBs in the image | The BB size is known approximately | :class:`~pylinac.core.metrics.GlobalSizedDiskLocator` | -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ -| Find the location of a square field in an image | The field size is known approximately | :class:`~pylinac.core.metrics.GlobalSizedFieldLocator` | -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ -| Find the locations of *N* square fields in an image | The field size is not known | :class:`~pylinac.core.metrics.GlobalFieldLocator` | -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ -| Find the location of a circular field in an image | The field size and location are known approximately | :class:`~pylinac.core.metrics.SizedDiskLocator` (``invert=False``) | -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ -| Find the ROI properties of a circular field in an image | The field size and location are known approximately | :class:`~pylinac.core.metrics.SizedDiskRegion` (``invert=False``) | -+---------------------------------------------------------+-----------------------------------------------------+--------------------------------------------------------------------+ ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ +| Use Case | Constraint | Class | ++=========================================================+=====================================================+=====================================================================+ +| Find the location of a BB in the image | The BB size and location is known approximately | :class:`~pylinac.metrics.image.SizedDiskLocator` | ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ +| Find the ROI properties of a BB in the image | The BB size and location is known approximately | :class:`~pylinac.metrics.image.SizedDiskRegion` | ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ +| Find the location of *N* BBs in the image | The BB size is known approximately | :class:`~pylinac.metrics.image.GlobalSizedDiskLocator` | ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ +| Find the location of a square field in an image | The field size is known approximately | :class:`~pylinac.metrics.image.GlobalSizedFieldLocator` | ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ +| Find the locations of *N* square fields in an image | The field size is not known | :class:`~pylinac.metrics.image.GlobalFieldLocator` | ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ +| Find the location of a circular field in an image | The field size and location are known approximately | :class:`~pylinac.metrics.image.SizedDiskLocator` (``invert=False``) | ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ +| Find the ROI properties of a circular field in an image | The field size and location are known approximately | :class:`~pylinac.metrics.image.SizedDiskRegion` (``invert=False``) | ++---------------------------------------------------------+-----------------------------------------------------+---------------------------------------------------------------------+ Basic Usage ----------- @@ -47,7 +47,7 @@ To calculate metrics on an image, simply pass the metric(s) to the ``compute`` m .. code-block:: python from pylinac.core.image import DicomImage - from pylinac.core.metrics import DiskLocator, DiskRegion + from pylinac.metrics.image import DiskLocator, DiskRegion img = DicomImage("my_image.dcm") metric = img.compute( @@ -65,7 +65,7 @@ You may compute multiple metrics by passing a list of metrics: .. code-block:: python from pylinac.core.image import DicomImage - from pylinac.core.metrics import DiskLocator, DiskRegion + from pylinac.metrics.image import DiskLocator, DiskRegion img = DicomImage("my_image.dcm") metrics = img.compute( @@ -93,7 +93,7 @@ Metrics might have something to plot on the image. If so, the ``plot`` method of .. code-block:: python from pylinac.core.image import DicomImage - from pylinac.core.metrics import DiskLocator, DiskRegion + from pylinac.metrics.image import DiskLocator, DiskRegion img = DicomImage("my_image.dcm") metrics = img.compute( @@ -127,13 +127,13 @@ Sized Disk Locator The values provided below are in pixels. The following sections show how variants of how to use the metrics using physical units and relative to the center of the image. -Here's an example of using the :class:`~pylinac.core.metrics.SizedDiskLocator`: +Here's an example of using the :class:`~pylinac.metrics.image.SizedDiskLocator`: .. code-block:: python :caption: Search for a disk 100 pixels right and 100 pixels down from the top left of the image from pylinac.core.image import DicomImage - from pylinac.core.metrics import DiskLocator, DiskRegion + from pylinac.metrics.image import DiskLocator, DiskRegion img = DicomImage("my_image.dcm") img.compute( @@ -163,7 +163,7 @@ To perform the same Disk/BB location using mm instead of pixels: :caption: Search for a disk 30mm right and 30mm down from the top left of the image from pylinac.core.image import DicomImage - from pylinac.core.metrics import DiskLocator, DiskRegion + from pylinac.metrics.image import DiskLocator, DiskRegion img = DicomImage("my_image.dcm") img.compute( @@ -195,7 +195,7 @@ This will look for the disk/BB 30 pixels right and 30 pixels down from the cente :caption: Relative to center using pixels from pylinac.core.image import DicomImage - from pylinac.core.metrics import DiskLocator, DiskRegion + from pylinac.metrics.image import DiskLocator, DiskRegion img = DicomImage("my_image.dcm") img.compute( @@ -233,18 +233,18 @@ Sized Disk Region ^^^^^^^^^^^^^^^^^ -The :class:`~pylinac.core.metrics.SizedDiskRegion` metric is the same as the :class:`~pylinac.core.metrics.SizedDiskLocator`, but instead of returning the location, it returns a +The :class:`~pylinac.metrics.image.SizedDiskRegion` metric is the same as the :class:`~pylinac.metrics.image.SizedDiskLocator`, but instead of returning the location, it returns a `scikit-image regionprops `__ object that is the region of the disk. This allows one to then calculate things like the weighted centroid, area, etc. -It also supports the same class methods as the :class:`~pylinac.core.metrics.SizedDiskLocator` metric. +It also supports the same class methods as the :class:`~pylinac.metrics.image.SizedDiskLocator` metric. Global Sized Disk Locator ^^^^^^^^^^^^^^^^^^^^^^^^^ .. versionadded:: 3.17 -The :class:`~pylinac.core.metrics.GlobalSizedDiskLocator` metric is similar to the :class:`~pylinac.core.metrics.SizedDiskLocator` metric +The :class:`~pylinac.metrics.image.GlobalSizedDiskLocator` metric is similar to the :class:`~pylinac.metrics.image.SizedDiskLocator` 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, even if the locations are known. @@ -254,7 +254,7 @@ For example, here is an example analysis of an MPC image: .. code-block:: python from pylinac.core.image import XIM - from pylinac.core.metrics import GlobalDiskLocator + from pylinac.metrics.image import GlobalDiskLocator img = XIM("my_image.xim") bbs = img.compute( @@ -279,7 +279,7 @@ Global Sized Field Locator .. versionadded:: 3.17 -The :class:`~pylinac.core.metrics.GlobalSizedFieldLocator` metric is similar to the :class:`~pylinac.core.metrics.GlobalSizedDiskLocator` metric. +The :class:`~pylinac.metrics.image.GlobalSizedFieldLocator` metric is similar to the :class:`~pylinac.metrics.image.GlobalSizedDiskLocator` metric. 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. @@ -371,7 +371,7 @@ as well as reuse them where needed. To write a custom plugin, you must -* Inherit from the :class:`~pylinac.core.metrics.MetricBase` class +* Inherit from the :class:`~pylinac.metrics.image.MetricBase` class * Specify a ``name`` attribute. * Implement the ``calculate`` method. * (Optional) Implement the ``plot`` method if you want the metric to plot on the image. @@ -386,7 +386,7 @@ For example, let's built a simple plugin that finds and plots an "X" at the cent from pylinac.core.image_generator import AS1000Image, FilteredFieldLayer, GaussianFilterLayer from pylinac.core.image import DicomImage - from pylinac.core.metrics import MetricBase + from pylinac.metrics.image import MetricBase class ImageCenterMetric(MetricBase): name = "Image Center" @@ -567,26 +567,26 @@ Here is the plot of the final image with the BB location and threshold boundary API --- -.. autoclass:: pylinac.core.metrics.MetricBase +.. autoclass:: pylinac.metrics.image.MetricBase :inherited-members: :members: -.. autoclass:: pylinac.core.metrics.SizedDiskLocator +.. autoclass:: pylinac.metrics.image.SizedDiskLocator :inherited-members: :members: -.. autoclass:: pylinac.core.metrics.SizedDiskRegion +.. autoclass:: pylinac.metrics.image.SizedDiskRegion :inherited-members: :members: -.. autoclass:: pylinac.core.metrics.GlobalSizedDiskLocator +.. autoclass:: pylinac.metrics.image.GlobalSizedDiskLocator :inherited-members: :members: -.. autoclass:: pylinac.core.metrics.GlobalSizedFieldLocator +.. autoclass:: pylinac.metrics.image.GlobalSizedFieldLocator :inherited-members: :members: -.. autoclass:: pylinac.core.metrics.GlobalFieldLocator +.. autoclass:: pylinac.metrics.image.GlobalFieldLocator :inherited-members: :members: diff --git a/docs/source/topics/profiles.rst b/docs/source/topics/profiles.rst index e51fe343..f556ea79 100644 --- a/docs/source/topics/profiles.rst +++ b/docs/source/topics/profiles.rst @@ -1,7 +1,7 @@ .. _profiles: -Profiles -======== +Profiles & 1D Metrics +===================== Profiles, in the context of pylinac, are 1D arrays of data that contain a single radiation field. Colloquially, these are what @@ -122,7 +122,7 @@ a penumbra plugin could be written that calculates the penumbra of the profile. Several plugins are provided out of the box, and writing new plugins is straightforward. -Metrics are calculated by passing it as a list to the ``compute`` method. +Metrics are calculated by passing it as a list to the ``calculate`` method. The examples below show its usage. .. _profile_builtin_plugins: @@ -676,31 +676,31 @@ API :inherited-members: :members: -.. autoclass:: pylinac.core.profile.PenumbraRightMetric +.. autoclass:: pylinac.metrics.profile.PenumbraRightMetric :inherited-members: :members: -.. autoclass:: pylinac.core.profile.PenumbraLeftMetric +.. autoclass:: pylinac.metrics.profile.PenumbraLeftMetric :inherited-members: :members: -.. autoclass:: pylinac.core.profile.SymmetryPointDifferenceMetric +.. autoclass:: pylinac.metrics.profile.SymmetryPointDifferenceMetric :inherited-members: :members: -.. autoclass:: pylinac.core.profile.SymmetryPointDifferenceQuotientMetric +.. autoclass:: pylinac.metrics.profile.SymmetryPointDifferenceQuotientMetric :inherited-members: :members: -.. autoclass:: pylinac.core.profile.TopDistanceMetric +.. autoclass:: pylinac.metrics.profile.TopDistanceMetric :inherited-members: :members: -.. autoclass:: pylinac.core.profile.FlatnessRatioMetric +.. autoclass:: pylinac.metrics.profile.FlatnessRatioMetric :inherited-members: :members: -.. autoclass:: pylinac.core.profile.FlatnessDifferenceMetric +.. autoclass:: pylinac.metrics.profile.FlatnessDifferenceMetric :inherited-members: :members: diff --git a/docs/source/winston_lutz.rst b/docs/source/winston_lutz.rst index 4fd74862..2a7e5af8 100644 --- a/docs/source/winston_lutz.rst +++ b/docs/source/winston_lutz.rst @@ -532,10 +532,7 @@ The algorithm works like such: the threshold is associated with the open field. The pixels are converted to black & white and the center of mass of the pixels is assumed to be the field CAX. -* **Find the BB** -- The image is converted to binary based on pixel values *both* above the 50% threshold as above, - and below the upper threshold. The upper threshold is an iterative value, starting at the image maximum value, - that is lowered slightly when the BB is not found. If the binary image has a reasonably circular ROI, the BB is - considered found and the pixel-weighted center of mass of the BB is considered the BB location. +* **Find the BB** -- The algorithm for finding the BB can be found here: :ref:`image_metric_algorithm`. .. note:: Strictly speaking, the following aren't required analyses, but are explained for fullness and clarity. diff --git a/pylinac/core/contrast.py b/pylinac/core/contrast.py index e5e14e07..0b9c31f0 100644 --- a/pylinac/core/contrast.py +++ b/pylinac/core/contrast.py @@ -2,7 +2,7 @@ import numpy as np -from pylinac.core.utilities import OptionListMixin +from ..core.utilities import OptionListMixin class Contrast(OptionListMixin): diff --git a/pylinac/core/image.py b/pylinac/core/image.py index 1345e493..623ae828 100644 --- a/pylinac/core/image.py +++ b/pylinac/core/image.py @@ -31,6 +31,7 @@ from skimage.draw import disk from skimage.transform import rotate +from ..metrics.image import MetricBase from ..settings import PATH_TRUNCATION_LENGTH, get_dicom_cmap from .array_utils import ( bit_invert, @@ -49,7 +50,6 @@ retrieve_dicom_file, retrieve_filenames, ) -from .metrics import MetricBase from .profile import stretch as stretcharray from .scale import wrap360 from .utilities import decode_binary, is_close, simple_round diff --git a/pylinac/core/metrics.py b/pylinac/core/metrics.py index 897dc09f..8ed5f87d 100644 --- a/pylinac/core/metrics.py +++ b/pylinac/core/metrics.py @@ -1,853 +1,7 @@ -from __future__ import annotations +from ..metrics.features import * # noqa:F403 +from ..metrics.image import * # noqa:F403 +from ..metrics.utils import * # noqa:F403 -import math -from abc import ABC, abstractmethod -from collections.abc import Callable -from typing import TYPE_CHECKING, Any - -import matplotlib.pyplot as plt -import numpy as np -from skimage import measure, segmentation -from skimage.measure._regionprops import RegionProperties -from skimage.segmentation import find_boundaries - -from pylinac.core.array_utils import invert, stretch -from pylinac.core.geometry import Point - -if TYPE_CHECKING: - from pylinac.core.image import BaseImage - - -def is_symmetric(region: RegionProperties, *args, **kwargs) -> bool: - """Whether the binary object's dimensions are symmetric, i.e. a perfect circle. Used to find the BB.""" - ymin, xmin, ymax, xmax = region.bbox - y = abs(ymax - ymin) - x = abs(xmax - xmin) - if x > max(y * 1.05, y + 3) or x < min(y * 0.95, y - 3): - return False - return True - - -def is_near_center(region: RegionProperties, *args, **kwargs) -> bool: - """Whether the bb is <2cm from the center of the field""" - dpmm = kwargs["dpmm"] - shape = kwargs["shape"] - extent_limit_mm = 20 - bottom, left, top, right = region.bbox - bb_center_x = left + (right - left) / 2 - bb_center_y = bottom + (top - bottom) / 2 - x_lo_limit = shape[1] / 2 - dpmm * extent_limit_mm - x_hi_limit = shape[1] / 2 + dpmm * extent_limit_mm - is_bb_x_centered = x_lo_limit < bb_center_x < x_hi_limit - y_lo_limit = shape[0] / 2 - dpmm * extent_limit_mm - y_hi_limit = shape[0] / 2 + dpmm * extent_limit_mm - is_bb_y_centered = y_lo_limit < bb_center_y < y_hi_limit - return is_bb_x_centered and is_bb_y_centered - - -def is_right_size_bb(region: RegionProperties, *args, **kwargs) -> bool: - """Decide whether the ROI is roughly the size of a BB; not noise and not an artifact. Used to find the BB.""" - bb_area = region.area_filled / (kwargs["dpmm"] ** 2) - bb_size = kwargs["bb_size"] # radius in mm - tolerance = kwargs["tolerance"] # radius tolerance in mm - # A = pi * r^2 - larger_bb_area = np.pi * (bb_size + tolerance) ** 2 - smaller_bb_area = max( - (np.pi * (bb_size - tolerance) ** 2, 2) - ) # set a min of 2 to avoid a lower bound of 0 when radius<=2. Having a small lower bound is much more likely to find noise in a block. - # this is actually really important. A lower bound of 1 will catch SIGNIFICANT noise and produce erroneous results. - return smaller_bb_area < bb_area < larger_bb_area - - -def is_solid(region: RegionProperties, *args, **kwargs) -> bool: - """Whether the ROI is spiculated. We want nice, round ROIs, - and this will drop such ROIs. Generally, these spiculations are noise or a BB rod. - """ - return region.solidity > 0.9 - - -def is_round(region: RegionProperties, *args, **kwargs) -> bool: - """Decide if the ROI is circular in nature by testing the filled area vs bounding box. Used to find the BB.""" - expected_fill_ratio = np.pi / 4 # area of a circle inside a square - actual_fill_ratio = region.filled_area / region.bbox_area - return expected_fill_ratio * 1.2 > actual_fill_ratio > expected_fill_ratio * 0.8 - - -def is_right_circumference(region: RegionProperties, *args, **kwargs) -> bool: - """Test the regionprop's perimeter attr to see if it matches - that of an equivalent circle""" - upper_circumference = 2 * np.pi * (kwargs["bb_size"] + kwargs["tolerance"]) - lower_circumference = 2 * np.pi * (kwargs["bb_size"] - kwargs["tolerance"]) - actual_perimeter = region.perimeter / kwargs["dpmm"] - return upper_circumference > actual_perimeter > lower_circumference - - -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. 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 = 1.20 * 2 * ( - kwargs["field_width_mm"] + kwargs["field_tolerance_mm"] - ) + 2 * (kwargs["field_height_mm"] + kwargs["field_tolerance_mm"]) - lower_perimeter = 2 * ( - kwargs["field_width_mm"] - kwargs["field_tolerance_mm"] - ) + 2 * (kwargs["field_height_mm"] - kwargs["field_tolerance_mm"]) - return upper_perimeter > actual_perimeter > lower_perimeter - - -def is_square(region: RegionProperties, *args, **kwargs) -> bool: - """Decide if the ROI is square in nature by testing the filled area vs bounding box. Used to find the BB.""" - actual_fill_ratio = region.filled_area / region.bbox_area - return actual_fill_ratio > 0.8 - - -def is_right_area_square(region: RegionProperties, *args, **kwargs) -> bool: - """Decide if the ROI is square in nature by testing the filled area vs bounding box. Used to find the BB.""" - field_area = region.area_filled / (kwargs["dpmm"] ** 2) - low_bound_expected_area = ( - kwargs["field_width_mm"] - kwargs["field_tolerance_mm"] - ) * (kwargs["field_height_mm"] - kwargs["field_tolerance_mm"]) - high_bound_expected_area = ( - kwargs["field_width_mm"] + kwargs["field_tolerance_mm"] - ) * (kwargs["field_height_mm"] + kwargs["field_tolerance_mm"]) - return low_bound_expected_area < field_area < high_bound_expected_area - - -def deduplicate_points_and_boundaries( - original_points: list[Point], - new_points: list[Point], - min_separation_px: float, - original_boundaries: list[np.ndarray], - new_boundaries: list[np.ndarray], -) -> (list[Point], list[np.ndarray]): - """Deduplicate points that are too close together. The original points should be the - starting point. The new point's x, y, and z values are compared to the existing points. - If the new point is too close to the original point, it's dropped. If it's sufficiently - far away, it is added. Will return a new combined list of points. - - We assume the original points are already deduplicated. When used in a loop starting from an empty list - this is true.""" - combined_points = original_points - combined_boundaries = original_boundaries - for new_point, new_boundary in zip(new_points, new_boundaries): - for original_point in original_points: - if new_point.distance_to(original_point) < min_separation_px: - break - else: - combined_points.append(new_point) - combined_boundaries.append(new_boundary) - return combined_points, combined_boundaries - - -class MetricBase(ABC): - """Base class for any 2D metric. This class is abstract and should not be instantiated. - - The subclass should implement the ``calculate`` method and the ``name`` attribute. - - As a best practice, the ``image_compatibility`` attribute should be set to a list of image classes that the metric - is compatible with. Image types that are not in the list will raise an error. This allows - compatibility to be explicit. However, by default this is None and no compatibility checking is done. - """ - - unit: str = "" - image: BaseImage - image_compatibility: list[BaseImage] | None = None - name: str - - def inject_image(self, image: BaseImage): - """Inject the image into the metric.""" - if self.image_compatibility is not None and not isinstance( - image, self.image_compatibility - ): - raise TypeError(f"Image must be one of {self.image_compatibility}") - self.image = image - - def context_calculate(self) -> Any: - """Calculate the metric. This also checks the image hash to attempt to ensure no changes were made.""" - img_hash = hash(self.image.array.tobytes()) - calculation = self.calculate() - # check no modifications - if hash(self.image.array.tobytes()) != img_hash: - raise RuntimeError( - "A metric modified an image. This is not allowed as this could affect other, downstream metrics. Change" - "the calculate method to not modify the underlying image." - ) - return calculation - - @abstractmethod - def calculate(self) -> Any: - """Calculate the metric. Can return anything""" - pass - - def plot(self, axis: plt.Axes, **kwargs) -> None: - """Plot the metric""" - pass - - def additional_plots(self) -> list[plt.figure]: - """Plot additional information on a separate figure as needed. - - This should NOT show the figure. The figure will be shown - via the ``metric_plots`` method. Calling show here would - block other metrics from plotting their own separate metrics. - """ - pass - - -class GlobalSizedDiskLocator(MetricBase): - name: str - points: list[Point] - y_boundaries: list[np.ndarray] - x_boundaries: list[np.ndarray] - - def __init__( - self, - radius_mm: float, - radius_tolerance_mm: float, - detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( - is_round, - is_right_size_bb, - is_right_circumference, - ), - min_number: int = 1, - max_number: int | None = None, - min_separation_mm: float = 5, - name="Global Disk Locator", - ): - """Finds BBs globally within an image. - - Parameters - ---------- - radius_mm : float - The radius of the BB in mm. - radius_tolerance_mm : float - The tolerance of the BB radius in mm. - detection_conditions : list[callable] - A list of functions that take a regionprops object and return a boolean. - The functions should be used to determine whether the regionprops object - is a BB. - min_number : int - The minimum number of BBs to find. If not found, an error is raised. - max_number : int, None - The maximum number of BBs to find. If None, no maximum is set. - min_separation_mm : float - The minimum distance between BBs in mm. If BBs are found that are closer than this, - they are deduplicated. - name : str - The name of the metric. - """ - self.radius = radius_mm - self.radius_tolerance = radius_tolerance_mm - self.detection_conditions = detection_conditions - self.name = name - self.min_number = min_number - self.max_number = max_number or 1e3 - self.min_separation_mm = min_separation_mm - - def _calculate_sample( - self, sample: np.ndarray, top_offset: int, left_offset: int - ) -> (list[Point], list[np.ndarray], list[RegionProperties]): - """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.""" - - # The implementation difference here from the original isn't large, - # But we need to detect MULTIPLE bbs instead of just one. - bbs = [] - boundaries = [] - detected_regions = {} - # uses the same algo as original WL; this is better than a percentile method as the percentile method - # can often be thrown off at the very ends of the distribution. It's more linear and faster to use the simple - # spread of min/max. - sample = stretch(sample, min=0, max=1) - imin, imax = sample.min(), sample.max() - spread = imax - imin - step_size = ( - spread / 50 - ) # move in 1/50 increments; maximum of 50 passes per image - cutoff = ( - imin + step_size - ) # start at the min + 1 step; we know the min cutoff will be a blank, full image - while cutoff <= imax and len(bbs) < self.max_number: - try: - binary_array = sample > cutoff - labeled_arr = measure.label(binary_array, connectivity=1) - regions = measure.regionprops(labeled_arr, intensity_image=sample) - detected_regions = {i: r for i, r in enumerate(regions)} - for condition in self.detection_conditions: - to_pop = [] - for key, region in sorted( - detected_regions.items(), - key=lambda item: item[1].filled_area, - reverse=True, - ): - if not condition( - region, - dpmm=self.image.dpmm, - bb_size=self.radius, - tolerance=self.radius_tolerance, - shape=binary_array.shape, - ): - to_pop.append(key) - detected_regions = { - key: region - for key, region in detected_regions.items() - if key not in to_pop - } - if len(detected_regions) == 0: - raise ValueError - else: - points = [ - Point(region.weighted_centroid[1], region.weighted_centroid[0]) - for region in detected_regions.values() - ] - new_boundaries = [ - get_boundary( - detected_region, - top_offset=top_offset, - left_offset=left_offset, - ) - for detected_region in detected_regions.values() - ] - bbs, boundaries = deduplicate_points_and_boundaries( - original_points=bbs, - new_points=points, - min_separation_px=self.min_separation_mm * self.image.dpmm, - original_boundaries=boundaries, - new_boundaries=new_boundaries, - ) - except (IndexError, ValueError): - pass - finally: - cutoff += step_size - if len(bbs) < self.min_number: - # didn't find the number we needed - raise ValueError( - f"Couldn't find the minimum number of disks in the image. Found {len(bbs)}; required: {self.min_number}" - ) - return bbs, boundaries, list(detected_regions.values()) - - 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.""" - sample = invert(self.image.array) - self.points, boundaries, _ = self._calculate_sample( - sample, top_offset=0, left_offset=0 - ) - self.y_boundaries = [] - self.x_boundaries = [] - for boundary in boundaries: - boundary_y, boundary_x = np.nonzero(boundary) - self.y_boundaries.append(boundary_y) - self.x_boundaries.append(boundary_x) - return self.points - - 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""" - for point in self.points: - axis.plot(point.x, point.y, "o", color=color) - if show_boundaries: - for boundary_y, boundary_x in zip(self.y_boundaries, self.x_boundaries): - axis.scatter( - boundary_x, - boundary_y, - c=color, - marker="s", - alpha=alpha, - s=markersize, - ) - - -class SizedDiskRegion(GlobalSizedDiskLocator): - """A metric to find a disk/BB in an image where the BB is near an expected position and size. - This will calculate the scikit-image regionprops of the BB.""" - - x_offset: float - y_offset: float - is_from_physical: bool - is_from_center: bool - max_number = 1 - min_number = 1 - min_separation_mm = 1e4 - - def __init__( - self, - expected_position: Point | tuple[float, float], - search_window: tuple[float, float], - radius: float, - radius_tolerance: float, - detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( - is_right_size_bb, - is_round, - is_right_circumference, - is_symmetric, - is_solid, - ), - invert: bool = True, - name: str = "Disk Region", - ): - # purposely avoid super call as parent defaults to mm. We set the values ourselves. - self.expected_position = Point(expected_position) - self.radius = radius - self.radius_tolerance = radius_tolerance - self.search_window = search_window - self.detection_conditions = detection_conditions - self.name = name - self.invert = invert - self.is_from_center = False - self.is_from_physical = False - - @classmethod - def from_physical( - cls, - expected_position_mm: Point | tuple[float, float], - search_window_mm: tuple[float, float], - radius_mm: float, - radius_tolerance_mm: float, - detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( - is_right_size_bb, - is_round, - is_right_circumference, - is_symmetric, - is_solid, - ), - invert: bool = True, - name="Disk Region", - ): - """Create a DiskRegion using physical dimensions.""" - # We set a flag so we know to convert from physical sizes to pixels later. - # We don't have the image/dpmm yet so we can't do it now. - instance = cls( - expected_position=expected_position_mm, - search_window=search_window_mm, - radius=radius_mm, - radius_tolerance=radius_tolerance_mm, - detection_conditions=detection_conditions, - name=name, - invert=invert, - ) - instance.is_from_physical = True - return instance - - @classmethod - def from_center( - cls, - expected_position: Point | tuple[float, float], - search_window: tuple[float, float], - radius: float, - radius_tolerance: float, - detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( - is_right_size_bb, - is_round, - is_right_circumference, - is_symmetric, - is_solid, - ), - invert: bool = True, - name="Disk Region", - ): - """Create a DiskRegion from a center point.""" - # We set a flag so we know to convert from image edge to center later. - # We don't have the image/dpmm yet so we can't do it now - instance = cls( - expected_position=expected_position, - search_window=search_window, - radius=radius, - radius_tolerance=radius_tolerance, - detection_conditions=detection_conditions, - name=name, - invert=invert, - ) - instance.is_from_center = True - return instance - - @classmethod - def from_center_physical( - cls, - expected_position_mm: Point | tuple[float, float], - search_window_mm: tuple[float, float], - radius_mm: float, - radius_tolerance_mm: float = 0.25, - detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( - is_right_size_bb, - is_round, - is_right_circumference, - is_symmetric, - is_solid, - ), - invert: bool = True, - name="Disk Region", - ): - """Create a DiskRegion using physical dimensions from the center point.""" - # We set a flag so we know to convert from physical sizes to pixels later. - # We don't have the image/dpmm yet so we can't do it now - instance = cls( - expected_position=expected_position_mm, - search_window=search_window_mm, - radius=radius_mm, - radius_tolerance=radius_tolerance_mm, - detection_conditions=detection_conditions, - name=name, - invert=invert, - ) - instance.is_from_physical = True - instance.is_from_center = True - return instance - - def calculate(self) -> RegionProperties: - """Find the scikit-image regiongprops of the BB. - - This will apply a high-pass filter to the image iteratively. - The filter starts at a very low percentile and increases until - a region is found that meets the detection conditions. - """ - if self.is_from_physical: - # convert from physical sizes to pixels - self.expected_position * self.image.dpmm - self.search_window = np.asarray(self.search_window) * self.image.dpmm - else: - # convert from pixels to physical sizes - # I know, it's weird. The functions - # for detection historically have expected - # sizes in physical dimensions - self.radius /= self.image.dpmm - self.radius_tolerance /= self.image.dpmm - if self.is_from_center: - # convert from image edge to center - self.expected_position.x += self.image.shape[1] / 2 - self.expected_position.y += self.image.shape[0] / 2 - # sample the image in the search window; need to convert to mm - left = math.floor(self.expected_position.x - self.search_window[0] / 2) - right = math.ceil(self.expected_position.x + self.search_window[0] / 2) - top = math.floor(self.expected_position.y - self.search_window[1] / 2) - bottom = math.ceil(self.expected_position.y + self.search_window[1] / 2) - sample = self.image[top:bottom, left:right] - # we might need to invert the image so that the BB pixel intensity is higher than the background - if self.invert: - sample = invert(sample) - points, boundaries, regions = self._calculate_sample( - sample, top_offset=top, left_offset=left - ) - self.x_offset = left - self.y_offset = top - y_boundary, x_boundary = np.nonzero(boundaries[0]) - self.y_boundaries = [y_boundary] - self.x_boundaries = [x_boundary] - return regions[0] - - def plot( - self, - axis: plt.Axes, - show_boundaries: bool = True, - color: str = "red", - markersize: float = 3, - alpha: float = 0.25, - ) -> None: - """Plot the BB center""" - if show_boundaries: - axis.scatter( - self.x_boundaries[0], - self.y_boundaries[0], - c=color, - marker="s", - alpha=alpha, - s=markersize, - ) - - -class SizedDiskLocator(SizedDiskRegion): - """Calculates the weighted centroid of a disk/BB as a Point in an image where the disk is near an expected position and size.""" - - point: Point - - def calculate(self) -> Point: - """Get the weighted centroid of the region prop of the BB.""" - region = super().calculate() - self.point = Point( - region.weighted_centroid[1] + self.x_offset, - region.weighted_centroid[0] + self.y_offset, - ) - return self.point - - def plot( - self, - axis: plt.Axes, - show_boundaries: bool = True, - color: str = "red", - markersize: float = 3, - alpha: float = 0.25, - ) -> None: - """Plot the BB center""" - super().plot( - axis, - show_boundaries=show_boundaries, - color=color, - markersize=markersize, - alpha=alpha, - ) - axis.plot(self.point.x, self.point.y, "o", color=color, markersize=10) - - -class GlobalSizedFieldLocator(MetricBase): - fields: list[Point] - boundaries: list[np.ndarray] - is_from_physical: bool = False - - def __init__( - self, - field_width_px: float, - field_height_px: float, - field_tolerance_px: 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, - ), - ): - """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. - """ - 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 - - @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, - ), - ): - """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. - """ - 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, - ) - instance.is_from_physical = True - return instance - - def calculate(self) -> list[Point]: - """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. - imin, imax = sample.min(), sample.max() - spread = imax - imin - step_size = ( - spread / 50 - ) # move in 1/50 increments; maximum of 50 passes per image - cutoff = imin + step_size * 5 # start at 10% height - while cutoff <= imax and len(fields) < self.max_number: - try: - binary_array = sample > cutoff - binary_array = segmentation.clear_border(binary_array, buffer_size=3) - labeled_arr = measure.label(binary_array) - regions = measure.regionprops(labeled_arr, intensity_image=sample) - conditions_met = [ - all( - condition( - region, - dpmm=self.image.dpmm, - field_width_mm=self.field_width_mm, - field_height_mm=self.field_height_mm, - field_tolerance_mm=self.field_tolerance_mm, - shape=binary_array.shape, - ) - for condition in self.detection_conditions - ) - for region in regions - ] - if not any(conditions_met): - raise ValueError - else: - fields_regions = [ - regions[idx] - for idx, value in enumerate(conditions_met) - if value - ] - points = [ - 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. - new_boundaries = [ - get_boundary(region, top_offset=0, left_offset=0) - for region in fields_regions - ] - # the separation is the minimum value + field size - fields, boundaries = deduplicate_points_and_boundaries( - original_points=fields, - new_points=points, - min_separation_px=max( - r.equivalent_diameter_area for r in fields_regions - ) - / self.image.dpmm, - original_boundaries=boundaries, - new_boundaries=new_boundaries, - ) - except (IndexError, ValueError): - pass - finally: - cutoff += step_size - if len(fields) < self.min_number: - # didn't find the number we needed - raise ValueError( - f"Couldn't find the minimum number of fields in the image. Found {len(fields)}; required: {self.min_number}" - ) - self.fields = fields - self.boundaries = boundaries - return fields - - 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, 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, - ) - - -class GlobalFieldLocator(GlobalSizedFieldLocator): - def __init__( - self, - 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, - ), - ): - """Finds fields globally within an image, irrespective of size.""" - # we override to set the width/height/tolerance to be very large - # in this case we are more likely to get noise since the size is unconstrained. - super().__init__( - field_width_px=1e4, - field_height_px=1e4, - field_tolerance_px=1e4, - min_number=min_number, - max_number=max_number, - name=name, - detection_conditions=detection_conditions, - ) - - @classmethod - def from_physical( - cls, - *args, - **kwargs, - ): - raise NotImplementedError( - "This method is not implemented for global field-finding. Use the " - "standard initializer instead." - ) - - -def get_boundary( - region: RegionProperties, top_offset: int, left_offset: int -) -> np.ndarray: - """Find the boundary of the region as the absolute position in the image. - This will calculate the outline of the region. Mostly used for plotting.""" - # we pad the boundaries to offset the ROI to the right and down to the absolute position - # on the image. - return 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] + top_offset - 1, 0), (region.bbox[1] + left_offset - 1, 0)), - mode="constant", - constant_values=0, - ) +raise DeprecationWarning( + "This module has been moved to pylinac.metrics. Please import from there in the future." +) diff --git a/pylinac/core/profile.py b/pylinac/core/profile.py index 25e765d2..4adb38ad 100644 --- a/pylinac/core/profile.py +++ b/pylinac/core/profile.py @@ -1,5 +1,4 @@ -"""Module of objects that resemble or contain a profile, i.e. a 1 or 2-D f(x) representation.""" -from __future__ import annotations +from __future__ import annotations # noqa:I001 import enum import math @@ -12,22 +11,27 @@ import matplotlib.pyplot as plt import numpy as np from matplotlib.patches import Circle as mpl_Circle -from matplotlib.patches import Rectangle from scipy import ndimage, signal from scipy.interpolate import UnivariateSpline, interp1d from scipy.ndimage import gaussian_filter1d, zoom from scipy.optimize import OptimizeWarning, minimize from scipy.stats import linregress +# noqa:I001 +from ..metrics.profile import FlatnessDifferenceMetric # noqa:F401 +from ..metrics.profile import FlatnessRatioMetric # noqa:F401 +from ..metrics.profile import PenumbraLeftMetric # noqa:F401 +from ..metrics.profile import PenumbraRightMetric # noqa:F401 +from ..metrics.profile import SymmetryAreaMetric # noqa:F401 +from ..metrics.profile import SymmetryPointDifferenceMetric # noqa:F401 +from ..metrics.profile import SymmetryPointDifferenceQuotientMetric # noqa:F401 +from ..metrics.profile import LEFT, RIGHT, ProfileMetric from . import array_utils as utils from . import validators from .geometry import Circle, Point from .hill import Hill from .utilities import convert_to_enum -RIGHT = "right" -LEFT = "left" - # for Hill fits of 2D device data the # of points can be small. # This results in optimization warnings about the variance of the fit (the variance isn't of concern for us for that particular item) warnings.simplefilter("ignore", OptimizeWarning) @@ -149,331 +153,6 @@ def stretch( return stretched_array -class ProfileMetric(ABC): - """Abstract base class for profile metrics. A profile metric is a value that can be calculated from a profile - and potentially has plot features associated with it. - Examples include penumbra, flatness, and symmetry""" - - name: str - unit: str = "" - profile: ProfileBase | PhysicalProfileMixin - - def __init__(self, color: str | None = None, linestyle: str | None = None): - self.color = color - self.linestyle = linestyle - - def inject_profile(self, profile: ProfileBase) -> None: - """Inject the profile into the metric class. - We can't do this at instantiation because we don't have - the profile yet. We also don't want to force the user - to have to save it manually as they might forget. - Finally, we want to have it around for any method we might use.""" - self.profile = profile - - def plot(self, axis: plt.Axes): - """Plot the metric on the given axis.""" - pass - - @abstractmethod - def calculate(self) -> Any: - """Calculate the metric on the given profile.""" - pass - - -class FlatnessDifferenceMetric(ProfileMetric): - """Flatness as defined by IAEA Rad Onc Handbook pg 196: https://www-pub.iaea.org/MTCD/Publications/PDF/Pub1196_web.pdf""" - - name = "Flatness (Difference)" - unit = "%" - - def __init__(self, in_field_ratio: float = 0.8, color="g", linestyle="-."): - self.in_field_ratio = in_field_ratio - super().__init__(color=color, linestyle=linestyle) - - def calculate(self) -> float: - """Calculate the flatness ratio of the profile.""" - return ( - 100 - * (self.profile.field_values().max() - self.profile.field_values().min()) - / (self.profile.field_values().max() + self.profile.field_values().min()) - ) - - def plot(self, axis: plt.Axes) -> None: - """Plot the points of largest flattness difference as well as the search bounding box.""" - data = self.profile.field_values() - left, _, width = self.profile.field_indices(in_field_ratio=self.in_field_ratio) - # plot the search bounding box - axis.add_patch( - Rectangle( - (left, np.min(data)), - width, - np.max(data) - np.min(data), - fill=False, - color=self.color, - label=self.label + " Bounding box", - ) - ) - # plot the max and min values - axis.plot( - [np.argmax(data) + left, np.argmin(data) + left], - [np.max(data), np.min(data)], - "o", - color=self.color, - label=self.name, - ) - - -class FlatnessRatioMetric(FlatnessDifferenceMetric): - """Flatness as (apparently) defined by IEC.""" - - name = "Flatness (Ratio)" - - def calculate(self) -> float: - """Calculate the flatness ratio of the profile.""" - return ( - 100 * self.profile.field_values().max() / self.profile.field_values().min() - ) - - -class SymmetryPointDifferenceMetric(ProfileMetric): - """Symmetry using the point difference method.""" - - unit = "%" - name = "Point Difference Symmetry" - - def __init__( - self, - in_field_ratio: float = 0.8, - color="magenta", - linestyle="--", - max_sym_range: float = 2, - min_sym_range: float = -2, - ): - self.in_field_ratio = in_field_ratio - self.max_sym = max_sym_range - self.min_sym = min_sym_range - super().__init__(color=color, linestyle=linestyle) - - @staticmethod - def _calc_point(lt: float, rt: float, cax: float) -> float: - return 100 * (lt - rt) / cax - - @cached_property - def symmetry_values(self) -> list[float]: - field_values = self.profile.field_values(in_field_ratio=self.in_field_ratio) - cax_value = self.profile.y_at_x(self.profile.center_idx) - return [ - self._calc_point(lt, rt, cax_value) - for lt, rt in zip(field_values, field_values[::-1]) - ] - - def calculate(self) -> float: - """Calculate the symmetry ratio of the profile.""" - max_sym_idx = np.argmax(np.abs(self.symmetry_values)) - return self.symmetry_values[max_sym_idx] - - def plot(self, axis: plt.Axes, markers: (str, str) = ("^", "v")) -> None: - idx = np.argmax(self.symmetry_values) - left_edge, right_edge, _ = self.profile.field_indices( - in_field_ratio=self.in_field_ratio - ) - # plot max sym value - max_x = self.profile.x_at_x_idx(self.profile.x_idx_at_x(left_edge) + idx) - axis.plot( - max_x, - self.profile.y_at_x(max_x), - markers[0], - color=self.color, - label=self.name, - ) - # plot min sym value - min_x = self.profile.x_at_x_idx(self.profile.x_idx_at_x(right_edge) - idx) - axis.plot( - min_x, - self.profile.y_at_x(min_x), - markers[1], - color=self.color, - ) - - # plot the symmetry on a secondary axis - sec_ax = axis.twinx() - sec_ax.set_ylabel(self.name) - - # plot the symmetry on the secondary axis - # add some vertical padding and/or use the minimum/maximum symmetry values - ylim_top = max((max(self.symmetry_values) + 0.5, self.max_sym + 0.5)) - ylim_bottom = min((min(self.symmetry_values) - 0.5, self.min_sym - 0.5)) - sec_ax.set_ylim(ylim_bottom, ylim_top) - sec_ax.plot( - self.profile.field_x_values(self.in_field_ratio), - self.symmetry_values, - color=self.color, - linestyle=self.linestyle, - ) - - -class SymmetryPointDifferenceQuotientMetric(SymmetryPointDifferenceMetric): - """Symmetry as defined by IEC.""" - - name = "Point Difference Quotient Symmetry" - - def __init__( - self, - in_field_ratio: float = 0.8, - color="magenta", - linestyle="--", - max_sym_range: float = 2, - min_sym_range: float = 0, - ): - super().__init__(in_field_ratio, color, linestyle, max_sym_range, min_sym_range) - - @staticmethod - def _calc_point(lt: float, rt: float, cax: float) -> float: - """Calculate an individual point's symmetry.""" - return 100 * max((lt / rt), (rt / lt)) - - def plot(self, axis: plt.Axes, markers: (str, str) = ("x", "x")) -> None: - super().plot(axis, markers) - - -class SymmetryAreaMetric(ProfileMetric): - """The symmetry using ratios of the areas of the left and right sides of the profile.""" - - name = "Symmetry (Area)" - - def __init__( - self, - in_field_ratio: float = 0.8, - ): - self.in_field_ratio = in_field_ratio - - def calculate(self) -> float: - """Calculate the symmetry ratio of the profile using the area of the left side vs the right side.""" - _, _, width = self.profile.field_indices(in_field_ratio=self.in_field_ratio) - area_left = np.sum( - self.profile.field_values(self.in_field_ratio)[: math.floor(width / 2) + 1] - ) - area_right = np.sum( - self.profile.field_values(self.in_field_ratio)[math.ceil(width / 2) :] - ) - return 100 * (area_left - area_right) / (area_left + area_right) - - def plot(self, axis: plt.Axes): - """Plot the symmetry by shading the left and right areas""" - field_values = self.profile.field_values(self.in_field_ratio) - x_values = self.profile.field_x_values(self.in_field_ratio) - split = math.floor(len(field_values) / 2) - left_data = field_values[: split + 1] - left_x = x_values[: split + 1] - right_data = field_values[split:] - right_x = x_values[split:] - axis.fill_between(left_x, left_data, alpha=0.2, label="Left Area") - axis.fill_between(right_x, right_data, alpha=0.2, label="Right Area") - - -class PenumbraLeftMetric(ProfileMetric): - unit = "%" - name = "Left Penumbra" - side = LEFT - - def __init__(self, lower: float = 20, upper: float = 80, color="pink", ls="-."): - self.lower = lower - self.upper = upper - super().__init__(color=color, linestyle=ls) - - def calculate(self) -> float: - """Calculate the left penumbra in mm. - We first find the edge point and then return the - distance from the lower penumbra value to upper penumbra value. - The trick is that wherever the field edge is, is assumed to be 50% - height. It's okay if it's not actually (like for FFF). - """ - left_edge = self.profile.field_edge_idx(side=self.side) - left_edge_value = self.profile.y_at_x(left_edge) - lower_search_value = left_edge_value * 2 * self.lower / 100 - lower_index = self.profile.x_at_y(y=lower_search_value, side=self.side) - upper_search_value = left_edge_value * 2 * self.upper / 100 - upper_index = self.profile.x_at_y(y=upper_search_value, side=self.side) - self.lower_index = lower_index - self.upper_index = upper_index - return abs(upper_index - lower_index) / self.profile.dpmm - - def plot(self, axis: plt.Axes): - axis.vlines( - x=[self.lower_index, self.upper_index], - ymin=self.profile.values.min(), - ymax=self.profile.values.max(), - color=self.color, - linestyle=self.linestyle, - label=self.name, - ) - - -class PenumbraRightMetric(PenumbraLeftMetric): - side = RIGHT - name = "Right Penumbra" - - -class TopDistanceMetric(ProfileMetric): - """The distance from an FFF beam's "top" to the center of the field. Similar, although - not 100% faithful to NCS-33. The NCS report uses the middle 5cm but we use a field ratio. - In practice, this shouldn't make a difference.""" - - name = "Top Distance" - unit = "mm" - - def __init__(self, top_region_ratio: float = 0.2, color="orange"): - self.top_region_ratio = top_region_ratio - super().__init__(color=color) - - def calculate(self) -> float: - """Calculate the distance from the top to the field center. Positive means the top is to the right, - negative means the top is to the left.""" - values = self.profile.field_values(in_field_ratio=self.top_region_ratio) - left, right, _ = self.profile.field_indices( - in_field_ratio=self.top_region_ratio - ) - fit_params = np.polyfit( - range(left, right + 1), - values, - deg=2, - ) - - # minimize the polynomial function - min_f = minimize( - lambda x: -np.polyval( - fit_params, x - ), # return the negative since we're MINIMIZING and want the top value - method="Nelder-Mead", - x0=self.profile.center_idx, - bounds=((left, right),), - ) - top_idx = min_f.x[0] - self.top_idx = top_idx - self.top_values = np.polyval(fit_params, range(left, right + 1)) - return (top_idx - self.profile.center_idx) / self.profile.dpmm - - def plot(self, axis: plt.Axes): - """Plot the top point and the fitted curve.""" - axis.plot( - self.top_idx, - self.profile.y_at_x(self.top_idx), - "o", - color=self.color, - label=self.name, - ) - left, right, _ = self.profile.field_indices( - in_field_ratio=self.top_region_ratio - ) - axis.plot( - range(left, right + 1), - self.top_values, - color=self.color, - linestyle=self.linestyle, - label=self.name + " Fit", - ) - - class ProfileMixin: """A mixin to provide various manipulations of 1D profile data.""" diff --git a/pylinac/metrics/__init__.py b/pylinac/metrics/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/pylinac/metrics/features.py b/pylinac/metrics/features.py new file mode 100644 index 00000000..cfd445c7 --- /dev/null +++ b/pylinac/metrics/features.py @@ -0,0 +1,101 @@ +from __future__ import annotations + +import numpy as np +from skimage.measure._regionprops import RegionProperties + + +def is_symmetric(region: RegionProperties, *args, **kwargs) -> bool: + """Whether the binary object's dimensions are symmetric, i.e. a perfect circle. Used to find the BB.""" + ymin, xmin, ymax, xmax = region.bbox + y = abs(ymax - ymin) + x = abs(xmax - xmin) + if x > max(y * 1.05, y + 3) or x < min(y * 0.95, y - 3): + return False + return True + + +def is_near_center(region: RegionProperties, *args, **kwargs) -> bool: + """Whether the bb is <2cm from the center of the field""" + dpmm = kwargs["dpmm"] + shape = kwargs["shape"] + extent_limit_mm = 20 + bottom, left, top, right = region.bbox + bb_center_x = left + (right - left) / 2 + bb_center_y = bottom + (top - bottom) / 2 + x_lo_limit = shape[1] / 2 - dpmm * extent_limit_mm + x_hi_limit = shape[1] / 2 + dpmm * extent_limit_mm + is_bb_x_centered = x_lo_limit < bb_center_x < x_hi_limit + y_lo_limit = shape[0] / 2 - dpmm * extent_limit_mm + y_hi_limit = shape[0] / 2 + dpmm * extent_limit_mm + is_bb_y_centered = y_lo_limit < bb_center_y < y_hi_limit + return is_bb_x_centered and is_bb_y_centered + + +def is_right_size_bb(region: RegionProperties, *args, **kwargs) -> bool: + """Decide whether the ROI is roughly the size of a BB; not noise and not an artifact. Used to find the BB.""" + bb_area = region.area_filled / (kwargs["dpmm"] ** 2) + bb_size = kwargs["bb_size"] # radius in mm + tolerance = kwargs["tolerance"] # radius tolerance in mm + # A = pi * r^2 + larger_bb_area = np.pi * (bb_size + tolerance) ** 2 + smaller_bb_area = max( + (np.pi * (bb_size - tolerance) ** 2, 2) + ) # set a min of 2 to avoid a lower bound of 0 when radius<=2. Having a small lower bound is much more likely to find noise in a block. + # this is actually really important. A lower bound of 1 will catch SIGNIFICANT noise and produce erroneous results. + return smaller_bb_area < bb_area < larger_bb_area + + +def is_solid(region: RegionProperties, *args, **kwargs) -> bool: + """Whether the ROI is spiculated. We want nice, round ROIs, + and this will drop such ROIs. Generally, these spiculations are noise or a BB rod. + """ + return region.solidity > 0.9 + + +def is_round(region: RegionProperties, *args, **kwargs) -> bool: + """Decide if the ROI is circular in nature by testing the filled area vs bounding box. Used to find the BB.""" + expected_fill_ratio = np.pi / 4 # area of a circle inside a square + actual_fill_ratio = region.filled_area / region.bbox_area + return expected_fill_ratio * 1.2 > actual_fill_ratio > expected_fill_ratio * 0.8 + + +def is_right_circumference(region: RegionProperties, *args, **kwargs) -> bool: + """Test the regionprop's perimeter attr to see if it matches + that of an equivalent circle""" + upper_circumference = 2 * np.pi * (kwargs["bb_size"] + kwargs["tolerance"]) + lower_circumference = 2 * np.pi * (kwargs["bb_size"] - kwargs["tolerance"]) + actual_perimeter = region.perimeter / kwargs["dpmm"] + return upper_circumference > actual_perimeter > lower_circumference + + +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. 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 = 1.20 * 2 * ( + kwargs["field_width_mm"] + kwargs["field_tolerance_mm"] + ) + 2 * (kwargs["field_height_mm"] + kwargs["field_tolerance_mm"]) + lower_perimeter = 2 * ( + kwargs["field_width_mm"] - kwargs["field_tolerance_mm"] + ) + 2 * (kwargs["field_height_mm"] - kwargs["field_tolerance_mm"]) + return upper_perimeter > actual_perimeter > lower_perimeter + + +def is_square(region: RegionProperties, *args, **kwargs) -> bool: + """Decide if the ROI is square in nature by testing the filled area vs bounding box. Used to find the BB.""" + actual_fill_ratio = region.filled_area / region.bbox_area + return actual_fill_ratio > 0.8 + + +def is_right_area_square(region: RegionProperties, *args, **kwargs) -> bool: + """Decide if the ROI is square in nature by testing the filled area vs bounding box. Used to find the BB.""" + field_area = region.area_filled / (kwargs["dpmm"] ** 2) + low_bound_expected_area = ( + kwargs["field_width_mm"] - kwargs["field_tolerance_mm"] + ) * (kwargs["field_height_mm"] - kwargs["field_tolerance_mm"]) + high_bound_expected_area = ( + kwargs["field_width_mm"] + kwargs["field_tolerance_mm"] + ) * (kwargs["field_height_mm"] + kwargs["field_tolerance_mm"]) + return low_bound_expected_area < field_area < high_bound_expected_area diff --git a/pylinac/metrics/image.py b/pylinac/metrics/image.py new file mode 100644 index 00000000..62135f5c --- /dev/null +++ b/pylinac/metrics/image.py @@ -0,0 +1,713 @@ +from __future__ import annotations + +import math +import typing +from abc import ABC, abstractmethod +from typing import Any, Callable + +import numpy as np +from matplotlib import pyplot as plt +from skimage import measure, segmentation +from skimage.measure._regionprops import RegionProperties + +from ..core.array_utils import invert, stretch +from ..core.geometry import Point +from ..metrics.features import ( + is_right_area_square, + is_right_circumference, + is_right_size_bb, + is_right_square_perimeter, + is_round, + is_solid, + is_symmetric, +) +from ..metrics.utils import deduplicate_points_and_boundaries, get_boundary + +if typing.TYPE_CHECKING: + from ..core.image import BaseImage + + +class MetricBase(ABC): + """Base class for any 2D metric. This class is abstract and should not be instantiated. + + The subclass should implement the ``calculate`` method and the ``name`` attribute. + + As a best practice, the ``image_compatibility`` attribute should be set to a list of image classes that the metric + is compatible with. Image types that are not in the list will raise an error. This allows + compatibility to be explicit. However, by default this is None and no compatibility checking is done. + """ + + unit: str = "" + image: BaseImage + image_compatibility: list[BaseImage] | None = None + name: str + + def inject_image(self, image: BaseImage): + """Inject the image into the metric.""" + if self.image_compatibility is not None and not isinstance( + image, self.image_compatibility + ): + raise TypeError(f"Image must be one of {self.image_compatibility}") + self.image = image + + def context_calculate(self) -> Any: + """Calculate the metric. This also checks the image hash to attempt to ensure no changes were made.""" + img_hash = hash(self.image.array.tobytes()) + calculation = self.calculate() + # check no modifications + if hash(self.image.array.tobytes()) != img_hash: + raise RuntimeError( + "A metric modified an image. This is not allowed as this could affect other, downstream metrics. Change" + "the calculate method to not modify the underlying image." + ) + return calculation + + @abstractmethod + def calculate(self) -> Any: + """Calculate the metric. Can return anything""" + pass + + def plot(self, axis: plt.Axes, **kwargs) -> None: + """Plot the metric""" + pass + + def additional_plots(self) -> list[plt.figure]: + """Plot additional information on a separate figure as needed. + + This should NOT show the figure. The figure will be shown + via the ``metric_plots`` method. Calling show here would + block other metrics from plotting their own separate metrics. + """ + pass + + +class GlobalSizedDiskLocator(MetricBase): + name: str + points: list[Point] + y_boundaries: list[np.ndarray] + x_boundaries: list[np.ndarray] + + def __init__( + self, + radius_mm: float, + radius_tolerance_mm: float, + detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( + is_round, + is_right_size_bb, + is_right_circumference, + ), + min_number: int = 1, + max_number: int | None = None, + min_separation_mm: float = 5, + name="Global Disk Locator", + ): + """Finds BBs globally within an image. + + Parameters + ---------- + radius_mm : float + The radius of the BB in mm. + radius_tolerance_mm : float + The tolerance of the BB radius in mm. + detection_conditions : list[callable] + A list of functions that take a regionprops object and return a boolean. + The functions should be used to determine whether the regionprops object + is a BB. + min_number : int + The minimum number of BBs to find. If not found, an error is raised. + max_number : int, None + The maximum number of BBs to find. If None, no maximum is set. + min_separation_mm : float + The minimum distance between BBs in mm. If BBs are found that are closer than this, + they are deduplicated. + name : str + The name of the metric. + """ + self.radius = radius_mm + self.radius_tolerance = radius_tolerance_mm + self.detection_conditions = detection_conditions + self.name = name + self.min_number = min_number + self.max_number = max_number or 1e3 + self.min_separation_mm = min_separation_mm + + def _calculate_sample( + self, sample: np.ndarray, top_offset: int, left_offset: int + ) -> (list[Point], list[np.ndarray], list[RegionProperties]): + """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.""" + + # The implementation difference here from the original isn't large, + # But we need to detect MULTIPLE bbs instead of just one. + bbs = [] + boundaries = [] + detected_regions = {} + # uses the same algo as original WL; this is better than a percentile method as the percentile method + # can often be thrown off at the very ends of the distribution. It's more linear and faster to use the simple + # spread of min/max. + sample = stretch(sample, min=0, max=1) + imin, imax = sample.min(), sample.max() + spread = imax - imin + step_size = ( + spread / 50 + ) # move in 1/50 increments; maximum of 50 passes per image + cutoff = ( + imin + step_size + ) # start at the min + 1 step; we know the min cutoff will be a blank, full image + while cutoff <= imax and len(bbs) < self.max_number: + try: + binary_array = sample > cutoff + labeled_arr = measure.label(binary_array, connectivity=1) + regions = measure.regionprops(labeled_arr, intensity_image=sample) + detected_regions = {i: r for i, r in enumerate(regions)} + for condition in self.detection_conditions: + to_pop = [] + for key, region in sorted( + detected_regions.items(), + key=lambda item: item[1].filled_area, + reverse=True, + ): + if not condition( + region, + dpmm=self.image.dpmm, + bb_size=self.radius, + tolerance=self.radius_tolerance, + shape=binary_array.shape, + ): + to_pop.append(key) + detected_regions = { + key: region + for key, region in detected_regions.items() + if key not in to_pop + } + if len(detected_regions) == 0: + raise ValueError + else: + points = [ + Point(region.weighted_centroid[1], region.weighted_centroid[0]) + for region in detected_regions.values() + ] + new_boundaries = [ + get_boundary( + detected_region, + top_offset=top_offset, + left_offset=left_offset, + ) + for detected_region in detected_regions.values() + ] + bbs, boundaries = deduplicate_points_and_boundaries( + original_points=bbs, + new_points=points, + min_separation_px=self.min_separation_mm * self.image.dpmm, + original_boundaries=boundaries, + new_boundaries=new_boundaries, + ) + except (IndexError, ValueError): + pass + finally: + cutoff += step_size + if len(bbs) < self.min_number: + # didn't find the number we needed + raise ValueError( + f"Couldn't find the minimum number of disks in the image. Found {len(bbs)}; required: {self.min_number}" + ) + return bbs, boundaries, list(detected_regions.values()) + + 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.""" + sample = invert(self.image.array) + self.points, boundaries, _ = self._calculate_sample( + sample, top_offset=0, left_offset=0 + ) + self.y_boundaries = [] + self.x_boundaries = [] + for boundary in boundaries: + boundary_y, boundary_x = np.nonzero(boundary) + self.y_boundaries.append(boundary_y) + self.x_boundaries.append(boundary_x) + return self.points + + 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""" + for point in self.points: + axis.plot(point.x, point.y, "o", color=color) + if show_boundaries: + for boundary_y, boundary_x in zip(self.y_boundaries, self.x_boundaries): + axis.scatter( + boundary_x, + boundary_y, + c=color, + marker="s", + alpha=alpha, + s=markersize, + ) + + +class SizedDiskRegion(GlobalSizedDiskLocator): + """A metric to find a disk/BB in an image where the BB is near an expected position and size. + This will calculate the scikit-image regionprops of the BB.""" + + x_offset: float + y_offset: float + is_from_physical: bool + is_from_center: bool + max_number = 1 + min_number = 1 + min_separation_mm = 1e4 + + def __init__( + self, + expected_position: Point | tuple[float, float], + search_window: tuple[float, float], + radius: float, + radius_tolerance: float, + detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( + is_right_size_bb, + is_round, + is_right_circumference, + is_symmetric, + is_solid, + ), + invert: bool = True, + name: str = "Disk Region", + ): + # purposely avoid super call as parent defaults to mm. We set the values ourselves. + self.expected_position = Point(expected_position) + self.radius = radius + self.radius_tolerance = radius_tolerance + self.search_window = search_window + self.detection_conditions = detection_conditions + self.name = name + self.invert = invert + self.is_from_center = False + self.is_from_physical = False + + @classmethod + def from_physical( + cls, + expected_position_mm: Point | tuple[float, float], + search_window_mm: tuple[float, float], + radius_mm: float, + radius_tolerance_mm: float, + detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( + is_right_size_bb, + is_round, + is_right_circumference, + is_symmetric, + is_solid, + ), + invert: bool = True, + name="Disk Region", + ): + """Create a DiskRegion using physical dimensions.""" + # We set a flag so we know to convert from physical sizes to pixels later. + # We don't have the image/dpmm yet so we can't do it now. + instance = cls( + expected_position=expected_position_mm, + search_window=search_window_mm, + radius=radius_mm, + radius_tolerance=radius_tolerance_mm, + detection_conditions=detection_conditions, + name=name, + invert=invert, + ) + instance.is_from_physical = True + return instance + + @classmethod + def from_center( + cls, + expected_position: Point | tuple[float, float], + search_window: tuple[float, float], + radius: float, + radius_tolerance: float, + detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( + is_right_size_bb, + is_round, + is_right_circumference, + is_symmetric, + is_solid, + ), + invert: bool = True, + name="Disk Region", + ): + """Create a DiskRegion from a center point.""" + # We set a flag so we know to convert from image edge to center later. + # We don't have the image/dpmm yet so we can't do it now + instance = cls( + expected_position=expected_position, + search_window=search_window, + radius=radius, + radius_tolerance=radius_tolerance, + detection_conditions=detection_conditions, + name=name, + invert=invert, + ) + instance.is_from_center = True + return instance + + @classmethod + def from_center_physical( + cls, + expected_position_mm: Point | tuple[float, float], + search_window_mm: tuple[float, float], + radius_mm: float, + radius_tolerance_mm: float = 0.25, + detection_conditions: list[Callable[[RegionProperties, ...], bool]] = ( + is_right_size_bb, + is_round, + is_right_circumference, + is_symmetric, + is_solid, + ), + invert: bool = True, + name="Disk Region", + ): + """Create a DiskRegion using physical dimensions from the center point.""" + # We set a flag so we know to convert from physical sizes to pixels later. + # We don't have the image/dpmm yet so we can't do it now + instance = cls( + expected_position=expected_position_mm, + search_window=search_window_mm, + radius=radius_mm, + radius_tolerance=radius_tolerance_mm, + detection_conditions=detection_conditions, + name=name, + invert=invert, + ) + instance.is_from_physical = True + instance.is_from_center = True + return instance + + def calculate(self) -> RegionProperties: + """Find the scikit-image regiongprops of the BB. + + This will apply a high-pass filter to the image iteratively. + The filter starts at a very low percentile and increases until + a region is found that meets the detection conditions. + """ + if self.is_from_physical: + # convert from physical sizes to pixels + self.expected_position * self.image.dpmm + self.search_window = np.asarray(self.search_window) * self.image.dpmm + else: + # convert from pixels to physical sizes + # I know, it's weird. The functions + # for detection historically have expected + # sizes in physical dimensions + self.radius /= self.image.dpmm + self.radius_tolerance /= self.image.dpmm + if self.is_from_center: + # convert from image edge to center + self.expected_position.x += self.image.shape[1] / 2 + self.expected_position.y += self.image.shape[0] / 2 + # sample the image in the search window; need to convert to mm + left = math.floor(self.expected_position.x - self.search_window[0] / 2) + right = math.ceil(self.expected_position.x + self.search_window[0] / 2) + top = math.floor(self.expected_position.y - self.search_window[1] / 2) + bottom = math.ceil(self.expected_position.y + self.search_window[1] / 2) + sample = self.image[top:bottom, left:right] + # we might need to invert the image so that the BB pixel intensity is higher than the background + if self.invert: + sample = invert(sample) + points, boundaries, regions = self._calculate_sample( + sample, top_offset=top, left_offset=left + ) + self.x_offset = left + self.y_offset = top + y_boundary, x_boundary = np.nonzero(boundaries[0]) + self.y_boundaries = [y_boundary] + self.x_boundaries = [x_boundary] + return regions[0] + + def plot( + self, + axis: plt.Axes, + show_boundaries: bool = True, + color: str = "red", + markersize: float = 3, + alpha: float = 0.25, + ) -> None: + """Plot the BB center""" + if show_boundaries: + axis.scatter( + self.x_boundaries[0], + self.y_boundaries[0], + c=color, + marker="s", + alpha=alpha, + s=markersize, + ) + + +class SizedDiskLocator(SizedDiskRegion): + """Calculates the weighted centroid of a disk/BB as a Point in an image where the disk is near an expected position and size.""" + + point: Point + + def calculate(self) -> Point: + """Get the weighted centroid of the region prop of the BB.""" + region = super().calculate() + self.point = Point( + region.weighted_centroid[1] + self.x_offset, + region.weighted_centroid[0] + self.y_offset, + ) + return self.point + + def plot( + self, + axis: plt.Axes, + show_boundaries: bool = True, + color: str = "red", + markersize: float = 3, + alpha: float = 0.25, + ) -> None: + """Plot the BB center""" + super().plot( + axis, + show_boundaries=show_boundaries, + color=color, + markersize=markersize, + alpha=alpha, + ) + axis.plot(self.point.x, self.point.y, "o", color=color, markersize=10) + + +class GlobalSizedFieldLocator(MetricBase): + fields: list[Point] + boundaries: list[np.ndarray] + is_from_physical: bool = False + + def __init__( + self, + field_width_px: float, + field_height_px: float, + field_tolerance_px: 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, + ), + ): + """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. + """ + 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 + + @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, + ), + ): + """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. + """ + 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, + ) + instance.is_from_physical = True + return instance + + def calculate(self) -> list[Point]: + """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. + imin, imax = sample.min(), sample.max() + spread = imax - imin + step_size = ( + spread / 50 + ) # move in 1/50 increments; maximum of 50 passes per image + cutoff = imin + step_size * 5 # start at 10% height + while cutoff <= imax and len(fields) < self.max_number: + try: + binary_array = sample > cutoff + binary_array = segmentation.clear_border(binary_array, buffer_size=3) + labeled_arr = measure.label(binary_array) + regions = measure.regionprops(labeled_arr, intensity_image=sample) + conditions_met = [ + all( + condition( + region, + dpmm=self.image.dpmm, + field_width_mm=self.field_width_mm, + field_height_mm=self.field_height_mm, + field_tolerance_mm=self.field_tolerance_mm, + shape=binary_array.shape, + ) + for condition in self.detection_conditions + ) + for region in regions + ] + if not any(conditions_met): + raise ValueError + else: + fields_regions = [ + regions[idx] + for idx, value in enumerate(conditions_met) + if value + ] + points = [ + 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. + new_boundaries = [ + get_boundary(region, top_offset=0, left_offset=0) + for region in fields_regions + ] + # the separation is the minimum value + field size + fields, boundaries = deduplicate_points_and_boundaries( + original_points=fields, + new_points=points, + min_separation_px=max( + r.equivalent_diameter_area for r in fields_regions + ) + / self.image.dpmm, + original_boundaries=boundaries, + new_boundaries=new_boundaries, + ) + except (IndexError, ValueError): + pass + finally: + cutoff += step_size + if len(fields) < self.min_number: + # didn't find the number we needed + raise ValueError( + f"Couldn't find the minimum number of fields in the image. Found {len(fields)}; required: {self.min_number}" + ) + self.fields = fields + self.boundaries = boundaries + return fields + + 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, 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, + ) + + +class GlobalFieldLocator(GlobalSizedFieldLocator): + def __init__( + self, + 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, + ), + ): + """Finds fields globally within an image, irrespective of size.""" + # we override to set the width/height/tolerance to be very large + # in this case we are more likely to get noise since the size is unconstrained. + super().__init__( + field_width_px=1e4, + field_height_px=1e4, + field_tolerance_px=1e4, + min_number=min_number, + max_number=max_number, + name=name, + detection_conditions=detection_conditions, + ) + + @classmethod + def from_physical( + cls, + *args, + **kwargs, + ): + raise NotImplementedError( + "This method is not implemented for global field-finding. Use the " + "standard initializer instead." + ) diff --git a/pylinac/metrics/profile.py b/pylinac/metrics/profile.py new file mode 100644 index 00000000..80216292 --- /dev/null +++ b/pylinac/metrics/profile.py @@ -0,0 +1,344 @@ +from __future__ import annotations + +import math +import typing +from abc import ABC, abstractmethod +from functools import cached_property +from typing import Any + +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.patches import Rectangle +from scipy.optimize import minimize + +if typing.TYPE_CHECKING: + from ..core.profile import PhysicalProfileMixin, ProfileBase + + +RIGHT = "right" +LEFT = "left" + + +class ProfileMetric(ABC): + """Abstract base class for profile metrics. A profile metric is a value that can be calculated from a profile + and potentially has plot features associated with it. + Examples include penumbra, flatness, and symmetry""" + + name: str + unit: str = "" + profile: ProfileBase | PhysicalProfileMixin + + def __init__(self, color: str | None = None, linestyle: str | None = None): + self.color = color + self.linestyle = linestyle + + def inject_profile(self, profile: ProfileBase) -> None: + """Inject the profile into the metric class. + We can't do this at instantiation because we don't have + the profile yet. We also don't want to force the user + to have to save it manually as they might forget. + Finally, we want to have it around for any method we might use.""" + self.profile = profile + + def plot(self, axis: plt.Axes): + """Plot the metric on the given axis.""" + pass + + @abstractmethod + def calculate(self) -> Any: + """Calculate the metric on the given profile.""" + pass + + +class FlatnessDifferenceMetric(ProfileMetric): + """Flatness as defined by IAEA Rad Onc Handbook pg 196: https://www-pub.iaea.org/MTCD/Publications/PDF/Pub1196_web.pdf""" + + name = "Flatness (Difference)" + unit = "%" + + def __init__(self, in_field_ratio: float = 0.8, color="g", linestyle="-."): + self.in_field_ratio = in_field_ratio + super().__init__(color=color, linestyle=linestyle) + + def calculate(self) -> float: + """Calculate the flatness ratio of the profile.""" + return ( + 100 + * (self.profile.field_values().max() - self.profile.field_values().min()) + / (self.profile.field_values().max() + self.profile.field_values().min()) + ) + + def plot(self, axis: plt.Axes) -> None: + """Plot the points of largest flattness difference as well as the search bounding box.""" + data = self.profile.field_values() + left, _, width = self.profile.field_indices(in_field_ratio=self.in_field_ratio) + # plot the search bounding box + axis.add_patch( + Rectangle( + (left, np.min(data)), + width, + np.max(data) - np.min(data), + fill=False, + color=self.color, + label=self.label + " Bounding box", + ) + ) + # plot the max and min values + axis.plot( + [np.argmax(data) + left, np.argmin(data) + left], + [np.max(data), np.min(data)], + "o", + color=self.color, + label=self.name, + ) + + +class FlatnessRatioMetric(FlatnessDifferenceMetric): + """Flatness as (apparently) defined by IEC.""" + + name = "Flatness (Ratio)" + + def calculate(self) -> float: + """Calculate the flatness ratio of the profile.""" + return ( + 100 * self.profile.field_values().max() / self.profile.field_values().min() + ) + + +class SymmetryPointDifferenceMetric(ProfileMetric): + """Symmetry using the point difference method.""" + + unit = "%" + name = "Point Difference Symmetry" + + def __init__( + self, + in_field_ratio: float = 0.8, + color="magenta", + linestyle="--", + max_sym_range: float = 2, + min_sym_range: float = -2, + ): + self.in_field_ratio = in_field_ratio + self.max_sym = max_sym_range + self.min_sym = min_sym_range + super().__init__(color=color, linestyle=linestyle) + + @staticmethod + def _calc_point(lt: float, rt: float, cax: float) -> float: + return 100 * (lt - rt) / cax + + @cached_property + def symmetry_values(self) -> list[float]: + field_values = self.profile.field_values(in_field_ratio=self.in_field_ratio) + cax_value = self.profile.y_at_x(self.profile.center_idx) + return [ + self._calc_point(lt, rt, cax_value) + for lt, rt in zip(field_values, field_values[::-1]) + ] + + def calculate(self) -> float: + """Calculate the symmetry ratio of the profile.""" + max_sym_idx = np.argmax(np.abs(self.symmetry_values)) + return self.symmetry_values[max_sym_idx] + + def plot(self, axis: plt.Axes, markers: (str, str) = ("^", "v")) -> None: + idx = np.argmax(self.symmetry_values) + left_edge, right_edge, _ = self.profile.field_indices( + in_field_ratio=self.in_field_ratio + ) + # plot max sym value + max_x = self.profile.x_at_x_idx(self.profile.x_idx_at_x(left_edge) + idx) + axis.plot( + max_x, + self.profile.y_at_x(max_x), + markers[0], + color=self.color, + label=self.name, + ) + # plot min sym value + min_x = self.profile.x_at_x_idx(self.profile.x_idx_at_x(right_edge) - idx) + axis.plot( + min_x, + self.profile.y_at_x(min_x), + markers[1], + color=self.color, + ) + + # plot the symmetry on a secondary axis + sec_ax = axis.twinx() + sec_ax.set_ylabel(self.name) + + # plot the symmetry on the secondary axis + # add some vertical padding and/or use the minimum/maximum symmetry values + ylim_top = max((max(self.symmetry_values) + 0.5, self.max_sym + 0.5)) + ylim_bottom = min((min(self.symmetry_values) - 0.5, self.min_sym - 0.5)) + sec_ax.set_ylim(ylim_bottom, ylim_top) + sec_ax.plot( + self.profile.field_x_values(self.in_field_ratio), + self.symmetry_values, + color=self.color, + linestyle=self.linestyle, + ) + + +class SymmetryPointDifferenceQuotientMetric(SymmetryPointDifferenceMetric): + """Symmetry as defined by IEC.""" + + name = "Point Difference Quotient Symmetry" + + def __init__( + self, + in_field_ratio: float = 0.8, + color="magenta", + linestyle="--", + max_sym_range: float = 2, + min_sym_range: float = 0, + ): + super().__init__(in_field_ratio, color, linestyle, max_sym_range, min_sym_range) + + @staticmethod + def _calc_point(lt: float, rt: float, cax: float) -> float: + """Calculate an individual point's symmetry.""" + return 100 * max((lt / rt), (rt / lt)) + + def plot(self, axis: plt.Axes, markers: (str, str) = ("x", "x")) -> None: + super().plot(axis, markers) + + +class SymmetryAreaMetric(ProfileMetric): + """The symmetry using ratios of the areas of the left and right sides of the profile.""" + + name = "Symmetry (Area)" + + def __init__( + self, + in_field_ratio: float = 0.8, + ): + self.in_field_ratio = in_field_ratio + + def calculate(self) -> float: + """Calculate the symmetry ratio of the profile using the area of the left side vs the right side.""" + _, _, width = self.profile.field_indices(in_field_ratio=self.in_field_ratio) + area_left = np.sum( + self.profile.field_values(self.in_field_ratio)[: math.floor(width / 2) + 1] + ) + area_right = np.sum( + self.profile.field_values(self.in_field_ratio)[math.ceil(width / 2) :] + ) + return 100 * (area_left - area_right) / (area_left + area_right) + + def plot(self, axis: plt.Axes): + """Plot the symmetry by shading the left and right areas""" + field_values = self.profile.field_values(self.in_field_ratio) + x_values = self.profile.field_x_values(self.in_field_ratio) + split = math.floor(len(field_values) / 2) + left_data = field_values[: split + 1] + left_x = x_values[: split + 1] + right_data = field_values[split:] + right_x = x_values[split:] + axis.fill_between(left_x, left_data, alpha=0.2, label="Left Area") + axis.fill_between(right_x, right_data, alpha=0.2, label="Right Area") + + +class PenumbraLeftMetric(ProfileMetric): + unit = "%" + name = "Left Penumbra" + side = LEFT + + def __init__(self, lower: float = 20, upper: float = 80, color="pink", ls="-."): + self.lower = lower + self.upper = upper + super().__init__(color=color, linestyle=ls) + + def calculate(self) -> float: + """Calculate the left penumbra in mm. + We first find the edge point and then return the + distance from the lower penumbra value to upper penumbra value. + The trick is that wherever the field edge is, is assumed to be 50% + height. It's okay if it's not actually (like for FFF). + """ + left_edge = self.profile.field_edge_idx(side=self.side) + left_edge_value = self.profile.y_at_x(left_edge) + lower_search_value = left_edge_value * 2 * self.lower / 100 + lower_index = self.profile.x_at_y(y=lower_search_value, side=self.side) + upper_search_value = left_edge_value * 2 * self.upper / 100 + upper_index = self.profile.x_at_y(y=upper_search_value, side=self.side) + self.lower_index = lower_index + self.upper_index = upper_index + return abs(upper_index - lower_index) / self.profile.dpmm + + def plot(self, axis: plt.Axes): + axis.vlines( + x=[self.lower_index, self.upper_index], + ymin=self.profile.values.min(), + ymax=self.profile.values.max(), + color=self.color, + linestyle=self.linestyle, + label=self.name, + ) + + +class PenumbraRightMetric(PenumbraLeftMetric): + side = RIGHT + name = "Right Penumbra" + + +class TopDistanceMetric(ProfileMetric): + """The distance from an FFF beam's "top" to the center of the field. Similar, although + not 100% faithful to NCS-33. The NCS report uses the middle 5cm but we use a field ratio. + In practice, this shouldn't make a difference.""" + + name = "Top Distance" + unit = "mm" + + def __init__(self, top_region_ratio: float = 0.2, color="orange"): + self.top_region_ratio = top_region_ratio + super().__init__(color=color) + + def calculate(self) -> float: + """Calculate the distance from the top to the field center. Positive means the top is to the right, + negative means the top is to the left.""" + values = self.profile.field_values(in_field_ratio=self.top_region_ratio) + left, right, _ = self.profile.field_indices( + in_field_ratio=self.top_region_ratio + ) + fit_params = np.polyfit( + range(left, right + 1), + values, + deg=2, + ) + + # minimize the polynomial function + min_f = minimize( + lambda x: -np.polyval( + fit_params, x + ), # return the negative since we're MINIMIZING and want the top value + method="Nelder-Mead", + x0=self.profile.center_idx, + bounds=((left, right),), + ) + top_idx = min_f.x[0] + self.top_idx = top_idx + self.top_values = np.polyval(fit_params, range(left, right + 1)) + return (top_idx - self.profile.center_idx) / self.profile.dpmm + + def plot(self, axis: plt.Axes): + """Plot the top point and the fitted curve.""" + axis.plot( + self.top_idx, + self.profile.y_at_x(self.top_idx), + "o", + color=self.color, + label=self.name, + ) + left, right, _ = self.profile.field_indices( + in_field_ratio=self.top_region_ratio + ) + axis.plot( + range(left, right + 1), + self.top_values, + color=self.color, + linestyle=self.linestyle, + label=self.name + " Fit", + ) diff --git a/pylinac/metrics/utils.py b/pylinac/metrics/utils.py new file mode 100644 index 00000000..2e6f36ab --- /dev/null +++ b/pylinac/metrics/utils.py @@ -0,0 +1,59 @@ +from __future__ import annotations + +import numpy as np +from skimage.measure._regionprops import RegionProperties +from skimage.segmentation import find_boundaries + +from ..core.geometry import Point + + +def deduplicate_points_and_boundaries( + original_points: list[Point], + new_points: list[Point], + min_separation_px: float, + original_boundaries: list[np.ndarray], + new_boundaries: list[np.ndarray], +) -> (list[Point], list[np.ndarray]): + """Deduplicate points that are too close together. The original points should be the + starting point. The new point's x, y, and z values are compared to the existing points. + If the new point is too close to the original point, it's dropped. If it's sufficiently + far away, it is added. Will return a new combined list of points. + + We assume the original points are already deduplicated. When used in a loop starting from an empty list + this is true.""" + combined_points = original_points + combined_boundaries = original_boundaries + for new_point, new_boundary in zip(new_points, new_boundaries): + for original_point in original_points: + if new_point.distance_to(original_point) < min_separation_px: + break + else: + combined_points.append(new_point) + combined_boundaries.append(new_boundary) + return combined_points, combined_boundaries + + +def get_boundary( + region: RegionProperties, top_offset: int, left_offset: int +) -> np.ndarray: + """Find the boundary of the region as the absolute position in the image. + This will calculate the outline of the region. Mostly used for plotting.""" + # we pad the boundaries to offset the ROI to the right and down to the absolute position + # on the image. + return 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] + top_offset - 1, 0), (region.bbox[1] + left_offset - 1, 0)), + mode="constant", + constant_values=0, + ) diff --git a/pylinac/winston_lutz.py b/pylinac/winston_lutz.py index 3a596429..66b10138 100644 --- a/pylinac/winston_lutz.py +++ b/pylinac/winston_lutz.py @@ -53,16 +53,16 @@ from .core.image import DicomImageStack, LinacDicomImage, is_image, tiff_to_dicom from .core.io import TemporaryZipDirectory, get_url, retrieve_demo_file from .core.mask import bounding_box -from .core.metrics import ( - SizedDiskLocator, +from .core.scale import MachineScale, convert +from .core.utilities import ResultBase, convert_to_enum, is_close +from .metrics.features import ( is_right_circumference, is_right_size_bb, is_round, is_solid, is_symmetric, ) -from .core.scale import MachineScale, convert -from .core.utilities import ResultBase, convert_to_enum, is_close +from .metrics.image import SizedDiskLocator BB_ERROR_MESSAGE = ( "Unable to locate the BB. Make sure the field edges do not obscure the BB, that there are no artifacts in the images, that the 'bb_size' parameter is close to reality, " diff --git a/tests_basic/core/test_profile.py b/tests_basic/core/test_profile.py index 2c5154d4..61723b53 100644 --- a/tests_basic/core/test_profile.py +++ b/tests_basic/core/test_profile.py @@ -17,8 +17,6 @@ from pylinac.core.profile import ( CircleProfile, CollapsedCircleProfile, - FlatnessDifferenceMetric, - FlatnessRatioMetric, FWXMProfile, FWXMProfilePhysical, HillProfile, @@ -28,15 +26,19 @@ Interpolation, MultiProfile, Normalization, + SingleProfile, + gamma_1d, + stretch, +) +from pylinac.metrics.profile import ( + FlatnessDifferenceMetric, + FlatnessRatioMetric, PenumbraLeftMetric, PenumbraRightMetric, - SingleProfile, SymmetryAreaMetric, SymmetryPointDifferenceMetric, SymmetryPointDifferenceQuotientMetric, TopDistanceMetric, - gamma_1d, - stretch, ) from tests_basic.utils import get_file_from_cloud_test_repo diff --git a/tests_basic/core/test_profile_metrics.py b/tests_basic/core/test_profile_metrics.py index eb3cde86..f8b8226f 100644 --- a/tests_basic/core/test_profile_metrics.py +++ b/tests_basic/core/test_profile_metrics.py @@ -13,13 +13,13 @@ GaussianFilterLayer, Layer, ) -from pylinac.core.metrics import ( +from pylinac.metrics.image import ( GlobalFieldLocator, GlobalSizedDiskLocator, GlobalSizedFieldLocator, SizedDiskLocator, - deduplicate_points_and_boundaries, ) +from pylinac.metrics.utils import deduplicate_points_and_boundaries from tests_basic.utils import get_file_from_cloud_test_repo diff --git a/tests_basic/test_winstonlutz.py b/tests_basic/test_winstonlutz.py index 84865e21..63ce602d 100644 --- a/tests_basic/test_winstonlutz.py +++ b/tests_basic/test_winstonlutz.py @@ -14,8 +14,8 @@ from pylinac.core.array_utils import create_dicom_files_from_3d_array from pylinac.core.geometry import Vector, vector_is_close from pylinac.core.io import TemporaryZipDirectory -from pylinac.core.metrics import is_round from pylinac.core.scale import MachineScale +from pylinac.metrics.features import is_round from pylinac.winston_lutz import ( Axis, WinstonLutz2D,