Add benchmarking test framework and improve featurization performance

docs/src/benchmarking.md

@@ -0,0 +1,58 @@
+# Benchmarking and Accuracy Locks
+
+ZEDprofiler performance work should be treated as a score-improvement loop:
+make one optimization, keep the same feature outputs, then compare the scorecard.
+
+The test suite now has two layers for this.
+
+## Accuracy locks
+
+Accuracy locks run with normal pytest:
+
+```bash
+uv run pytest tests/test_benchmark_contracts.py
+```
+
+These tests run deterministic synthetic image-set cases across the major feature families and compare each output dataframe to a current-result fingerprint.
+If an optimization changes a feature value, column name, row identity, or missing-value pattern, the fingerprint changes and the test fails.
+
+Only update these fingerprints when the project intentionally changes the
+scientific result contract.
+
+## Benchmark scorecard
+
+The benchmark scorecard is opt-in:
+
+```bash
+uv run pytest tests/test_benchmark_contracts.py --run-benchmarks -s
+```
+
+or:
+
+```bash
+ZEDPROFILER_RUN_BENCHMARKS=1 uv run pytest tests/test_benchmark_contracts.py -s
+```
+
+The scorecard prints elapsed seconds, row count, column count, and the same
+result fingerprint for each feature family. It includes the small accuracy-lock
+level, a larger many-object synthetic scaling level, and a representative
+real-world data level from checked-in tutorial image data. Timing is
+intentionally not a hard gate yet because local machines and shared compute
+environments vary. Use the scorecard to compare before and after runs on the
+same environment.
+
+## Game loop
+
+Use this loop for each performance pass:
+
+1. Run the accuracy locks.
+1. Run the benchmark scorecard and keep the printed output.
+1. Optimize one hotspot.
+1. Re-run the accuracy locks.
+1. Re-run the benchmark scorecard in the same environment.
+1. Keep the optimization only when fingerprints stay fixed and scorecard timing
+   improves enough to matter.
+
+For larger-scale runs, pair this with workflow-level characterization. Local
+scorecards explain feature-kernel cost, while workflow traces explain scheduling,
+batching, queue, and filesystem behavior.

docs/src/index.md

@@ -49,4 +49,5 @@ api
 :maxdepth: 2
 :caption: A Note on Scalability
 scalability
+benchmarking
 ```

src/zedprofiler/contracts.py

@@ -17,7 +17,6 @@
 from typing import Any
 
 import numpy as np
-import pandas as pd
 import pandera.pandas as pa
 from beartype import beartype
 from pydantic import (
@@ -567,43 +566,38 @@ def validate_column_name_schema(
                 f"underscores, got {len(parts)} parts in '{column_name}'",
             )
         return True
-    feature_components = pd.DataFrame(
-        [
-            {
-                "compartment": parts[0],
-                "channel": parts[1],
-                "feature": parts[2],
-            },
-        ],
-    )
-
-    feature_component_schema = pa.DataFrameSchema(
-        {
-            "compartment": pa.Column(
-                str,
-                checks=pa.Check.isin(expected_values.get("compartments", [])),
-                nullable=False,
-                coerce=True,
-            ),
-            "channel": pa.Column(
-                str,
-                checks=pa.Check.isin(expected_values.get("channels", [])),
-                nullable=False,
-                coerce=True,
-            ),
-            "feature": pa.Column(
-                str,
-                checks=pa.Check.isin(expected_values.get("features", [])),
-                nullable=False,
-                coerce=True,
-            ),
-        },
-        strict=True,
-    )
-
-    try:
-        feature_component_schema.validate(feature_components)
-    except (pa.errors.SchemaError, pa.errors.SchemaErrors) as e:
-        raise ContractError(f"Column name schema validation failed: {e}") from e
+    # Validate the non-metadata column's compartment/channel/feature components.
+    #
+    # This is equivalent to the prior pandera DataFrameSchema check
+    # (Check.isin with coerce=True, nullable=False, strict=True applied to
+    # parts[0..2]) but expressed as a plain set-membership test. ``parts`` are
+    # always strings produced by ``column_name.split("_")``, so the pandera
+    # ``coerce`` (str cast), ``nullable`` (no NaN), and ``strict`` (exactly
+    # three components) considerations are no-ops here, and pass/fail semantics
+    # are identical.
+    #
+    # The motivation is overhead: this function is called once per output
+    # column by every feature module's tail validation loop, and rebuilding a
+    # pandas DataFrame plus a pandera DataFrameSchema and running ``.validate``
+    # on a one-row frame per column dominated feature extraction time (e.g.
+    # ~67% of texture runtime). A direct membership check removes that cost
+    # without changing any feature value, column name, or row identity.
+    feature_components = {
+        "compartment": parts[0],
+        "channel": parts[1],
+        "feature": parts[2],
+    }
+    for component, allowed_key in (
+        ("compartment", "compartments"),
+        ("channel", "channels"),
+        ("feature", "features"),
+    ):
+        allowed = expected_values.get(allowed_key, [])
+        if feature_components[component] not in allowed:
+            raise ContractError(
+                f"Column name schema validation failed: "
+                f"{component} '{feature_components[component]}' in "
+                f"'{column_name}' is not in allowed {allowed_key} {list(allowed)}.",
+            )
 
     return True

src/zedprofiler/featurization/colocalization.py

@@ -598,15 +598,28 @@ def compute_colocalization(  # noqa: C901, PLR0912
     if channel1 is None or channel2 is None:
         raise ValueError("channel1 and channel2 must be provided for feature naming.")
     list_of_dfs = []
-    for object_id in two_object_loader.object_ids:
-        cropped_image1, cropped_image2 = prepare_two_images_for_colocalization(
-            label_object1=two_object_loader.label_image,
-            label_object2=two_object_loader.label_image,
-            image_object1=two_object_loader.image1,
-            image_object2=two_object_loader.image2,
-            object_id1=object_id,
-            object_id2=object_id,
+    props = skimage.measure.regionprops_table(
+        two_object_loader.label_image,
+        properties=["label", "bbox"],
+    )
+    label_to_bbox = {
+        int(label): (
+            int(props["bbox-0"][index]),
+            int(props["bbox-1"][index]),
+            int(props["bbox-2"][index]),
+            int(props["bbox-3"][index]),
+            int(props["bbox-4"][index]),
+            int(props["bbox-5"][index]),
         )
+        for index, label in enumerate(props.get("label", []))
+    }
+
+    for object_id in two_object_loader.object_ids:
+        bbox = label_to_bbox.get(int(object_id))
+        if bbox is None:
+            continue
+        cropped_image1 = crop_3D_image(two_object_loader.image1, bbox)
+        cropped_image2 = crop_3D_image(two_object_loader.image2, bbox)
         colocalization_features = calculate_colocalization(
             cropped_image_1=cropped_image1,
             cropped_image_2=cropped_image2,

src/zedprofiler/featurization/granularity.py

@@ -327,8 +327,15 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
 
     if nobjects > 0:
         # CellProfiler: self.labels[~im.mask] = 0
-        masked_labels = original_labels.copy()
-        masked_labels[~original_mask] = 0
+        # When the mask covers the whole image (the common unmasked case,
+        # image_mask=None) no labels get zeroed, so skip the full-array copy
+        # and reuse original_labels directly. scipy.ndimage.mean does not
+        # mutate its label input, so this is safe.
+        if original_mask.all():
+            masked_labels = original_labels
+        else:
+            masked_labels = original_labels.copy()
+            masked_labels[~original_mask] = 0
 
         if numpy.any(masked_labels > 0):
             per_object_current_mean = _fix_scipy_ndimage_result(
@@ -348,7 +355,17 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
     # CellProfiler computes startmean AFTER background subtraction but
     # BEFORE zeroing pixels outside mask (zeroing is implicit via indexing).
     # ------------------------------------------------------------------
-    startmean = numpy.mean(pixels[mask]) if mask.any() else 0.0
+    # Whether the (possibly subsampled) mask covers the whole image. In the
+    # common unmasked case (image_mask=None) this is True, and several per-scale
+    # operations below become no-ops or can avoid full-array copies: the masked
+    # erosion zeroing, the rec[pixels[mask]] mean, and the startmean mean. Compute
+    # it once here instead of re-evaluating mask.any()/mask.all() each scale.
+    mask_all_true = bool(mask.all())
+    startmean = (
+        pixels.mean()
+        if mask_all_true
+        else (numpy.mean(pixels[mask]) if mask.any() else 0.0)
+    )
     ero = pixels.copy()
     # Mask the test image so masked pixels have no effect during reconstruction
     ero[~mask] = 0
@@ -364,19 +381,49 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
             f"Spectrum length: {granular_spectrum_length}",
         )
 
+    # Precompute the upsample coordinate grids once before the spectrum loop.
+    # They depend only on the fixed subsampled/original shapes (not on the
+    # per-scale ``rec``), so rebuilding three full-resolution float arrays via
+    # ``numpy.mgrid`` on every scale (as ``_upsample_3d`` does internally) is
+    # wasted work. Reuse the same coordinate tuple for every
+    # ``map_coordinates`` call below. Only needed when subsampling is active
+    # and there are objects to measure.
+    upsample_coords: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray] | None = None
+    if subsample_size < 1.0 and nobjects > 0:
+        k, i, j = numpy.mgrid[
+            0 : original_shape[0],
+            0 : original_shape[1],
+            0 : original_shape[2],
+        ].astype(float)
+        if original_shape[0] > 1:
+            k *= float(new_shape[0] - 1) / float(original_shape[0] - 1)
+        if original_shape[1] > 1:
+            i *= float(new_shape[1] - 1) / float(original_shape[1] - 1)
+        if original_shape[2] > 1:
+            j *= float(new_shape[2] - 1) / float(original_shape[2] - 1)
+        upsample_coords = (k, i, j)
+
     for scale in range(1, granular_spectrum_length + 1):
         prevmean = currentmean
 
-        # Masked erosion
-        ero_masked = numpy.zeros_like(ero)
-        ero_masked[mask] = ero[mask]
-        ero = skimage.morphology.erosion(ero_masked, footprint=footprint)
+        # Masked erosion: zero pixels outside the mask before eroding. When the
+        # mask is all-True this is a no-op (ero is already zeroed at ~mask from
+        # the prior iteration), so skip the full-array ``numpy.where`` copy and
+        # erode ero directly.
+        ero_marker = ero if mask_all_true else numpy.where(mask, ero, 0)
+        ero = skimage.morphology.erosion(ero_marker, footprint=footprint)
 
         # Reconstruction
         rec = skimage.morphology.reconstruction(ero, pixels, footprint=footprint)
 
-        # Image-level granularity
-        currentmean = numpy.mean(rec[mask]) if mask.any() else 0.0
+        # Image-level granularity. When the mask is all-True, rec[mask] is all
+        # of rec, so rec.mean() avoids the boolean-index copy that rec[mask]
+        # would perform.
+        currentmean = (
+            rec.mean()
+            if mask_all_true
+            else (numpy.mean(rec[mask]) if mask.any() else 0.0)
+        )
         gs = (prevmean - currentmean) * 100 / startmean if startmean > 0 else 0.0
 
         if verbose and scale == 1:
@@ -387,11 +434,11 @@ def compute_granularity(  # noqa: C901, PLR0912, PLR0913, PLR0915
         # then compute per-label means using masked_labels.
         # ----------------------------------------------------------
         if nobjects > 0:
-            if subsample_size < 1.0:
-                rec_full = _upsample_3d(
+            if upsample_coords is not None:
+                rec_full = scipy.ndimage.map_coordinates(
                     rec,
-                    subsampled_shape=new_shape,
-                    original_shape=original_shape,
+                    upsample_coords,
+                    order=1,
                 )
             else:
                 rec_full = rec