add fix to anisotropy bugs

pyproject.toml

@@ -77,6 +77,10 @@ root = "."
 [tool.ruff]
 target-version = "py311"
 line-length = 88
+# ignore all linting in the theoretical_explainations directory, which is used for testing and development of new features
+extend-exclude = [
+  "theorectical_explainations/notebooks/interpolation_edge_demo.ipynb",
+]  # Ignore `E501` (line too long) for Jupyter notebooks
 fix = true
 lint.select = [
   # flake8-builtins
@@ -101,9 +105,10 @@ lint.select = [
   "SIM",
   "W",
 ]
-# Ignore `E501` (line too long) for Jupyter notebooks
 lint.per-file-ignores."*.ipynb" = [ "C408", "E402", "E501" ]
 lint.per-file-ignores."src/**/__init__.py" = [ "F401" ]
+lint.per-file-ignores."src/zedprofiler/featurization/neighbors.py" = [ "PLR0913", "PLR0917" ]
+lint.per-file-ignores."src/zedprofiler/featurization/texture.py" = [ "RUF046" ]
 
 [tool.codespell]
 ignore-words-list = "doesnot,inout"

src/zedprofiler/IO/loading_classes.py

@@ -11,7 +11,10 @@
 import numpy
 from beartype import beartype
 
-from zedprofiler.contracts import ImageArrayModel
+from zedprofiler.contracts import (
+    ImageArrayModel,
+    validate_anisotropy_factor_with_pydantic,
+)
 from zedprofiler.identifiers import build_image_id
 
 logging.basicConfig(level=logging.INFO)
@@ -216,7 +219,11 @@ def __init__(  # noqa: PLR0913
         channel_tokens = [str(value) for value in channel_mapping.values()]
         self._label_key_names = list(config.label_key_name or [])
         self.anisotropy_spacing = anisotropy_spacing
-        self.anisotropy_factor = self.anisotropy_spacing[0] / self.anisotropy_spacing[1]
+        self.anisotropy_factor = validate_anisotropy_factor_with_pydantic(
+            self.anisotropy_spacing[0] / self.anisotropy_spacing[1],
+            y_spacing=self.anisotropy_spacing[1],
+            x_spacing=self.anisotropy_spacing[2],
+        ).anisotropy_factor
         self.image_set_name = config.image_set_name
         self.label_set_path = label_set_path
         # Deterministic imaging identifier for the warehouse join key
@@ -298,7 +305,11 @@ def from_image_dict(  # noqa: PLR0913
             self.image_set_dict[key] = ImageArrayModel(array=array).array
         self._label_key_names = list(label_key_names or [])
         self.anisotropy_spacing = anisotropy_spacing
-        self.anisotropy_factor = self.anisotropy_spacing[0] / self.anisotropy_spacing[1]
+        self.anisotropy_factor = validate_anisotropy_factor_with_pydantic(
+            self.anisotropy_spacing[0] / self.anisotropy_spacing[1],
+            y_spacing=self.anisotropy_spacing[1],
+            x_spacing=self.anisotropy_spacing[2],
+        ).anisotropy_factor
         self.image_set_name = image_set_name
         self.label_set_path = None
         config = ImageSetConfig(
@@ -538,7 +549,11 @@ def get_anisotropy(self) -> float:
             Ratio of z-spacing to y-spacing.
 
         """
-        return self.anisotropy_spacing[0] / self.anisotropy_spacing[1]
+        return validate_anisotropy_factor_with_pydantic(
+            self.anisotropy_spacing[0] / self.anisotropy_spacing[1],
+            y_spacing=self.anisotropy_spacing[1],
+            x_spacing=self.anisotropy_spacing[2],
+        ).anisotropy_factor
 
 
 class ObjectLoader:

src/zedprofiler/contracts.py

@@ -14,6 +14,7 @@
 
 from __future__ import annotations
 
+import math
 from typing import Any
 
 import numpy as np
@@ -29,6 +30,7 @@
 
 from zedprofiler.exceptions import ContractError
 
+MIN_ANISOTROPY_FACTOR = 1
 EXPECTED_SPATIAL_DIMS = 3
 TWO_DIMENSIONAL = 2
 FOUR_DIMENSIONAL = 4
@@ -128,6 +130,52 @@ def validate_array_dtype_and_shape(_cls, arr: np.ndarray) -> np.ndarray:
         return arr
 
 
+class AnisotropyFactorModel(BaseModel):
+    """Pydantic model for validating the anisotropy factor.
+
+    Feature modules (e.g. texture, neighbors) assume the z-spacing is never
+    finer than the x/y-spacing, so ``anisotropy_factor`` (z_spacing /
+    y_spacing) must be 1 or greater. ``anisotropy_factor`` collapses to a
+    single z/y ratio, so it is only meaningful when the transverse (y, x)
+    spacings are equal; ``y_spacing``/``x_spacing`` are optional and, when
+    both are provided, are checked for equality.
+    """
+
+    anisotropy_factor: float
+    y_spacing: float | None = None
+    x_spacing: float | None = None
+
+    @field_validator("anisotropy_factor", mode="after")
+    @classmethod
+    def validate_at_least_one(_cls, value: float) -> float:
+        """Ensure the anisotropy factor is finite and 1 or greater."""
+        if not math.isfinite(value):
+            raise ValueError(
+                f"Anisotropy factor must be a finite number, got {value}.",
+            )
+        if value < MIN_ANISOTROPY_FACTOR:
+            raise ValueError(
+                f"Anisotropy factor must be {MIN_ANISOTROPY_FACTOR} or greater, "
+                f"got {value}.",
+            )
+        return value
+
+    @model_validator(mode="after")
+    def validate_transverse_spacing_equal(self) -> AnisotropyFactorModel:
+        """Ensure y/x spacing are equal, when both are provided."""
+        if (
+            self.y_spacing is not None
+            and self.x_spacing is not None
+            and self.y_spacing != self.x_spacing
+        ):
+            raise ValueError(
+                "anisotropy_factor is a single z/y ratio and requires equal "
+                f"transverse spacing, got y_spacing={self.y_spacing} and "
+                f"x_spacing={self.x_spacing}.",
+            )
+        return self
+
+
 class FeatureDictModel(BaseModel):
     """Pydantic model for validating feature dictionaries."""
 
@@ -371,6 +419,51 @@ def validate_return_with_pydantic(
         raise ContractError(msg)
 
 
+def validate_anisotropy_factor_with_pydantic(
+    anisotropy_factor: float,
+    y_spacing: float | None = None,
+    x_spacing: float | None = None,
+) -> AnisotropyFactorModel:
+    """Validate the anisotropy factor using a Pydantic model.
+
+    Parameters
+    ----------
+    anisotropy_factor : float
+        Ratio of z-spacing to y-spacing to validate.
+    y_spacing : float | None, optional
+        Y (transverse) spacing. When provided along with ``x_spacing``, they
+        are checked for equality since ``anisotropy_factor`` only captures a
+        single z/y ratio.
+    x_spacing : float | None, optional
+        X (transverse) spacing. See ``y_spacing``.
+
+    Returns
+    -------
+    AnisotropyFactorModel
+        Validated anisotropy factor model.
+
+    Raises
+    ------
+    ContractError
+        If the anisotropy factor is less than ``MIN_ANISOTROPY_FACTOR``, is
+        not finite, or if ``y_spacing`` and ``x_spacing`` are unequal.
+
+    """
+    try:
+        return AnisotropyFactorModel(
+            anisotropy_factor=anisotropy_factor,
+            y_spacing=y_spacing,
+            x_spacing=x_spacing,
+        )
+    except Exception as e:
+        msg = (
+            "Anisotropy factor validation failed. Please ensure that the "
+            "z-spacing is not finer than the y-spacing and that the y- and "
+            f"x-spacing are equal: {e}"
+        )
+        raise ContractError(msg)
+
+
 def validate_image_with_pydantic(arr: np.ndarray) -> ImageArrayModel:
     """Validate the input image array using Pydantic model.
 

src/zedprofiler/featurization/granularity.py

@@ -2,6 +2,8 @@
 
 from __future__ import annotations
 
+import math
+
 import numpy
 import pandas
 import scipy.ndimage
@@ -12,6 +14,56 @@
 from zedprofiler.IO.loading_classes import ObjectLoader
 
 
+def anisotropic_ball(
+    radius: int,
+    spacing: tuple[float, float, float] | None = None,
+) -> numpy.ndarray:
+    """Build a spherical structuring element that is physically isotropic.
+
+    Parameters
+    ----------
+    radius : int
+        Radius of the structuring element in voxel units.
+    spacing : tuple[float, float, float] or None
+        Physical spacing of the image in (z, y, x) order.
+        If None, the structuring element is isotropic in voxel space.
+        If provided, the structuring element will be isotropic
+        in physical space, taking into account the anisotropy
+        of the voxel spacing.
+
+    Returns
+    -------
+    numpy.ndarray
+        A boolean array representing the structuring element,
+        where True values indicate the presence of the struct
+        during element and False values indicate the absence.
+
+    ...
+    """
+    if spacing is None:
+        return skimage.morphology.ball(radius, dtype=bool)
+
+    z_spacing, y_spacing, x_spacing = spacing
+    if z_spacing == y_spacing == x_spacing:
+        return skimage.morphology.ball(radius, dtype=bool)
+
+    min_spacing = min(z_spacing, y_spacing, x_spacing)
+    physical_radius = radius * min_spacing
+
+    # Largest voxel offset on each axis that can still land within the
+    # physical radius. floor() (not round()) so a coarse axis correctly
+    # collapses to 0 when even one voxel step overshoots the radius.
+    rz = math.floor(physical_radius / z_spacing)
+    ry = math.floor(physical_radius / y_spacing)
+    rx = math.floor(physical_radius / x_spacing)
+
+    zz, yy, xx = numpy.ogrid[-rz : rz + 1, -ry : ry + 1, -rx : rx + 1]
+    physical_dist_sq = (
+        (zz * z_spacing) ** 2 + (yy * y_spacing) ** 2 + (xx * x_spacing) ** 2
+    )
+    return physical_dist_sq <= physical_radius**2
+
+
 def _fix_scipy_ndimage_result(result: float | list | numpy.ndarray) -> numpy.ndarray:
     """Convert scipy.ndimage aggregation results to a consistent array.
 
@@ -156,11 +208,17 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
 ) -> pandas.DataFrame:
     """Calculate the granularity spectrum of a 3D image.
 
-    Follows the CellProfiler MeasureGranularity algorithm exactly for 3D:
+    Based on the CellProfiler MeasureGranularity algorithm, generalized to 3D:
     1. Subsample the image uniformly (same factor for Z, Y, X).
     2. Further subsample for background tophat removal.
-    3. Iteratively erode with ball(1) and reconstruct, measuring
-    signal lost at each scale as image-level and per-object values.
+    3. Iteratively erode with a spherical structuring element and
+    reconstruct, measuring signal lost at each scale as image-level and
+    per-object values.
+
+    The structuring elements used for background removal and the erosion
+    spectrum are physically isotropic spheres (see ``anisotropic_ball``),
+    not raw voxel-space spheres, so results are correct rather than biased
+    when z-spacing differs from x/y-spacing.
 
     Parameters
     ----------
@@ -220,6 +278,7 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
     original_pixels = object_loader.image
     original_labels = object_loader.label_image
     original_shape = original_pixels.shape
+    spacing = object_loader.image_set_loader.anisotropy_spacing
 
     # Mask: CellProfiler uses im.mask (typically all-True for unmasked images)
     if image_mask is None:
@@ -265,24 +324,24 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
     # ------------------------------------------------------------------
     # Step 2: Background removal via tophat filter
     #
-    # CellProfiler 3D BUG (replicated for compatibility):
-    #   The 3D branch uses new_shape for grid bounds and subsample_size
-    #   for coordinate division, instead of back_shape and
-    #   image_sample_size as the 2D branch does. This means:
-    #   - back_pixels has the SAME shape as pixels (not smaller)
-    #   - Many coordinates are out of bounds → map_coordinates returns 0
-    #   We replicate this exactly to match CellProfiler output.
+    # Downsample the (already subsampled) image and mask to back_shape
+    # for the background estimate, exactly mirroring CellProfiler's 2D
+    # branch: grid bounds are back_shape, and coordinates are divided by
+    # image_sample_size to map back into the new_shape-sized `pixels`
+    # array. CellProfiler's actual 3D implementation uses new_shape /
+    # subsample_size here instead, a bug that leaves back_pixels the same
+    # size as pixels (mostly zero-filled from out-of-bounds sampling) and
+    # applies the tophat radius at the wrong scale. We intentionally do
+    # not replicate that 3D bug.
     # ------------------------------------------------------------------
     if image_sample_size < 1.0:
         back_shape = new_shape * image_sample_size
 
-        # CellProfiler 3D: mgrid[0:new_shape] / subsample_size
-        # (NOT mgrid[0:back_shape] / image_sample_size as 2D does)
         k, i, j = (
-            numpy.mgrid[0 : new_shape[0], 0 : new_shape[1], 0 : new_shape[2]].astype(
+            numpy.mgrid[0 : back_shape[0], 0 : back_shape[1], 0 : back_shape[2]].astype(
                 float,
             )
-            / subsample_size
+            / image_sample_size
         )
         back_pixels = scipy.ndimage.map_coordinates(pixels, (k, i, j), order=1)
         back_mask = (
@@ -302,7 +361,7 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
         back_shape = new_shape
 
     # Tophat filter: masked erosion + masked dilation
-    footprint_bg = skimage.morphology.ball(radius, dtype=bool)
+    footprint_bg = anisotropic_ball(radius, spacing)
 
     back_pixels_masked = numpy.zeros_like(back_pixels)
     back_pixels_masked[back_mask] = back_pixels[back_mask]
@@ -315,10 +374,10 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
         footprint=footprint_bg,
     )
 
-    # Upsample background back to subsampled image size
+    # Upsample background back to subsampled image size: grid over
+    # new_shape, with coordinates scaled by (back_shape - 1) / (new_shape - 1)
+    # to map back into the back_shape-sized back_pixels array.
     if image_sample_size < 1.0:
-        # CellProfiler 3D: mgrid[0:new_shape] with coords scaled by
-        # (back_shape - 1) / (new_shape - 1)
         k, i, j = numpy.mgrid[
             0 : new_shape[0],
             0 : new_shape[1],
@@ -400,8 +459,9 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
     ero[~mask] = 0
     currentmean = startmean
 
-    # CellProfiler uses ball(1) for the iterative erosion/reconstruction loop
-    footprint = skimage.morphology.ball(1, dtype=bool)
+    # Physically-isotropic radius-1 structuring element for the iterative
+    # erosion/reconstruction loop (see anisotropic_ball).
+    footprint = anisotropic_ball(1, spacing)
 
     if verbose:
         print(

src/zedprofiler/featurization/intensity.py

@@ -30,10 +30,7 @@ def get_outline(mask: numpy.ndarray) -> numpy.ndarray:
         The outline of the mask.
 
     """
-    outline = numpy.zeros_like(mask)
-    for z in range(mask.shape[0]):
-        outline[z] = skimage.segmentation.find_boundaries(mask[z], mode="inner")
-    return outline
+    return skimage.segmentation.find_boundaries(mask, mode="inner")
 
 
 def compute_intensity(  # noqa: C901, PLR0915
@@ -189,9 +186,16 @@ def compute_intensity(  # noqa: C901, PLR0915
         else:
             cmi_x = cmi_y = cmi_z = numpy.nan
         # calculate the center of mass distance
-        diff_x = cm_x - cmi_x
-        diff_y = cm_y - cmi_y
-        diff_z = cm_z - cmi_z
+        # Scale each axis by its physical voxel spacing before combining into
+        # one Euclidean distance -- diff_x/diff_y/diff_z are raw voxel-index
+        # offsets, and mixing them unscaled would bias the result whenever
+        # z-spacing differs from x/y-spacing.
+        z_spacing, y_spacing, x_spacing = (
+            object_loader.image_set_loader.anisotropy_spacing
+        )
+        diff_x = (cm_x - cmi_x) * x_spacing
+        diff_y = (cm_y - cmi_y) * y_spacing
+        diff_z = (cm_z - cmi_z) * z_spacing
         # mass displacement
         mass_displacement = numpy.sqrt(diff_x**2 + diff_y**2 + diff_z**2)
         # mean absolute deviation