feat(pyramid): pyramids on images and label images; mantis-v2 adds them after tracking and QC (#359)

.claude/skills/reconstruct-with-nextflow/SKILL.md

@@ -3,10 +3,10 @@ name: reconstruct-with-nextflow
 description: >-
   Run a mantis-v2 microscope dataset through the biahub Nextflow reconstruction
   pipeline (nextflow/mantis-v2.nf: flat-field → deskew → reconstruct →
-  virtual-stain → assemble → project → segment → track → QC) on the Bruno HPC
-  cluster. Assemble, project, segment, track and QC are optional and selected per
+  virtual-stain → assemble → project → segment → track → pyramid → QC) on the Bruno
+  HPC cluster. Assemble, project, segment, track and QC are optional and selected per
   dataset family: A549 runs all of them, neuromast runs assemble + QC and no
-  segmentation or tracking. Locates the raw
+  segmentation or tracking; pyramids are on whenever assemble runs. Locates the raw
   acquisition under /hpc/instruments/cm.mantis, picks the output project
   directory, scaffolds configs, launches the run in a tmux session in the
   foreground so it can be watched, monitors it unattended every 30 minutes
@@ -358,8 +358,14 @@ Do not run anything yet. Show the user:
 
    | family | steps |
    |---|---|
-   | A549 / cell-line / organelle | flat-field → deskew → reconstruct → virtual-stain → **assemble → project → segment → track → QC** |
-   | neuromast / zebrafish / dynatrack | flat-field → deskew → reconstruct → virtual-stain → **assemble → QC** (no projection, segmentation or tracking) |
+   | A549 / cell-line / organelle | flat-field → deskew → reconstruct → virtual-stain → **assemble → project → segment → track → pyramid → QC** |
+   | neuromast / zebrafish / dynatrack | flat-field → deskew → reconstruct → virtual-stain → **assemble → pyramid → QC** (no projection, segmentation or tracking) |
+
+   Pyramid is on by default whenever assemble runs: it adds two coarser levels
+   (3 in total, DCA v0.2) to every image and label image of the assembled and
+   projected plates, in place, once segmentation and tracking have finished
+   writing into them. QC is the very last step and waits for the pyramids.
+   `--pyramid_levels 1` skips them.
 
    Write the list out in the plan rather than naming the family, and say which
    optional steps are being skipped and why. If the user asks to skip anything

.claude/skills/reconstruct-with-nextflow/references/caveats.md

@@ -265,30 +265,30 @@ Assemble is one large SLURM job — the single longest step.
 `biahub concatenate --resume` (passed by the Nextflow task) makes a late kill
 recoverable per (t, c) unit.
 
-## 7. Shard geometry and `shards_ratio`
-
-`shards_ratio` is a ratio over the **chunk** shape
-(`shards = chunks * shards_ratio`, elementwise TCZYX), not an absolute shard
-shape — changing T without recomputing ZYX silently changes the spatial extent.
-Larger T shards mean a larger RAM request and a worse blast radius when one
-tears; leave the T ratio alone without a measured reason. Channels cannot be
-sharded — keep C at 1. The same two keys exist in `project.yml`: the A549
-template uses one-frame `256 x 256` chunks (the assembled and tracking grid)
-and `shards_ratio: [16, 1, 1, 7, 5]`, 16-timepoint full-frame shards; its YX
-entries are `ceil(Y/256), ceil(X/256)` and must be recomputed for a different
-frame size, exactly as for `concatenate.yml`. Auto-sizing is
-[iohub#458](https://github.com/czbiohub-sf/iohub/issues/458); most of this
-section goes away when it lands.
-
-**A T ratio above 1 needs [iohub#460](https://github.com/czbiohub-sf/iohub/pull/460),
+## 7. Shard geometry and `shards`
+
+`shards` in `concatenate.yml` and `project.yml` is a **size target**
+(`shards: 5GB`, `shards: 200MB`): iohub builds the largest chunk-aligned shard
+under it, filling X, Y and Z first and then timepoints, and never spans
+channels ([iohub#461](https://github.com/czbiohub-sf/iohub/pull/461)). Nothing
+is recomputed per dataset: the same `5GB` gives 5 timepoints by the full volume
+on the A549 plate and 9 on the dynatrack one. `"XY"` / `"XYZ"` (one shard per
+plane / volume) and an explicit TCZYX shape in array elements are also accepted.
+The old `shards_ratio` (a multiplier over the chunk shape) is refused with a
+message naming the new key. A larger target means more timepoints per shard,
+hence a larger RAM request (concatenate reads the T extent back from the plate it
+initialised to size its workers) and a worse blast radius when one tears; leave
+the templates' values alone without a measured reason.
+
+**Shards spanning several timepoints need [iohub#460](https://github.com/czbiohub-sf/iohub/pull/460),
 which shipped in iohub 0.3.11** — the floor `pyproject.toml` already requires, so a
 synced checkout has it.
 
 Before that fix, `concatenate` aborted with a message that points nowhere near the
 cause — `numpy ... Unable to allocate 651. TiB for an array with shape
 (1048576, 86, 1664, 1193)`, whose first axis is not T, not a shard, not
 anything. It is **not** evidence that the shard buffer is too large, and
-unsetting `shards_ratio` "fixes" it for the wrong reason. What happens: a
+unsetting `shards` "fixes" it for the wrong reason. What happens: a
 shard-aligned batch writes a whole shard's worth of timepoints in one call,
 iohub drops the blank ones (`Skipping t=1, c=0 due to all zeros or nans` in the
 step's `.command.out`), and the surviving indices are left with a *gap* — which
@@ -440,3 +440,17 @@ Consequences worth knowing:
   differ" means nothing here; `corr` and max abs diff are the measures that do. A
   real virtual-stain regression moves the distribution, not the last bit. Deskew,
   reconstruct and assemble are fully deterministic.
+
+## 15. Pyramids follow tracking and precede QC, and are not a resume unit
+
+`biahub pyramid` rewrites each position's `zarr.json` (the multiscales list and
+its `biahub-pyramid` record), as do tracking and imaging-qc, so the pipeline
+runs it after tracking and holds QC, the very last step, until it has finished. A run that fails after the pyramid step
+left the plates with 3 levels; a relaunch with `-resume` reuses the cached
+pyramid tasks, and a relaunch that recomputes earlier steps rewrites level 0
+only — rerun `biahub pyramid -i <plate>/*/*/* --levels 3` by hand in that
+case, since the coarser levels are not recomputed by anything else. A plate
+that already carries a pyramid of another depth is refused at init
+(`already has a pyramid of a different depth`); the fix is iohub's
+`delete_pyramid` on each position or matching `--pyramid_levels`.
+

.claude/skills/reconstruct-with-nextflow/templates/run_mantis_v2.sh

@@ -204,6 +204,11 @@ unset CLAUDECODE
 # --qc_track_config was removed: qc_project.yaml counts the labels tracking writes
 # into the projected plate, so passing the old flag errors at launch.
 #
+# Pyramids are not a flag here: the assembled and projected plates get 3
+# resolution levels (images mean-pooled, labels max-pooled; DCA v0.2, biahub#359)
+# after tracking by default, and QC runs last, after them. Add
+# `--pyramid_levels 1` to skip them.
+#
 # TO SKIP A STEP, DELETE ITS LINE. Do not comment it out: `#` inside a
 # backslash-continued command does not start a comment line — the continuation
 # swallows it, every flag below it is dropped including `-resume`, and bash then

README.md

@@ -99,7 +99,7 @@ Some steps include an `estimate-*` command that writes a YAML of parameters you
 | | |
 | `flat-field` | Apply flat-field correction to selected channels |
 | `flip` | Flip images in a dataset |
-| `pyramid` | Create multiscale pyramid levels |
+| `pyramid` | Compute multiscale pyramids of images and label images, in place |
 | `process-with-config` | Apply arbitrary YAML-defined functions to a dataset |
 | | |
 | `estimate-deskew` | Estimate deskewing parameters |

biahub/cli/main.py

@@ -93,7 +93,7 @@ def format_options(self, ctx, formatter):
     {
         "name": "pyramid",
         "import_path": "biahub.pyramid.pyramid_cli",
-        "help": "Create pyramid levels for a dataset",
+        "help": "Compute multiscale pyramids of a plate's images and label images, in place",
     },
     {
         "name": "apply-transform",

biahub/concatenate.py

@@ -651,7 +651,7 @@ def _resolve_concatenate_inputs(
     output_metadata = {
         "shape": (len(input_time_indices), len(all_channel_names)) + tuple(cropped_shape_zyx),
         "chunks": chunk_size,
-        "shards_ratio": settings.shards_ratio,
+        "shards": settings.shards,
         "version": resolve_ome_zarr_version(
             all_data_paths[0], settings.output_ome_zarr_version
         ),
@@ -793,7 +793,12 @@ def concatenate(
     # 16-worker tasks took 7-26 min; 0.15 min/volume budgets 60 min here.
     T_out, C_out, _, _, _ = prep["output_metadata"]["shape"]
     _, _, Z, Y, X = prep["shape"]
-    batch_size = settings.shards_ratio[0] if settings.shards_ratio else 1
+    # iohub writes one shard's timepoints per unit, so the shard's T extent is the
+    # batch a worker holds. Read it back from the plate init just created rather than
+    # predicting it: the config may give the shards as a size target.
+    with open_ome_zarr(output_dirpath, layout="hcs", mode="r") as output_plate:
+        _, first_position = next(output_plate.positions())
+        batch_size = (first_position.data.shards or first_position.data.chunks)[0]
     time_minutes, num_cpus, gb_ram_per_cpu = estimate_resources(
         shape=(max(T_out // batch_size, 1), C_out, Z, Y, X),
         ram_multiplier=8 * batch_size,

biahub/project.py

@@ -221,7 +221,9 @@ def _init_output_plate(
         ),
         dtype=dtype,
         chunks=(1, *settings.chunks_czyx) if settings.chunks_czyx else None,
-        shards_ratio=tuple(settings.shards_ratio) if settings.shards_ratio else None,
+        shards=tuple(settings.shards)
+        if isinstance(settings.shards, list)
+        else settings.shards,
         metadata_sources=input_plate,
         metadata_keys=PROVENANCE_METADATA_KEYS,
         extra_metadata={PROVENANCE_KEY: settings.model_dump(mode="json")},