Skip to content

Commit a2249d5

Browse files
lguerardclaude
andcommitted
fix(tiling): πŸ› drop the halo on a one-voxel-thick axis
`tile_shape: "auto"` with `do_3D: false` produces z=1 tiles, so the shipped `overlap: [4, 30, 30]` expanded every tile to 9 planes to keep 1 -- a 9x read amplification. Worse, the 2-D guard added in v2.0.0 then refused the block, because Cellpose without do_3D would read those planes as channels. So `auto` plus a 2-D method could not run at all. There is no context to gather along an axis the tile does not span, so normalize_overlap now takes the tile shape and drops the halo where the tile is a single voxel thick -- the behaviour the tiling guide already described ("z-tiles of size 1 get depth=0 in z even if you pass overlap=20") but that stage_tile never implemented. Explicit tiles with room keep their halo, and omitting tile_shape keeps the previous behaviour. Also fixes the amplification prepare reports: it did not clip the halo to the image the way stage_tile does, so a tile spanning the whole image claimed 5x while reading exactly 1x. That figure is meant to be tuned from. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
1 parent 58b68ea commit a2249d5

3 files changed

Lines changed: 90 additions & 6 deletions

File tree

β€Žsrc/patchworks/_distributed.pyβ€Ž

Lines changed: 33 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
from __future__ import annotations
1313

1414
import itertools
15+
import logging
1516
from pathlib import Path
1617
from typing import Callable, Sequence, Union
1718

@@ -21,10 +22,16 @@
2122
from ._io import zarr_compressor_kwargs
2223
from ._relabel import relabel_sequential_array
2324

25+
logger = logging.getLogger(__name__)
26+
2427
Overlap = Union[int, Sequence[int]]
2528

2629

27-
def normalize_overlap(overlap: Overlap, ndim: int) -> tuple[int, ...]:
30+
def normalize_overlap(
31+
overlap: Overlap,
32+
ndim: int,
33+
tile_shape: "Sequence[int] | None" = None,
34+
) -> tuple[int, ...]:
2835
"""Expand an overlap spec to one halo width per axis.
2936
3037
A scalar applies the same halo to every axis (the historical behaviour).
@@ -33,12 +40,21 @@ def normalize_overlap(overlap: Overlap, ndim: int) -> tuple[int, ...]:
3340
``76 x 1084 x 1084`` to keep ``16 x 1024 x 1024`` -- 5.3x more voxels than
3441
it uses, nearly all of it in z.
3542
43+
With *tile_shape*, an axis only one voxel thick gets **no** halo. There is
44+
no context to gather along an axis the tile does not span, and a 2-D
45+
method handed the extra planes would read them as channels. This is the
46+
``tile_shape: "auto"`` + ``do_3D: false`` case, where tiles come out one
47+
plane thick: a z-overlap of 4 would otherwise read 9 planes per tile to
48+
keep 1.
49+
3650
Parameters
3751
----------
3852
overlap : int or sequence of int
3953
Halo width, shared or per-axis.
4054
ndim : int
4155
Number of axes the halo is applied to.
56+
tile_shape : sequence of int, optional
57+
Tile extent per axis. Used to drop halos an axis has no room for.
4258
4359
Returns
4460
-------
@@ -55,6 +71,21 @@ def normalize_overlap(overlap: Overlap, ndim: int) -> tuple[int, ...]:
5571
)
5672
if any(o < 0 for o in values):
5773
raise ValueError(f"overlap must be non-negative, got {values}")
74+
75+
if tile_shape is not None:
76+
clipped = tuple(
77+
0 if int(t) <= 1 else o for o, t in zip(values, tile_shape)
78+
)
79+
if clipped != values:
80+
dropped = [
81+
i for i, (a, b) in enumerate(zip(values, clipped)) if a != b
82+
]
83+
logger.info(
84+
"dropping the halo on axis %s: the tile is 1 voxel thick "
85+
"there, so there is no context to read.",
86+
dropped,
87+
)
88+
values = clipped
5889
return values
5990

6091

@@ -173,7 +204,7 @@ def stage_tile(
173204
"""
174205
shape = image.shape
175206
sl = spatial_tiles(shape, tile_shape)[index]
176-
halo = normalize_overlap(overlap, len(sl))
207+
halo = normalize_overlap(overlap, len(sl), tile_shape=tile_shape)
177208
expanded, trims = [], []
178209
for s, dim, ov in zip(sl, shape, halo):
179210
lo = max(0, s.start - ov)

β€Žtests/test_distributed.pyβ€Ž

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,53 @@ def test_normalize_overlap_scalar_and_per_axis():
4141
normalize_overlap([-1, 0, 0], 3)
4242

4343

44+
def test_no_halo_on_a_one_voxel_thick_axis():
45+
"""A 1-plane tile must get no halo on that axis.
46+
47+
`tile_shape: "auto"` with `do_3D: false` produces z=1 tiles, so a z-overlap
48+
of 4 would read 9 planes to keep 1 (9x) -- and hand a 2-D method a stack it
49+
would read as channels.
50+
"""
51+
assert normalize_overlap([4, 30, 30], 3, tile_shape=(1, 8000, 8000)) == (
52+
0,
53+
30,
54+
30,
55+
)
56+
# A tile with room keeps its halo.
57+
assert normalize_overlap([4, 30, 30], 3, tile_shape=(16, 1024, 1024)) == (
58+
4,
59+
30,
60+
30,
61+
)
62+
# Without a tile_shape nothing is clipped (the historical behaviour).
63+
assert normalize_overlap([4, 30, 30], 3) == (4, 30, 30)
64+
65+
66+
def test_single_plane_tiles_segment_without_a_z_halo(tmp_path):
67+
"""End to end: z=1 tiles must not pull neighbouring planes into the tile."""
68+
img = np.zeros((4, 32), "uint16")
69+
img[1, 4:10] = 500 # only plane 1 has signal
70+
tile = (1, 16)
71+
72+
stage = str(tmp_path / "stage.zarr")
73+
create_stage(stage, img.shape, tile)
74+
seen = []
75+
76+
def _recording(block):
77+
seen.append(block.shape)
78+
return _fn(block)
79+
80+
for i in range(len(spatial_tiles(img.shape, tile))):
81+
stage_tile(img, _recording, stage, i, tile_shape=tile, overlap=[3, 2])
82+
83+
assert all(s[0] == 1 for s in seen), (
84+
f"a z=1 tile was handed multiple planes: {seen}"
85+
)
86+
written = np.asarray(zarr.open_group(stage, mode="r")["staged"])
87+
assert (written[0] == 0).all(), "plane 0 has no signal and must stay empty"
88+
assert (written[1, 4:10] > 0).all()
89+
90+
4491
def test_auto_overlap_anisotropy_shrinks_z():
4592
"""A 20x coarser z step needs 20x fewer planes for the same distance."""
4693
assert auto_overlap(30) == 30 # unchanged without voxel_size

β€Žworkflow/scripts/prepare_tiles.pyβ€Ž

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -85,17 +85,23 @@
8585
# A halo wider than the tile itself is not boundary context -- it means every
8686
# tile re-reads and re-segments its neighbours. Catch it here rather than
8787
# paying 5x the GPU time for results that get trimmed away.
88-
overlap = normalize_overlap(cfg.get("overlap", 0), len(tile_shape))
88+
overlap = normalize_overlap(
89+
cfg.get("overlap", 0), len(tile_shape), tile_shape=tile_shape
90+
)
8991
for axis, (ov, extent) in enumerate(zip(overlap, tile_shape)):
9092
if ov >= extent:
9193
raise ValueError(
9294
f"overlap[{axis}]={ov} >= tile_shape[{axis}]={extent}: each tile "
9395
f"would read past its neighbours. Use a per-axis overlap, e.g. "
9496
f"overlap: {list(max(1, t // 4) for t in tile_shape)}"
9597
)
96-
amplification = np.prod(
97-
[(t + 2 * o) for t, o in zip(tile_shape, overlap)]
98-
) / np.prod(tile_shape)
98+
# Clip the halo to the image the way stage_tile does, or a tile that already
99+
# spans an axis reports a halo it will never actually read -- a tile covering
100+
# the whole image would claim 5x while reading exactly 1x.
101+
read = np.prod(
102+
[min(t + 2 * o, s) for t, o, s in zip(tile_shape, overlap, image.shape)]
103+
)
104+
amplification = read / np.prod(tile_shape)
99105
print(f"[patchworks] halo read amplification: {amplification:.2f}x")
100106

101107
tiles = spatial_tiles(image.shape, tile_shape)

0 commit comments

Comments
Β (0)