From b647d4344a0c96418356f318f62e726e189fd967 Mon Sep 17 00:00:00 2001 From: Tristan Simas Date: Sun, 4 Oct 2026 09:01:23 -0400 Subject: [PATCH] Measure 3D shape on prepared label cohorts Refs #624. Derive Lewiner geometry during existing family preparation and share the complete feature mapping between cohort measurements and the numerical reference path. --- .../processing/backends/cellprofiler/shape.py | 608 +++++++++++------- tests/unit/test_cellprofiler_shape_hotpath.py | 97 ++- 2 files changed, 475 insertions(+), 230 deletions(-) diff --git a/openhcs/processing/backends/cellprofiler/shape.py b/openhcs/processing/backends/cellprofiler/shape.py index 46a797c83d..d810c8bf50 100644 --- a/openhcs/processing/backends/cellprofiler/shape.py +++ b/openhcs/processing/backends/cellprofiler/shape.py @@ -1,7 +1,7 @@ """Shape-measurement backends for CellProfiler-compatible processing.""" from __future__ import annotations -from typing import Annotated, TYPE_CHECKING, TypeAlias +from typing import Annotated, ClassVar, TYPE_CHECKING, TypeAlias from openhcs.interop.cellprofiler.settings_binder import ( SettingToKeywordBinding, parse_cellprofiler_bool, @@ -603,45 +603,6 @@ def row_mappings(self) -> tuple[Mapping[str, object], ...]: return self._rows -@dataclass(frozen=True, slots=True) -class SurfaceArea3DRegions: - """Object labels and bounded 3-D regions for surface-area measurement.""" - - label_ids: np.ndarray - bounds_zyxzyx: np.ndarray - - def __post_init__(self) -> None: - label_ids = np.ascontiguousarray(self.label_ids, dtype=np.int64) - bounds = np.ascontiguousarray(self.bounds_zyxzyx, dtype=np.int64) - if bounds.shape != (label_ids.size, 6): - raise ValueError( - f"3-D surface-area bounds must have shape ({label_ids.size}, 6), got {bounds.shape!r}." - ) - object.__setattr__(self, "label_ids", label_ids) - object.__setattr__(self, "bounds_zyxzyx", bounds) - - @classmethod - def from_label_array( - cls, labels: np.ndarray, label_ids: np.ndarray - ) -> "SurfaceArea3DRegions": - label_id_array = np.asarray(label_ids, dtype=np.int64) - if label_id_array.size == 0: - return cls(label_id_array, np.zeros((0, 6), dtype=np.int64)) - label_array = np.asarray(labels) - bounds = np.zeros((label_id_array.size, 6), dtype=np.int64) - for index, label_id in enumerate(label_id_array): - positions = np.argwhere(label_array == int(label_id)) - if positions.size == 0: - continue - minimum = positions.min(axis=0) - maximum = positions.max(axis=0) + 1 - bounds[index, :3] = minimum - bounds[index, 3:] = maximum - return cls( - label_id_array, _expanded_surface_area_bounds(bounds, label_array.shape) - ) - - @dataclass(frozen=True, slots=True) class ObjectSizeShapeFeatureArrayOwner(ABC, metaclass=AutoRegisterMeta): """Shared AreaShape feature-array invocation policy for backend owners.""" @@ -908,8 +869,6 @@ def _feature_arrays_2d(self, labels: np.ndarray) -> ShapeFeatureArrays: return (features, measured_labels) def _feature_arrays_3d(self, labels: np.ndarray) -> ShapeFeatureArrays: - import skimage.measure - total_started_at = time.perf_counter() capture_array_fixture( "measure_object_size_shape_3d_input", @@ -917,140 +876,21 @@ def _feature_arrays_3d(self, labels: np.ndarray) -> ShapeFeatureArrays: calculate_advanced=np.asarray(self.calculate_advanced), calculate_zernikes=np.asarray(self.calculate_zernikes), ) - phase_started_at = time.perf_counter() shape_backend = ShapeMeasurementBackendStrategy.for_memory_type( backend_provider=self.shape_backend_provider ) - runtime_profiler.log( - "moss_backend_resolution_3d", - time.perf_counter() - phase_started_at, - function="measure_object_size_shape", - ) - phase_started_at = time.perf_counter() - regions = skimage.measure.regionprops(labels, cache=True) - object_count = len(regions) - measured_labels = np.empty(object_count, dtype=np.int64) - area = np.empty(object_count, dtype=np.float64) - centroid = np.empty((object_count, 3), dtype=np.float64) - bounds = np.empty((object_count, 6), dtype=np.int64) - bounding_box_volume = np.empty(object_count, dtype=np.float64) - inertia_tensor_eigenvalues = np.empty((object_count, 3), dtype=np.float64) - extent = np.empty(object_count, dtype=np.float64) - equivalent_diameter = np.empty(object_count, dtype=np.float64) - euler_number = np.empty(object_count, dtype=np.int64) - solidity = ( - np.empty(object_count, dtype=np.float64) - if self.calculate_advanced - else None - ) - for object_index, region in enumerate(regions): - measured_labels[object_index] = region.label - area[object_index] = region.area - centroid[object_index] = region.centroid - bounds[object_index] = region.bbox - bounding_box_volume[object_index] = region.area_bbox - inertia_tensor_eigenvalues[object_index] = region.inertia_tensor_eigvals - extent[object_index] = region.extent - equivalent_diameter[object_index] = region.equivalent_diameter_area - euler_number[object_index] = region.euler_number - if solidity is not None: - solidity[object_index] = region.solidity - runtime_profiler.log( - "moss_regionprops_table_3d", - time.perf_counter() - phase_started_at, - function="measure_object_size_shape", - ) - surface_regions = SurfaceArea3DRegions( - measured_labels, - _expanded_surface_area_bounds(bounds, labels.shape), - ) - phase_started_at = time.perf_counter() - major_axis_length, minor_axis_length = _cellprofiler_3d_axis_lengths( - inertia_tensor_eigenvalues - ) - runtime_profiler.log( - "moss_axis_lengths_3d", - time.perf_counter() - phase_started_at, - function="measure_object_size_shape", - objects=int(measured_labels.size), - ) - phase_started_at = time.perf_counter() - surface_areas = shape_backend.surface_areas_3d( + features, measured_labels = shape_backend.feature_arrays_3d( labels, - surface_regions, + calculate_advanced=self.calculate_advanced, spacing=self.feature_source_voxel_spacing.spacing_for_ndim(labels.ndim), ) - runtime_profiler.log( - "moss_surface_areas_3d", - time.perf_counter() - phase_started_at, - function="measure_object_size_shape", - objects=int(measured_labels.size), - ) - features = { - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.VOLUME - ): area, - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.SURFACE_AREA - ): surface_areas, - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.MAJOR_AXIS_LENGTH - ): major_axis_length, - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.MINOR_AXIS_LENGTH - ): minor_axis_length, - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.CENTER_X - ): centroid[:, 2], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.CENTER_Y - ): centroid[:, 1], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.CENTER_Z - ): centroid[:, 0], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_VOLUME - ): bounding_box_volume, - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MINIMUM_X - ): bounds[:, 2], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MAXIMUM_X - ): bounds[:, 5], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MINIMUM_Y - ): bounds[:, 1], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MAXIMUM_Y - ): bounds[:, 4], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MINIMUM_Z - ): bounds[:, 0], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MAXIMUM_Z - ): bounds[:, 3], - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.EXTENT - ): extent, - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.EULER_NUMBER - ): euler_number, - _shape_feature( - MeasureObjectSizeShapeModule.MeasurementFeature.EQUIVALENT_DIAMETER - ): equivalent_diameter, - } - if self.calculate_advanced: - assert solidity is not None - features[ - _shape_feature(MeasureObjectSizeShapeModule.MeasurementFeature.SOLIDITY) - ] = solidity runtime_profiler.log( "moss_features_3d_total", time.perf_counter() - total_started_at, function="measure_object_size_shape", objects=int(measured_labels.size), ) - return (features, measured_labels) + return features, measured_labels def measure_object_size_shape_feature_arrays( @@ -1404,15 +1244,133 @@ def feret_diameters( ) -> tuple[np.ndarray, np.ndarray]: """Return minimum and maximum Feret diameters.""" - def surface_areas_3d( + _unit_surface_triangles: ClassVar[np.ndarray | None] = None + _cube_euler_coefficients: ClassVar[np.ndarray | None] = None + + def feature_arrays_3d( self, labels: np.ndarray, - regions: SurfaceArea3DRegions, *, - spacing: tuple[float, ...] | None = None, - ) -> np.ndarray: - """Return Lewiner marching-cubes surface areas for 3-D dense labels.""" - return _surface_areas_3d_from_regions(labels, regions, spacing=spacing) + calculate_advanced: bool, + spacing: tuple[float, ...], + ) -> ShapeFeatureArrays: + """Reference numerical owner for unprepared calls and oversized moments.""" + import skimage.measure + + regions = skimage.measure.regionprops(labels, cache=True) + count = len(regions) + measured_labels = np.empty(count, dtype=np.int64) + area = np.empty(count, dtype=np.float64) + centroid = np.empty((count, 3), dtype=np.float64) + bounds = np.empty((count, 6), dtype=np.int64) + eigenvalues = np.empty((count, 3), dtype=np.float64) + euler_number = np.empty(count, dtype=np.int64) + surface_areas = np.empty(count, dtype=np.float64) + solidity = np.empty(count, dtype=np.float64) if calculate_advanced else None + for index, region in enumerate(regions): + measured_labels[index] = region.label + area[index] = region.area + centroid[index] = region.centroid + bounds[index] = region.bbox + eigenvalues[index] = region.inertia_tensor_eigvals + euler_number[index] = region.euler_number + lower = np.maximum(bounds[index, :3] - 1, 0) + upper = np.minimum(bounds[index, 3:] + 1, labels.shape) + binary = ( + labels[tuple(slice(a, b) for a, b in zip(lower, upper))] == region.label + ) + surface_areas[index] = _surface_area(binary, spacing) + if solidity is not None: + solidity[index] = region.solidity + return self._feature_arrays_from_3d_statistics( + measured_labels, + area, + centroid, + bounds, + eigenvalues, + euler_number, + surface_areas, + solidity, + ) + + @staticmethod + def _feature_arrays_from_3d_statistics( + measured_labels: np.ndarray, + area: np.ndarray, + centroid: np.ndarray, + bounds: np.ndarray, + inertia_tensor_eigenvalues: np.ndarray, + euler_number: np.ndarray, + surface_areas: np.ndarray, + solidity: np.ndarray | None, + ) -> ShapeFeatureArrays: + """Retain the single declared AreaShape mapping for both algorithms.""" + bounding_box_volume = np.prod(bounds[:, 3:] - bounds[:, :3], axis=1).astype( + np.float64 + ) + extent = area / bounding_box_volume + equivalent_diameter = (6.0 * area / np.pi) ** (1.0 / 3.0) + major_axis_length, minor_axis_length = _cellprofiler_3d_axis_lengths( + inertia_tensor_eigenvalues + ) + features = { + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.VOLUME + ): area, + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.SURFACE_AREA + ): surface_areas, + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.MAJOR_AXIS_LENGTH + ): major_axis_length, + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.MINOR_AXIS_LENGTH + ): minor_axis_length, + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.CENTER_X + ): centroid[:, 2], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.CENTER_Y + ): centroid[:, 1], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.CENTER_Z + ): centroid[:, 0], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_VOLUME + ): bounding_box_volume, + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MINIMUM_X + ): bounds[:, 2], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MAXIMUM_X + ): bounds[:, 5], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MINIMUM_Y + ): bounds[:, 1], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MAXIMUM_Y + ): bounds[:, 4], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MINIMUM_Z + ): bounds[:, 0], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.BOUNDING_BOX_MAXIMUM_Z + ): bounds[:, 3], + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.EXTENT + ): extent, + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.EULER_NUMBER + ): euler_number, + _shape_feature( + MeasureObjectSizeShapeModule.MeasurementFeature.EQUIVALENT_DIAMETER + ): equivalent_diameter, + } + if solidity is not None: + features[ + _shape_feature(MeasureObjectSizeShapeModule.MeasurementFeature.SOLIDITY) + ] = solidity + return features, measured_labels @abstractmethod def distance_to_edge(self, labels: np.ndarray) -> np.ndarray: @@ -1437,6 +1395,154 @@ def color_labels(self, labels: np.ndarray) -> np.ndarray: class NumbaShapeMeasurementMixin(ABC): """Shared Numba-backed shape leaves reused by concrete backend policies.""" + def prepare_3d_shape_measurements(self) -> None: + """Derive binary Lewiner geometry and warm both supported input layouts.""" + if ShapeMeasurementBackendStrategy._unit_surface_triangles is None: + import skimage.measure + from skimage.measure._regionprops_utils import EULER_COEFS3D_26 + + triangles = [] + euler_coefficients = np.empty(256, dtype=np.int64) + for case in range(256): + cube = np.asarray( + [(case >> corner) & 1 for corner in range(8)], + dtype=np.float32, + ).reshape((2, 2, 2)) + if case in (0, 255): + triangles.append(np.empty((0, 3, 3), dtype=np.float32)) + else: + vertices, faces, _normals, _values = skimage.measure.marching_cubes( + cube, + method="lewiner", + level=0, + ) + triangles.append(vertices[faces]) + # Match the existing Euler convolution's reversed Z/Y/X + # neighborhood and its X/Y bit ordering; retain its coefficients. + euler_case = sum( + int(value) << corner + for corner, value in enumerate( + cube[::-1, ::-1, ::-1].transpose((0, 2, 1)).ravel() + ) + ) + euler_coefficients[case] = EULER_COEFS3D_26[euler_case] + geometry = np.zeros( + (len(triangles), max(map(len, triangles)), 3, 3), dtype=np.float32 + ) + for case, values in enumerate(triangles): + geometry[case, : len(values)] = values + geometry.setflags(write=False) + euler_coefficients.setflags(write=False) + ShapeMeasurementBackendStrategy._cube_euler_coefficients = ( + euler_coefficients + ) + ShapeMeasurementBackendStrategy._unit_surface_triangles = geometry + + labels = np.zeros((3, 3, 3), dtype=np.int32) + labels[1, 1, 1] = 1 + for writable in (True, False): + labels.setflags(write=writable) + self.feature_arrays_3d( + labels, calculate_advanced=False, spacing=(1.0, 1.0, 1.0) + ) + + def feature_arrays_3d( + self, + labels: np.ndarray, + *, + calculate_advanced: bool, + spacing: tuple[float, ...], + ) -> ShapeFeatureArrays: + """Measure the full label cohort without per-object binary meshes.""" + geometry = ShapeMeasurementBackendStrategy._unit_surface_triangles + if geometry is None: + # Direct public calls retain reference behavior before library READY. + return super().feature_arrays_3d( + labels, calculate_advanced=calculate_advanced, spacing=spacing + ) + label_array = np.ascontiguousarray(labels, dtype=np.int32) + label_ids = np.unique(label_array) + label_ids = label_ids[label_ids > 0] + counts, bounds = _label_bounds_counts_3d_numba(label_array, label_ids) + if not self._local_moments_fit_int64(counts, bounds): + return super().feature_arrays_3d( + label_array, calculate_advanced=calculate_advanced, spacing=spacing + ) + coefficients = ShapeMeasurementBackendStrategy._cube_euler_coefficients + assert coefficients is not None + sums, products, surface_cases, euler_scaled = ( + _label_local_moments_surface_euler_3d_numba( + label_array, label_ids, bounds, coefficients + ) + ) + area = counts.astype(np.float64) + centroid = sums / area[:, None] + bounds[:, :3] + inertia = np.zeros((label_ids.size, 3, 3), dtype=np.float64) + for index, count in enumerate(counts): + # Integer central numerators avoid subtracting large float moments. + # Python integers also avoid overflow in the derived count products. + n = int(count) + local_sums = tuple(int(value) for value in sums[index]) + central = tuple( + (n * int(products[index, product]) - local_sums[a] * local_sums[b]) + / (n * n) + for product, (a, b) in enumerate( + ((0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)) + ) + ) + inertia[index, 0, 0] = central[1] + central[2] + inertia[index, 1, 1] = central[0] + central[2] + inertia[index, 2, 2] = central[0] + central[1] + inertia[index, 0, 1] = inertia[index, 1, 0] = -central[3] + inertia[index, 0, 2] = inertia[index, 2, 0] = -central[4] + inertia[index, 1, 2] = inertia[index, 2, 1] = -central[5] + eigenvalues = np.maximum(np.linalg.eigvalsh(inertia)[:, ::-1], 0.0) + if spacing == (1.0, 1.0, 1.0): + vertices = geometry + else: + vertices = geometry * np.asarray(spacing, dtype=np.float64) + edge_a = vertices[:, :, 0] - vertices[:, :, 1] + edge_b = vertices[:, :, 0] - vertices[:, :, 2] + case_areas = ( + np.sqrt((np.cross(edge_a, edge_b) ** 2).sum(axis=2)).sum(axis=1) / 2.0 + ) + surface_areas = surface_cases @ case_areas + solidity = None + if calculate_advanced: + from skimage.morphology import convex_hull_image + + convex_volumes = np.empty(label_ids.size, dtype=np.int64) + for index, label_id in enumerate(label_ids): + lower, upper = bounds[index, :3], bounds[index, 3:] + binary = ( + label_array[tuple(slice(a, b) for a, b in zip(lower, upper))] + == label_id + ) + convex_volumes[index] = np.count_nonzero(convex_hull_image(binary)) + solidity = area / convex_volumes + return self._feature_arrays_from_3d_statistics( + label_ids.astype(np.int64), + area, + centroid, + bounds, + eigenvalues, + euler_scaled // 8, + surface_areas, + solidity, + ) + + @staticmethod + def _local_moments_fit_int64(counts: np.ndarray, bounds: np.ndarray) -> bool: + """Check one conservative bound before accumulating local moments.""" + capacity = np.iinfo(np.int64).max + for count, region in zip(counts, bounds): + extent = max( + int(region[axis + 3]) - int(region[axis]) - 1 for axis in range(3) + ) + if int(count) * extent * extent > capacity: + return False + return True + def prepare_numba_shape_leaves(self) -> None: labels = np.array([[0, 1, 1], [0, 1, 0], [2, 2, 0]], dtype=np.int32) image = np.arange(9, dtype=np.float64).reshape((3, 3)) @@ -1453,6 +1559,7 @@ def prepare_numba_shape_leaves(self) -> None: immutable_ids.setflags(write=False) self.maximum_position_of_labels(ordering_image, labels, immutable_ids) self.color_labels(labels) + self.prepare_3d_shape_measurements() def form_factor_values( self, labels: np.ndarray, label_ids: np.ndarray @@ -1528,12 +1635,6 @@ class NumbaNumpyShapeMeasurementBackendStrategy( def prepare_backend(self) -> None: self.prepare_numba_shape_leaves() - label_ids = np.array([1, 2], dtype=np.int32) - labels_3d = np.zeros((3, 3, 3), dtype=np.int32) - labels_3d[0:2, 0:2, 0:2] = 1 - labels_3d[1:3, 1:3, 1:3] = 2 - regions_3d = SurfaceArea3DRegions.from_label_array(labels_3d, label_ids) - self.surface_areas_3d(labels_3d, regions_3d) def _distance_to_edge_planewise( @@ -1737,6 +1838,108 @@ def _zernike_features( } +@njit(cache=True) +def _label_bounds_counts_3d_numba( + labels: np.ndarray, label_ids: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """First cohort pass captures bounds for stable local moment coordinates.""" + nz, ny, nx = labels.shape + counts = np.zeros(label_ids.size, dtype=np.int64) + bounds = np.zeros((label_ids.size, 6), dtype=np.int64) + bounds[:, 0] = nz + bounds[:, 1] = ny + bounds[:, 2] = nx + for z in range(nz): + for y in range(ny): + for x in range(nx): + label = labels[z, y, x] + if label <= 0: + continue + index = np.searchsorted(label_ids, label) + counts[index] += 1 + bounds[index, 0] = min(bounds[index, 0], z) + bounds[index, 1] = min(bounds[index, 1], y) + bounds[index, 2] = min(bounds[index, 2], x) + bounds[index, 3] = max(bounds[index, 3], z + 1) + bounds[index, 4] = max(bounds[index, 4], y + 1) + bounds[index, 5] = max(bounds[index, 5], x + 1) + return counts, bounds + + +@njit(cache=True) +def _label_local_moments_surface_euler_3d_numba( + labels: np.ndarray, + label_ids: np.ndarray, + bounds: np.ndarray, + euler_coefficients: np.ndarray, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Second cohort pass collects local moments and correlated cube counts.""" + count = label_ids.size + sums = np.zeros((count, 3), dtype=np.int64) + products = np.zeros((count, 6), dtype=np.int64) + surface_cases = np.zeros((count, 256), dtype=np.int64) + euler_scaled = np.zeros(count, dtype=np.int64) + corners = np.zeros(8, dtype=np.int32) + nz, ny, nx = labels.shape + for z in range(-1, nz): + for y in range(-1, ny): + for x in range(-1, nx): + for corner in range(8): + zz = z + (corner >> 2) + yy = y + ((corner >> 1) & 1) + xx = x + (corner & 1) + if 0 <= zz < nz and 0 <= yy < ny and 0 <= xx < nx: + corners[corner] = labels[zz, yy, xx] + else: + corners[corner] = 0 + if corners[0] > 0: + index = np.searchsorted(label_ids, corners[0]) + local_z = z - bounds[index, 0] + local_y = y - bounds[index, 1] + local_x = x - bounds[index, 2] + sums[index, 0] += local_z + sums[index, 1] += local_y + sums[index, 2] += local_x + products[index, 0] += local_z * local_z + products[index, 1] += local_y * local_y + products[index, 2] += local_x * local_x + products[index, 3] += local_z * local_y + products[index, 4] += local_z * local_x + products[index, 5] += local_y * local_x + if ( + corners[0] + == corners[1] + == corners[2] + == corners[3] + == corners[4] + == corners[5] + == corners[6] + == corners[7] + ): + continue + for corner in range(8): + label = corners[corner] + if label <= 0: + continue + seen = False + for previous in range(corner): + if corners[previous] == label: + seen = True + if seen: + continue + case = 0 + for other in range(8): + if corners[other] == label: + case |= 1 << other + index = np.searchsorted(label_ids, label) + euler_scaled[index] += euler_coefficients[case] + # Euler includes an exterior zero border; marching cubes + # retains the original clipped domain, including open edges. + if 0 <= z < nz - 1 and 0 <= y < ny - 1 and 0 <= x < nx - 1: + surface_cases[index, case] += 1 + return sums, products, surface_cases, euler_scaled + + def _surface_area( volume: np.ndarray, spacing: tuple[float, ...] | None = None ) -> float: @@ -1757,59 +1960,6 @@ def _surface_area( return float(((np.cross(edge_a, edge_b) ** 2).sum(axis=1) ** 0.5).sum() / 2.0) -def _expanded_surface_area_bounds( - bounds_zyxzyx: np.ndarray, labels_shape: tuple[int, ...] -) -> np.ndarray: - bounds = np.ascontiguousarray(bounds_zyxzyx, dtype=np.int64) - if bounds.shape[1:] != (6,): - raise ValueError( - f"3-D surface-area bounds must have six columns (z0, y0, x0, z1, y1, x1), got {bounds.shape!r}." - ) - if len(labels_shape) != 3: - raise ValueError( - f"3-D surface-area labels require a 3-D shape, got {labels_shape!r}." - ) - expanded = bounds.copy() - shape = np.asarray(labels_shape, dtype=np.int64) - expanded[:, :3] = np.maximum(expanded[:, :3] - 1, 0) - expanded[:, 3:] = np.minimum(expanded[:, 3:] + 1, shape) - return np.ascontiguousarray(expanded, dtype=np.int64) - - -def _surface_areas_3d_from_labels( - labels: np.ndarray, - label_ids: np.ndarray, - *, - spacing: tuple[float, ...] | None = None, -) -> np.ndarray: - return _surface_areas_3d_from_regions( - labels, - SurfaceArea3DRegions.from_label_array(labels, label_ids), - spacing=spacing, - ) - - -def _surface_areas_3d_from_regions( - labels: np.ndarray, - regions: SurfaceArea3DRegions, - *, - spacing: tuple[float, ...] | None = None, -) -> np.ndarray: - if regions.label_ids.size == 0: - return np.zeros(0, dtype=np.float64) - label_array = np.asarray(labels) - if label_array.ndim != 3: - raise ValueError("3-D surface-area measurement requires a 3-D label array.") - surface_areas = np.zeros(regions.label_ids.size, dtype=np.float64) - for object_index, label_id in enumerate(regions.label_ids): - z0, y0, x0, z1, y1, x1 = ( - int(value) for value in regions.bounds_zyxzyx[object_index] - ) - volume = label_array[z0:z1, y0:y1, x0:x1] == int(label_id) - surface_areas[object_index] = _surface_area(volume, spacing=spacing) - return surface_areas - - def _form_factor_values_from_labels( labels: np.ndarray, label_ids: np.ndarray ) -> np.ndarray: diff --git a/tests/unit/test_cellprofiler_shape_hotpath.py b/tests/unit/test_cellprofiler_shape_hotpath.py index cff04c3dca..3a5a1bae7b 100644 --- a/tests/unit/test_cellprofiler_shape_hotpath.py +++ b/tests/unit/test_cellprofiler_shape_hotpath.py @@ -22,12 +22,18 @@ from openhcs.core.runtime_tabular_values import MeasurementObjectRowIdentity from openhcs.processing.backends.cellprofiler.shape import ( MeasureObjectSizeShapeModule, + ShapeMeasurementBackendStrategy, ShapeObjectFeatureValueTable, measure_object_size_shape, measure_object_size_shape_feature_arrays, ) +@pytest.fixture(scope="module", autouse=True) +def prepared_shape_backends() -> None: + ShapeMeasurementBackendStrategy.prepare_registered_family() + + def _assert_rows_strict( actual_rows: Sequence[Mapping[str, object]], expected_rows: Sequence[Mapping[str, object]], @@ -184,7 +190,9 @@ def test_nonempty_stacked_shape_vectors_and_surface_areas_match_oracles() -> Non np.testing.assert_allclose( [row["SurfaceArea"] for row in actual_rows], expected_surface_areas, - rtol=0.0, + # Binary cube aggregation uses the approved CellProfiler tolerance. + # Float32 mesh translation/reduction can differ by a few micro-units. + rtol=1e-6, atol=1e-6, ) @@ -252,3 +260,90 @@ def test_sparse_high_id_shape_preserves_label_domain_and_row_identity() -> None: assert ( rows[0][MeasureObjectSizeShapeModule.MeasurementFeature.CENTER_Y.value] == 2.5 ) + + +@pytest.mark.parametrize("shape", ((1, 6, 9), (5, 9, 11), (5, 9, 4096))) +@pytest.mark.parametrize("readonly", (False, True)) +def test_prepared_3d_cohort_matches_reference_at_boundaries_and_offsets( + shape: tuple[int, int, int], + readonly: bool, +) -> None: + labels = np.zeros(shape, dtype=np.int32) + labels[:, 1:5, -6:-1] = 892 + labels[0, :2, :2] = 19 + labels[:2, 2:3, -4:-3] = 0 + original = labels.copy() + labels.setflags(write=not readonly) + backend = ShapeMeasurementBackendStrategy.for_memory_type() + spacing = (2.75, 0.65, 0.65) + expected, expected_ids = ShapeMeasurementBackendStrategy.feature_arrays_3d( + backend, labels, calculate_advanced=False, spacing=spacing + ) + actual, actual_ids = backend.feature_arrays_3d( + labels, calculate_advanced=False, spacing=spacing + ) + assert set(actual) == set(expected) + assert actual_ids.dtype == expected_ids.dtype + np.testing.assert_array_equal(actual_ids, expected_ids) + for key in expected: + np.testing.assert_allclose(actual[key], expected[key], atol=1e-6, rtol=1e-6) + np.testing.assert_array_equal(labels, original) + + +def test_direct_unprepared_3d_call_retains_reference_behavior(monkeypatch) -> None: + labels = np.zeros((4, 6, 8), dtype=np.int32) + labels[1:3, 1:4, 2:6] = 5 + backend = ShapeMeasurementBackendStrategy.for_memory_type() + expected, expected_ids = ShapeMeasurementBackendStrategy.feature_arrays_3d( + backend, labels, calculate_advanced=False, spacing=(1.0, 1.0, 1.0) + ) + monkeypatch.setattr( + ShapeMeasurementBackendStrategy, "_unit_surface_triangles", None + ) + actual, actual_ids = measure_object_size_shape_feature_arrays( + labels, calculate_advanced=False, calculate_zernikes=False + ) + np.testing.assert_array_equal(actual_ids, expected_ids) + for key in expected: + np.testing.assert_array_equal(actual[key], expected[key]) + + +def test_local_moment_capacity_depends_on_extent_not_global_offset() -> None: + backend = ShapeMeasurementBackendStrategy.for_memory_type() + distant = np.asarray( + [[10**12, 10**12, 10**12, 10**12 + 5, 10**12 + 7, 10**12 + 9]], dtype=np.int64 + ) + assert backend._local_moments_fit_int64(np.asarray([100]), distant) + oversized = np.asarray([[0, 0, 0, 2_000_000_001, 1, 1]], dtype=np.int64) + assert backend._local_moments_fit_int64(np.asarray([2]), oversized) + assert not backend._local_moments_fit_int64(np.asarray([3]), oversized) + + +def test_prepared_3d_advanced_solidity_and_kernel_signatures_match_reference() -> None: + from openhcs.processing.backends.cellprofiler.shape import ( + _label_bounds_counts_3d_numba, + _label_local_moments_surface_euler_3d_numba, + ) + + kernels = ( + _label_bounds_counts_3d_numba, + _label_local_moments_surface_euler_3d_numba, + ) + signatures = tuple(tuple(kernel.signatures) for kernel in kernels) + labels = np.zeros((7, 18, 20), dtype=np.int32) + labels[1:5, 2:9, 3:11] = 1 + labels[2:6, 10:16, 12:18] = 3 + labels[2:4, 12:14, 12:15] = 0 + backend = ShapeMeasurementBackendStrategy.for_memory_type() + for writable in (True, False): + labels.setflags(write=writable) + expected, expected_ids = ShapeMeasurementBackendStrategy.feature_arrays_3d( + backend, labels, calculate_advanced=True, spacing=(1.0, 1.0, 1.0) + ) + actual, actual_ids = backend.feature_arrays_3d( + labels, calculate_advanced=True, spacing=(1.0, 1.0, 1.0) + ) + np.testing.assert_array_equal(actual_ids, expected_ids) + for key in expected: + np.testing.assert_allclose(actual[key], expected[key], atol=1e-6, rtol=1e-6) + assert signatures == tuple(tuple(kernel.signatures) for kernel in kernels)