diff --git a/src/squidpy/experimental/im/_calculate_image_features.py b/src/squidpy/experimental/im/_calculate_image_features.py index 4ae5de1f7..ee6f7c291 100644 --- a/src/squidpy/experimental/im/_calculate_image_features.py +++ b/src/squidpy/experimental/im/_calculate_image_features.py @@ -31,7 +31,6 @@ _run_tiled, build_tile_specs, compute_cell_info, - compute_cell_info_multiscale, compute_cell_info_tiled, extract_labels_tile_lazy, extract_tile_lazy, @@ -107,6 +106,14 @@ "cp_measure:correlation_rwc": {"correlation_rwc": True}, } +# cp_measure granularity samples each crop on a 1/4 then 1/16 grid anchored at the crop +# origin and estimates background with an opening reaching ~320 px. Tiles that request it +# are aligned to that grid (1 / (subsample_size * image_sample_size) at cp_measure's +# defaults) and padded for context; on real images this makes the lower granularity bins +# match an untiled run closely. +_GRANULARITY_ALIGN = 16 +_GRANULARITY_PAD = 256 + # cp_measure correlation features need >=2 channels (they correlate channel pairs). _CP_CORRELATION_KEYS = frozenset(_CPMEASURE_FLAGS["cp_measure:correlation"]) @@ -315,6 +322,35 @@ def _build_cp_config(cp_flags: dict[str, bool], channel_names: list[str]) -> dic # Per-tile dispatcher # --------------------------------------------------------------------------- +# Coordinate-valued feature columns (prefixes) by the axis they measure. Tiles are +# featurized on crops, so these come out crop-local and are shifted by the crop origin. +_Y_POSITION_COLS = ( + "Center_Y", + "BoundingBoxMinimum_Y", + "BoundingBoxMaximum_Y", + "Location_CenterMassIntensity_Y__", + "Location_MaxIntensity_Y__", + "centroid-0", +) +_X_POSITION_COLS = ( + "Center_X", + "BoundingBoxMinimum_X", + "BoundingBoxMaximum_X", + "Location_CenterMassIntensity_X__", + "Location_MaxIntensity_X__", + "centroid-1", +) + + +def _shift_positions(df: pd.DataFrame, dy: int, dx: int) -> pd.DataFrame: + """Shift coordinate-valued feature columns by ``(dy, dx)`` pixels (in place).""" + for col in df.columns: + if col.startswith(_Y_POSITION_COLS): + df[col] += dy + elif col.startswith(_X_POSITION_COLS): + df[col] += dx + return df + def _featurize_tile( tile_image: np.ndarray | None, @@ -639,12 +675,13 @@ def _align_to_image_grid( image_da: xr.DataArray, labels_da: xr.DataArray, align_mode: Literal["strict", "rasterize"], -) -> tuple[xr.DataArray, xr.DataArray]: +) -> tuple[xr.DataArray, xr.DataArray, tuple[int, int]]: """Crop image and labels to their pixel-grid overlap, honoring transforms. Cells falling outside the overlap rectangle are dropped (logged). Under ``align_mode='strict'`` a non-pixel-aligned relative transform raises; under - ``'rasterize'`` the labels are resampled onto the image grid. + ``'rasterize'`` the labels are resampled onto the image grid. Also returns + the ``(y, x)`` origin of the labels crop in the labels' pixel grid. """ cs = _shared_coordinate_system(sdata, image_key, labels_key) affine = _relative_affine(sdata, image_key, labels_key, cs) @@ -706,7 +743,7 @@ def _align_to_image_grid( f"Dropped {cells_outside} cell(s) fully and {len(partial_ids)} cell(s) partially outside the image extent." ) - return image_crop, labels_crop + return image_crop, labels_crop, (lbl_y0, lbl_x0) # --------------------------------------------------------------------------- @@ -762,13 +799,15 @@ def _prepare_lazy( scale: str | None, channels: list[str] | None, align_mode: Literal["strict", "rasterize"], -) -> tuple[xr.DataArray | None, xr.DataArray, list[str]]: - """Return lazy image and labels DataArrays, plus channel names. +) -> tuple[xr.DataArray | None, xr.DataArray, list[str], tuple[int, int]]: + """Return lazy image and labels DataArrays, channel names, and the labels origin. ``image_da`` is ``None`` (and ``channel_names`` empty) for a morphology-only - run with no ``image_key``. Does NOT call ``.compute()`` - arrays stay lazy - for on-demand tile reads. For the shapes->labels path, labels are - materialized but wrapped in a DataArray for a uniform interface. + run with no ``image_key``. The origin is the ``(y, x)`` offset of + ``labels_da`` in the labels' pixel grid (non-zero only when alignment + crops). Does NOT call ``.compute()`` - arrays stay lazy for on-demand + tile reads. For the shapes->labels path, labels are materialized but + wrapped in a DataArray for a uniform interface. """ _validate_inputs(sdata, image_key, labels_key, shapes_key, scale) @@ -798,11 +837,14 @@ def _prepare_lazy( # Align labels to the image pixel grid via SpatialData transformations. # Only meaningful with a real labels element + an image; the shapes->labels # path already rasterized onto the image grid (identity transform -> no-op). + origin = (0, 0) if image_da is not None and labels_key is not None: - image_da, labels_da = _align_to_image_grid(sdata, image_key, labels_key, image_da, labels_da, align_mode) + image_da, labels_da, origin = _align_to_image_grid( + sdata, image_key, labels_key, image_da, labels_da, align_mode + ) if image_da is None: - return image_da, labels_da, [] + return image_da, labels_da, [], origin # Resolve channel names through spatialdata's canonical accessor so we # honor c_coords set at parse time. Always cast to str. @@ -830,29 +872,11 @@ def _prepare_lazy( else: ch_names = all_ch - return image_da, labels_da, ch_names + return image_da, labels_da, ch_names, origin -def _compute_centroids( - sdata: SpatialData, - labels_key: str | None, - labels_da: xr.DataArray, - scale: str | None, -) -> dict[int, CellInfo]: - """Compute cell centroids using the most efficient strategy available.""" - # Multiscale: the coarse-scale fast path is only valid when alignment did not - # crop labels_da; after a crop, recompute from it so centroids and tiling - # share one frame. - if labels_key is not None and isinstance(sdata.labels[labels_key], xr.DataTree): - full = _select_scale_array(sdata.labels[labels_key], scale) - full_grid = (full.sizes.get("y"), full.sizes.get("x")) - cur_grid = (labels_da.sizes.get("y"), labels_da.sizes.get("x")) - if cur_grid == full_grid: - logg.info("Computing centroids from coarse scale.") - return compute_cell_info_multiscale(sdata.labels[labels_key], target_scale=scale or "scale0") - logg.info("Computing centroids in tiled mode (aligned multiscale labels).") - return compute_cell_info_tiled(labels_da) - +def _compute_centroids(labels_da: xr.DataArray) -> dict[int, CellInfo]: + """Compute cell centroids and bounding boxes on the featurized labels grid.""" # Small enough to fit in memory - direct regionprops n_pixels = labels_da.sizes.get("y", 1) * labels_da.sizes.get("x", 1) if n_pixels <= 4096 * 4096: @@ -861,11 +885,49 @@ def _compute_centroids( lbl_np = lbl_np.squeeze() return compute_cell_info(lbl_np) - # Large single-scale - tiled centroid computation - logg.info("Computing centroids in tiled mode (large single-scale labels).") + # Large - tiled centroid computation + logg.info("Computing centroids in tiled mode (large labels).") return compute_cell_info_tiled(labels_da) +def _stack_tiles(tile_dfs: list[pd.DataFrame], drop_constant: bool) -> tuple[np.ndarray, np.ndarray, list[str]]: + """Stack per-tile float32 features into one label-sorted matrix; return ``(labels, X, columns)``. + + Fills a single preallocated array instead of concat -> sort -> cast, which + copied the full table three times. With ``drop_constant``, zero-variance + columns (all-NaN, or NaN-free with a single value) are found from per-tile + min/max and never copied. Skipped for a single cell, where every column is + trivially constant. + """ + # Peak is tiles + output (~2 copies); freeing tiles while filling would halve it if ever needed. + columns = list(dict.fromkeys(c for df in tile_dfs for c in df.columns)) + blocks = [(df if list(df.columns) == columns else df.reindex(columns=columns)).to_numpy() for df in tile_dfs] + labels = np.concatenate([df.index.to_numpy() for df in tile_dfs]) + + keep = np.ones(len(columns), dtype=bool) + if drop_constant and len(labels) > 1: + # min/max propagate NaN (so a NaN-mixed column never compares equal); + # nanmax is NaN only where a column is all-NaN. + col_min = np.min([b.min(axis=0) for b in blocks], axis=0) + col_max = np.max([b.max(axis=0) for b in blocks], axis=0) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) # all-NaN slices + all_nan = np.isnan(np.nanmax([np.nanmax(b, axis=0) for b in blocks], axis=0)) + keep = ~(all_nan | (col_min == col_max)) + if not keep.all(): + logg.warning(f"Dropped {int((~keep).sum())} constant feature(s) with no variance across cells.") + + order = np.argsort(labels, kind="stable") + rows = np.empty_like(order) + rows[order] = np.arange(len(order)) # output row of each stacked cell + X = np.empty((len(labels), int(keep.sum())), dtype=np.float32) + start = 0 + for b in blocks: + X[rows[start : start + len(b)]] = b[:, keep] + start += len(b) + return labels[order], X, [c for c, k in zip(columns, keep, strict=True) if k] + + # --------------------------------------------------------------------------- # Main function # --------------------------------------------------------------------------- @@ -949,7 +1011,9 @@ def calculate_image_features( the skimage-only props (``centroid_local``, ``feret_diameter_max``) are kept. cp_measure computes its groups all-or-nothing, so it wins. tile_size - Side length of the tiling grid (pixels). + Side length of the tiling grid (pixels). With ``"cp_measure:granularity"`` + each tile is padded by 256 px of image context, so prefer + ``tile_size >= 2048`` there to keep the extra reads small. align_mode How to handle image/labels whose pixel grids do not match (via their SpatialData transformations). @@ -994,6 +1058,16 @@ def calculate_image_features( constant features removed by ``drop_constant_features`` are logged at WARNING level. + Positional features (cp_measure ``Center_*``, ``BoundingBox*``, + ``Location_*``; skimage ``centroid-*``) are in pixel units of the labels + grid at ``scale`` (the image grid when labels are resampled via + ``align_mode="rasterize"`` or ``shapes_key``), independent of ``tile_size``. + ``"cp_measure:granularity"`` depends on the image around each cell. Tiles + that compute it are aligned to its sampling grid and padded, which keeps + the lower (fine-granule) bins close to an untiled run; the higher bins + depend on image context far beyond any tile (as they do across CellProfiler + fields of view) and still vary with ``tile_size``. + With ``n_jobs > 1`` a ``LocalCluster`` is started, which spawns worker processes. On macOS/Windows (spawn start method) the calling code must be guarded by ``if __name__ == "__main__":`` (the standard Python multiprocessing @@ -1049,7 +1123,7 @@ def calculate_image_features( if channels is not None: raise ValueError("`channels` selection requires `image_key`.") - image_da, labels_da, channel_names = _prepare_lazy( + image_da, labels_da, channel_names, origin = _prepare_lazy( sdata, image_key, labels_key, shapes_key, scale, channels, align_mode ) @@ -1074,16 +1148,24 @@ def calculate_image_features( cp_config = _build_cp_config(parsed.cp_flags, channel_names) if parsed.cp_flags is not None else None # --- Warmup: compute centroids without materializing full arrays --- - cell_info = _compute_centroids(sdata, labels_key, labels_da, scale) + cell_info = _compute_centroids(labels_da) if not cell_info: raise ValueError("No cells found in labels (all zeros).") H, W = yx_size(labels_da) # --- Tile --- - # overlap_margin="auto" derives the minimum safe margin from the largest cell; - # not exposed -- any manual value either truncates boundary cells or wastes reads. - specs = build_tile_specs((H, W), cell_info, tile_size=tile_size, overlap_margin="auto") + # overlap_margin="auto" crops each tile to its owned cells' bounding boxes (+1 px); + # not exposed -- a fixed margin either truncates boundary cells or wastes reads. + granularity = parsed.cp_flags is not None and (not parsed.cp_flags or parsed.cp_flags.get("granularity", False)) + specs = build_tile_specs( + (H, W), + cell_info, + tile_size=tile_size, + overlap_margin="auto", + pad=_GRANULARITY_PAD if granularity else 1, + align=_GRANULARITY_ALIGN if granularity else 1, + ) total_tiles = len(specs) logg.info(f"Tiling input into {total_tiles} tile(s) of size {tile_size} px.") @@ -1092,10 +1174,13 @@ def calculate_image_features( def _process_one(spec, image_da, labels_da): with threadpool_limits(limits=1): if image_da is None: - tile_lbl = extract_labels_tile_lazy(labels_da, spec) - return _featurize_tile(None, tile_lbl, parsed, channel_names, cp_config=cp_config) - tile_img, tile_lbl = extract_tile_lazy(image_da, labels_da, spec) - return _featurize_tile(tile_img, tile_lbl, parsed, channel_names, cp_config=cp_config) + tile_img, tile_lbl = None, extract_labels_tile_lazy(labels_da, spec) + else: + tile_img, tile_lbl = extract_tile_lazy(image_da, labels_da, spec) + df = _featurize_tile(tile_img, tile_lbl, parsed, channel_names, cp_config=cp_config) + # Report positions in the labels' pixel grid, not the tile crop's. The output + # is float32, so cast here: half the memory and transfer per tile. + return _shift_positions(df, origin[0] + spec.crop[0], origin[1] + spec.crop[1]).astype(np.float32) # cp_measure is GIL-bound, so kind="processes" (an active Client wins if set). results = _run_tiled( @@ -1107,28 +1192,16 @@ def _process_one(spec, image_da, labels_da): if not tile_dfs: raise ValueError("No features computed for any tile.") - # Sort by cell label for deterministic output. inf/NaN handling happens - # in one numpy pass below to avoid two extra full-table allocations. - combined = pd.concat(tile_dfs, axis=0).sort_index() - - # Drop zero-variance features (nunique(dropna=False) treats an all-NaN column - # as constant too). Skipped for a single cell, where every column is trivially - # constant and the filter would drop everything. - if drop_constant_features and len(combined) > 1: - constant_cols = list(combined.columns[combined.nunique(dropna=False) <= 1]) - if constant_cols: - logg.warning(f"Dropped {len(constant_cols)} constant feature(s) with no variance across cells.") - combined = combined.drop(columns=constant_cols) + labels, arr, columns = _stack_tiles(tile_dfs, drop_constant_features) # --- Build AnnData --- # Exactly one of labels_key / shapes_key is set (enforced in _validate_inputs). region_key_value = labels_key or shapes_key - arr = combined.to_numpy(dtype=np.float32, copy=True) if invalid_as_zero: np.nan_to_num(arr, copy=False, nan=0.0, posinf=0.0, neginf=0.0) adata = ad.AnnData(X=arr) - adata.var_names = list(combined.columns) + adata.var_names = columns adata.uns["spatialdata_attrs"] = { "region": region_key_value, @@ -1140,7 +1213,7 @@ def _process_one(spec, image_da, labels_da): if shapes_key is not None and len(sdata.shapes[shapes_key]) == len(adata): adata.obs["label_id"] = sdata.shapes[shapes_key].index.values else: - adata.obs["label_id"] = combined.index.values + adata.obs["label_id"] = labels # obs_names are the cell's label-image ID (the label_id), as str for AnnData. adata.obs_names = adata.obs["label_id"].astype(str).values diff --git a/src/squidpy/experimental/im/_tiling.py b/src/squidpy/experimental/im/_tiling.py index d23720990..187d830f2 100644 --- a/src/squidpy/experimental/im/_tiling.py +++ b/src/squidpy/experimental/im/_tiling.py @@ -69,8 +69,8 @@ class TileSpec: The non-overlapping region ``(y0, x0, y1, x1)`` used for centroid ownership. Tiles partition the image into a grid of base regions. crop - The extended region ``(y0, x0, y1, x1)`` that includes the overlap - margin. This is the actual slice extracted from the image/labels. + The region ``(y0, x0, y1, x1)`` that contains every owned cell. + This is the actual slice extracted from the image/labels. owned_ids Label IDs whose centroid falls inside ``base``. Only these labels are kept in the tile's mask; all others are zeroed out. @@ -112,48 +112,6 @@ def compute_cell_info(labels: np.ndarray) -> dict[int, CellInfo]: return info -def compute_cell_info_multiscale( - labels_node: xr.DataTree, - target_scale: str = "scale0", -) -> dict[int, CellInfo]: - """Compute centroids using the coarsest scale of a multiscale label pyramid. - - Reads only the smallest resolution, then scales coordinates to *target_scale*. - """ - available = list(labels_node.keys()) - if not available: - return {} - - def _spatial_size(k: str) -> int: - h, w = yx_size(labels_node[k].ds["image"]) - return h * w - - coarsest = min(available, key=_spatial_size) - coarse_labels = np.asarray(labels_node[coarsest].ds["image"].values).squeeze() - - if coarse_labels.ndim != 2: - raise ValueError(f"Expected 2-D labels at scale {coarsest}, got shape {coarse_labels.shape}") - - target_h, target_w = yx_size(labels_node[target_scale].ds["image"]) - coarse_h, coarse_w = coarse_labels.shape - scale_y = target_h / coarse_h - scale_x = target_w / coarse_w - - props = regionprops(coarse_labels) - return { - p.label: CellInfo( - label=p.label, - centroid_y=p.centroid[0] * scale_y, - centroid_x=p.centroid[1] * scale_x, - bbox_h=int(np.ceil((p.bbox[2] - p.bbox[0]) * scale_y)), - bbox_w=int(np.ceil((p.bbox[3] - p.bbox[1]) * scale_x)), - bbox_y0=int(np.floor(p.bbox[0] * scale_y)), - bbox_x0=int(np.floor(p.bbox[1] * scale_x)), - ) - for p in props - } - - @dataclass class _Accum: """Per-label running totals while streaming chunks (a cell may span chunks).""" @@ -219,21 +177,13 @@ def compute_cell_info_tiled( # Tile spec building -def _auto_margin(cell_info: dict[int, CellInfo]) -> int: - """Compute the minimum margin that covers the largest cell's half-extent.""" - if not cell_info: - return 0 - max_extent = max(max(c.bbox_h, c.bbox_w) for c in cell_info.values()) - # Centroid can be at most half a bbox away from the cell's edge. - # Add 1 pixel for safety (rounding / off-by-one). - return int(np.ceil(max_extent / 2)) + 1 - - def build_tile_specs( grid_shape: tuple[int, int], cell_info: dict[int, CellInfo], tile_size: int = 2048, overlap_margin: int | Literal["auto"] = "auto", + pad: int = 1, + align: int = 1, ) -> list[TileSpec]: """Build tile specifications from pre-computed centroids. @@ -244,13 +194,21 @@ def build_tile_specs( grid_shape ``(height, width)`` of the full-resolution labels grid. cell_info - Pre-computed centroids from :func:`compute_cell_info`, - :func:`compute_cell_info_multiscale`, or :func:`compute_cell_info_tiled`. + Pre-computed centroids from :func:`compute_cell_info` or + :func:`compute_cell_info_tiled`. tile_size Side length of the non-overlapping base grid cells. overlap_margin - Pixel margin added around each base region. ``"auto"`` computes the - minimum margin from the largest cell's bounding box. + ``"auto"`` crops each tile to the union of its owned cells' bounding + boxes plus ``pad`` pixels, so every owned cell is whole and has + background around its boundary (edge and radial features need 1 px). + An integer instead adds that fixed margin around the base region. + pad + Context in pixels around the owned cells for ``overlap_margin="auto"``. + align + Snap each crop's origin down, and its size up, to multiples of + ``align`` pixels, for features that sample the crop on a grid anchored + at its origin. Returns ------- @@ -261,8 +219,8 @@ def build_tile_specs( if tile_size <= 0: raise ValueError(f"tile_size must be positive, got {tile_size}") - margin = _auto_margin(cell_info) if overlap_margin == "auto" else int(overlap_margin) - if margin < 0: + margin = None if overlap_margin == "auto" else int(overlap_margin) + if margin is not None and margin < 0: raise ValueError(f"overlap_margin must be non-negative, got {margin}") cell_to_tile: dict[int, tuple[int, int]] = {} @@ -282,10 +240,22 @@ def build_tile_specs( by1 = min(by0 + tile_size, height) bx1 = min(bx0 + tile_size, width) - cy0 = max(by0 - margin, 0) - cx0 = max(bx0 - margin, 0) - cy1 = min(by1 + margin, height) - cx1 = min(bx1 + margin, width) + if margin is None: + cells = [cell_info[lid] for lid in owned] + cy0 = max(min(c.bbox_y0 for c in cells) - pad, 0) + cx0 = max(min(c.bbox_x0 for c in cells) - pad, 0) + cy1 = min(max(c.bbox_y0 + c.bbox_h for c in cells) + pad, height) + cx1 = min(max(c.bbox_x0 + c.bbox_w for c in cells) + pad, width) + else: + cy0 = max(by0 - margin, 0) + cx0 = max(bx0 - margin, 0) + cy1 = min(by1 + margin, height) + cx1 = min(bx1 + margin, width) + + if align > 1: + cy0, cx0 = cy0 - cy0 % align, cx0 - cx0 % align + cy1 = min(cy1 + (cy0 - cy1) % align, height) + cx1 = min(cx1 + (cx0 - cx1) % align, width) specs.append( TileSpec( diff --git a/src/squidpy/experimental/tl/_tiling_qc.py b/src/squidpy/experimental/tl/_tiling_qc.py index a8feb43c9..01581158a 100644 --- a/src/squidpy/experimental/tl/_tiling_qc.py +++ b/src/squidpy/experimental/tl/_tiling_qc.py @@ -46,7 +46,6 @@ _run_tiled, build_tile_specs, compute_cell_info, - compute_cell_info_multiscale, compute_cell_info_tiled, extract_labels_tile_lazy, ) @@ -401,17 +400,8 @@ def _score_tile( # Centroid computation (shared logic with _feature.py) -def _compute_centroids_for_labels( - sdata: sd.SpatialData, - labels_key: str, - labels_da: xr.DataArray, - scale: str | None, -) -> dict: - """Compute cell centroids using the most efficient strategy available.""" - if isinstance(sdata.labels[labels_key], xr.DataTree): - logg.info("Computing centroids from coarse scale.") - return compute_cell_info_multiscale(sdata.labels[labels_key], target_scale=scale or "scale0") - +def _compute_centroids_for_labels(labels_da: xr.DataArray) -> dict: + """Compute cell centroids and bounding boxes on the scored labels grid.""" n_pixels = labels_da.sizes.get("y", 1) * labels_da.sizes.get("x", 1) if n_pixels <= 4096 * 4096: lbl_np = labels_da.values @@ -419,7 +409,7 @@ def _compute_centroids_for_labels( lbl_np = lbl_np.squeeze() return compute_cell_info(lbl_np) - logg.info("Computing centroids in tiled mode (large single-scale labels).") + logg.info("Computing centroids in tiled mode (large labels).") return compute_cell_info_tiled(labels_da) @@ -469,8 +459,9 @@ def calculate_tiling_qc( tile_size Side length of the tiling grid (pixels). overlap_margin - Overlap around each tile. ``"auto"`` computes the minimum from - the largest cell's bounding box. + Overlap around each tile. ``"auto"`` crops each tile to the + bounding boxes of the cells it owns (plus 1 pixel), so no cell is + truncated; an integer adds that fixed margin around the tile. downsample Factor by which to downsample each cell's bounding-box crop before contour extraction. Straightness is scale-invariant, @@ -561,7 +552,7 @@ def calculate_tiling_qc( labels_da = resolve_labels_array(sdata, labels_key, scale) - cell_info = _compute_centroids_for_labels(sdata, labels_key, labels_da, scale) + cell_info = _compute_centroids_for_labels(labels_da) if not cell_info: raise ValueError("No cells found in labels (all zeros).") diff --git a/tests/experimental/test_calculate_image_features.py b/tests/experimental/test_calculate_image_features.py index e67212321..f8dbf2ce6 100644 --- a/tests/experimental/test_calculate_image_features.py +++ b/tests/experimental/test_calculate_image_features.py @@ -438,19 +438,28 @@ def test_channel_selection_invalid(self, sdata_synthetic): # --- Tiled vs non-tiled equivalence --- - def test_tiled_vs_single_tile_equivalence(self, sdata_synthetic): - """Tile-invariant features should be identical whether we tile or not. + @pytest.mark.parametrize( + "features", + [ + ["skimage:morphology", "squidpy:summary"], + ["cp_measure:sizeshape", "cp_measure:intensity", "cp_measure:radial"], + ], + ) + def test_tiled_vs_single_tile_equivalence(self, sdata_synthetic, features): + """Features are identical whether we tile or not. - Position-dependent features (centroid, perimeter_crofton) are expected - to differ across tile boundaries, so we test with ``area`` and - ``squidpy:summary`` which depend only on the cell's pixel values. + Covers positional features (centroids, bounding boxes, intensity + locations), which must be reported in the labels' pixel frame rather + than the tile's, and edge/radial features, which need background + around the cell boundary inside the tile crop. """ kw = { "image_key": "test_img", "labels_key": "test_labels", - "features": ["skimage:morphology:area", "squidpy:summary"], + "features": features, "inplace": False, "invalid_as_zero": True, + "drop_constant_features": False, } # Single tile (tile_size >= image -> no tiling) result_single = sq.experimental.im.calculate_image_features(sdata_synthetic, tile_size=1000, **kw) @@ -473,6 +482,41 @@ def test_tiled_vs_single_tile_equivalence(self, sdata_synthetic): np.testing.assert_array_equal(df_single.index, df_tiled.index) np.testing.assert_allclose(df_single.values, df_tiled.values, rtol=1e-5, atol=1e-5) + def test_granularity_independent_of_tile_size_within_context(self, sdata_synthetic): + """Granularity needs image context; when the tile padding covers the image, tiling changes nothing.""" + kw = { + "image_key": "test_img", + "labels_key": "test_labels", + "features": ["cp_measure:granularity"], + "inplace": False, + "drop_constant_features": False, + } + single = sq.experimental.im.calculate_image_features(sdata_synthetic, tile_size=1000, **kw).to_df() + tiled = sq.experimental.im.calculate_image_features(sdata_synthetic, tile_size=100, **kw).to_df() + pd.testing.assert_frame_equal(tiled.loc[single.index], single) + + def test_asymmetric_cell_not_truncated_by_tiling(self): + """A cell reaching far from its centroid (blob + long process) stays whole when tiled.""" + labels = np.zeros((200, 200), dtype=np.int32) + yy, xx = np.ogrid[:200, :200] + labels[(yy - 50) ** 2 + (xx - 60) ** 2 <= 20**2] = 1 + labels[49:52, 60:198] = 1 # the process pulls the centroid to x~77, still in tile 0 + labels[140:160, 140:160] = 2 + sdata = SpatialData( + images={"img": Image2DModel.parse(np.ones((1, 200, 200), dtype=np.uint8), dims=("c", "y", "x"))}, + labels={"lbl": Labels2DModel.parse(labels, dims=("y", "x"))}, + ) + result = sq.experimental.im.calculate_image_features( + sdata, + image_key="img", + labels_key="lbl", + features=["skimage:morphology:area"], + tile_size=100, + inplace=False, + drop_constant_features=False, + ) + np.testing.assert_array_equal(result[["1", "2"], "area"].X.ravel(), [(labels == 1).sum(), 400]) + # --- Parallelization --- def test_n_jobs_produces_same_result(self, sdata_synthetic): @@ -917,6 +961,24 @@ def test_multiscale_featurized_with_scale(self): assert set(adata.obs["label_id"].astype(int)) == set(range(1, 17)) np.testing.assert_array_equal(adata[:, "area"].X.ravel(), np.full(16, 900.0)) + def test_multiscale_keeps_cells_absent_from_coarse_scales(self): + """Cells too small to survive downsampling are still featurized at the requested scale.""" + labels = np.zeros((128, 128), dtype=np.int32) + labels[8:40, 8:40] = 1 + tiny = [(y, x) for y in range(50, 120, 9) for x in range(50, 120, 9)] + for lid, (y, x) in enumerate(tiny, start=2): + labels[y, x] = lid # 1-px cells: gone at scale1/scale2 + sdata = SpatialData(labels={"lbl": Labels2DModel.parse(labels, dims=("y", "x"), scale_factors=[2, 2])}) + adata = sq.experimental.im.calculate_image_features( + sdata, + labels_key="lbl", + scale="scale0", + features=["skimage:morphology:area"], + inplace=False, + drop_constant_features=False, + ) + assert set(adata.obs["label_id"].astype(int)) == set(range(1, len(tiny) + 2)) + def test_invalid_scale_name(self): sdata = _multiscale_sdata(multiscale_image=True, multiscale_labels=True) with pytest.raises(ValueError, match="Scale 'scale9' not found"): @@ -1072,6 +1134,24 @@ def test_translation_drops_outside_and_partial_cells(self): assert 0 < result.n_obs < n_cells np.testing.assert_array_equal(result[:, "area"].X.ravel(), 625.0) + def test_positions_in_labels_pixel_frame(self): + """Centroids are reported in the labels' own pixel grid, not the cropped overlap.""" + # Labels pixel (30, 30) lands on image pixel (0, 0): the overlap crop starts at 30. + sdata = _toy_sdata(labels_translation=(-30, -30)) + labels = sdata.labels["lbl"].values + result = sq.experimental.im.calculate_image_features( + sdata, + image_key="img", + labels_key="lbl", + features=["skimage:morphology:centroid"], + inplace=False, + drop_constant_features=False, + ) + assert result.n_obs > 0 + for lid in result.obs["label_id"].astype(int): + ys, xs = np.nonzero(labels == lid) + np.testing.assert_allclose(result[str(lid), ["centroid-0", "centroid-1"]].X.ravel(), [ys.mean(), xs.mean()]) + def test_multiscale_rasterize_raises(self): sdata = _toy_sdata(labels_scale=(1.3, 1.3), multiscale=True) with pytest.raises(ValueError, match="not supported for multiscale"): diff --git a/tests/experimental/test_tiling.py b/tests/experimental/test_tiling.py index f3728f199..fb44d370d 100644 --- a/tests/experimental/test_tiling.py +++ b/tests/experimental/test_tiling.py @@ -19,7 +19,6 @@ _zero_non_owned, build_tile_specs, compute_cell_info, - compute_cell_info_multiscale, compute_cell_info_tiled, extract_tile_lazy, ) @@ -255,6 +254,21 @@ def test_crop_contains_owned_cells_fully(self, brick_labels): f"Cell {lid} x-range [{cell_x0:.0f},{cell_x1:.0f}] not in crop x-range [{cx0},{cx1}]" ) + def test_aligned_crops_on_grid(self, brick_labels): + """With ``align``, crops start and span on the alignment grid and still cover the unaligned crop.""" + labels, _ = brick_labels + H, W = labels.shape + cell_info = compute_cell_info(labels) + plain = build_tile_specs(labels.shape, cell_info, tile_size=_TILE_SIZE) + aligned = build_tile_specs(labels.shape, cell_info, tile_size=_TILE_SIZE, align=16) + for p, a in zip(plain, aligned, strict=True): + (py0, px0, py1, px1), (ay0, ax0, ay1, ax1) = p.crop, a.crop + assert a.owned_ids == p.owned_ids + assert ay0 % 16 == 0 and ax0 % 16 == 0 + assert (ay1 - ay0) % 16 == 0 or ay1 == H + assert (ax1 - ax0) % 16 == 0 or ax1 == W + assert ay0 <= py0 and ax0 <= px0 and ay1 >= py1 and ax1 >= px1 + class TestBuildTileSpecsEdgeCases: def test_empty_labels(self): @@ -404,42 +418,7 @@ def _plot_tile_assignment(labels, specs, title=""): ax.set_ylabel("y") -# Lazy / multiscale helpers - - -def _make_multiscale_tree(labels: np.ndarray, n_scales: int = 3) -> xr.DataTree: - """Build a tiny multiscale DataTree by integer-downsampling.""" - scales: dict[str, xr.DataTree] = {} - for i in range(n_scales): - step = 2**i - sub = labels[::step, ::step] - ds = xr.Dataset({"image": xr.DataArray(sub, dims=("y", "x"))}) - scales[f"scale{i}"] = xr.DataTree(ds) - return xr.DataTree.from_dict(scales) - - -class TestComputeCellInfoMultiscale: - def test_target_is_coarsest_matches_eager(self): - labels, _ = _make_brick_labels(gap=10) - tree = _make_multiscale_tree(labels, n_scales=3) - # scale2 is coarsest. Target it -> use that scale directly. - info_ms = compute_cell_info_multiscale(tree, target_scale="scale2") - info_eager = compute_cell_info(tree["scale2"].ds["image"].values) - assert set(info_ms.keys()) == set(info_eager.keys()) - for lid in info_ms: - assert info_ms[lid].centroid_y == pytest.approx(info_eager[lid].centroid_y, abs=0.5) - assert info_ms[lid].centroid_x == pytest.approx(info_eager[lid].centroid_x, abs=0.5) - - def test_rescale_to_finer(self): - labels, _ = _make_brick_labels(gap=10) - tree = _make_multiscale_tree(labels, n_scales=3) - info_ms = compute_cell_info_multiscale(tree, target_scale="scale0") - info_eager = compute_cell_info(labels) - # Centroids should be close (within ~1 px due to coarse-scale quantization) - assert set(info_ms.keys()) == set(info_eager.keys()) - for lid in info_ms: - assert info_ms[lid].centroid_y == pytest.approx(info_eager[lid].centroid_y, abs=4.0) - assert info_ms[lid].centroid_x == pytest.approx(info_eager[lid].centroid_x, abs=4.0) +# Lazy helpers class TestComputeCellInfoTiled: diff --git a/tests/experimental/test_tiling_qc.py b/tests/experimental/test_tiling_qc.py index b0ec4adcb..2e2b53fb7 100644 --- a/tests/experimental/test_tiling_qc.py +++ b/tests/experimental/test_tiling_qc.py @@ -160,6 +160,19 @@ def test_few_cells_below_k(self): for col in ["smoothed_cut_score", "is_outlier", "nhood_outlier_fraction"]: assert col in adata.obs.columns + def test_multiscale_scores_cells_absent_from_coarse_scales(self): + """Thin cells lost by downsampling are still scored at the requested scale.""" + from spatialdata import SpatialData + from spatialdata.models import Labels2DModel + + labels = np.zeros((128, 128), dtype=np.int32) + rows = range(9, 120, 8) + for lid, y in enumerate(rows, start=1): + labels[y, 10:118] = lid # 1-px-wide cells: gone at scale1/scale2 + sdata = SpatialData(labels={"labels": Labels2DModel.parse(labels, dims=("y", "x"), scale_factors=[2, 2])}) + adata = sq.experimental.tl.calculate_tiling_qc(sdata, labels_key="labels", scale="scale0", inplace=False) + assert set(adata.obs["label_id"]) == set(range(1, len(rows) + 1)) + def test_both_gates_disabled_raises(self, sdata_tile_boundary): sdata, _ = sdata_tile_boundary with pytest.raises(ValueError, match="At least one outlier gate"):