From f57b5fb1d4bc90a06f18773ac53d7fcc65f3e7dc Mon Sep 17 00:00:00 2001
From: anon <tim.treis@helmholtz-munich.de>
Date: Mon, 21 Sep 2026 23:21:01 +0200
Subject: [PATCH 1/3] fix: tile-invariant positions, no truncated or dropped
 cells in tiling

- Positional features (Center_*, BoundingBox*, Location_*, centroid-*)
  were relative to each tile's crop; report them in the labels' pixel
  grid, including the alignment crop offset.
- The auto margin assumed a cell reaches at most half its bbox from its
  centroid, truncating asymmetric cells at tile edges. Crop each tile to
  the union of its owned cells' bounding boxes plus 1 px instead.
- Multiscale labels computed centroids on the coarsest scale only, so
  cells that vanish there were never featurized or QC'd. Scan the
  requested scale.
---
 .../im/_calculate_image_features.py           | 109 +++++++++++-------
 src/squidpy/experimental/im/_tiling.py        |  85 ++++----------
 src/squidpy/experimental/tl/_tiling_qc.py     |  23 ++--
 .../test_calculate_image_features.py          |  79 ++++++++++++-
 tests/experimental/test_tiling.py             |  38 +-----
 tests/experimental/test_tiling_qc.py          |  13 +++
 6 files changed, 182 insertions(+), 165 deletions(-)

diff --git a/src/squidpy/experimental/im/_calculate_image_features.py b/src/squidpy/experimental/im/_calculate_image_features.py
index 4ae5de1f7..2a2170a42 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,
@@ -315,6 +314,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 +667,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 +735,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 +791,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 +829,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 +864,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,8 +877,8 @@ 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)
 
 
@@ -994,6 +1010,13 @@ 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, so its
+    values 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 +1072,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,15 +1097,15 @@ 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.
+    # 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.
     specs = build_tile_specs((H, W), cell_info, tile_size=tile_size, overlap_margin="auto")
     total_tiles = len(specs)
     logg.info(f"Tiling input into {total_tiles} tile(s) of size {tile_size} px.")
@@ -1092,10 +1115,12 @@ 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.
+        return _shift_positions(df, origin[0] + spec.crop[0], origin[1] + spec.crop[1])
 
     # cp_measure is GIL-bound, so kind="processes" (an active Client wins if set).
     results = _run_tiled(
diff --git a/src/squidpy/experimental/im/_tiling.py b/src/squidpy/experimental/im/_tiling.py
index d23720990..51aa0e73e 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,16 +177,6 @@ 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],
@@ -244,13 +192,15 @@ 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 1 pixel, so every owned cell is whole and has background
+        around its boundary (edge and radial features need it).  An integer
+        instead adds that fixed margin around the base region.
 
     Returns
     -------
@@ -261,8 +211,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 +232,17 @@ 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) - 1, 0)
+            cx0 = max(min(c.bbox_x0 for c in cells) - 1, 0)
+            cy1 = min(max(c.bbox_y0 + c.bbox_h for c in cells) + 1, height)
+            cx1 = min(max(c.bbox_x0 + c.bbox_w for c in cells) + 1, width)
+        else:
+            cy0 = max(by0 - margin, 0)
+            cx0 = max(bx0 - margin, 0)
+            cy1 = min(by1 + margin, height)
+            cx1 = min(bx1 + margin, 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..ad903a4dc 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,28 @@ 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_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 +948,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 +1121,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..882c789ef 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,
 )
@@ -404,42 +403,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"):

From 62513f063786cefd546ca125eeab90de83e2924d Mon Sep 17 00:00:00 2001
From: anon <tim.treis@helmholtz-munich.de>
Date: Mon, 21 Sep 2026 23:21:01 +0200
Subject: [PATCH 2/3] perf: assemble featurization output without full-table
 copies

Cast each tile to float32 and fill one preallocated, label-sorted array
instead of concat -> sort_index -> nunique -> cast. Constant columns are
found from per-tile min/max (same semantics) and never copied. ~3x lower
peak memory and ~27x faster on 1e8 values; output is bit-identical.
---
 .../im/_calculate_image_features.py           | 61 +++++++++++++------
 1 file changed, 44 insertions(+), 17 deletions(-)

diff --git a/src/squidpy/experimental/im/_calculate_image_features.py b/src/squidpy/experimental/im/_calculate_image_features.py
index 2a2170a42..9c4799a3b 100644
--- a/src/squidpy/experimental/im/_calculate_image_features.py
+++ b/src/squidpy/experimental/im/_calculate_image_features.py
@@ -882,6 +882,44 @@ def _compute_centroids(labels_da: xr.DataArray) -> dict[int, CellInfo]:
     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
 # ---------------------------------------------------------------------------
@@ -1119,8 +1157,9 @@ def _process_one(spec, image_da, labels_da):
             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.
-        return _shift_positions(df, origin[0] + spec.crop[0], origin[1] + spec.crop[1])
+        # 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(
@@ -1132,28 +1171,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,
@@ -1165,7 +1192,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
 

From ba04760a6d303b4117cf622eec73cd2ea80d7fa6 Mon Sep 17 00:00:00 2001
From: anon <tim.treis@helmholtz-munich.de>
Date: Mon, 21 Sep 2026 23:26:44 +0200
Subject: [PATCH 3/3] fix: align and pad tiles for cp_measure granularity

Granularity samples each crop on a 1/4 then 1/16 grid anchored at the
crop origin and removes background with an opening reaching ~320 px, so
tiled values drifted from an untiled run. When granularity is requested,
snap tile crops to its 16 px grid and pad them by 256 px. On real Xenium
DAPI nuclei, tiled-vs-untiled Spearman for bins 1-5 rises from
0.77-0.90 to 0.94-0.997; higher bins remain context-dependent
(documented).
---
 .../im/_calculate_image_features.py           | 29 ++++++++++++++++---
 src/squidpy/experimental/im/_tiling.py        | 27 ++++++++++++-----
 .../test_calculate_image_features.py          | 13 +++++++++
 tests/experimental/test_tiling.py             | 15 ++++++++++
 4 files changed, 73 insertions(+), 11 deletions(-)

diff --git a/src/squidpy/experimental/im/_calculate_image_features.py b/src/squidpy/experimental/im/_calculate_image_features.py
index 9c4799a3b..ee6f7c291 100644
--- a/src/squidpy/experimental/im/_calculate_image_features.py
+++ b/src/squidpy/experimental/im/_calculate_image_features.py
@@ -106,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"])
 
@@ -1003,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).
@@ -1052,8 +1062,11 @@ def calculate_image_features(
     ``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, so its
-    values vary with ``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
@@ -1144,7 +1157,15 @@ def calculate_image_features(
     # --- Tile ---
     # 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.
-    specs = build_tile_specs((H, W), cell_info, tile_size=tile_size, overlap_margin="auto")
+    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.")
 
diff --git a/src/squidpy/experimental/im/_tiling.py b/src/squidpy/experimental/im/_tiling.py
index 51aa0e73e..187d830f2 100644
--- a/src/squidpy/experimental/im/_tiling.py
+++ b/src/squidpy/experimental/im/_tiling.py
@@ -182,6 +182,8 @@ def build_tile_specs(
     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.
 
@@ -198,9 +200,15 @@ def build_tile_specs(
         Side length of the non-overlapping base grid cells.
     overlap_margin
         ``"auto"`` crops each tile to the union of its owned cells' bounding
-        boxes plus 1 pixel, so every owned cell is whole and has background
-        around its boundary (edge and radial features need it).  An integer
-        instead adds that fixed margin around the base region.
+        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
     -------
@@ -234,16 +242,21 @@ def build_tile_specs(
 
         if margin is None:
             cells = [cell_info[lid] for lid in owned]
-            cy0 = max(min(c.bbox_y0 for c in cells) - 1, 0)
-            cx0 = max(min(c.bbox_x0 for c in cells) - 1, 0)
-            cy1 = min(max(c.bbox_y0 + c.bbox_h for c in cells) + 1, height)
-            cx1 = min(max(c.bbox_x0 + c.bbox_w for c in cells) + 1, width)
+            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(
                 base=(by0, bx0, by1, bx1),
diff --git a/tests/experimental/test_calculate_image_features.py b/tests/experimental/test_calculate_image_features.py
index ad903a4dc..f8dbf2ce6 100644
--- a/tests/experimental/test_calculate_image_features.py
+++ b/tests/experimental/test_calculate_image_features.py
@@ -482,6 +482,19 @@ def test_tiled_vs_single_tile_equivalence(self, sdata_synthetic, features):
         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)
diff --git a/tests/experimental/test_tiling.py b/tests/experimental/test_tiling.py
index 882c789ef..fb44d370d 100644
--- a/tests/experimental/test_tiling.py
+++ b/tests/experimental/test_tiling.py
@@ -254,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):