docs: 📝 document the pipeline changes and ship the optimized defaults

Brings the documentation in line with the eight preceding commits, and turns
the optimizations on by default rather than leaving them opt-in.

Docs:
- skip_empty: separates the fast-but-approximate estimate_empty_tiles preview
  from the exact build_occupancy_map/tile_occupancy pair, and says plainly
  that the centred window covers 6.25% of a [16, 1024, 1024] tile.
- tiling: per-axis overlap with the amplification arithmetic, and
  auto_overlap(voxel_size=...) for deriving it from the physical anisotropy.
- gpu_distributed: corrects the "at most half the available VRAM" claim (only
  auto_tile_shape does that; the Cellpose estimator uses its own model),
  documents the headroom, the CUDA_VISIBLE_DEVICES-aware device pick, and the
  shared-GPU OOM retry.
- merging: the algorithm is now five steps, with the offsets-from-counts step
  and the sequential renumber folded into the LUT.
- performance: "machine-aware" was the bug — it is allocation-aware now.
- snakemake: tiles_per_job, concurrent multi-runs, retries/attempt-scaled
  memory, and the merge writing level 0 in place.
- api: pages for the newly public helpers.

Defaults: tiles_per_job: 4 in every shipped config, so the model load is
amortized out of the box; 1 restores one job per tile.

Verified: pytest, ruff, markdownlint and `mkdocs build --strict` all clean,
plus an end-to-end concurrent multi-config run.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: lguerard

README.md

@@ -286,9 +286,12 @@ tiles where the dask-image approach stalls.
 | ------------------------------ | ----------------------------------------- | ------------------------------------------------------------- |
 | In-process Dask client         | `FutureCancelledError: lost dependencies` | Detected at startup, raises immediately with fix instructions |
 | 3-4× fn recompute during merge | Cellpose runs 3× per tile                 | Staging writes labels once, merge reads from disk             |
-| O(n²) sequential relabelling   | Graph construction hangs at 1000+ tiles   | Linear post-pass O(voxels) via `np.unique` + LUT              |
+| O(n²) sequential relabelling   | Graph construction hangs at 1000+ tiles   | Folded into the merge's own LUT — no extra pass over the volume |
 | Wrong overlap boundary         | Output shape mismatch                     | Always uses `boundary="none"`                                 |
 | Persisting large arrays        | Worker OOM                                | Never persists; keeps dask graph lazy and streams             |
+| Sizing work to the whole node  | Job OOM-killed on a shared cluster node   | Reads the SLURM/cgroup allocation, not `os.cpu_count()`       |
+| Rechunking a pyramid level     | Threaded scheduler stockpiles intermediates | Levels stream one source chunk per task, bounded by construction |
+| Isotropic halo on flat tiles   | 5× the voxels read and segmented, then trimmed | `overlap` takes one width per axis                       |
 
 ---
 

docs/api/chunks.md

@@ -3,3 +3,17 @@
 ::: patchworks.auto_tile_shape
 
 ::: patchworks.auto_tile_shape_cellpose
+
+::: patchworks.auto_overlap
+
+::: patchworks.normalize_overlap
+
+## Cluster resource detection
+
+On a shared node the machine's core count and free RAM say nothing about what
+this job was granted. These read the allocation instead, and everything that
+sizes a worker pool goes through them.
+
+::: patchworks.cpu_allocation
+
+::: patchworks.safe_worker_count

docs/api/io.md

@@ -2,4 +2,18 @@
 
 ::: patchworks.load_ome_zarr
 
+## Deciding which tiles hold signal
+
+`estimate_empty_tiles` is a fast **preview** — it samples a centred window per
+tile, so it can miss signal at a tile's edge. `build_occupancy_map` +
+`tile_occupancy` are **exact**: a brick maximum exceeds the threshold exactly
+when some voxel in that brick does. Use the latter pair when the result is
+used as a skip list. See [Skipping empty tiles](../guide/skip_empty.md).
+
 ::: patchworks.estimate_empty_tiles
+
+::: patchworks.build_occupancy_map
+
+::: patchworks.tile_occupancy
+
+::: patchworks.auto_empty_threshold

docs/guide/gpu_distributed.md

@@ -73,12 +73,35 @@ FutureCancelledError: lost dependencies
 ## GPU memory sizing
 
 When `use_gpu=True`, patchworks queries free GPU VRAM via `nvidia-ml-py`
-(install: `pip install "patchworks[gpu]"`). Tile size is set so each tile
-uses at most half the available VRAM.
+(install: `pip install "patchworks[gpu]"`) and keeps 20% of it as headroom,
+because `info.free` is a point-in-time reading of a device you usually do not
+own outright. `auto_tile_shape` then sizes each tile to at most half of that
+budget; `auto_tile_shape_cellpose` uses Cellpose's own memory model instead
+(roughly 20× the raw tile bytes, plus ~2 GiB for the model).
+
+The device is resolved from `CUDA_VISIBLE_DEVICES`. This matters on
+multi-GPU nodes: NVML enumerates **every** GPU regardless of `--gres=gpu:1`,
+so querying index 0 unconditionally would read a different card's free memory
+than the one your job was granted.
 
 Without `nvidia-ml-py`, a conservative 8 GiB default is used with a warning.
 Install it for accurate sizing on large-VRAM cards (A100, H100):
 
 ```bash
 pip install "patchworks[gpu]"
 ```
+
+!!! tip "Sizing for a GPU you can't see"
+    The Snakemake workflow plans tiles in `prepare`, which runs on a CPU
+    node, so no GPU can be queried there. Set `gpu_memory_gb` in the config
+    to the segment GPU's VRAM (e.g. `24` for an RTX 4090) instead of relying
+    on the fallback.
+
+## Surviving a shared GPU
+
+An out-of-memory error on a shared device is often transient — a co-tenant
+job grew, not a tile that doesn't fit. patchworks retries on the GPU with a
+backoff rather than falling back to the CPU, where a single tile can take
+over an hour. Before each backoff it releases its own device memory
+(including any cached Cellpose model), since holding that is exactly what
+starves the other job. Both torch and cupy OOM errors are recognised.

docs/guide/merging.md

@@ -38,22 +38,40 @@ tile_process calls fn once per tile → staged zarr
                          merge reads from staged zarr (no fn calls)
 ```
 
-### Step 2: boundary scan
+### Step 2: make the ids globally unique
+
+Tiles write local `1..n`, which collide, so the boundary scan could not
+otherwise tell two different objects apart. If each tile's label **count** is
+known, this is just an exclusive cumulative sum — global id is
+`offset[tile] + local`, computed in `O(n_tiles)` with no read of the volume
+at all. `stage_tile` returns that count for exactly this purpose; pass the
+counts as `label_counts=`.
+
+Without counts, the merge falls back to streaming every chunk and renumbering
+it in place — correct, but a full read **and write** of the volume.
+
+### Step 3: boundary scan
 
 Only the two voxels on either side of each tile boundary are read. For any
-pair of touching non-zero labels `(a, b)`, they must be the same object.
+pair of touching non-zero labels `(a, b)`, they must be the same object. The
+per-tile offsets are applied here, on the fly.
 
 I/O cost: `O(n_boundaries × face_area)`, not `O(full_volume)`.
 
-### Step 3: connected components
+### Step 4: connected components
 
 scipy sparse connected components on the touching pairs produces a relabeling
 lookup table. All labels that transitively touch each other are mapped to the
 same canonical label.
 
 Cost: `O(n_touching_pairs)`.
 
-### Step 4: parallel relabel
+With `sequential_labels=True` the contiguous renumbering is folded into this
+same LUT. Because the id domain is dense by construction, the surviving ids
+are exactly the distinct LUT values — a `np.unique` over an array the length
+of the object count, with no scan of the volume.
+
+### Step 5: parallel relabel
 
 The LUT is applied to every tile in parallel via `multiprocessing.Pool`. The
 LUT is shared via process initializer to avoid re-pickling it for every chunk
@@ -86,17 +104,20 @@ merged = merge_tile_labels(
 
 ## Sequential label numbering
 
-By default, merged labels are globally unique but may be **gappy**
-(block-encoded IDs like 1, 2, 500001, 500002, …). This is fine for
-counting, `regionprops`, and measurement — the IDs just aren't consecutive.
+By default, merged labels are globally unique but may be **gappy** — boundary
+merging fuses ids, leaving holes where the absorbed ones were. This is fine
+for counting, `regionprops`, and measurement — the IDs just aren't
+consecutive.
 
 For contiguous 1..N numbering, use `sequential_labels=True`:
 
 ```python
 tile_process("image.zarr", fn, write_to="labels.zarr", sequential_labels=True)
 ```
 
-This runs a cheap linear post-pass: `np.unique` + lookup-table remap, O(voxels).
+This is free: it composes into the relabel LUT the merge already applies, so
+it costs a `np.unique` over the object count rather than another pass over
+the volume.
 
 !!! warning "Do not use dask's built-in sequential relabel"
     `dask_image.ndmeasure.merge_labels_across_chunk_boundaries` has a

docs/guide/performance.md

@@ -3,19 +3,27 @@
 `tile_process` is built so a run **adapts to whatever machine it lands on** and
 can't run out of RAM/VRAM or freeze the box — without you tuning anything.
 
-## Automatic, machine-aware concurrency
+## Automatic, allocation-aware concurrency
 
 The staging step (running your `fn` once per tile to a temp store) and the
-merge step are sized to the host automatically:
+merge step are sized automatically:
 
 - **GPU** (`use_gpu=True`) → **one tile at a time**, so concurrent evaluations
   can never exhaust VRAM.
 - **CPU** → as many tiles in flight as fit **80 % of available RAM** (estimated
   from the tile size), and always **leaving one core free** so the machine
   stays responsive — it never pins every core.
 
-The RAM figure is read live via `psutil`; without it, a conservative default is
-used instead of guessing high.
+"Available" means available **to this process**, not to the machine. On a
+shared cluster node those differ wildly — a 32-core, 128 GB job on a 128-core,
+512 GB node would otherwise size itself for the whole box and get OOM-killed
+while its own accounting said it had room. patchworks takes the smallest of
+`SLURM_MEM_PER_NODE`, `SLURM_MEM_PER_CPU × cpus`, the cgroup limit (the one
+that actually triggers the kill) and `psutil`'s free RAM, and reads the core
+count from `SLURM_CPUS_PER_TASK` or the process' CPU affinity mask.
+
+Without any of those signals, a conservative default is used instead of
+guessing high.
 
 ## Live progress dashboard (GPU runs)
 

docs/guide/skip_empty.md

@@ -37,7 +37,9 @@ tile_process(
 )
 ```
 
-## How `estimate_empty_tiles` works
+## Two ways to decide
+
+### `estimate_empty_tiles` — a fast preview
 
 For each tile in the grid, only a small centred **sample window** is read
 (default: 24×256×256 voxels). If the maximum value in that window exceeds
@@ -47,11 +49,37 @@ This is **bounded I/O**: the total data read is `n_tiles × sample_window`,
 not the full image. For a 2200-tile image with the default window, this reads
 ≈ 30 MB instead of 250 GB — and it runs in seconds.
 
-!!! warning "Approximate"
-    `estimate_empty_tiles` inspects only the tile centre. Signal confined to
-    a tile's edge can be missed. The actual `tile_process` run always inspects
-    the **full tile** max inline — so no objects are ever dropped in the real
-    run, only in the preview.
+!!! warning "Approximate — a preview, not a skip list"
+    Only the tile centre is inspected. On a `(16, 1024, 1024)` tile the
+    default window covers **6.25% of the tile's area**, so an object in the
+    outer ring is invisible. `tile_process` re-tests the **full tile** max
+    inline, so nothing is dropped in a real run — but anything that uses this
+    result *as* the skip list would drop those tiles for good.
+
+### `build_occupancy_map` + `tile_occupancy` — exact
+
+Reduces every `block`-sized brick of the image to its **maximum** and stores
+that beside `image.zarr` (about 1/16384 of the image at the default 128 px
+block). `tile_occupancy` then reduces the map over each tile's full footprint.
+
+```python
+from patchworks import auto_empty_threshold, build_occupancy_map, tile_occupancy
+
+build_occupancy_map("image.zarr")             # once per image, all channels
+info = tile_occupancy(
+    "image.zarr", TILE, channel=0, threshold=auto_empty_threshold(img, 0, 0)
+)
+```
+
+This is **exact, not approximate**: `block_max > threshold` is true exactly
+when some voxel in that block exceeds the threshold, so comparing pooled
+maxima against a threshold derived from raw voxels answers the same question
+as scanning every voxel. Max-pooling cannot lose a bright voxel.
+
+The map is built once and shared by every segmentation reading that image, so
+a three-config run pays for one pooling pass instead of three sampling passes.
+This is what the Snakemake workflow uses, and it builds the map on first use
+for stores converted before the map existed.
 
 ## Threshold selection