From 4eeed3abcb8b1d242b756b883e313a7896f2b4ff Mon Sep 17 00:00:00 2001 From: Vince Reuter Date: Tue, 19 Mar 2024 03:04:19 +0100 Subject: [PATCH] first-pass at full flexibility with prelim tests (accordance with previous settings) passing --- poetry.lock | 22 +- pyproject.toml | 4 +- spotfishing/__init__.py | 2 + spotfishing/_exceptions.py | 4 + spotfishing/_numeric_types.py | 2 +- spotfishing/deprecated.py | 153 ------------- spotfishing/detection_result.py | 27 +-- spotfishing/detectors.py | 119 +++++----- spotfishing/dog_transform.py | 216 ++++++++++++++++++ .../original__difference_of_gaussians.json | 7 + tests/helpers.py | 19 ++ tests/test_accord_with_original_settings.py | 58 ++++- tests/test_old_new_equivalence.py | 45 ---- tests/test_spot_detectors.py | 11 +- 14 files changed, 399 insertions(+), 290 deletions(-) delete mode 100644 spotfishing/deprecated.py create mode 100644 spotfishing/dog_transform.py create mode 100644 spotfishing/examples/original__difference_of_gaussians.json delete mode 100644 tests/test_old_new_equivalence.py diff --git a/poetry.lock b/poetry.lock index 35bc664..018082e 100644 --- a/poetry.lock +++ b/poetry.lock @@ -508,6 +508,24 @@ files = [ {file = "numpy-1.26.4.tar.gz", hash = "sha256:2a02aba9ed12e4ac4eb3ea9421c420301a0c6460d9830d74a9df87efa4912010"}, ] +[[package]] +name = "numpydoc-decorator" +version = "0.0.0" +description = "" +optional = false +python-versions = ">=3.9,<3.13" +files = [] +develop = false + +[package.dependencies] +typing_extensions = "*" + +[package.source] +type = "git" +url = "https://github.com/alimanfoo/numpydoc_decorator.git" +reference = "main" +resolved_reference = "a9cbb0a29544788bbfd83ed47082e3442d9e5364" + [[package]] name = "packaging" version = "24.0" @@ -1162,5 +1180,5 @@ files = [ [metadata] lock-version = "2.0" -python-versions = "^3.10.0" -content-hash = "5aff3846f7e83945d2c678ae0c6f098e378ab0155a851509f0eff486448dce57" +python-versions = ">= 3.10.0, < 3.13" +content-hash = "75aeae2c872a6a04b97e25bb1cab53e44491f267e909808db0407e0cd125f4bc" diff --git a/pyproject.toml b/pyproject.toml index abebe17..354f338 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,14 +17,16 @@ classifiers = [ "Topic :: Scientific/Engineering :: Bio-Informatics", "Typing :: Typed", ] +include = ["examples"] [tool.poetry.dependencies] # These are the main runtime dependencies. -python = "^3.10.0" +python = ">= 3.10.0, < 3.13" numpy = [ {version = "^1.24.2", python = "<3.12"}, {version = "^1.25", python = ">=3.12"} ] +numpydoc_decorator = { git = "https://github.com/alimanfoo/numpydoc_decorator.git", branch = "main" } pandas = "^1.5.3" scikit-image = "^0.20.0" scipy = "^1.10.1" diff --git a/spotfishing/__init__.py b/spotfishing/__init__.py index fff9066..f02b8ae 100644 --- a/spotfishing/__init__.py +++ b/spotfishing/__init__.py @@ -4,12 +4,14 @@ from ._exceptions import * from .detection_result import DetectionResult from .detectors import detect_spots_dog, detect_spots_int +from .dog_transform import DifferenceOfGaussiansTransformation __author__ = "Vince Reuter" __credits__ = ["Vince Reuter", "Kai Sandoval Beckwith"] __all__ = [ "ROI_MEAN_INTENSITY_KEY", + "DifferenceOfGaussiansTransformation", "DimensionalityError", "DetectionResult", "detect_spots_dog", diff --git a/spotfishing/_exceptions.py b/spotfishing/_exceptions.py index 133665c..cef551f 100644 --- a/spotfishing/_exceptions.py +++ b/spotfishing/_exceptions.py @@ -3,6 +3,10 @@ __author__ = "Vince Reuter" __credits__ = ["Vince Reuter"] +__all__ = [ + "DimensionalityError", +] + class DimensionalityError(Exception): """Error subtype when dimensionality of something isn't as expected""" diff --git a/spotfishing/_numeric_types.py b/spotfishing/_numeric_types.py index f5ba6ff..411225a 100644 --- a/spotfishing/_numeric_types.py +++ b/spotfishing/_numeric_types.py @@ -9,6 +9,6 @@ __all__ = ["NumpyInt", "NumpyFloat", "PixelValue"] -NumpyInt = Union[np.int8, np.int16, np.int32, np.int64] NumpyFloat = Union[np.float16, np.float32, np.float64] +NumpyInt = Union[np.int8, np.int16, np.int32, np.int64] PixelValue = Union[np.int8, np.int16] diff --git a/spotfishing/deprecated.py b/spotfishing/deprecated.py deleted file mode 100644 index e68b005..0000000 --- a/spotfishing/deprecated.py +++ /dev/null @@ -1,153 +0,0 @@ -from typing import Union - -import numpy as np -import pandas as pd -from scipy import ndimage as ndi -from skimage.filters import gaussian -from skimage.measure import regionprops_table -from skimage.morphology import ball, remove_small_objects, white_tophat -from skimage.segmentation import expand_labels - -Numeric = Union[int, float] - -CENTROID_COLUMNS_REMAPPING = { - "centroid_weighted-0": "zc", - "centroid_weighted-1": "yc", - "centroid_weighted-2": "xc", -} - - -def detect_spots_dog_old(input_img, spot_threshold: Numeric, expand_px: int = 10): - """Spot detection by difference of Gaussians filter - - Arguments - --------- - img : ndarray - Input 3D image - spot_threshold : NumberLike - Threshold to use for spots - expand_px : int - Number of pixels by which to expand contiguous subregion, - up to point of overlap with neighboring subregion of image - - Returns - ------- - pd.DataFrame, np.ndarray, np.ndarray: - The centroids and roi_IDs of the spots found, - the image used for spot detection, and - numpy array with only sufficiently large regions - retained (bigger than threshold number of pixels), - and dilated by expansion amount (possibly) - """ - # TODO: Do not use hard-coded sigma values (second parameter to gaussian(...)). - # See: https://github.com/gerlichlab/looptrace/issues/124 - - img = white_tophat(image=input_img, footprint=ball(2)) - img = gaussian(img, 0.8) - gaussian(img, 1.3) - img = img / gaussian(input_img, 3) - img = (img - np.mean(img)) / np.std(img) - labels, _ = ndi.label(img > spot_threshold) - labels = expand_labels(labels, expand_px) - - # Make a DataFrame with the ROI info. - spot_props = _reindex_to_roi_id( - pd.DataFrame( - regionprops_table( - label_image=labels, - intensity_image=input_img, - properties=("label", "centroid_weighted", "intensity_mean"), - ) - ) - .drop(["label"], axis=1) - .rename(columns=CENTROID_COLUMNS_REMAPPING) - ) - - return spot_props, img, labels - - -def detect_spots_int_old(input_img, spot_threshold: Numeric, expand_px: int = 1): - """Spot detection by intensity filter - - Arguments - --------- - img : ndarray - Input 3D image - spot_threshold : NumberLike - Threshold to use for spots - expand_px : int - Number of pixels by which to expand contiguous subregion, - up to point of overlap with neighboring subregion of image - - Returns - ------- - pd.DataFrame, np.ndarray, np.ndarray: - The centroids and roi_IDs of the spots found, - the image used for spot detection, and - numpy array with only sufficiently large regions - retained (bigger than threshold number of pixels), - and dilated by expansion amount (possibly) - """ - # TODO: enforce that output column names don't vary with code path walked. - # See: https://github.com/gerlichlab/looptrace/issues/125 - - binary = input_img > spot_threshold - binary = ndi.binary_fill_holes(binary) - struct = ndi.generate_binary_structure(input_img.ndim, 2) - labels, n_obj = ndi.label(binary, structure=struct) - if n_obj > 1: # Do not need this with area filtering below - labels = remove_small_objects(labels, min_size=5) - if expand_px > 0: - labels = expand_labels(labels, expand_px) - if np.all(labels == 0): # No substructures (ROIs) exist after filtering. - spot_props = pd.DataFrame( - columns=[ - "label", - "z_min", - "y_min", - "x_min", - "z_max", - "y_max", - "x_max", - "area", - "zc", - "yc", - "xc", - "intensity_mean", - ] - ) - else: - spot_props = _reindex_to_roi_id( - pd.DataFrame( - regionprops_table( - labels, - input_img, - properties=( - "label", - "bbox", - "area", - "centroid_weighted", - "intensity_mean", - ), - ) - ).rename( - columns={ - **CENTROID_COLUMNS_REMAPPING, - "bbox-0": "z_min", - "bbox-1": "y_min", - "bbox-2": "x_min", - "bbox-3": "z_max", - "bbox-4": "y_max", - "bbox-5": "x_max", - } - ) - ) - - return spot_props, input_img, labels - - -def _index_as_roi_id(props_table: pd.DataFrame) -> pd.DataFrame: - return props_table.rename(columns={"index": "roi_id"}) - - -def _reindex_to_roi_id(props_table: pd.DataFrame) -> pd.DataFrame: - return _index_as_roi_id(props_table.reset_index(drop=True)) diff --git a/spotfishing/detection_result.py b/spotfishing/detection_result.py index 1e79fed..5b5042b 100644 --- a/spotfishing/detection_result.py +++ b/spotfishing/detection_result.py @@ -1,8 +1,10 @@ """Abstraction over the result of application of a spot detection procedure to an input image""" -import dataclasses +from dataclasses import dataclass from typing import TYPE_CHECKING, Iterable +from numpydoc_decorator import doc # type: ignore[import-untyped] + if TYPE_CHECKING: import numpy.typing as npt import pandas as pd @@ -39,21 +41,16 @@ ] -@dataclasses.dataclass(frozen=True) +@doc( + summary="The result of applying spot detection to an input image", + parameters=dict( + table="A table of detected spot coordinates and measurements", + image="The image (after possibly some preprocessing) that was actually used in detection", + labels="Array of the same size/shape as the image, with nonnegative integer entries; a zero indicates that the pixel isn't in an ROI, and a nonzero indicates the index of the ROI to which the pixel's been assigned", + ), +) +@dataclass(frozen=True, kw_only=True) class DetectionResult: - """ - The result of applying spot detection to an input image - - Parameters - ---------- - table : pd.DataFrame - The table of detected spot coordinates and measurements - image : np.ndarray - The image (after possibly some preprocessing) that was actually used in detection - labels : np.ndarray - Array of the same size/shape as the image, with nonnegative integer entries; a zero indicates that the pixel isn't in an ROI, and a nonzero indicates the index of the ROI to which the pixel's been assigned - """ - table: "pd.DataFrame" image: "npt.NDArray[PixelValue]" labels: "npt.NDArray[NumpyInt]" diff --git a/spotfishing/detectors.py b/spotfishing/detectors.py index e1af79a..de68436 100644 --- a/spotfishing/detectors.py +++ b/spotfishing/detectors.py @@ -1,15 +1,17 @@ """Different spot detection implementations""" +from dataclasses import dataclass from typing import Optional, Tuple, Union import numpy as np import numpy.typing as npt import pandas as pd +from numpydoc_decorator import doc # type: ignore[import-untyped] from scipy import ndimage as ndi -from skimage.filters import gaussian # type: ignore[import-untyped] from skimage.measure import regionprops_table -from skimage.morphology import ball, remove_small_objects, white_tophat +from skimage.morphology import remove_small_objects from skimage.segmentation import expand_labels # type: ignore[import-untyped] +from typing_extensions import Annotated, Doc from ._constants import * from ._exceptions import DimensionalityError @@ -20,40 +22,63 @@ SKIMAGE_REGIONPROPS_TABLE_COLUMNS_EXPANDED, DetectionResult, ) +from .dog_transform import DifferenceOfGaussiansTransformation __author__ = "Vince Reuter" __credits__ = ["Vince Reuter", "Kai Sandoval Beckwith"] __all__ = ["detect_spots_dog", "detect_spots_int"] +Image = npt.NDArray[PixelValue] Numeric = Union[int, float] +@doc(summary="Parameter descriptions common to various spot detection procedures") +@dataclass +class detection_signature: + image = Annotated[npt.NDArray[PixelValue], Doc("3D image in which to detect spots")] + threshold = Annotated[ + Numeric, + Doc( + "The minimum (after any transformations) pixel value required to regard a pixel as part of a spot" + ), + ] + expand_px = Annotated[ + Optional[Numeric], + Doc("The number of pixels by which to expand a detected and defined region"), + ] + result = Annotated[ + DetectionResult, + Doc( + "Bundle of table of ROI coordinates and data, the image used for spot detection, and region labels array" + ), + ] + + +@doc( + summary="Detect spots by difference of Gaussians filter.", + parameters=dict( + transform="The subtraction-after-smoothing parameterisation that defined DoG", + ), + raises=dict( + TypeError="If the given `transform` isn't specifically a `DifferenceOfGaussiansTransformation`", + ), +) def detect_spots_dog( - input_image: npt.NDArray[PixelValue], + input_image: detection_signature.image, *, - spot_threshold: Numeric, - expand_px: Optional[Numeric], -) -> DetectionResult: - """Spot detection by difference of Gaussians filter - - Arguments - --------- - input_image : ndarray - 3D image in which to detect spots - spot_threshold : int or float - Minimum peak value after the DoG transformation to call a peak/spot - expand_px : float or int or NoneType - Number of pixels by which to expand contiguous subregion, - up to point of overlap with neighboring subregion of image - - Returns - ------- - spotfishing.DetectionResult - Bundle of table of ROI coordinates and data, the image used for spot detection, and region labels array - """ + spot_threshold: detection_signature.threshold, + expand_px: detection_signature.expand_px, + transform: DifferenceOfGaussiansTransformation, +) -> detection_signature.result: + # TODO: consider replacing by something from scikit-image. + # See: https://github.com/gerlichlab/spotfishing/issues/5 _check_input_image(input_image) - img = _preprocess_for_difference_of_gaussians(input_image) + if not isinstance(transform, DifferenceOfGaussiansTransformation): + raise TypeError( + f"For DoG-based detection, the transformation must be of type {DifferenceOfGaussiansTransformation.__name__}; got {type(transform).__name__}" + ) + img = transform(input_image) labels, _ = ndi.label(img > spot_threshold) # type: ignore[no-untyped-call] spot_props, labels = _build_props_table( labels=labels, input_image=input_image, expand_px=expand_px @@ -61,38 +86,24 @@ def detect_spots_dog( return DetectionResult(table=spot_props, image=img, labels=labels) +@doc( + summary="Detect spots by a simply pixel value threshold.", + raises=dict( + TypeError="If the given `transform` isn't specifically a `DifferenceOfGaussiansTransformation`", + ), +) def detect_spots_int( - input_image: npt.NDArray[PixelValue], + input_image: detection_signature.image, *, - spot_threshold: Numeric, - expand_px: Optional[Numeric], -) -> DetectionResult: - """Spot detection by intensity filter - - Arguments - --------- - input_image : ndarray - 3D image in which to detect spots - spot_threshold : NumberLike - Minimum intensity value in a pixel to consider it as part of a spot region - expand_px : float or int or NoneType - Number of pixels by which to expand contiguous subregion, - up to point of overlap with neighboring subregion of image - - Returns - ------- - spotfishing.DetectionResult - Bundle of table of ROI coordinates and data, the image used for spot detection, and region labels array - """ - # TODO: enforce that output column names don't vary with code path walked. - # See: https://github.com/gerlichlab/looptrace/issues/125 + spot_threshold: detection_signature.threshold, + expand_px: detection_signature.expand_px, +) -> detection_signature.result: _check_input_image(input_image) binary = input_image > spot_threshold binary = ndi.binary_fill_holes(binary) # type: ignore[no-untyped-call] struct = ndi.generate_binary_structure(input_image.ndim, 2) # type: ignore[no-untyped-call] labels, num_obj = ndi.label(binary, structure=struct) # type: ignore[no-untyped-call] - if num_obj > 1: - labels = remove_small_objects(labels, min_size=5) + labels = remove_small_objects(labels, min_size=5) if num_obj > 1 else labels spot_props, labels = _build_props_table( labels=labels, input_image=input_image, expand_px=expand_px ) @@ -135,13 +146,3 @@ def _check_input_image(img: npt.NDArray[PixelValue]) -> None: raise DimensionalityError( f"Expected 3D input image but got {img.ndim}-dimensional" ) - - -def _preprocess_for_difference_of_gaussians( - input_image: npt.NDArray[PixelValue], -) -> npt.NDArray[PixelValue]: - img = white_tophat(image=input_image, footprint=ball(2)) - img = gaussian(img, 0.8) - gaussian(img, 1.3) - img = img / gaussian(input_image, 3) - img = (img - np.mean(img)) / np.std(img) - return img # type: ignore[no-any-return] diff --git a/spotfishing/dog_transform.py b/spotfishing/dog_transform.py new file mode 100644 index 0000000..927aefc --- /dev/null +++ b/spotfishing/dog_transform.py @@ -0,0 +1,216 @@ +"""Image processing related to spot detection, but to be run first""" + +import json +from dataclasses import dataclass +from operator import itemgetter +from pathlib import Path +from typing import Callable, Dict, Optional, Union + +import numpy as np +import numpy.typing as npt +from numpydoc_decorator import doc # type: ignore[import-untyped] +from skimage.filters import gaussian as gaussian_filter # type: ignore[import-untyped] +from skimage.morphology import ball, white_tophat +from typing_extensions import Annotated, Doc + +from ._numeric_types import NumpyFloat, NumpyInt, PixelValue + +__author__ = "Vince Reuter" +__credits__ = ["Vince Reuter"] + +__all__ = ["DifferenceOfGaussiansTransformation"] + +Image = npt.NDArray[PixelValue] +Numeric = Union[float, int, NumpyFloat, NumpyInt] + + +@doc( + summary="Function that, when parameterised, takes an image and gives another image of same shape.", + parameters=dict( + _full_func="The function to which to pass an input image and parameters", + parameters="Parameters to pass to `_full_func`, besides the input image argument", + ), +) +@dataclass(eq=False) +class ImageEndomorphism: # type: ignore[misc] + _full_func: Callable[..., Image] # type: ignore[misc] + parameters: Dict[str, object] + + def __call__(self, img: Image) -> Image: + return self._full_func(img, **self.parameters) + + def __eq__(self, other: "ImageEndomorphism") -> bool: # type: ignore[override] + # Handle fact that a parameter may be an array (e.g. ball(2) for white_tophat). + if type(other) is not type(self) or self._full_func != other._full_func: + return False + kvs_this = sorted(self.parameters.items(), key=itemgetter(0)) + kvs_that = sorted(other.parameters.items(), key=itemgetter(0)) + return all( + k1 == k2 + and ( + np.array_equal(v1, v2) # type: ignore[arg-type] + if self._is_arr(v1) and self._is_arr(v2) + else v1 == v2 + ) + for (k1, v1), (k2, v2) in zip(kvs_this, kvs_that) + ) + + @staticmethod + def _is_arr(obj: object) -> bool: + return isinstance(obj, np.ndarray) + + +@doc( + summary="Function that, when parameterised, takes an two images and gives an image of same shape as the second.", + parameters=dict( + _full_func="The function to which to pass an input images and parameters", + parameters="Parameters to pass to `_full_func`, besides the input image arguments", + ), +) +@dataclass +class ImageMorphism21: # type: ignore[misc] + _full_func: Callable[..., Image] # type: ignore[misc] + parameters: Dict[str, object] + + def __call__(self, *, old_img: Image, new_img: Image) -> Image: + return self._full_func(old_img=old_img, new_img=new_img, **self.parameters) + + +class common_params: + sigma_narrow = Annotated[ + str, + Doc( + "Standard deviation of the more concentrated of the two smoothing distributions" + ), + ] + sigma_wide = Annotated[ + str, + Doc( + "Standard deviation of the more diffuse of the two smoothing distributions" + ), + ] + standardise = Annotated[str, Doc("Standardise (mean-0 and spread-1) the image")] + + +@doc( + summary="Bundle of parameters defining preprocessing for Difference of Gaussians (DoG) detector", + extended_summary=""" + Chain together a series of operations that to prepare an image for spot detection by + difference of Gaussians. This includes an optional step to do before taking the difference + of Gaussian-smoothed images, the spreads of the Gaussians to use for the smoothing, and + what (if anything) to do after subtracting the smoothed versions of the image, + """, + parameters=dict( + pre_diff="Transformation to apply to image before taking difference of smoothed versions", + sigma_narrow=common_params.sigma_narrow, + sigma_wide=common_params.sigma_wide, + post_diff="What (if anything) to do to the transformed image after difference but before standardisation; note that if non-null, this will be fed the *original* image as the first argument, and the *transformed image as the second argument", + standardise=common_params.standardise, + ), + raises=dict( + TypeError="If either of the standard deviations is non-numeric.", + ValueError="If the narrower Gaussian's standard deviation isn't less than the wider Gaussian's.", + ), + returns="A structure of the same shape as the input, just with all transformations applied", +) +@dataclass(frozen=True, kw_only=True) +class DifferenceOfGaussiansTransformation: + pre_diff: Optional[ImageEndomorphism] + sigma_narrow: Numeric + sigma_wide: Numeric + post_diff: Optional[ImageMorphism21] + standardise: bool + + def __post_init__(self) -> None: + # skimage raises errors for negative sigma, but we raise these. + if not (is_numeric(self.sigma_narrow) and is_numeric(self.sigma_wide)): + raise TypeError( + f"At least one sigma is non-numeric! narrow={self.sigma_narrow} ({type(self.sigma_narrow).__name__}), wide={self.sigma_wide} ({type(self.sigma_wide).__name__}" + ) + if self.sigma_narrow >= self.sigma_wide: + raise ValueError( + f"sigma for narrow Gaussian must be strictly less than sigma for wide Gaussian, but {self.sigma_narrow} >= {self.sigma_wide}" + ) + + @doc( + summary="Apply the sequence of transformations in this instance to given image.", + extended_summary=""" + The construction of an instance of this preprocessor parameterises several + operations that, when run in sequence, comprise a preprocessing workflow often + used before running spot detection. + """, + parameters=dict(input_image="The image (array of pixel values) to preprocess"), + returns="The array of values after transformations' application", + ) + def __call__(self, input_image: Image) -> Image: + img = self.pre_diff(input_image) if self.pre_diff is not None else input_image + img_narrow = gaussian_filter(img, self.sigma_narrow) + img_wide = gaussian_filter(img, self.sigma_wide) + img = img_narrow - img_wide + img = ( + self.post_diff(old_img=input_image, new_img=img) + if self.post_diff is not None + else img + ) + return (img - np.mean(img)) / np.std(img) if self.standardise else img # type: ignore[no-any-return] + + # @doc( + # summary="Build an instance with optional white tophat filter as preprocessing, and optional division by a Gaussian blur as postprocessing.", + # parameters=dict( + # apply_white_tophat="Whether to apply a white tophat transformation before taking the difference of the Gaussians", + # sigma_narrow=common_params.sigma_narrow, + # sigma_wide=common_params.sigma_wide, + # sigma_post_divide="The standard deviation of the Gaussian blur to apply after differencing but before potential standardisation", + # standardise=common_params.standardise, + # ), + # returns="A parameterised instance built according to the specification of the arguments provided here", + # ) + @classmethod + def build( + cls, + *, + apply_white_tophat: bool, + sigma_narrow: Numeric, + sigma_wide: Numeric, + sigma_post_divide: Optional[Numeric], + standardise: bool, + ) -> "DifferenceOfGaussiansTransformation": + # https://git.embl.de/grp-ellenberg/looptrace/-/blob/master/looptrace/image_processing_functions.py?ref_type=heads#L252 + return cls( + pre_diff=ImageEndomorphism(white_tophat, dict(footprint=ball(2))) + if apply_white_tophat + else None, + sigma_narrow=sigma_narrow, + sigma_wide=sigma_wide, + post_diff=None + if sigma_post_divide is None + else ImageMorphism21(div_by_gauss, dict(sigma=3)), + standardise=standardise, + ) + + # @doc( + # summary="Build an instance from parameters in a JSON file.", + # parameters=dict( + # fp="Path to the file to parse for parameters", + # ), + # returns="A newly constructed instance, parameterised according to the data in the given config file", + # see_also=":meth:`DifferenceOfGaussiansTransformation.build`", + # ) + @classmethod + def from_json_file(cls, fp: Path) -> "DifferenceOfGaussiansTransformation": + with open(fp, "r") as fh: + data = json.load(fh) + return cls.build(**data) + + +def div_by_gauss(*, old_img: Image, new_img: Image, sigma: Numeric) -> Image: + # https://git.embl.de/grp-ellenberg/looptrace/-/blob/master/looptrace/image_processing_functions.py?ref_type=heads#L252 + return new_img / gaussian_filter(old_img, sigma) # type: ignore[no-any-return] + + +@doc( + summary="Test the given object for membership in a numeric type", + parameters=dict(obj="Object to test for membership in a numeric type"), +) +def is_numeric(obj: object) -> bool: + return isinstance(obj, Numeric) # type: ignore diff --git a/spotfishing/examples/original__difference_of_gaussians.json b/spotfishing/examples/original__difference_of_gaussians.json new file mode 100644 index 0000000..c1118b0 --- /dev/null +++ b/spotfishing/examples/original__difference_of_gaussians.json @@ -0,0 +1,7 @@ +{ + "apply_white_tophat": true, + "sigma_narrow": 0.8, + "sigma_wide": 1.3, + "sigma_post_divide": 3, + "standardise": true +} \ No newline at end of file diff --git a/tests/helpers.py b/tests/helpers.py index ed8c9b1..55193b3 100644 --- a/tests/helpers.py +++ b/tests/helpers.py @@ -2,16 +2,35 @@ import os from pathlib import Path +from typing import Union import numpy as np +from skimage.morphology import ball, white_tophat + +from spotfishing import DifferenceOfGaussiansTransformation as DogTransform +from spotfishing.dog_transform import ImageEndomorphism, ImageMorphism21, div_by_gauss __author__ = "Vince Reuter" __credits__ = ["Vince Reuter"] __all__ = [ + "ORIGINAL_TRANSFORM", + "get_img_data_file", "load_image_file", ] +Numeric = Union[float, int] + + +# original parameterisation of the transformation for detection with DoG +ORIGINAL_TRANSFORM = DogTransform( + pre_diff=ImageEndomorphism(white_tophat, dict(footprint=ball(2))), + sigma_narrow=0.8, + sigma_wide=1.3, + post_diff=ImageMorphism21(div_by_gauss, dict(sigma=3)), + standardise=True, +) + def get_img_data_file(fn: str) -> Path: """Get the path to an input file (image data).""" diff --git a/tests/test_accord_with_original_settings.py b/tests/test_accord_with_original_settings.py index 65225fa..611b9d0 100644 --- a/tests/test_accord_with_original_settings.py +++ b/tests/test_accord_with_original_settings.py @@ -1,20 +1,40 @@ """Tests for maintaining accordance with original, fixed settings, as flexibility is added""" +import importlib.resources as importlib_resources +from functools import partial from pathlib import Path +from typing import Callable, Optional import pandas as pd import pytest -from helpers import load_image_file +from helpers import ORIGINAL_TRANSFORM, load_image_file from pandas.testing import assert_frame_equal +import spotfishing from spotfishing import ROI_MEAN_INTENSITY_KEY, detect_spots_dog, detect_spots_int -from spotfishing.detection_result import ROI_CENTER_KEYS +from spotfishing.detection_result import ROI_CENTER_KEYS, DetectionResult +from spotfishing.detectors import Image, Numeric +from spotfishing.dog_transform import DifferenceOfGaussiansTransformation __author__ = "Vince Reuter" __credits__ = ["Vince Reuter"] -COLUMNS_OF_INTEREST = ROI_CENTER_KEYS + [ROI_MEAN_INTENSITY_KEY] +# original pixel expansion setting for detection with DoG +ORIGINAL_DOG_EXPAND_PX = 10 + +# original pixel expansion setting for detection with simple intensity threshold +ORIGINAL_DOG_EXPAND_PX = 1 + + +def test_original_settings_from_json(): + data_file = ( + importlib_resources.files(spotfishing) + / "examples" + / "original__difference_of_gaussians.json" + ) + obs = DifferenceOfGaussiansTransformation.from_json_file(data_file) + assert obs == ORIGINAL_TRANSFORM @pytest.mark.parametrize( @@ -28,7 +48,13 @@ id=f"{func_type_name}__{threshold}", ) for detect, func_type_name, threshold, expand_px in [ - (detect_spots_dog, "diff_gauss", t, 10) for t in [15, 10] + ( + partial(detect_spots_dog, transform=ORIGINAL_TRANSFORM), + "diff_gauss", + t, + 10, + ) + for t in [15, 10] ] + [(detect_spots_int, "intensity", t, 1) for t in [300, 200]] ], @@ -42,15 +68,25 @@ def test_output_is_correct_with_original_settings( ) exp_path = Path(__file__).parent / "data" / "outputs__original_settings" / exp_name exp_table = pd.read_csv(exp_path, index_col=0) - - arr_name = f"img__{data_name}__smaller.npy" - input_image = load_image_file(arr_name) - obs = detect(input_image, spot_threshold=threshold, expand_px=expand_px) - obs_table = obs.table[COLUMNS_OF_INTEREST] - + obs_table = get_table_result( + detect=detect, threshold=threshold, expand_px=expand_px, data_name=data_name + ) print("EXPECTED (below):") print(exp_table) print("OBSERVED (below):") print(obs_table) - assert_frame_equal(obs_table, exp_table) + + +def get_table_result( + *, + detect: Callable[[Image, Numeric, Optional[Numeric]], DetectionResult], + threshold: Numeric, + expand_px: Optional[Numeric], + data_name: str, +) -> pd.DataFrame: + arr_name = f"img__{data_name}__smaller.npy" + input_image = load_image_file(arr_name) + obs = detect(input_image, spot_threshold=threshold, expand_px=expand_px) # type: ignore[call-arg] + cols = ROI_CENTER_KEYS + [ROI_MEAN_INTENSITY_KEY] + return obs.table[cols] diff --git a/tests/test_old_new_equivalence.py b/tests/test_old_new_equivalence.py deleted file mode 100644 index 2f06b17..0000000 --- a/tests/test_old_new_equivalence.py +++ /dev/null @@ -1,45 +0,0 @@ -import os -from pathlib import Path - -import numpy as np -import pytest -from helpers import load_image_file - -from spotfishing import detect_spots_dog, detect_spots_int -from spotfishing.deprecated import detect_spots_dog_old, detect_spots_int_old - -OLD_PIXEL_EXPANSION_INT = 1 -OLD_PIXEL_EXPANSION_DOG = 10 - - -def get_img_data_file(fn: str) -> Path: - return Path(os.path.dirname(__file__)) / "data" / "inputs" / fn - - -@pytest.mark.parametrize( - "input_image", - [ - load_image_file(fn) - for fn in ("img__p0_t57_c0__smaller.npy", "img__p13_t57_c0__smaller.npy") - ], -) -@pytest.mark.parametrize( - ["old_fun", "new_fun", "threshold", "expand_px"], - [ - (detect_spots_int_old, detect_spots_int, threshold, OLD_PIXEL_EXPANSION_INT) - for threshold in (500, 300, 100) - ] - + [ - (detect_spots_dog_old, detect_spots_dog, threshold, OLD_PIXEL_EXPANSION_DOG) - for threshold in (20, 15, 10) - ], -) -def test_eqv(input_image, old_fun, new_fun, threshold, expand_px): - old_table, _, _ = old_fun(input_image, threshold) - new_res = new_fun(input_image, spot_threshold=threshold, expand_px=expand_px) - new_table = new_res.table - assert np.all(old_table.index == new_table.index) - cols = ["zc", "yc", "xc", "area", "intensity_mean"] - sub_cols = [c for c in cols if c != "area"] # wasn't in DoG before - assert list(new_table.columns) == cols - assert np.all(old_table[sub_cols].to_numpy() == new_table[sub_cols].to_numpy()) diff --git a/tests/test_spot_detectors.py b/tests/test_spot_detectors.py index 2d5d7d3..fa10799 100644 --- a/tests/test_spot_detectors.py +++ b/tests/test_spot_detectors.py @@ -1,10 +1,12 @@ """Tests about column names of spot tables""" +from functools import partial + import hypothesis as hyp import numpy as np import pandas as pd import pytest -from helpers import load_image_file +from helpers import ORIGINAL_TRANSFORM, load_image_file from spotfishing import DimensionalityError, detect_spots_dog, detect_spots_int from spotfishing.detection_result import DETECTION_RESULT_TABLE_COLUMNS @@ -25,7 +27,7 @@ # detection function, threshold setting, and pixel expansion setting BASE_INTENSITY_BUNDLE = (detect_spots_int, 300, 1) -BASE_DOG_BUNDLE = (detect_spots_dog, 15, 10) +BASE_DOG_BUNDLE = (partial(detect_spots_dog, transform=ORIGINAL_TRANSFORM), 15, 10) # @pytest.mark.parametrize("detect", [detect_spots_dog, detect_spots_int]) @@ -58,7 +60,10 @@ def test_simple_intensity_detector_result_always_contains_original_image(input_i assert np.all(result.image == input_image) -@pytest.mark.parametrize("detect", [detect_spots_dog, detect_spots_int]) +@pytest.mark.parametrize( + "detect", + [partial(detect_spots_dog, transform=ORIGINAL_TRANSFORM), detect_spots_int], +) @pytest.mark.parametrize( ["input_image", "expected_error_type", "expected_message"], [