diff --git a/docs/src/benchmarking.md b/docs/src/benchmarking.md new file mode 100644 index 0000000..f9ab22a --- /dev/null +++ b/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. diff --git a/docs/src/index.md b/docs/src/index.md index 0ce86e8..63a17aa 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -49,4 +49,5 @@ api :maxdepth: 2 :caption: A Note on Scalability scalability +benchmarking ``` diff --git a/src/zedprofiler/contracts.py b/src/zedprofiler/contracts.py index c4bfaf8..b84fd9d 100644 --- a/src/zedprofiler/contracts.py +++ b/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 diff --git a/src/zedprofiler/featurization/colocalization.py b/src/zedprofiler/featurization/colocalization.py index 02f1bcb..138083d 100644 --- a/src/zedprofiler/featurization/colocalization.py +++ b/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, diff --git a/src/zedprofiler/featurization/granularity.py b/src/zedprofiler/featurization/granularity.py index 72c1f78..efd9840 100644 --- a/src/zedprofiler/featurization/granularity.py +++ b/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 diff --git a/src/zedprofiler/featurization/intensity.py b/src/zedprofiler/featurization/intensity.py index d153e2f..eaf4ed6 100644 --- a/src/zedprofiler/featurization/intensity.py +++ b/src/zedprofiler/featurization/intensity.py @@ -8,6 +8,7 @@ import numpy import pandas import scipy.ndimage +import skimage.measure import skimage.segmentation from zedprofiler.contracts import validate_column_name_schema @@ -35,7 +36,7 @@ def get_outline(mask: numpy.ndarray) -> numpy.ndarray: return outline -def compute_intensity( # noqa: PLR0915 +def compute_intensity( # noqa: C901, PLR0915 object_loader: ObjectLoader, ) -> pandas.DataFrame: """Measure the intensity of objects in a 3D image. @@ -65,58 +66,82 @@ def compute_intensity( # noqa: PLR0915 "compartment": [], "value": [], } - # loop through each object and calculate measurements - for index, label in enumerate(labels): - selected_label_object = label_object.copy() - selected_image_object = image_object.copy() - selected_label_object[selected_label_object != label] = 0 - selected_label_object[selected_label_object > 0] = ( - 1 # binarize the label for volume calcs + props = skimage.measure.regionprops_table( + label_object, + 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]), ) - selected_image_object[selected_label_object != 1] = 0 - non_zero_pixels_object = selected_image_object[selected_image_object > 0] - if non_zero_pixels_object.size == 0: - non_zero_pixels_object = numpy.array([0], dtype=numpy.float32) - mask_outlines = get_outline(selected_label_object) - - # Extract only coordinates where object exists - z_indices, y_indices, x_indices = numpy.where(selected_label_object > 0) - bbox_min_z, bbox_max_z = numpy.min(z_indices), numpy.max(z_indices) - bbox_min_y, bbox_max_y = numpy.min(y_indices), numpy.max(y_indices) - bbox_min_x, bbox_max_x = numpy.min(x_indices), numpy.max(x_indices) - - # Crop to bounding box for efficiency - cropped_label = selected_label_object[ - bbox_min_z : bbox_max_z + 1, - bbox_min_y : bbox_max_y + 1, - bbox_min_x : bbox_max_x + 1, + for index, label in enumerate(props.get("label", [])) + } + + # loop through each object and calculate measurements + for label in labels: + bbox = label_to_bbox.get(int(label)) + if bbox is None: + continue + bbox_min_z, bbox_min_y, bbox_min_x, bbox_max_z, bbox_max_y, bbox_max_x = bbox + # regionprops bbox max coords are exclusive (half-open [min, max)), + # so the slices below need no +1. (The old code computed the bbox with + # numpy.max, which is an inclusive last index and required +1; the two + # produce the identical cropped region.) + cropped_label_values = label_object[ + bbox_min_z:bbox_max_z, + bbox_min_y:bbox_max_y, + bbox_min_x:bbox_max_x, ] - cropped_image = selected_image_object[ - bbox_min_z : bbox_max_z + 1, - bbox_min_y : bbox_max_y + 1, - bbox_min_x : bbox_max_x + 1, + cropped_image_values = image_object[ + bbox_min_z:bbox_max_z, + bbox_min_y:bbox_max_y, + bbox_min_x:bbox_max_x, ] + # The bbox always bounds at least one voxel of ``label`` because it was + # derived from regionprops for that exact label, so ``object_mask`` is + # never empty here. Kept as a defensive guard; excluded from coverage. + object_mask = cropped_label_values == label + if not numpy.any(object_mask): # pragma: no cover + continue # pragma: no cover - # Create coordinate grids for the bounding box + object_pixels = cropped_image_values[object_mask] + non_zero_pixels_object = object_pixels[object_pixels > 0] + if non_zero_pixels_object.size == 0: + non_zero_pixels_object = numpy.array([0], dtype=numpy.float32) + + cropped_label = object_mask.astype(numpy.uint8) + cropped_image = numpy.where(object_mask, cropped_image_values, 0) + + padded_label = numpy.pad(cropped_label, pad_width=1, mode="constant") + mask_outlines = get_outline(padded_label)[1:-1, 1:-1, 1:-1] + + # Create coordinate grids for the bounding box. Same exclusive + # bbox convention as the crops above: regionprops max coords are + # half-open [min, max), so no +1 is needed here either. mesh_z, mesh_y, mesh_x = numpy.mgrid[ - bbox_min_z : bbox_max_z + 1, # + 1 to include the max index - bbox_min_y : bbox_max_y + 1, - bbox_min_x : bbox_max_x + 1, + bbox_min_z:bbox_max_z, + bbox_min_y:bbox_max_y, + bbox_min_x:bbox_max_x, ] # calculate the integrated intensity - integrated_intensity = scipy.ndimage.sum( - selected_image_object, - selected_label_object, - index=1, - ) + integrated_intensity = numpy.sum(object_pixels) # calculate the volume - volume = numpy.sum(selected_label_object) + volume = numpy.sum(object_mask) - # Skip if volume is zero to avoid division by zero - if volume == 0: - continue + # Skip if volume is zero to avoid division by zero. Unreachable in + # practice because ``object_mask`` (label voxels within the bbox) + # always has at least one True voxel — the bbox is derived from + # regionprops for this exact label. Kept as a defensive guard; + # excluded from coverage. + if volume == 0: # pragma: no cover + continue # pragma: no cover # calculate the mean intensity mean_intensity = integrated_intensity / volume @@ -133,12 +158,15 @@ def compute_intensity( # noqa: PLR0915 # median intensity median_intensity = numpy.median(non_zero_pixels_object) # location of maximum intensity pixel (z, y, x) - max_intensity_z, max_intensity_y, max_intensity_x = ( - scipy.ndimage.maximum_position(selected_image_object) + max_position = numpy.unravel_index( + numpy.argmax(cropped_image), + cropped_image.shape, ) + max_intensity_z = bbox_min_z + max_position[0] + max_intensity_y = bbox_min_y + max_position[1] + max_intensity_x = bbox_min_x + max_position[2] # Calculate center of mass (geometric center) using cropped arrays - object_mask = cropped_label > 0 cm_x = numpy.mean(mesh_x[object_mask]) cm_y = numpy.mean(mesh_y[object_mask]) cm_z = numpy.mean(mesh_z[object_mask]) @@ -169,11 +197,12 @@ def compute_intensity( # noqa: PLR0915 # mean absolute deviation mad_intensity = numpy.mean(numpy.abs(non_zero_pixels_object - mean_intensity)) edge_count = scipy.ndimage.sum(mask_outlines) - integrated_intensity_edge = numpy.sum(selected_image_object[mask_outlines > 0]) + edge_pixels = cropped_image[mask_outlines > 0] + integrated_intensity_edge = numpy.sum(edge_pixels) mean_intensity_edge = integrated_intensity_edge / edge_count - std_intensity_edge = numpy.std(selected_image_object[mask_outlines > 0]) - min_intensity_edge = numpy.min(selected_image_object[mask_outlines > 0]) - max_intensity_edge = numpy.max(selected_image_object[mask_outlines > 0]) + std_intensity_edge = numpy.std(edge_pixels) + min_intensity_edge = numpy.min(edge_pixels) + max_intensity_edge = numpy.max(edge_pixels) measurements_dict = { "IntegratedIntensity": integrated_intensity, "MeanIntensity": mean_intensity, diff --git a/src/zedprofiler/featurization/neighbors.py b/src/zedprofiler/featurization/neighbors.py index 9ab04ca..c1741ba 100644 --- a/src/zedprofiler/featurization/neighbors.py +++ b/src/zedprofiler/featurization/neighbors.py @@ -122,24 +122,35 @@ def compute_neighbors( "NeighborsCountAdjacent": [], f"NeighborsCountByDistance-{distance_threshold}": [], } - for index, label in enumerate(labels): - props_label = skimage.measure.regionprops_table( - (label_object == label).astype(numpy.uint8), - properties=["bbox"], + props = skimage.measure.regionprops_table( + label_object, + 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 label in labels: + bbox_label = label_to_bbox.get(int(label)) + if bbox_label is None: + continue # get the number of neighbors for each object distance_x_y = distance_threshold distance_z = numpy.ceil(distance_threshold / anisotropy_factor).astype(int) # find how many other indexes are within a specified distance of the object # first expand the mask image by a specified distance - z_min, y_min, x_min, z_max, y_max, x_max = ( - props_label["bbox-0"][0], - props_label["bbox-1"][0], - props_label["bbox-2"][0], - props_label["bbox-3"][0], - props_label["bbox-4"][0], - props_label["bbox-5"][0], - ) + # regionprops bbox returns all min coords first then all max coords, + # in axis order (z, y, x): bbox-0/1/2 = min_z/y/x, bbox-3/4/5 = + # max_z/y/x. So this unpacks to each dimension's own min and max, + # all sourced from this label's bounding box. + z_min, y_min, x_min, z_max, y_max, x_max = bbox_label new_z_min, new_z_max = neighbors_expand_box( min_coor=image_global_min_coord_z, max_coord=image_global_max_coord_z, @@ -163,9 +174,40 @@ def compute_neighbors( ) bbox = (new_z_min, new_y_min, new_x_min, new_z_max, new_y_max, new_x_max) croppped_neighbor_image = crop_3D_image(image=label_object, bbox=bbox) - binary_mask = label_object == label + + adjacent_z_min, adjacent_z_max = neighbors_expand_box( + min_coor=image_global_min_coord_z, + max_coord=image_global_max_coord_z, + current_min=z_min, + current_max=z_max, + expand_by=1, + ) + adjacent_y_min, adjacent_y_max = neighbors_expand_box( + min_coor=image_global_min_coord_y, + max_coord=image_global_max_coord_y, + current_min=y_min, + current_max=y_max, + expand_by=1, + ) + adjacent_x_min, adjacent_x_max = neighbors_expand_box( + min_coor=image_global_min_coord_x, + max_coord=image_global_max_coord_x, + current_min=x_min, + current_max=x_max, + expand_by=1, + ) + adjacent_bbox = ( + adjacent_z_min, + adjacent_y_min, + adjacent_x_min, + adjacent_z_max, + adjacent_y_max, + adjacent_x_max, + ) + adjacent_label_crop = crop_3D_image(image=label_object, bbox=adjacent_bbox) + binary_mask = adjacent_label_crop == label dilated_mask = skimage.morphology.dilation(binary_mask) - labels_in_dilation = label_object[dilated_mask] + labels_in_dilation = adjacent_label_crop[dilated_mask] adjacent_labels = numpy.unique(labels_in_dilation) n_neighbors_adjacent = int( numpy.sum((adjacent_labels != 0) & (adjacent_labels != label)) diff --git a/src/zedprofiler/featurization/volumesizeshape.py b/src/zedprofiler/featurization/volumesizeshape.py index 123b5be..8799f80 100644 --- a/src/zedprofiler/featurization/volumesizeshape.py +++ b/src/zedprofiler/featurization/volumesizeshape.py @@ -103,6 +103,7 @@ def calculate_surface_area( label_object: np.ndarray, props: dict[str, np.ndarray], spacing: tuple[float, float, float], + label: int | None = None, ) -> float: """Calculate surface area for one labeled object using marching cubes.""" measure = _get_skimage_measure() @@ -112,7 +113,7 @@ def calculate_surface_area( max(props["bbox-1"][0], 0) : min(props["bbox-4"][0], label_object.shape[1]), max(props["bbox-2"][0], 0) : min(props["bbox-5"][0], label_object.shape[2]), ] - volume_truths = volume > 0 + volume_truths = volume == label if label is not None else volume > 0 verts, faces, _normals, _values = measure.marching_cubes( volume_truths, method="lewiner", @@ -136,6 +137,7 @@ def measure_3D_volume_size_shape( features_to_record = _empty_feature_result() desired_properties = [ + "label", "area", # for 3D it is volume but skimage uses "area" naming for the property "bbox", "centroid", @@ -144,20 +146,26 @@ def measure_3D_volume_size_shape( "euler_number", "equivalent_diameter", ] + all_props = measure.regionprops_table( + label_object, + properties=desired_properties, + ) + label_to_index = { + int(label): index for index, label in enumerate(all_props.get("label", [])) + } + for label in unique_objects: # avoid the 0 index which is the background and not an object, if label == 0: continue - subset_lab_object = label_object.copy() - # subset here means zeroing out all other objects except the - # one we want to measure, so that we can use - # skimage's regionprops_table to compute - # features for that object - subset_lab_object[subset_lab_object != label] = 0 - props = measure.regionprops_table( - subset_lab_object, - properties=desired_properties, - ) + props_index = label_to_index.get(int(label)) + if props_index is None: + continue + props = { + prop_name: np.asarray([prop_values[props_index]]) + for prop_name, prop_values in all_props.items() + if prop_name != "label" + } features_to_record["Metadata_Object_ObjectID"].append(label) features_to_record["Volume"].append(props["area"].item()) @@ -180,9 +188,10 @@ def measure_3D_volume_size_shape( try: features_to_record["SurfaceArea"].append( calculate_surface_area( - label_object=subset_lab_object, + label_object=label_object, props=props, spacing=spacing, + label=int(label), ), ) except (RuntimeError, ValueError): diff --git a/tests/benchmarking.py b/tests/benchmarking.py new file mode 100644 index 0000000..2f50261 --- /dev/null +++ b/tests/benchmarking.py @@ -0,0 +1,434 @@ +"""Shared deterministic feature cases for accuracy and performance tests.""" + +from __future__ import annotations + +import hashlib +import json +import math +import time +from collections.abc import Callable, Iterable +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import pandas as pd + +from zedprofiler.featurization.colocalization import compute_colocalization +from zedprofiler.featurization.granularity import compute_granularity +from zedprofiler.featurization.intensity import compute_intensity +from zedprofiler.featurization.neighbors import compute_neighbors +from zedprofiler.featurization.texture import compute_texture +from zedprofiler.featurization.volumesizeshape import compute_volume_size_shape +from zedprofiler.IO.loading_classes import ( + ImageSetConfig, + ImageSetLoader, + ObjectLoader, + TwoObjectLoader, +) + +try: + import tifffile +except ImportError: # pragma: no cover - exercised only without optional test dep + tifffile = None + + +CELLPROFILER_TUTORIAL_ROOT = ( + Path(__file__).resolve().parent + / "data" + / "CP_tutorial_3D_noise_nuclei_segmentation" +) + + +@dataclass +class BenchmarkImageSet: + """Minimal image-set loader surface used by feature functions.""" + + image_set_name: str = "benchmark-level" + anisotropy_spacing: tuple[float, float, float] = (10.0, 1.0, 1.0) + + +@dataclass +class BenchmarkObjectLoader: + """Minimal object loader for single-channel feature functions.""" + + image_set_loader: BenchmarkImageSet + image: np.ndarray + label_image: np.ndarray + object_ids: list[int] + compartment: str = "Nuclei" + channel: str = "DNA" + + +@dataclass +class BenchmarkTwoObjectLoader: + """Minimal object loader for paired-channel feature functions.""" + + image_set_loader: BenchmarkImageSet + image1: np.ndarray + image2: np.ndarray + label_image: np.ndarray + object_ids: list[int] + compartment: str = "Nuclei" + + +FeatureCase = tuple[str, Callable[[], pd.DataFrame]] + + +def make_benchmark_loaders() -> tuple[BenchmarkObjectLoader, BenchmarkTwoObjectLoader]: + """Create a deterministic two-object image set for benchmark contracts.""" + shape = (16, 32, 32) + z, y, x = np.indices(shape) + image1 = ((z * 11 + y * 7 + x * 5) % 251 + 1).astype(np.uint16) + image2 = ((image1.astype(np.uint32) * 3 + z * 13 + x * 17) % 251 + 1).astype( + np.uint16, + ) + + labels = np.zeros(shape, dtype=np.int32) + labels[2:7, 3:8, 4:9] = 1 + labels[9:14, 20:26, 18:24] = 2 + object_ids = [1, 2] + + image_set = BenchmarkImageSet() + object_loader = BenchmarkObjectLoader( + image_set_loader=image_set, + image=image1, + label_image=labels, + object_ids=object_ids, + ) + two_object_loader = BenchmarkTwoObjectLoader( + image_set_loader=image_set, + image1=image1, + image2=image2, + label_image=labels, + object_ids=object_ids, + ) + return object_loader, two_object_loader + + +def make_grid_benchmark_loaders( + *, + image_set_name: str = "benchmark-grid", + shape: tuple[int, int, int] = (32, 128, 128), + object_count: int = 32, + cube_size: int = 4, +) -> tuple[BenchmarkObjectLoader, BenchmarkTwoObjectLoader]: + """Create a deterministic many-object image set for scaling scorecards.""" + z, y, x = np.indices(shape) + image1 = ((z * 11 + y * 7 + x * 5) % 4093 + 1).astype(np.uint16) + image2 = ((image1.astype(np.uint32) * 3 + z * 13 + x * 17) % 4093 + 1).astype( + np.uint16, + ) + + labels = np.zeros(shape, dtype=np.int32) + object_ids: list[int] = [] + object_id = 1 + for z_start in range(1, shape[0] - cube_size, cube_size + 2): + for y_start in range(1, shape[1] - cube_size, cube_size + 4): + for x_start in range(1, shape[2] - cube_size, cube_size + 4): + labels[ + z_start : z_start + cube_size, + y_start : y_start + cube_size, + x_start : x_start + cube_size, + ] = object_id + object_ids.append(object_id) + object_id += 1 + if len(object_ids) >= object_count: + image_set = BenchmarkImageSet(image_set_name=image_set_name) + return ( + BenchmarkObjectLoader( + image_set_loader=image_set, + image=image1, + label_image=labels, + object_ids=object_ids, + ), + BenchmarkTwoObjectLoader( + image_set_loader=image_set, + image1=image1, + image2=image2, + label_image=labels, + object_ids=object_ids, + ), + ) + + raise ValueError( + f"Could not place {object_count} objects in shape {shape} " + f"with cube size {cube_size}.", + ) + + +def feature_cases() -> list[FeatureCase]: + """Return feature computations that should remain result-stable.""" + object_loader, two_object_loader = make_benchmark_loaders() + return [ + ("intensity", lambda: compute_intensity(object_loader)), + ( + "volume_size_shape", + lambda: compute_volume_size_shape( + image_set_loader=object_loader.image_set_loader, + object_loader=object_loader, + ), + ), + ( + "neighbors", + lambda: compute_neighbors( + object_loader=object_loader, + distance_threshold=10, + anisotropy_factor=10, + ), + ), + ("texture", lambda: compute_texture(object_loader, distance=1, grayscale=256)), + ( + "granularity", + lambda: compute_granularity( + object_loader, + # The production default is 16, but granularity cost scales + # linearly with spectrum length (one morphology pass per + # scale). Kept small here so the accuracy lock stays fast and + # its fingerprint stable; the scaling scorecard below uses the + # realistic default of 16 for representative timing. + granular_spectrum_length=3, + subsample_size=0.5, + image_sample_size=0.5, + radius=2, + ), + ), + ( + "colocalization", + lambda: compute_colocalization( + two_object_loader, + fast_costes="Faster", + channel1="DNA", + channel2="GFP", + ), + ), + ] + + +def scaling_feature_cases() -> list[FeatureCase]: + """Return a many-object benchmark level for opt-in scorecards.""" + object_loader, two_object_loader = make_grid_benchmark_loaders() + return [ + ("scaling_intensity", lambda: compute_intensity(object_loader)), + ( + "scaling_volume_size_shape", + lambda: compute_volume_size_shape( + image_set_loader=object_loader.image_set_loader, + object_loader=object_loader, + ), + ), + ( + "scaling_neighbors", + lambda: compute_neighbors( + object_loader=object_loader, + distance_threshold=10, + anisotropy_factor=10, + ), + ), + ( + "scaling_texture", + # CellProfiler's default Haralick distance is 3; use it here so + # the scorecard reflects representative cost. The anisotropy spacing + # is (10, 1, 1), so 3 voxels in z is 30 physical units while 3 in + # x/y is 3 — the GLCM distance is in voxel space and does not + # adjust for anisotropy. + lambda: compute_texture(object_loader, distance=3, grayscale=256), + ), + ( + "scaling_granularity", + # Use the production default of 16 so the scaling scorecard + # reflects representative granularity cost, which scales with + # spectrum length. (No locked fingerprint for scaling cases.) + lambda: compute_granularity( + object_loader, + granular_spectrum_length=16, + subsample_size=0.5, + image_sample_size=0.5, + radius=2, + ), + ), + ( + "scaling_colocalization", + lambda: compute_colocalization( + two_object_loader, + fast_costes="Faster", + channel1="DNA", + channel2="GFP", + ), + ), + ] + + +def _load_real_world_object_loader( + *, + image_name: str = "nuclei1_out_c00_dr90_image", +) -> ObjectLoader: + """Load a representative real-world image/mask pair for scorecards.""" + if tifffile is None: + raise ModuleNotFoundError("tifffile is required for real-world benchmarks.") + + image = tifffile.imread(CELLPROFILER_TUTORIAL_ROOT / "input" / f"{image_name}.tif") + label = tifffile.imread( + CELLPROFILER_TUTORIAL_ROOT + / "output" + / "masks" + / f"{image_name}SegmentationMask.tiff", + ) + image_set_loader = ImageSetLoader( + image_set_path=None, + label_set_path=None, + image_set_array=image, + label_set_array=label, + anisotropy_spacing=(1.0, 1.0, 1.0), + channel_mapping={ + "DNA": image_name, + "Nuclei": "SegmentationMask", + }, + config=ImageSetConfig( + image_set_name=image_name, + label_key_name=["Nuclei"], + raw_image_key_name=["DNA"], + ), + ) + return ObjectLoader( + image_set_loader=image_set_loader, + channel_name="DNA", + compartment_name="Nuclei", + ) + + +def _load_real_world_two_object_loader() -> TwoObjectLoader: + """Load a representative paired-channel real-world case for scorecards.""" + if tifffile is None: + raise ModuleNotFoundError("tifffile is required for real-world benchmarks.") + + first_image_name = "nuclei1_out_c00_dr90_image" + second_image_name = "nuclei2_out_c90_dr90_image" + label = tifffile.imread( + CELLPROFILER_TUTORIAL_ROOT + / "output" + / "masks" + / f"{first_image_name}SegmentationMask.tiff", + ) + object_ids = [int(x) for x in np.unique(label) if x != 0] + image_set_loader = ImageSetLoader.__new__(ImageSetLoader) + image_set_loader.image_set_name = "real-world-dr90-c00-c90" + image_set_loader.image_set_dict = { + "DNA1": tifffile.imread( + CELLPROFILER_TUTORIAL_ROOT / "input" / f"{first_image_name}.tif", + ), + "DNA2": tifffile.imread( + CELLPROFILER_TUTORIAL_ROOT / "input" / f"{second_image_name}.tif", + ), + "Nuclei": label, + } + image_set_loader.unique_compartment_objects = {"Nuclei": object_ids} + return TwoObjectLoader( + image_set_loader=image_set_loader, + compartment="Nuclei", + channel1="DNA1", + channel2="DNA2", + ) + + +def real_world_feature_cases() -> list[FeatureCase]: + """Return real-world benchmark cases backed by checked-in tutorial data.""" + object_loader = _load_real_world_object_loader() + two_object_loader = _load_real_world_two_object_loader() + return [ + ("real_world_intensity", lambda: compute_intensity(object_loader)), + ( + "real_world_volume_size_shape", + lambda: compute_volume_size_shape( + image_set_loader=object_loader.image_set_loader, + object_loader=object_loader, + ), + ), + ( + "real_world_neighbors", + lambda: compute_neighbors( + object_loader=object_loader, + distance_threshold=50, + anisotropy_factor=1, + ), + ), + ( + "real_world_texture", + lambda: compute_texture(object_loader, distance=1, grayscale=256), + ), + ( + "real_world_granularity", + lambda: compute_granularity( + object_loader, + radius=1, + granular_spectrum_length=2, + subsample_size=1.0, + image_sample_size=1.0, + ), + ), + ( + "real_world_colocalization", + lambda: compute_colocalization( + two_object_loader, + fast_costes="Faster", + channel1="DNA1", + channel2="DNA2", + ), + ), + ] + + +def _normalize_scalar(value: object, precision: int) -> object: + """Normalize dataframe scalar values for stable JSON fingerprints.""" + if pd.isna(value): + return "NaN" + if isinstance(value, np.generic): + value = value.item() + if isinstance(value, float): + if math.isinf(value): + return str(value) + rounded = round(value, precision) + return 0.0 if rounded == 0 else rounded + return value + + +def canonical_records(dataframe: pd.DataFrame, precision: int = 6) -> list[dict]: + """Return deterministic records independent of dataframe column order.""" + canonical = dataframe.copy() + if "Metadata_Object_ObjectID" in canonical.columns: + canonical = canonical.sort_values("Metadata_Object_ObjectID") + canonical = canonical.reindex(sorted(canonical.columns), axis=1).reset_index( + drop=True, + ) + return [ + {column: _normalize_scalar(value, precision) for column, value in row.items()} + for row in canonical.to_dict(orient="records") + ] + + +def dataframe_signature(dataframe: pd.DataFrame, precision: int = 6) -> str: + """Hash a dataframe after deterministic normalization.""" + payload = json.dumps( + canonical_records(dataframe, precision=precision), + sort_keys=True, + separators=(",", ":"), + ) + return hashlib.sha256(payload.encode()).hexdigest() + + +def time_feature_cases(cases: Iterable[FeatureCase]) -> list[dict[str, object]]: + """Run feature cases once and return a compact scorecard.""" + scorecard: list[dict[str, object]] = [] + for name, run_case in cases: + start = time.perf_counter() + dataframe = run_case() + elapsed = time.perf_counter() - start + scorecard.append( + { + "feature": name, + "seconds": round(elapsed, 6), + "rows": int(dataframe.shape[0]), + "columns": int(dataframe.shape[1]), + "signature": dataframe_signature(dataframe), + }, + ) + return scorecard diff --git a/tests/conftest.py b/tests/conftest.py index 349e378..9fbcb67 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -6,6 +6,8 @@ from __future__ import annotations +import os + import numpy as np import pytest @@ -13,6 +15,43 @@ from test_data_profiles import Profile +def pytest_addoption(parser: pytest.Parser) -> None: + """Add opt-in benchmark execution flag.""" + parser.addoption( + "--run-benchmarks", + action="store_true", + default=False, + help="Run opt-in benchmark scorecard tests.", + ) + + +def pytest_configure(config: pytest.Config) -> None: + """Register local test markers.""" + config.addinivalue_line( + "markers", + "benchmark: opt-in performance scorecard tests", + ) + + +def pytest_collection_modifyitems( + config: pytest.Config, + items: list[pytest.Item], +) -> None: + """Skip benchmark scorecards unless explicitly requested.""" + should_run = config.getoption("--run-benchmarks") or ( + os.environ.get("ZEDPROFILER_RUN_BENCHMARKS") == "1" + ) + if should_run: + return + + skip_benchmark = pytest.mark.skip( + reason="benchmark scorecards require --run-benchmarks", + ) + for item in items: + if "benchmark" in item.keywords: + item.add_marker(skip_benchmark) + + @pytest.fixture def my_data() -> str: """Provide a basic string fixture. diff --git a/tests/featurization/test_colocalization.py b/tests/featurization/test_colocalization.py index 7b34252..59a9f39 100644 --- a/tests/featurization/test_colocalization.py +++ b/tests/featurization/test_colocalization.py @@ -344,3 +344,33 @@ def test_calculate_colocalization_identical_images() -> None: # Manders should be non-negative assert results["MandersCoeffM1"] >= 0.0 assert results["MandersCoeffM2"] >= 0.0 + + +@pytest.mark.parametrize("shape,center", [((7, 7, 7), (3, 3, 3))]) +def test_compute_colocalization_skips_phantom_object_id_without_bbox( + shape: tuple[int, int, int], + center: tuple[int, int, int], +) -> None: + """Object ids absent from the label image are skipped via the bbox guard. + + ``compute_colocalization`` looks up each requested object id in the + ``regionprops`` bbox table and skips ids with no matching region via the + ``if bbox is None: continue`` guard. Requesting a phantom id (99) + alongside a real one (1) exercises that guard: only object 1 should appear + in the output, with no crash. + """ + imgset = ImageSetLoaderModel() + label, im1, im2 = make_pair(shape, center) + loader = TwoObjectLoaderModel( + image_set_loader=imgset, + compartment="Cell", + image1=im1, + image2=im2, + label_image=label, + object_ids=[1, 99], + ) + + df = compute_colocalization(loader, channel1="A", channel2="B") + + returned_ids = sorted(int(x) for x in df["Metadata_Object_ObjectID"].tolist()) + assert returned_ids == [1] diff --git a/tests/featurization/test_intensity.py b/tests/featurization/test_intensity.py index a2593ea..85916f6 100644 --- a/tests/featurization/test_intensity.py +++ b/tests/featurization/test_intensity.py @@ -152,3 +152,71 @@ def test_compute_intensity_basic( df = compute_intensity(loader) assert isinstance(df, pd.DataFrame) assert "Metadata_Object_ObjectID" in df.columns + + +def test_compute_intensity_skips_phantom_object_id_without_bbox() -> None: + """Object ids absent from the label image are skipped, not crashed on. + + ``compute_intensity`` looks up each requested object id in the + ``regionprops`` bbox table. When an id has no matching label region the + bbox lookup returns ``None`` and the object is skipped via the + ``if bbox is None: continue`` guard. This test exercises that guard by + requesting a phantom id (99) alongside a real one (1) and asserting only + the real object appears in the output. + """ + shape = (10, 10, 10) + image = np.zeros(shape, dtype=np.float32) + label = np.zeros(shape, dtype=np.int32) + image[3:7, 3:7, 3:7] = 100.0 + label[3:7, 3:7, 3:7] = 1 + + imgset = ImageSetLoaderModel() + loader = ObjectLoaderModel( + image=image, + label_image=label, + object_ids=[1, 99], + image_set_loader=imgset, + ) + + df = compute_intensity(loader) + + returned_ids = sorted(int(x) for x in df["Metadata_Object_ObjectID"].tolist()) + assert returned_ids == [1] + + +def test_compute_intensity_handles_all_zero_intensity_object() -> None: + """An object with no nonzero-intensity voxels yields NaN mass displacement. + + When every voxel of an object has zero intensity, the non-zero pixel + selection is empty and the ``if non_zero_pixels_object.size == 0`` guard + substitutes a single zero so downstream statistics are well-defined. With + zero integrated intensity the intensity-weighted center of mass (and thus + mass displacement) is undefined, so those features are reported as NaN + rather than a misleading 0. This test exercises the zero-pixel fallback + and the NaN mass-displacement branch together. + """ + shape = (10, 10, 10) + image = np.zeros(shape, dtype=np.float32) + label = np.zeros(shape, dtype=np.int32) + label[3:7, 3:7, 3:7] = 1 + + imgset = ImageSetLoaderModel() + loader = ObjectLoaderModel( + image=image, + label_image=label, + object_ids=[1], + image_set_loader=imgset, + ) + + df = compute_intensity(loader) + + row = df[df["Metadata_Object_ObjectID"] == 1] + assert len(row) == 1 + + mass_disp_col = next( + c for c in df.columns if c.endswith("Intensity_MassDisplacement") + ) + assert np.isnan(row[mass_disp_col].values[0]) + + cmi_col = next(c for c in df.columns if c.endswith("Intensity_CMI-X")) + assert np.isnan(row[cmi_col].values[0]) diff --git a/tests/featurization/test_neighbors.py b/tests/featurization/test_neighbors.py index 398d12f..859d539 100644 --- a/tests/featurization/test_neighbors.py +++ b/tests/featurization/test_neighbors.py @@ -139,6 +139,9 @@ def test_mahalanobis_small_and_regularized_and_singular() -> None: assert np.allclose(md_sing, 0.0) +EXPECTED_MULTIPLE_NEIGHBORS = 2 + + def test_neighbors_count_adjacent_detects_touching_cells() -> None: """NeighborsCountAdjacent must be > 0 for two directly touching objects (Bug 5). @@ -177,6 +180,39 @@ def test_neighbors_count_adjacent_detects_touching_cells() -> None: ) +def test_neighbors_count_adjacent_detects_multiple_neighbors() -> None: + """NeighborsCountAdjacent must count more than one neighbour. + + Object 1 is face-adjacent to two separate objects: object 2 along the z + axis and object 3 along the x axis. The 1-voxel dilation of object 1's + mask must reach both, so ``NeighborsCountAdjacent`` for object 1 is 2. + """ + shape = (12, 12, 12) + lab = np.zeros(shape, dtype=int) + # Object 1: z,y,x in {3,4,5}; touches object 2 at z=6 and object 3 at x=6. + lab[3:6, 3:6, 3:6] = 1 + lab[6:9, 3:6, 3:6] = 2 + lab[3:6, 3:6, 6:9] = 3 + + imgset = ImageSetLoaderModel() + loader = ObjectLoaderModel( + label_image=lab, + object_ids=[1, 2, 3], + image_set_loader=imgset, + ) + + df = compute_neighbors(loader, distance_threshold=5, anisotropy_factor=1) + obj1_row = df[df["Metadata_Object_ObjectID"] == 1] + adjacent_col = [c for c in df.columns if "NeighborsCountAdjacent" in c] + assert adjacent_col, "NeighborsCountAdjacent column not found in output" + + n_adj = int(obj1_row[adjacent_col[0]].values[0]) + assert n_adj == EXPECTED_MULTIPLE_NEIGHBORS, ( + f"Expected object 1 to have {EXPECTED_MULTIPLE_NEIGHBORS} adjacent " + f"neighbours, got {n_adj}" + ) + + def test_create_results_dataframe_and_errors_and_plots() -> None: # create a simple classification results dict results = { @@ -203,3 +239,29 @@ def test_create_results_dataframe_and_errors_and_plots() -> None: fig2 = plot_distance_distributions(results, n_shells=2) assert hasattr(fig2, "axes") + + +def test_compute_neighbors_skips_phantom_object_id_without_bbox() -> None: + """Object ids absent from the label image are skipped via the bbox guard. + + ``compute_neighbors`` looks up each requested object id in the + ``regionprops`` bbox table and skips ids with no matching region via the + ``if bbox_label is None: continue`` guard. Requesting a phantom id (99) + alongside a real one (1) exercises that guard: only object 1 should appear + in the output, with no crash and a finite neighbor count for the real + object. + """ + shape = (10, 10, 10) + centers = [(3, 3, 3), (6, 6, 6)] + lab = make_two_labels(shape, centers) + imgset = ImageSetLoaderModel() + loader = ObjectLoaderModel( + label_image=lab, + object_ids=[1, 99], + image_set_loader=imgset, + ) + + df = compute_neighbors(loader, distance_threshold=5, anisotropy_factor=1) + + returned_ids = sorted(int(x) for x in df["Metadata_Object_ObjectID"].tolist()) + assert returned_ids == [1] diff --git a/tests/featurization/test_real_world_data.py b/tests/featurization/test_real_world_data.py index 8625e65..5a4da47 100644 --- a/tests/featurization/test_real_world_data.py +++ b/tests/featurization/test_real_world_data.py @@ -111,6 +111,7 @@ class LoadedNucleiCase: @property def object_ids(self) -> list[int]: + """The segmented object ids loaded for this image case.""" return self.object_loader.object_ids @@ -124,6 +125,7 @@ class FeatureRunner: def _cellprofiler_tutorial_image_case(image_name: str) -> RealImageCase: + """Build a RealImageCase pointing at a named tutorial image and its mask.""" return RealImageCase( image_name=image_name, image_path=CELLPROFILER_TUTORIAL_ROOT / "input" / f"{image_name}.tif", @@ -188,14 +190,17 @@ def _cellprofiler_tutorial_image_case(image_name: str) -> RealImageCase: def _feature_runner_id(runner: FeatureRunner) -> str: + """pytest id for a FeatureRunner parametrization.""" return runner.name def _dataset_case_id(dataset_case: RealDatasetCase) -> str: + """pytest id for a RealDatasetCase parametrization.""" return dataset_case.name def _image_case_id(case: tuple[RealDatasetCase, RealImageCase]) -> str: + """pytest id for a (dataset, image) parametrization pair.""" dataset_case, image_case = case return f"{dataset_case.name}-{image_case.image_name}" @@ -203,6 +208,7 @@ def _image_case_id(case: tuple[RealDatasetCase, RealImageCase]) -> str: def _colocalization_case_id( case: tuple[RealDatasetCase, RealColocalizationCase], ) -> str: + """pytest id for a (dataset, colocalization) parametrization pair.""" dataset_case, colocalization_case = case return f"{dataset_case.name}-{colocalization_case.name}" @@ -211,6 +217,7 @@ def _load_nuclei_case( dataset_case: RealDatasetCase, image_case: RealImageCase, ) -> LoadedNucleiCase: + """Load a real image/mask pair into loaders from in-memory arrays.""" image = tifffile.imread(image_case.image_path) label = tifffile.imread(image_case.label_path) image_set_loader = ImageSetLoader( @@ -247,6 +254,7 @@ def _load_nuclei_case_from_paths( dataset_case: RealDatasetCase, image_case: RealImageCase, ) -> LoadedNucleiCase: + """Load a real image/mask pair into loaders from on-disk file paths.""" image_set_loader = ImageSetLoader( image_set_path=image_case.image_path.parent, label_set_path=image_case.label_path.parent, @@ -278,6 +286,7 @@ def _load_nuclei_case_from_paths( def _load_colocalization_case( colocalization_case: RealColocalizationCase, ) -> TwoObjectLoader: + """Pair two real single-channel images into a TwoObjectLoader for colocalization.""" label = tifffile.imread(colocalization_case.label_image_case.label_path) object_ids = [int(x) for x in np.unique(label) if x != 0] image_set_loader = ImageSetLoader.__new__(ImageSetLoader) @@ -304,6 +313,7 @@ def _load_colocalization_case( def _expected_volumes_from_label(label: np.ndarray) -> dict[int, int]: + """Ground-truth object id -> voxel count, computed directly from the label mask.""" return { int(object_id): int(np.count_nonzero(label == object_id)) for object_id in np.unique(label) @@ -316,6 +326,7 @@ def _assert_real_feature_frame_matches_objects( loaded_case: LoadedNucleiCase, expected_column_token: str, ) -> None: + """Assert a feature frame has expected objects, columns, and finite values.""" assert isinstance(df, pd.DataFrame) assert df.shape[0] == loaded_case.dataset_case.expected_object_count assert "Metadata_Object_ObjectID" in df.columns diff --git a/tests/featurization/test_volumesizeshape.py b/tests/featurization/test_volumesizeshape.py index 02637f0..5061199 100644 --- a/tests/featurization/test_volumesizeshape.py +++ b/tests/featurization/test_volumesizeshape.py @@ -75,3 +75,28 @@ def test_compute_volume_size_shape_returns_dataframe( # All object ids present returned_ids = sorted(int(x) for x in df["Metadata_Object_ObjectID"].tolist()) assert returned_ids == obj_ids + + +def test_compute_volume_size_shape_skips_phantom_object_id_without_props() -> None: + """Object ids absent from the label image are skipped via the props guard. + + ``measure_3D_volume_size_shape`` builds a label->index map from + ``regionprops_table`` and skips any requested id with no matching region + via the ``if props_index is None: continue`` guard. Requesting a phantom + id (99) alongside a real one (1) exercises that guard: only object 1 + should appear in the output, with no crash. + """ + imgset = ImageSetLoaderModel(anisotropy_spacing=(1.0, 1.0, 1.0)) + label = make_label_image((7, 7, 7), [(3, 3, 3)]) + loader = ObjectLoaderModel( + label_image=label, + object_ids=[1, 99], + image_set_loader=imgset, + compartment="Nucleus", + channel="DAPI", + ) + + df = compute_volume_size_shape(image_set_loader=imgset, object_loader=loader) + + returned_ids = sorted(int(x) for x in df["Metadata_Object_ObjectID"].tolist()) + assert returned_ids == [1] diff --git a/tests/test_benchmark_contracts.py b/tests/test_benchmark_contracts.py new file mode 100644 index 0000000..bc06b0b --- /dev/null +++ b/tests/test_benchmark_contracts.py @@ -0,0 +1,70 @@ +"""Accuracy locks and opt-in benchmark scorecards for feature extraction.""" + +from __future__ import annotations + +import json +from collections.abc import Callable + +import pandas as pd +import pytest +from benchmarking import ( + dataframe_signature, + feature_cases, + real_world_feature_cases, + scaling_feature_cases, + time_feature_cases, +) + +EXPECTED_SIGNATURES = { + "intensity": ("351f6508dcfc0978c8d5bfc1891847cfe0090c42d8243baf0acb0f528f3061c0"), + "volume_size_shape": ( + "1fc7482eb490256eca01cc6d54d4b96956d132cdf0a6a27a934eba37dbf83f39" + ), + "neighbors": ("8f2b18e6023d656ec6fad41ffa0ff9b802f7b3eadb98ea315ea11cfa400644ec"), + "texture": ("e6ae19f7b6bc9e635fb6e199dd45452fca1a17246e027dcb8385478a95e913fa"), + "granularity": ("b46cd8ae17d0d1e8b47d24821c480975dd0464159744bc3913e619616dc94295"), + "colocalization": ( + "8bf9495cabc617a5743614f8d57c3855283f6cc6aff9715d975fef416bc29d4a" + ), +} +EXPECTED_OBJECT_ROWS = 2 +EXPECTED_SCALING_ROWS = 32 +EXPECTED_REAL_WORLD_ROWS = 5 + + +@pytest.mark.parametrize(("feature_name", "run_case"), feature_cases()) +def test_feature_outputs_match_current_accuracy_lock( + feature_name: str, + run_case: Callable[[], pd.DataFrame], +) -> None: + """Feature refactors must preserve current deterministic outputs.""" + dataframe = run_case() + assert dataframe.shape[0] == EXPECTED_OBJECT_ROWS + assert "Metadata_Object_ObjectID" in dataframe.columns + assert dataframe_signature(dataframe) == EXPECTED_SIGNATURES[feature_name] + + +@pytest.mark.benchmark +def test_feature_benchmark_scorecard() -> None: + """Print an opt-in scorecard for comparing performance passes.""" + scorecard = time_feature_cases( + [*feature_cases(), *scaling_feature_cases(), *real_world_feature_cases()], + ) + print("\nZedProfiler feature benchmark scorecard") + print(json.dumps(scorecard, indent=2, sort_keys=True)) + + observed_features = {record["feature"] for record in scorecard} + assert set(EXPECTED_SIGNATURES).issubset(observed_features) + for record in scorecard: + expected_rows = ( + EXPECTED_SCALING_ROWS + if str(record["feature"]).startswith("scaling_") + else EXPECTED_OBJECT_ROWS + ) + if str(record["feature"]).startswith("real_world_"): + expected_rows = EXPECTED_REAL_WORLD_ROWS + assert record["rows"] == expected_rows + assert record["columns"] > 0 + assert record["seconds"] >= 0 + if record["feature"] in EXPECTED_SIGNATURES: + assert record["signature"] == EXPECTED_SIGNATURES[record["feature"]]