From ec7ebb7c0b789150ad46661029b78f090d293599 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 15:28:52 +0000 Subject: [PATCH 1/5] =?UTF-8?q?Find=20every=20image's=20masks=20and=20desc?= =?UTF-8?q?ribe=20how=20they=20relate=20(plan=20=C2=A712=20Q29)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Masks of any number of types are paired with the images: a masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, plus the label images of OME-Zarr images. `check` and `run` check them from their headers only (2D, the image's Y and X, integer or boolean pixels) and a run records them in the manifest's optional `masks` key. The optional mask step (`masks` command, `scan_masks`) reads the pixels one image at a time and writes mask_relationships.json: each type's kind (instances, grouped, one_value, empty), each ordered pair's relation (same_objects, grouped_by, inside, none) with the counts behind it, such as parents with 0, 1, 2 and 3+ children, and the hierarchy of direct parents (nucleoli inside nuclei inside cells). SciPy becomes a core dependency for counting pieces. Contract §13, plan Q29, user guide. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01U6vJfhGWQ7pTpsJhw48Zf7 --- AGENTS.md | 15 +- README.md | 8 +- docs/contract.md | 87 +++- docs/plan.md | 26 ++ docs/user_guide.md | 93 +++++ pyproject.toml | 2 + src/morphobricks_preprocessing/__init__.py | 23 ++ src/morphobricks_preprocessing/cli.py | 67 ++- src/morphobricks_preprocessing/manifest.py | 20 +- src/morphobricks_preprocessing/mask_scan.py | 409 +++++++++++++++++++ src/morphobricks_preprocessing/masks.py | 428 ++++++++++++++++++++ src/morphobricks_preprocessing/report.py | 97 ++++- src/morphobricks_preprocessing/run.py | 32 +- tests/test_mask_scan.py | 247 +++++++++++ tests/test_masks.py | 239 +++++++++++ 15 files changed, 1780 insertions(+), 13 deletions(-) create mode 100644 src/morphobricks_preprocessing/mask_scan.py create mode 100644 src/morphobricks_preprocessing/masks.py create mode 100644 tests/test_mask_scan.py create mode 100644 tests/test_masks.py diff --git a/AGENTS.md b/AGENTS.md index 14b9a4e5..519c9fab 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -89,7 +89,13 @@ src/morphobricks_preprocessing/ manifest.py RunManifest / RunSettings / FileEntry, JSON read/write, SHA-256 of the CSV run.py preflight(), run_experiment() (incremental), sample_planes() for previews report.py plain-language summaries shared by CLI and notebook, group statistics tables - cli.py, __main__ `morphobricks-preprocessing check|run` (also `python -m morphobricks_preprocessing`) + cli.py, __main__ `morphobricks-preprocessing check|run|masks` (also `python -m morphobricks_preprocessing`) + masks.py masks (Q29): find_masks pairs every image with its masks (masks folder with one + folder per type mirroring the experiment, OME-Zarr label images), header-only + checks (2D, image YX, integer/bool), MaskInventory (to_manifest), read_mask + mask_scan.py the optional mask step (Q29): piece_counts (scipy.ndimage), pair_counts, + scan_masks -> MaskReport (kinds, relations, direct-parent hierarchy, evidence), + write_mask_report (mask_relationships.json); fixed THRESHOLDS zarr_input.py OME-Zarr input: open_ngff (OME-NGFF 0.4/0.5 metadata -> level-0 array, axis letters and physical sizes from the scale and axis units), ngff_label_names/open_ngff_labels (label images, for later bricks), @@ -388,7 +394,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | Topic | Decision | |---|---| | Binning | Per-plane min–max, same bin count per experiment, pooled counts for groups (binning study) | -| Core dependencies | `numpy`, `tifffile` and `morphobricks-normalization` (git dependency pinned to its release tag v0.2.0; `[tool.hatch.metadata] allow-direct-references`). `ipywidgets` and `matplotlib` in the `notebook` extra; `zarr` in the `zarr` extra; `pytest` and `ruff` in `dev` | +| Core dependencies | `numpy`, `tifffile`, `scipy` (the mask step's piece counts, Q29) and `morphobricks-normalization` (git dependency pinned to its release tag v0.2.0; `[tool.hatch.metadata] allow-direct-references`). `ipywidgets` and `matplotlib` in the `notebook` extra; `zarr` in the `zarr` extra; `pytest` and `ruff` in `dev` | | OME-Zarr (Q12) | Input may be TIFF and/or OME-Zarr (folders ending in `.zarr`), mixed. OME-Zarr axes come from their metadata and are matched by role; level 0 of a single image only, several images per store is `multiple_series`. TIFF choice once per experiment (`convert_tiffs`): read directly (default) or convert to `/ome_zarr` (OME-NGFF 0.4) and compute rows from the copies; switching recomputes TIFF rows and removes unused copies. Reads OME-NGFF 0.4 and 0.5 (0.5 needs zarr 3). In the notebook: its own step-3 card "(Optional) Image conversion" after Settings (Q13): radio buttons (read directly = recommended), a disk-space warning (amber, red when the estimate exceeds free space, never blocking), and conversion settings (chunk size, compression + level, pyramid levels) shown only when converting. CLI `--convert-tiff` with `--chunk-size`, `--compression`, `--compression-level`, `--pyramid-levels`. No overlap setting: chunks never overlap | | System files (Q5) | Ignored silently; `single` check and sibling uniqueness are case-insensitive | | Too-deep folders (Q7) | Error | @@ -408,6 +414,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | Dataset consistency (Q21, P2) | Step 1 "Dataset consistency" panel: balance table (images and technical-repeat folders per condition x biological repeat), `unbalanced_images` (< 0.5 or > 2 x median), `repeat_missing_conditions`, `different_pixel_sizes`, `missing_pixel_size` (OME-Zarr scale included), `different_plane_counts`; `different_sizes` narrowed to width/height; consistency issues shown only in that panel, not in Messages. Warnings only; outputs unchanged | | Experiment design (Q22, P3) | Step 1 after the labels: pairing dropdowns, one yellow "Yes, they are comparable" question per repeat group (single whole-experiment question without biological-repeat folders), control checkboxes; `experiment_design.csv` written after every answer; reruns load it and fold the questions behind a Change button; step 3 gated. CLI `run` reads it or `--design`, else writes a template and exits 1. `run_experiment(design=...)` optional; manifest `experiment_design_sha256` (optional key, version 1). Two folders of one condition in a group block that group | | Position identity (Q20, D02) | `position_id` = image name without extension, in every CSV row and manifest file entry; FOV key = labels + `position_id`; duplicates in one leaf folder are the error `duplicate_position_id`. CSV columns and schema 2 unchanged; manifest version 1 with an optional key | +| Masks (Q29) | Optional masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, paired by leaf folder + `position_id`; OME-Zarr label images are masks too. Header-only errors (`masks_folder_not_found`, `masks_inside_experiment`, `mask_names_differ_in_case`, `duplicate_mask`, `unreadable_mask`, `mask_not_2d`, `mask_size_mismatch`, `mask_pixel_type`) and warnings (`mask_missing`, `mask_without_image`, `mask_outside_type_folder`) in `check`/`run`; manifest key `masks` (optional, version 1). Optional mask step (`masks` command, `scan_masks`): kinds, relations with evidence, direct-parent hierarchy, `mask_relationships.json` (report version 1); fixed thresholds; suggestions only, the user names and confirms. SciPy core dependency. Notebook step later. To confirm: 2D masks, integer/bool pixels, fixed thresholds | | Notebook style | Follow the lab's NucleiSky and LabChronicle notebooks (see §7) | | Widget widths | **99%**, never 100% | | Paths and filenames in widget HTML and notebook text | **`…` / `**…**`, never `` or inline code** | @@ -668,6 +675,7 @@ ruff format --check . pytest morphobricks-preprocessing check /data/exp1 morphobricks-preprocessing run /data/exp1 --axes ZCYX +morphobricks-preprocessing masks /data/exp1 --masks /data/exp1_masks jupyter lab notebooks/MorphoBricks_Preprocessing.ipynb ``` @@ -695,6 +703,9 @@ jupyter lab notebooks/MorphoBricks_Preprocessing.ipynb - **Released as v0.1.0** (27 October 2026), with everything above; the normalization brick is pinned to its v0.2.0 (question 28). The next tag, v0.2.0, follows follow-ups P4 and P5. +- **Masks (question 29, in review):** pairing every image with its masks, header checks, + the manifest's `masks` key, and the mask step with `mask_relationships.json`; library, + command line and documentation (the notebook step later). - **Next (not started):** follow-ups P4–P5 of the normalization plan (image type, intensity checks and plots). - **Not yet verified:** diff --git a/README.md b/README.md index 0bebd06d..b4dcaadd 100644 --- a/README.md +++ b/README.md @@ -34,8 +34,14 @@ Pre-alpha. Implemented and tested: are comparable, and marking controls, in `experiment_design.csv` (`design`) - Running an experiment end to end and updating earlier output incrementally, with a run manifest (`run`, `manifest`) +- Masks (optional): pairing every image with its masks of any number of kinds (a masks + folder with one folder per kind, and OME-Zarr label images), checking them from their + headers and recording them in the run manifest (`masks`); and a mask step that reads them + and tells how the kinds relate (inside, grouped by, same objects, with the counts behind + it, such as cells with two nuclei) and their hierarchy, in `mask_relationships.json` + (`mask_scan`, command `masks`) - A step-by-step notebook ([notebooks/MorphoBricks_Preprocessing.ipynb](notebooks/MorphoBricks_Preprocessing.ipynb)) - and a command line (`morphobricks-preprocessing check` / `run`) + and a command line (`morphobricks-preprocessing check` / `run` / `masks`) From the normalization brick (`morphobricks_normalization`, a required dependency, pinned to its release v0.2.0): per-plane histograms diff --git a/docs/contract.md b/docs/contract.md index 0a469a1d..0d4286aa 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -107,8 +107,8 @@ method is an offset and a scale, `(pixels − offset) / scale`. | Severity | Situations | Effect | |---|---|---| -| Error | Experiment design not confirmed (`design_not_confirmed`) or a repeat group holding two folders of one condition (`repeat_group_conflict`), when a design is used (§12); mixed dtypes; images with different numbers of channels (`different_channel_counts`, §2); unreadable TIFF (header or pixel data) or OME-Zarr folder in a leaf folder that is not excluded; OME-Zarr images, or converting TIFF images to OME-Zarr, without the `zarr` extra; OME-Zarr images whose declared axes do not fit the experiment's roles; every image excluded; more than one image series in a file; unsupported sample format or axis; missing/duplicate X or Y, duplicate roles, unassigned non-singleton axis, inconsistent mapping; plane that cannot be read within the memory limit, not even in blocks (for example a plane stored compressed in one piece); folder named `single` in another case, or beside other folders (§1); sibling folders differing only in case; images in one leaf folder with the same position identifier (§1); folder below the technical-repeat level; output folder inside the experiment; no usable planes | Processing is blocked | -| Warning | Files that are not images; empty folders; images outside leaf folders; images excluded by the user, or excluded paths that match no image; channel names differing between files; non-finite pixels; different image widths or heights; different numbers of depth slices, time points or other planes; physical pixel sizes that differ or are known for only some images; unbalanced numbers of images; biological repeats missing from some conditions (§11); planes without finite pixels while others remain; planes over the memory limit that are read in blocks | Processing continues; the item is reported | +| Error | Experiment design not confirmed (`design_not_confirmed`) or a repeat group holding two folders of one condition (`repeat_group_conflict`), when a design is used (§12); mixed dtypes; images with different numbers of channels (`different_channel_counts`, §2); unreadable TIFF (header or pixel data) or OME-Zarr folder in a leaf folder that is not excluded; OME-Zarr images, or converting TIFF images to OME-Zarr, without the `zarr` extra; OME-Zarr images whose declared axes do not fit the experiment's roles; every image excluded; more than one image series in a file; unsupported sample format or axis; missing/duplicate X or Y, duplicate roles, unassigned non-singleton axis, inconsistent mapping; plane that cannot be read within the memory limit, not even in blocks (for example a plane stored compressed in one piece); folder named `single` in another case, or beside other folders (§1); sibling folders differing only in case; images in one leaf folder with the same position identifier (§1); folder below the technical-repeat level; output folder inside the experiment; no usable planes; when masks are given, a masks folder missing or inside the experiment, mask names differing only in case, two masks of one type for an image, or a mask that is unreadable, not 2D, of another size than its image or not of integers or booleans (§13) | Processing is blocked | +| Warning | Files that are not images; empty folders; images outside leaf folders; images excluded by the user, or excluded paths that match no image; channel names differing between files; non-finite pixels; different image widths or heights; different numbers of depth slices, time points or other planes; physical pixel sizes that differ or are known for only some images; unbalanced numbers of images; biological repeats missing from some conditions (§11); planes without finite pixels while others remain; planes over the memory limit that are read in blocks; images without a mask of a type other images have, masks matching no analysed image, files outside the mask-type folders (§13) | Processing continues; the item is reported | Every message states what happened in plain language, lists the affected paths, and suggests an action. Technical details go to the diagnostic log. @@ -150,6 +150,7 @@ file's schema or statistical interpretation. | `files` | One entry per analysed file, sorted by path: `relative_path`, `size` (bytes), `mtime_ns`, the three labels, `axes`, `singleton_roles`, `plane_count` (its rows in the table: XYC planes), `channel_names` (list, null for unnamed channels; empty if none in metadata), `significant_bits` (null if not in metadata), `ome_zarr_path` (for a TIFF converted to OME-Zarr, the copy relative to the output folder; null otherwise or when the key is absent), `read_in_blocks` (true when the image's planes are over the memory limit and are read in blocks; information only, false when absent), `position_id` (the image's position identifier, §1; null when absent, in manifests written before it). For an OME-Zarr image, `size` is the total size of its files and `mtime_ns` the time of its newest file. | | `csv_sha256` | SHA-256 of the table written with this manifest. | | `experiment_design_sha256` | SHA-256 of the `experiment_design.csv` written with this run (§12), so later bricks can check that the design they read is the one used; null when the run had no design, or absent in manifests written before it. Optional key: the manifest stays version 1. | +| `masks` | The masks paired with the images (§13), when the run was given a masks folder: `folder` (absolute path of the masks folder, null when only OME-Zarr label images were used), `types` (sorted mask-type names) and `files` (one entry per image and mask type, sorted by image then type: `image`, `type`, `path`, `source` `folder` or `ome-zarr`). Absent when the run had no masks. Optional key: the manifest stays version 1, and the table does not depend on the masks. | | `versions` | Versions of `morphobricks-preprocessing`, `morphobricks-normalization` (which computes the histograms and writes the CSV; absent in manifests written before it), `numpy`, `tifffile` and Python. Information only: a different version never causes a rebuild. | Channel names and significant bits come from OME metadata and are for information only. @@ -239,3 +240,85 @@ folder. | Notebook | Step 1 asks once the labels are confirmed and writes the file after every answer (unconfirmed groups with `confirmed=false`). A design confirmed for the same folders is loaded with the questions folded away behind **Change the experiment design**. Step 3 runs only once every group is confirmed. | | Command line | `run` reads the file from the output folder, or `--design FILE`. Without a fully confirmed design it writes the design brought up to date with the folders (a template with `confirmed=false` for the groups still to confirm, keeping the answers already given) to the output folder, explains the columns, and stops with status 1. Confirmed rows without a date get the date of the run. | | Run | `run_experiment(..., design=...)` brings the design up to date with the folders and stops with `design_not_confirmed` or `repeat_group_conflict` if it is not fully confirmed. Otherwise the file is written with the table and published before the manifest, which records its SHA-256 (§9). Without a design the run is unchanged: no file is written, an existing one is left as it is, and `experiment_design_sha256` is null. `plane_histograms.csv` does not depend on the design. | + +## 13. Masks (optional) + +Masks are segmentations made elsewhere: label images of cells, nuclei, nucleoli, vesicles or +any other object, one per image and mask type. This brick finds them, checks them and +describes how they relate, for the bricks that measure objects. It never changes a mask, and +`plane_histograms.csv` does not depend on the masks. + +``` +masks/ masks folder, outside the experiment folder + Cells/ one folder per mask type, named by the user + control/bio_rep1/tech_rep1/img_01.tif + Nuclei/ + control/bio_rep1/tech_rep1/img_01.tif +``` + +| Rule | Detail | +|---|---| +| Layout | One subfolder per mask type, each mirroring the experiment's folders. Any number of mask types. The label images of an OME-Zarr image (its OME-NGFF `labels` group, with the `zarr` extra) are masks too, of the type named by the label image; they are used whenever masks are asked for, with or without a masks folder. | +| Pairing | A mask belongs to the image in the same leaf folder with the same `position_id` (§1), ignoring upper/lower case; its file is a TIFF (`.tif`, `.tiff`, optionally `.ome`, any case). Excluded images have no masks. | +| Names | A mask type is named by its folder or label image. Two names differing only in case are the error `mask_names_differ_in_case`; one image with two masks of one type (two extensions, or a file and a label image) is the error `duplicate_mask`. | +| Header checks | From the headers only, no pixels read (`check`, `run`, `preflight`): a mask is 2D once its length-1 axes are removed (`mask_not_2d`), has the Y and X sizes of its image (`mask_size_mismatch`), and integer or boolean pixels (`mask_pixel_type`); a mask that cannot be opened, or holds several series, is `unreadable_mask`. All four are errors. | +| Folder | The masks folder must exist (`masks_folder_not_found`) and be outside the experiment folder (`masks_inside_experiment`), where its folders would be read as labels. Both errors. | +| Warnings | `mask_missing` (images without a mask of a type that other images have: their objects of that type cannot be measured), `mask_without_image` (files matching no analysed image), `mask_outside_type_folder` (files directly in the masks folder). | +| Run | `run_experiment(..., masks_folder=...)` stops on a mask error and records the masks in the manifest's `masks` key (§9). Without masks, the run and its outputs are unchanged. | + +### 13.1 The mask step (`mask_relationships.json`, report version 1) + +The optional mask step reads every image's masks, one image at a time, and describes each +mask type and every ordered pair of types over the whole experiment. Pieces are counted with +4-connectivity; label 0 is background. Every conclusion is a suggestion with its evidence: +the user names the masks and confirms or corrects the relationships. + +| Kind of mask | Rule | +|---|---| +| `empty` | No foreground in any image. | +| `one_value` | Every image's mask has at most one foreground value: a binary mask, or a single object (only the user can tell which). | +| `grouped` | The median number of pieces per label, over every label of every image, is more than 2: each label is spread over many pieces, as in a mask of mitochondria whose pixels carry their cell's label. | +| `instances` | Otherwise: each label is one object. | + +For a child type and a parent type, the counts are summed over the images with both masks. +Each child object's parent is the parent object sharing most of its pixels. + +| Relation of child to parent | Rule (fixed thresholds, written in the report) | +|---|---| +| `same_objects` | At least 90% of the children and 90% of the parents are matched one to one at an intersection over union of at least 0.5: two segmentations of the same objects. | +| `grouped_by` | The child is `grouped`, at least 90% of its objects are inside as below, and at least 99% of its pixels carry the label of the parent object under them. | +| `inside` | At least 90% of the child objects have at least half of their pixels in one parent object, and that parent is larger than the child for at least 90% of them. | +| `none` | Otherwise. | + +| Evidence (per pair) | Meaning | +|---|---| +| `children`, `parents` | Objects of each type. | +| `inside`, `inside_parent_larger` | Children with at least half of their pixels in one parent; of those, children whose parent is larger. | +| `split_between_parents`, `outside_parents` | Children touching parents but with less than half of their pixels in any one; children touching no parent. | +| `parents_with_0_children`, `parents_with_1_child`, `parents_with_2_children`, `parents_with_3_or_more_children` | How many children each parent holds (two nuclei in one cell counts as 2). | +| `child_pixels`, `child_pixels_with_parent_label` | Child pixels; those whose label equals the parent label under them. | +| `children_matched_one_to_one`, `parents_matched_one_to_one` | Objects with a partner at an intersection over union of at least 0.5. | + +**Hierarchy.** A relation `inside` or `grouped_by` is **direct** unless another type lies +between: nucleoli inside nuclei inside cells have the direct parent nuclei, and the +relation nucleoli inside cells is kept in the report but not in the hierarchy. A type may +have several direct parents (cells inside an organoid and inside a region of the slide), and +two types that are the same objects are both kept. + +The report is UTF-8 JSON (two-space indent, LF), written atomically to the output folder: + +| Key | Meaning | +|---|---| +| `report_version` | `1`. | +| `thresholds` | The thresholds above. | +| `images` | Images scanned (every image with at least one mask). | +| `masks` | Per mask type: `kind`, `images`, `objects` (labels over all images) and `median_pieces_per_label`. | +| `hierarchy` | Per mask type, its direct parents (sorted; empty for none). | +| `same_objects` | Pairs of types that are the same objects, each pair sorted. | +| `relationships` | Every ordered pair with a relation other than `none`: `child`, `parent`, `relation`, `direct` and `evidence`. | +| `pairs` | Every ordered pair of types found together in at least one image, `none` included, in the same form. | +| `versions` | As in the manifest (§9); information only. | + +The step is the command `masks` and the functions `find_masks` and `scan_masks`; it reads +each image's masks once and holds only one image's masks in memory. A mask that is not 2D +when its pixels are read raises an error. diff --git a/docs/plan.md b/docs/plan.md index f9e283ef..52d6aef9 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -691,3 +691,29 @@ Raised during review; to be decided before the milestones that need them: depends on the normalization brick's release v0.2.0 (tag, not commit) instead of the Milestone 2b merge commit `0e53e67`; v0.2.0 adds only RGB and intensity checks, which this brick does not call yet. The next tag, v0.2.0, follows follow-ups P4 and P5.* +29. Can this brick find each image's masks and tell how they relate, so that a user who is + not a programmer only has to name the masks (MorphoFeature decision D78)? Bricks that + measure objects (MorphoObjects, MorphoRelate, CellPatcher, MorphoFeature's workflows) + need every image paired with its masks and the parent-child relationships between mask + types (two nuclei in one cell, nucleoli inside nuclei inside cells, mitochondria grouped + by cell), and today each asks the user again. + *Decided (October 2026), by the MorphoBricks owner: yes, for any number of mask types + (contract §13). Masks are in a separate masks folder, outside the experiment, one + subfolder per mask type mirroring the experiment's folders, paired with the images by + leaf folder and `position_id`; the label images of OME-Zarr images are masks too. The + checks of `check` and `run` read headers only; reading the pixels is an optional mask + step (`masks` command, `scan_masks`), which writes `mask_relationships.json` to the + output folder: each mask type's kind (`instances`, `grouped`, `one_value`, `empty`), + each ordered pair's relation (`same_objects`, `grouped_by`, `inside`, `none`) with the + counts behind it (including parents with 0, 1, 2 and 3 or more children), and the + hierarchy of direct parents. A run given masks records them in the manifest's optional + `masks` key (manifest version 1); the table never depends on them. The report suggests; + it decides nothing: users name the masks (by default the folder name, names unique) + and confirm or correct the relationships in the bricks that use them. Library, command + line and documentation first; the notebook step follows in a later change. SciPy + becomes a core dependency (the owner's agreement: "can use scipy if good reasons"): + counting the connected pieces of every label of one 2048 x 2048 HPMEC cell mask took + about 5 s with a NumPy union-find, and 0.3 to 0.9 s per HPMEC mask with `scipy.ndimage`. Proposed defaults, to confirm in + review: masks are 2D with the Y and X sizes of their image (3D masks later), their + pixels integers or booleans, and the thresholds fixed (written in every report), not + settings.* diff --git a/docs/user_guide.md b/docs/user_guide.md index 360d5f91..48e41aea 100644 --- a/docs/user_guide.md +++ b/docs/user_guide.md @@ -92,6 +92,36 @@ experiment/ To change a label, rename the folder and check again. +### Masks (optional) + +If you already have masks (segmentations of cells, nuclei, nucleoli, ...), keep them in a +**masks folder** next to the experiment folder, not inside it. Make one folder per kind of +mask, named as you like (`Cells`, `Nuclei`, ...), and inside each the same folders as the +experiment, with one mask per image under the image's name: + +``` +masks/ + Cells/ + control/bio_rep1/tech_rep1/img_01.tif + control/bio_rep1/tech_rep1/img_02.tif + Nuclei/ + control/bio_rep1/tech_rep1/img_01.tif + control/bio_rep1/tech_rep1/img_02.tif +``` + +- A mask is a 2D TIFF label image with the width and height of its image: background 0, + each object its own whole number (or a true/false mask). The extension may differ from + the image's (`img_01.tiff` for `img_01.ome.tif`). +- You may have any number of kinds of mask. Their folder names are their names: they must + differ by more than upper/lower case. +- The label images stored inside an OME-Zarr image (its `labels` group) are used as masks + too, named after the label image. +- An image without a mask of one kind is reported; its objects of that kind cannot be + measured. + +Masks do not change `plane_histograms.csv`. They are checked and recorded for the bricks +that measure objects; the optional `masks` command (§5) also tells how they relate. + ## 3. Supported images - **TIFF files:** `.tif` and `.tiff` (any case), including OME-TIFF and ImageJ hyperstacks. @@ -295,8 +325,38 @@ the design is confirmed. Other options: | `--compression lz4\|zstd\|none`, `--compression-level N` | With `--convert-tiff`: lossless compression (default lz4, level 5) | | `--pyramid-levels auto\|none\|N` | With `--convert-tiff`: smaller levels for viewing (default auto) | | `--log FILE` | Write technical details to a file | +| `--masks FOLDER` | The masks folder (§2). `check` lists the mask types found and checks every mask from its header; `run` stops on a mask error and records the masks in `run_manifest.json` | Both commands exit with status 1 when an error blocks processing. + +The optional `masks` command reads the masks and tells how they relate: + +```bash +morphobricks-preprocessing masks /data/exp1 --masks /data/exp1_masks +``` + +It reads every image's masks, one image at a time, and prints for each kind of mask +whether its labels are separate objects (`instances`), labels spread over many pieces +(`grouped`, such as mitochondria labelled with their cell's number), a single value per +image (a true/false mask or one object) or empty; then which kinds lie inside which, with +the counts behind it, for example how many cells hold 0, 1, 2 or 3 or more nuclei; and the +hierarchy of direct parents, for example: + +``` +Hierarchy (direct parents only): + Organoid + Cells + Nuclei + Nucleoli + Vesicles +``` + +Nucleoli inside nuclei inside cells have the nuclei as their direct parent; two kinds of +mask that segment the same objects are reported as such. These are suggestions from the +pixels, not decisions: you name each kind of mask and check the relationships in the bricks +that use them. The report is written to `mask_relationships.json` in the output folder +(§6). The command takes `--output`, `--exclude` and `--log` like the others, and exits with +status 1 when an error blocks it or when no masks are found. `python -m morphobricks_preprocessing` works the same way. ## 6. Output files @@ -338,6 +398,12 @@ The default output folder is next to the experiment folder: TIFF images directly removes the copies and recomputes their rows. Details in [contract.md](contract.md) §10. +- **`mask_relationships.json`** (only from the `masks` command): each kind of mask, the + relationship between every two kinds with its counts, and the hierarchy. The rules and + thresholds are in [contract.md](contract.md) §13. + +When masks are given, `run_manifest.json` also lists every image's masks (`masks`). + The table and manifest are both prepared before the earlier files are replaced. If a replacement fails, the earlier files are restored. This needs temporary disk space for the new outputs and copies of the earlier pair. If restoration also fails (for example, @@ -410,6 +476,16 @@ continues. | `non_finite_pixels` | Warning | NaN or infinite pixels are left out of the histograms and counted. | | `no_finite_pixels` | Warning | Some planes have a channel without finite pixels; those channels are left out of group statistics. | | `excluded_images`, `excluded_path_not_found` | Warning | Images excluded as requested, or an excluded path that matches no image. | +| `masks_folder_not_found` | Error | The masks folder does not exist. Select the folder that holds one folder per kind of mask. | +| `masks_inside_experiment` | Error | The masks folder is inside the experiment folder, where its folders would be read as conditions and repeats. Move it next to the experiment folder. | +| `mask_names_differ_in_case` | Error | Two kinds of mask have names that differ only in upper/lower case. Rename one. | +| `duplicate_mask` | Error | An image has two masks of one kind (two extensions, or a file and an OME-Zarr label image of the same name). Keep one. | +| `unreadable_mask` | Error | A mask cannot be opened, or holds several images. Save it again as a TIFF label image. | +| `mask_not_2d` | Error | A mask has more than one plane. Save one 2D label image per position. | +| `mask_size_mismatch` | Error | A mask does not have the width and height of its image. Check that it was made from that image at full resolution. | +| `mask_pixel_type` | Error | A mask's pixels are not whole numbers (for example decimals). Save it as an integer label image. | +| `mask_missing` | Warning | Some images have no mask of a kind that other images have; their objects of that kind cannot be measured. Add the masks, or exclude the images. | +| `mask_without_image`, `mask_outside_type_folder` | Warning | Files in the masks folder that match no analysed image (check the folder and file names), or that are not in a folder of their kind. They are not used. | | `zarr_not_installed` | Error | OME-Zarr images, or converting TIFF images, need the `zarr` extra. Install it, exclude the OME-Zarr images, or read the TIFF images directly. | ## 9. Python interface @@ -443,6 +519,23 @@ pairing = design.repeat_groups() # {(condition, biological_repeat): repeat_grou Later bricks read the file with `read_experiment_design(path)`. +Masks are found from their headers by `preflight` and read by `scan_masks`: + +```python +from morphobricks_preprocessing import preflight, scan_masks, write_mask_report + +check = preflight("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/data/exp1", masks_folder="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/data/exp1_masks") +masks = check.masks # MaskInventory: types, files, issues +cells = masks.for_image("control/bio_rep1/tech_rep1/img_01.tif")["Cells"] # MaskFile +report = scan_masks(masks, "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/data/exp1") # reads the pixels, one image at a time +report.hierarchy # {"Nuclei": ("Cells",), "Cells": (), ...}: direct parents +write_mask_report(report, "/data/exp1_MorphoBricks_Preprocessing/mask_relationships.json") +``` + +`preflight(..., ome_zarr_labels=True)` uses the label images of OME-Zarr images without a +masks folder; `read_mask` reads one mask's pixels; `pair_counts` and `piece_counts` give the +counts behind the report for one pair of masks or one mask. + Reading the table, group statistics and normalisation parameters come from the normalization brick (`morphobricks_normalization`, installed with this package): `read_plane_histograms`, `PooledHistogram`, `group_records`, `min_max`, diff --git a/pyproject.toml b/pyproject.toml index 29e167e4..7a85b695 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,8 @@ classifiers = [ dependencies = [ "numpy>=1.24", "tifffile>=2023.7.10", + # Counting the connected pieces of each mask label in the optional mask step. + "scipy>=1.10", # Per-plane histograms, the CSV table, pooling and normalization (moved there in P1). # Pinned to the brick's release v0.2.0, which writes plane_histograms.csv schema # version 3 with the channels in columns (follow-up P3b; P3b first pinned its diff --git a/src/morphobricks_preprocessing/__init__.py b/src/morphobricks_preprocessing/__init__.py index 5fe27b65..071e7514 100644 --- a/src/morphobricks_preprocessing/__init__.py +++ b/src/morphobricks_preprocessing/__init__.py @@ -66,6 +66,17 @@ read_manifest, write_manifest, ) +from .mask_scan import ( + MASK_REPORT_FILENAME, + MaskReport, + MaskTypeSummary, + PairRelation, + pair_counts, + piece_counts, + scan_masks, + write_mask_report, +) +from .masks import MaskFile, MaskInventory, find_masks, read_mask from .metadata import ( ImageMetadata, MetadataReport, @@ -141,6 +152,7 @@ def __dir__() -> list[str]: "DESIGN_FILENAME", "DESIGN_SCHEMA_VERSION", "MANIFEST_FILENAME", + "MASK_REPORT_FILENAME", "OUTPUT_FILENAME", "SINGLE", "AxisMapping", @@ -163,12 +175,17 @@ def __dir__() -> list[str]: "ImageMetadata", "LeafFolder", "ManifestFormatError", + "MaskFile", + "MaskInventory", + "MaskReport", + "MaskTypeSummary", "MetadataReport", "NgffError", "NgffImage", "NoValidPixelsError", "NormalizationParameters", "OmeZarrWriter", + "PairRelation", "PhysicalSizes", "Plane", "PlaneHistogram", @@ -208,6 +225,7 @@ def __dir__() -> list[str]: "discover_experiment", "estimate_conversion", "file_mapping_for", + "find_masks", "group_records", "group_summaries", "has_errors", @@ -226,7 +244,9 @@ def __dir__() -> list[str]: "ome_zarr_relative_path", "open_ngff", "open_ngff_labels", + "pair_counts", "percentile_range", + "piece_counts", "plane_layout", "pooled_statistics", "position_id", @@ -235,6 +255,7 @@ def __dir__() -> list[str]: "preflight", "read_experiment_design", "read_manifest", + "read_mask", "read_physical_sizes", "read_plane", "read_plane_histograms", @@ -245,8 +266,10 @@ def __dir__() -> list[str]: "robust_z_score", "run_experiment", "sample_planes", + "scan_masks", "write_experiment_design", "write_manifest", + "write_mask_report", "write_plane_histograms", "z_score", "zarr_available", diff --git a/src/morphobricks_preprocessing/cli.py b/src/morphobricks_preprocessing/cli.py index 4f940f31..c810d94a 100644 --- a/src/morphobricks_preprocessing/cli.py +++ b/src/morphobricks_preprocessing/cli.py @@ -1,4 +1,4 @@ -"""Command-line interface: ``morphobricks-preprocessing check`` and ``run``. +"""Command-line interface: ``morphobricks-preprocessing check``, ``run`` and ``masks``. ``check`` runs the preflight (folders and image headers, no pixels) and prints what will be analysed, the problems found and the suggested axis roles. @@ -6,7 +6,10 @@ are never taken from metadata without being stated. It also needs the experiment design confirmed in ``experiment_design.csv``: without it, ``run`` writes a template to the output folder for the user to fill in, and stops. -Both exit with status 1 when an error blocks processing. +``masks`` is the optional mask step: it pairs every image with its masks (``--masks``, and the +label images of OME-Zarr images), reads them, and writes ``mask_relationships.json`` with each +mask type and the relationships between them (contract §13). All three exit with status 1 when +an error blocks processing. """ from __future__ import annotations @@ -31,6 +34,7 @@ write_experiment_design, ) from .discovery import check_output_folder, default_output_folder +from .mask_scan import MASK_REPORT_FILENAME, scan_masks, write_mask_report from .planes import DEFAULT_MEMORY_LIMIT, format_bytes from .report import ( balance_summary, @@ -38,6 +42,8 @@ format_issues, label_tree, mapping_summary, + mask_report_summary, + mask_summary, preflight_summary, run_summary, ) @@ -188,6 +194,13 @@ def build_parser() -> argparse.ArgumentParser: f"none, or a number from 1 to {MAX_LEVELS - 1}", ) run.add_argument("--log", metavar="FILE", help="write technical details to this file") + + masks = commands.add_parser( + "masks", + help="find every image's masks, describe them and how they relate (reads the masks)", + ) + _common(masks) + masks.add_argument("--log", metavar="FILE", help="write technical details to this file") return parser @@ -197,7 +210,7 @@ def main(argv: Sequence[str] | None = None) -> int: args = build_parser().parse_args(argv) handler = _attach_log(args.log) if getattr(args, "log", None) else None try: - return _check(args) if args.command == "check" else _run(args) + return {"check": _check, "run": _run, "masks": _masks}[args.command](args) finally: if handler is not None: logging.getLogger("morphobricks_preprocessing").removeHandler(handler) @@ -218,14 +231,23 @@ def _common(parser: argparse.ArgumentParser) -> None: metavar="PATH", help="image to leave out, relative to the experiment folder (repeatable)", ) + parser.add_argument( + "--masks", + metavar="FOLDER", + help="masks folder: one subfolder per mask type, each mirroring the experiment's " + "folders, outside the experiment; the label images of OME-Zarr images are added", + ) def _check(args: argparse.Namespace) -> int: - pre = preflight(args.experiment, args.output, exclude=args.exclude) + pre = preflight(args.experiment, args.output, exclude=args.exclude, masks_folder=args.masks) print(preflight_summary(pre)) if pre.discovery.leaves: print() print(label_tree(pre.discovery)) + if pre.masks is not None: + print() + print(mask_summary(pre.masks)) if pre.dataset.balance.cells: print() print(balance_summary(pre.dataset.balance)) @@ -255,8 +277,42 @@ def _check(args: argparse.Namespace) -> int: return 0 if pre.ok else 1 +def _masks(args: argparse.Namespace) -> int: + pre = preflight( + args.experiment, + args.output, + exclude=args.exclude, + masks_folder=args.masks, + ome_zarr_labels=True, + ) + print(mask_summary(pre.masks)) + blocking = [*pre.discovery.issues, *pre.exclusion_issues, *pre.masks.issues] + if blocking: + print() + print(format_issues(blocking)) + if has_errors(blocking): + return 1 + if not pre.masks.types: + print("Give the masks folder with --masks, or use OME-Zarr images with label images.") + return 1 + + def progress(done: int, total: int, path: str) -> None: + if path: + print(f"[{done + 1}/{total}] {path}", flush=True) + + print() + report = scan_masks(pre.masks, pre.discovery.root, progress=progress) + output = Path(args.output) if args.output else default_output_folder(pre.discovery.root) + path = write_mask_report(report, output / MASK_REPORT_FILENAME) + print() + print(mask_report_summary(report)) + print() + print(f"Report written to {path}") + return 0 + + def _run(args: argparse.Namespace) -> int: - pre = preflight(args.experiment, args.output, exclude=args.exclude) + pre = preflight(args.experiment, args.output, exclude=args.exclude, masks_folder=args.masks) axes = reference_axes(pre.metadata.files) if axes is None: print(format_issues(pre.issues) or "No readable images were found.") @@ -316,6 +372,7 @@ def progress(done: int, total: int, path: str) -> None: convert_tiffs=args.convert_tiff, conversion=conversion, design=design, + masks_folder=args.masks, ) if result.issues: print(format_issues(result.issues)) diff --git a/src/morphobricks_preprocessing/manifest.py b/src/morphobricks_preprocessing/manifest.py index 35141d23..97fb6e18 100644 --- a/src/morphobricks_preprocessing/manifest.py +++ b/src/morphobricks_preprocessing/manifest.py @@ -165,6 +165,10 @@ class RunManifest: SHA-256 of the ``experiment_design.csv`` written with this run, or None when the run had no experiment design (or in manifests written before it was recorded). + masks + The masks of the images when masks were given (contract §13): the masks + folder, the mask types and one entry per mask (image, type, path, + source); None otherwise. The key is left out of the file when None. schema_version Version of the CSV format. manifest_version @@ -180,6 +184,7 @@ class RunManifest: versions: dict[str, str] = field(default_factory=dict) written_at: str = "" experiment_design_sha256: str | None = None + masks: dict | None = None schema_version: str = SCHEMA_VERSION manifest_version: str = MANIFEST_VERSION @@ -197,7 +202,7 @@ def file(self, relative_path: str) -> FileEntry | None: def to_dict(self) -> dict: """Return the manifest as JSON-compatible data.""" - return { + data = { "manifest_version": self.manifest_version, "schema_version": self.schema_version, "written_at": self.written_at, @@ -226,6 +231,9 @@ def to_dict(self) -> dict: "experiment_design_sha256": self.experiment_design_sha256, "versions": dict(self.versions), } + if self.masks is not None: + data["masks"] = self.masks + return data @classmethod def from_dict(cls, data: dict) -> RunManifest: @@ -265,6 +273,7 @@ def from_dict(cls, data: dict) -> RunManifest: if data.get("experiment_design_sha256") is None else str(data["experiment_design_sha256"]) ), + masks=_masks(data.get("masks")), schema_version=str(data["schema_version"]), ) except ManifestFormatError: @@ -380,6 +389,15 @@ def read_manifest(path: str | os.PathLike) -> RunManifest: return RunManifest.from_dict(data) +def _masks(data: object) -> dict | None: + """Check the optional ``masks`` entry of a manifest (absent in older manifests).""" + if data is None: + return None + if not isinstance(data, dict) or not isinstance(data.get("files"), list): + raise ValueError("its masks entry is not a mask list") + return data + + def _file_dict(entry: FileEntry) -> dict: return { "relative_path": entry.relative_path, diff --git a/src/morphobricks_preprocessing/mask_scan.py b/src/morphobricks_preprocessing/mask_scan.py new file mode 100644 index 00000000..9196c19d --- /dev/null +++ b/src/morphobricks_preprocessing/mask_scan.py @@ -0,0 +1,409 @@ +"""Describe every mask type and how the masks relate to each other (``docs/contract.md`` §13). + +The optional mask step: it reads the masks found by :func:`.masks.find_masks`, one image at a +time, and summarises each mask type and every ordered pair of types over the whole +experiment, with the evidence (counts) behind each conclusion: + +- **Kind of mask:** ``instances`` (each label is one object, in at most two pieces: a median + over all objects), ``grouped`` (each label is spread over more than two pieces, as in a + per-cell mask of mitochondria whose pixels carry their cell's label), ``one_value`` (every + image's mask has a single foreground value: a binary mask or one object, which only the user + can tell apart) or ``empty``. +- **Relationship of a child type to a parent type:** + + - ``same_objects``: one-to-one matches at intersection over union >= 0.5 for at least 90% of + the objects of both types (two segmentations of the same objects); + - ``grouped_by``: the child is ``grouped``, at least 99% of its pixels carry the label of the + parent object under them, and it is inside as below; + - ``inside``: at least 90% of the child objects have at least half of their pixels in one + parent object, which is larger than the child for at least 90% of them; + - ``none`` otherwise. + + Several children per parent are counted (parents with 0, 1, 2 and 3 or more children), as + are children outside every parent and children split between parents. +- **Hierarchy:** the direct parents of each type: a parent is left out when it is the parent + of another parent (nucleoli inside nuclei inside cells: the nucleoli's direct parent is the + nuclei). A type may have several direct parents (cells inside an organoid and inside a + region of the slide). + +Nothing is decided for the user: the report is a suggestion with its evidence; the user names +the masks and confirms or corrects the relationships. Pieces are counted with 4-connectivity +(``scipy.ndimage.label``). +""" + +from __future__ import annotations + +import json +import logging +import os +import tempfile +from collections.abc import Callable, Mapping +from dataclasses import dataclass, field +from pathlib import Path +from types import MappingProxyType + +import numpy as np +from scipy import ndimage + +from .manifest import package_versions +from .masks import MaskInventory, read_mask + +logger = logging.getLogger(__name__) + +MASK_REPORT_FILENAME = "mask_relationships.json" +"""Name of the report written to the output folder.""" + +MASK_REPORT_VERSION = "1" +"""Version of the report's format.""" + +THRESHOLDS = MappingProxyType( + { + "inside_pixel_share": 0.5, + "inside_object_share": 0.9, + "parent_larger_share": 0.9, + "same_objects_iou": 0.5, + "same_objects_share": 0.9, + "grouped_label_share": 0.99, + "grouped_median_pieces": 2, + } +) +"""Fixed thresholds of the relationships (contract §13); written in every report.""" + +_COUNTS = ( + "children", + "parents", + "inside", + "inside_parent_larger", + "split_between_parents", + "outside_parents", + "parents_with_0_children", + "parents_with_1_child", + "parents_with_2_children", + "parents_with_3_or_more_children", + "child_pixels", + "child_pixels_with_parent_label", + "children_matched_one_to_one", + "parents_matched_one_to_one", +) + + +@dataclass(frozen=True) +class MaskTypeSummary: + """One mask type over the experiment. + + Attributes + ---------- + name + The mask type. + kind + ``instances``, ``grouped``, ``one_value`` or ``empty`` (module docstring). + images + Images with a mask of this type. + objects + Distinct labels, summed over the images. + median_pieces_per_label + Median number of 4-connected pieces per label over all objects; 0 when empty. + """ + + name: str + kind: str + images: int + objects: int + median_pieces_per_label: float + + +@dataclass(frozen=True) +class PairRelation: + """How one mask type (the child) sits relative to another (the parent). + + Attributes + ---------- + child, parent + The two mask types. + relation + ``inside``, ``grouped_by``, ``same_objects`` or ``none`` (module docstring). + direct + True for a relationship kept in the hierarchy (not implied by another parent). + evidence + The counts behind the relation, summed over the images with both masks. + """ + + child: str + parent: str + relation: str + direct: bool + evidence: Mapping[str, int] + + +@dataclass(frozen=True) +class MaskReport: + """The kinds of the mask types and their relationships, with the evidence. + + Attributes + ---------- + images + Images scanned (those with at least one mask). + types + One summary per mask type, by name. + pairs + Every ordered pair of types seen together in at least one image, sorted. + """ + + images: int + types: Mapping[str, MaskTypeSummary] = field(default_factory=lambda: MappingProxyType({})) + pairs: tuple[PairRelation, ...] = () + + @property + def relationships(self) -> tuple[PairRelation, ...]: + """The pairs with a relation other than ``none``.""" + return tuple(pair for pair in self.pairs if pair.relation != "none") + + @property + def hierarchy(self) -> Mapping[str, tuple[str, ...]]: + """Each mask type's direct parents (empty for a type without a parent).""" + parents = {name: [] for name in self.types} + for pair in self.pairs: + if pair.direct: + parents[pair.child].append(pair.parent) + return MappingProxyType({name: tuple(sorted(found)) for name, found in parents.items()}) + + @property + def same_objects(self) -> tuple[tuple[str, str], ...]: + """Pairs of types that segment the same objects, each pair once.""" + found = { + tuple(sorted((pair.child, pair.parent))) + for pair in self.pairs + if pair.relation == "same_objects" + } + return tuple(sorted(found)) + + def to_dict(self) -> dict: + """Return the report as JSON-compatible data (contract §13).""" + return { + "report_version": MASK_REPORT_VERSION, + "thresholds": dict(THRESHOLDS), + "images": self.images, + "masks": { + name: { + "kind": summary.kind, + "images": summary.images, + "objects": summary.objects, + "median_pieces_per_label": summary.median_pieces_per_label, + } + for name, summary in sorted(self.types.items()) + }, + "hierarchy": {name: list(parents) for name, parents in self.hierarchy.items()}, + "same_objects": [list(pair) for pair in self.same_objects], + "relationships": [_pair_dict(pair) for pair in self.relationships], + "pairs": [_pair_dict(pair) for pair in self.pairs], + "versions": package_versions(), + } + + +def piece_counts(mask: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Return the labels of a 2D mask and the number of 4-connected pieces of each. + + Two pixels are in one piece when they carry the same label and touch along an edge. + + Returns + ------- + tuple of numpy.ndarray + The sorted nonzero labels, and the number of pieces of each. + """ + values, dense = np.unique(mask, return_inverse=True) + dense = dense.reshape(mask.shape) + if values.size and values[0] == 0: + values = values[1:] + else: + dense = dense + 1 # no background: every value is a label + counts = np.zeros(values.size, dtype=np.int64) + for index, box in enumerate(ndimage.find_objects(dense)): + if box is not None: + counts[index] = ndimage.label(dense[box] == index + 1)[1] + return values, counts + + +def pair_counts(child: np.ndarray, parent: np.ndarray) -> dict[str, int]: + """Return the counts relating the objects of one mask (child) to another (parent). + + Both are 2D label images of one position; label 0 is background. + """ + counts = dict.fromkeys(_COUNTS, 0) + on = child != 0 + child_labels, child_of = np.unique(child[on], return_inverse=True) + parent_labels, parent_sizes = np.unique(parent[parent != 0], return_counts=True) + counts["children"], counts["parents"] = child_labels.size, parent_labels.size + if child_labels.size == 0: + counts["parents_with_0_children"] = parent_labels.size + return counts + child_sizes = np.bincount(child_of, minlength=child_labels.size) + under = parent[on] + counts["child_pixels"] = int(on.sum()) + counts["child_pixels_with_parent_label"] = int(np.count_nonzero(child[on] == under)) + covered = under != 0 + slots = parent_labels.size + pairs, shared = np.unique( + child_of[covered] * slots + np.searchsorted(parent_labels, under[covered]), + return_counts=True, + ) + kid, host = pairs // slots, pairs % slots + if kid.size == 0: # no child pixel lies on a parent + counts["outside_parents"] = child_labels.size + counts["parents_with_0_children"] = slots + return counts + union = child_sizes[kid] + parent_sizes[host] - shared + matched = shared / union >= THRESHOLDS["same_objects_iou"] + counts["children_matched_one_to_one"] = int(np.unique(kid[matched]).size) + counts["parents_matched_one_to_one"] = int(np.unique(host[matched]).size) + best = np.lexsort((-shared, kid)) + lead = np.r_[True, kid[best][1:] != kid[best][:-1]] + best_kid, best_host, best_shared = kid[best][lead], host[best][lead], shared[best][lead] + inside = best_shared / child_sizes[best_kid] >= THRESHOLDS["inside_pixel_share"] + counts["inside"] = int(inside.sum()) + counts["inside_parent_larger"] = int( + np.count_nonzero(parent_sizes[best_host[inside]] > child_sizes[best_kid[inside]]) + ) + counts["split_between_parents"] = int(np.count_nonzero(~inside)) + counts["outside_parents"] = child_labels.size - best_kid.size + per_parent = np.bincount(best_host[inside], minlength=slots) + counts["parents_with_0_children"] = int(np.count_nonzero(per_parent == 0)) + counts["parents_with_1_child"] = int(np.count_nonzero(per_parent == 1)) + counts["parents_with_2_children"] = int(np.count_nonzero(per_parent == 2)) + counts["parents_with_3_or_more_children"] = int(np.count_nonzero(per_parent >= 3)) + return counts + + +def scan_masks( + inventory: MaskInventory, + experiment_root: str | os.PathLike, + *, + progress: Callable[[int, int, str], None] | None = None, +) -> MaskReport: + """Read every image's masks and describe the mask types and their relationships. + + Parameters + ---------- + inventory + The masks found by :func:`.masks.find_masks`; it should have no errors. + experiment_root + The experiment folder (for OME-Zarr label images). + progress + Called as ``progress(done, total, image)`` before each image is read. + + Returns + ------- + MaskReport + """ + root = Path(experiment_root).resolve() + images = inventory.images + pieces: dict[str, list[np.ndarray]] = {name: [] for name in inventory.types} + single: dict[str, list[bool]] = {name: [] for name in inventory.types} + totals: dict[tuple[str, str], dict[str, int]] = {} + for done, image in enumerate(images): + if progress is not None: + progress(done, len(images), image) + masks = { + name: read_mask(mask, inventory.folder, root) + for name, mask in inventory.for_image(image).items() + } + for name, pixels in masks.items(): + labels, counts = piece_counts(pixels) + pieces[name].append(counts) + single[name].append(labels.size <= 1) + for child in masks: + for parent in masks: + if child != parent: + found = pair_counts(masks[child], masks[parent]) + total = totals.setdefault((child, parent), dict.fromkeys(_COUNTS, 0)) + for key, value in found.items(): + total[key] += value + if progress is not None: + progress(len(images), len(images), "") + types = {name: _summary(name, pieces[name], single[name]) for name in inventory.types} + return MaskReport( + images=len(images), + types=MappingProxyType(types), + pairs=_relations(types, totals), + ) + + +def write_mask_report(report: MaskReport, path: str | os.PathLike) -> Path: + """Write the report as UTF-8 JSON (two-space indent, LF), atomically; return the path.""" + target = Path(path) + target.parent.mkdir(parents=True, exist_ok=True) + text = json.dumps(report.to_dict(), indent=2, ensure_ascii=False) + "\n" + handle, temporary = tempfile.mkstemp(prefix=".mask_report-", dir=target.parent) + try: + with os.fdopen(handle, "w", encoding="utf-8", newline="\n") as stream: + stream.write(text) + os.replace(temporary, target) + finally: + if os.path.exists(temporary): + os.remove(temporary) + return target + + +def _summary(name: str, pieces: list[np.ndarray], single: list[bool]) -> MaskTypeSummary: + """Return a mask type's summary from its pieces per label in every image.""" + every = np.concatenate(pieces) if pieces else np.zeros(0, dtype=np.int64) + median = float(np.median(every)) if every.size else 0.0 + if every.size == 0: + kind = "empty" + elif all(single): + kind = "one_value" + elif median > THRESHOLDS["grouped_median_pieces"]: + kind = "grouped" + else: + kind = "instances" + return MaskTypeSummary(name, kind, len(pieces), int(every.size), median) + + +def _relation(child: MaskTypeSummary, counts: Mapping[str, int]) -> str: + """Return the relation of a child type to a parent type from their counts.""" + children, parents = counts["children"], counts["parents"] + if not children or not parents: + return "none" + share = THRESHOLDS["same_objects_share"] + if ( + counts["children_matched_one_to_one"] >= share * children + and counts["parents_matched_one_to_one"] >= share * parents + ): + return "same_objects" + inside = counts["inside"] >= THRESHOLDS["inside_object_share"] * children + if ( + child.kind == "grouped" + and inside + and counts["child_pixels_with_parent_label"] + >= THRESHOLDS["grouped_label_share"] * counts["child_pixels"] + ): + return "grouped_by" + larger = counts["inside_parent_larger"] >= THRESHOLDS["parent_larger_share"] * counts["inside"] + return "inside" if inside and larger else "none" + + +def _relations( + types: Mapping[str, MaskTypeSummary], totals: Mapping[tuple[str, str], Mapping[str, int]] +) -> tuple[PairRelation, ...]: + """Return every pair's relation, marking the direct ones (module docstring).""" + found = { + (child, parent): _relation(types[child], counts) + for (child, parent), counts in totals.items() + } + nested = {pair for pair, relation in found.items() if relation in ("inside", "grouped_by")} + pairs = [] + for (child, parent), relation in sorted(found.items()): + between = any((child, other) in nested and (other, parent) in nested for other in types) + direct = (child, parent) in nested and not between + evidence = MappingProxyType(dict(totals[(child, parent)])) + pairs.append(PairRelation(child, parent, relation, direct, evidence)) + return tuple(pairs) + + +def _pair_dict(pair: PairRelation) -> dict: + return { + "child": pair.child, + "parent": pair.parent, + "relation": pair.relation, + "direct": pair.direct, + "evidence": dict(pair.evidence), + } diff --git a/src/morphobricks_preprocessing/masks.py b/src/morphobricks_preprocessing/masks.py new file mode 100644 index 00000000..75fc50e6 --- /dev/null +++ b/src/morphobricks_preprocessing/masks.py @@ -0,0 +1,428 @@ +"""Find the masks of every image and check them from their headers (``docs/contract.md`` §13). + +Masks (segmentations: label images of cells, nuclei, ...) are kept in a masks folder outside +the experiment, one subfolder per mask type, each mirroring the experiment's folders:: + + /////.tif + +A mask belongs to the image with the same position (same leaf folder, same ``position_id``, +contract §1); its file may be a TIFF (``.tif``, ``.tiff``, any case, with an optional +``.ome``). The label images of an OME-Zarr image (its OME-NGFF ``labels`` group) are masks too, +of the type named by the label image, when the ``zarr`` extra is installed. + +Any number of mask types is accepted; each is named by its folder (or label image), and two +names may not differ only in case. Only headers are read here, never pixels: a mask must be +2D (length-1 axes aside) with the image's Y and X sizes, and its pixels integers or booleans. +Reading the pixels and relating the masks to each other is :mod:`.mask_scan`. +""" + +from __future__ import annotations + +import logging +import os +from collections.abc import Iterable, Mapping +from dataclasses import dataclass +from pathlib import Path +from types import MappingProxyType + +import numpy as np +import tifffile + +from .discovery import DiscoveryResult, is_system_name, is_tiff_name, position_id +from .metadata import TiffMetadata +from .validation import ValidationIssue, error, has_errors, plural, warning +from .zarr_input import NgffError, ZarrUnavailableError, ngff_label_names, open_ngff_labels + +logger = logging.getLogger(__name__) + +FOLDER = "folder" +"""Source of a mask kept as a TIFF file in the masks folder.""" + +OME_ZARR = "ome-zarr" +"""Source of a mask kept as a label image inside an OME-Zarr image.""" + + +@dataclass(frozen=True) +class MaskFile: + """One mask of one image. + + Attributes + ---------- + image + The image's path relative to the experiment root. + mask_type + The mask type: its folder name in the masks folder, or the label image's name. + path + The mask's path: relative to the masks folder for a TIFF mask; for a label image of + an OME-Zarr image, ``/labels/`` relative to the experiment root. + source + :data:`FOLDER` or :data:`OME_ZARR`. + shape + Shape of the mask without its length-1 axes, ``(Y, X)`` when it is valid. + dtype + Pixel type, or None when the header could not be read. + """ + + image: str + mask_type: str + path: str + source: str + shape: tuple[int, ...] = () + dtype: np.dtype | None = None + + +@dataclass(frozen=True) +class MaskInventory: + """Every mask found for the images of an experiment. + + Attributes + ---------- + folder + The masks folder, as an absolute path; None when only OME-Zarr label images were used. + types + The mask types found, sorted. + files + One entry per mask, sorted by image then type. + issues + Errors and warnings about the masks. + """ + + folder: Path | None + types: tuple[str, ...] = () + files: tuple[MaskFile, ...] = () + issues: tuple[ValidationIssue, ...] = () + + @property + def ok(self) -> bool: + """True if no issue blocks using the masks.""" + return not has_errors(self.issues) + + @property + def images(self) -> tuple[str, ...]: + """Images with at least one mask, sorted.""" + return tuple(sorted({f.image for f in self.files})) + + def for_image(self, image: str) -> Mapping[str, MaskFile]: + """Return the masks of one image, by type (empty when it has none).""" + return MappingProxyType({f.mask_type: f for f in self.files if f.image == image}) + + def to_manifest(self) -> dict: + """Return the inventory as recorded in the run manifest (contract §9).""" + return { + "folder": None if self.folder is None else self.folder.as_posix(), + "types": list(self.types), + "files": [ + {"image": f.image, "type": f.mask_type, "path": f.path, "source": f.source} + for f in self.files + ], + } + + +def find_masks( + masks_folder: str | os.PathLike | None, + discovery: DiscoveryResult, + images: Iterable[TiffMetadata], + *, + ome_zarr_labels: bool = True, +) -> MaskInventory: + """Pair every image with its masks and check them from their headers, without pixels. + + Parameters + ---------- + masks_folder + The masks folder (contract §13), or None to use only the OME-Zarr images' own label + images. + discovery + The experiment's folders, which say where each image is. + images + Header information of the analysed images (excluded images left out); their Y and X + sizes are what each mask must match. + ome_zarr_labels + Also use the label images of OME-Zarr images (needs the ``zarr`` extra; skipped + silently without it). + + Returns + ------- + MaskInventory + The masks found and every issue about them. + """ + images = {meta.relative_path: meta for meta in images} + issues: list[ValidationIssue] = [] + files: list[MaskFile] = [] + folder = None + if masks_folder is not None: + folder = Path(masks_folder).resolve() + issues.extend(_folder_issues(folder, discovery.root)) + if not has_errors(issues): + found, extra = _folder_masks(folder, discovery, images) + files.extend(found) + issues.extend(extra) + if ome_zarr_labels: + files.extend(_label_masks(discovery.root, images)) + issues.extend(_name_issues(files)) + files.sort(key=lambda f: (f.image, f.mask_type)) + issues.extend(_duplicate_issues(files)) + issues.extend(_header_issues(files, images, folder, discovery.root)) + types = tuple(sorted({f.mask_type for f in files})) + issues.extend(_missing_issues(types, files, images)) + return MaskInventory( + folder=folder, + types=types, + files=tuple(files), + issues=tuple(sorted(issues, key=lambda issue: (not issue.is_error, issue.code))), + ) + + +def read_mask(mask: MaskFile, masks_folder: Path | None, experiment_root: Path) -> np.ndarray: + """Read one mask's pixels as a 2D array (its length-1 axes removed). + + Raises + ------ + ValueError + If the mask is not 2D once its length-1 axes are removed. + """ + if mask.source == FOLDER: + pixels = tifffile.imread(Path(masks_folder) / mask.path) + else: + name = mask.path.rsplit("/labels/", 1)[1] + image = experiment_root / mask.path.rsplit("/labels/", 1)[0] + pixels = np.asarray(open_ngff_labels(image, name).array) + pixels = np.squeeze(pixels) + if pixels.ndim != 2: + raise ValueError(f"the mask {mask.path} is not 2D ({pixels.shape})") + return pixels + + +def _folder_issues(folder: Path, experiment_root: Path) -> list[ValidationIssue]: + if not folder.is_dir(): + return [ + error( + "masks_folder_not_found", + f"The masks folder {folder} does not exist or is not a folder.", + action="Select the folder that holds one subfolder per mask type.", + ) + ] + if folder == experiment_root or experiment_root in folder.parents: + return [ + error( + "masks_inside_experiment", + "The masks folder is inside the experiment folder, where its folders would be " + "read as conditions and repeats.", + [folder.relative_to(experiment_root).as_posix() or "."], + "Move the masks folder next to the experiment folder and select it there.", + ) + ] + return [] + + +def _folder_masks( + folder: Path, discovery: DiscoveryResult, images: Mapping[str, TiffMetadata] +) -> tuple[list[MaskFile], list[ValidationIssue]]: + """Pair the TIFF files of every mask-type folder with the images.""" + wanted = {} + for leaf in discovery.leaves: + for image in leaf.images: + if image in images: + prefix = "" if leaf.relative_path == "." else leaf.relative_path + "/" + wanted[(prefix + position_id(image.rsplit("/", 1)[-1])).casefold()] = image + files, unpaired, loose = [], [], [] + for entry in sorted(folder.iterdir()): + if is_system_name(entry.name): + continue + if not entry.is_dir(): + loose.append(entry.name) + continue + for path in sorted(entry.rglob("*")): + relative = path.relative_to(entry) + if any(is_system_name(part) for part in relative.parts) or not path.is_file(): + continue + key = relative.with_name(position_id(path.name)).as_posix().casefold() + mask_path = path.relative_to(folder).as_posix() + if is_tiff_name(path.name) and key in wanted: + files.append(MaskFile(wanted[key], entry.name, mask_path, FOLDER)) + else: + unpaired.append(mask_path) + issues = [] + if loose: + issues.append( + warning( + "mask_outside_type_folder", + f"{plural(len(loose), 'file')} in the masks folder {_are(loose)} not in a " + "mask-type folder and will not be used.", + loose, + "Put each mask in the folder of its type, e.g. Nuclei//.tif.", + ) + ) + if unpaired: + issues.append( + warning( + "mask_without_image", + f"{plural(len(unpaired), 'file')} in the masks folder {_match(unpaired)} no " + "analysed image (same folders and image name) and will not be used.", + unpaired, + "Check the folder names and file names against the experiment, or remove " + "masks of excluded images.", + ) + ) + return files, issues + + +def _label_masks(root: Path, images: Mapping[str, TiffMetadata]) -> list[MaskFile]: + """Return the label images of every OME-Zarr image, as masks.""" + files = [] + for image, meta in images.items(): + if meta.kind != "ome-zarr" or not meta.readable: + continue + try: + names = ngff_label_names(root / image) + except (ZarrUnavailableError, NgffError) as exc: + logger.info("No label images read from %s: %s", image, exc) + continue + files.extend(MaskFile(image, name, f"{image}/labels/{name}", OME_ZARR) for name in names) + return files + + +def _name_issues(files: list[MaskFile]) -> list[ValidationIssue]: + by_folded: dict[str, set[str]] = {} + for f in files: + by_folded.setdefault(f.mask_type.casefold(), set()).add(f.mask_type) + clashes = sorted(name for names in by_folded.values() if len(names) > 1 for name in names) + if not clashes: + return [] + return [ + error( + "mask_names_differ_in_case", + "Some mask types have names that differ only in upper and lower case, so they " + "cannot be told apart: " + ", ".join(clashes) + ".", + clashes, + "Give every mask type its own name.", + ) + ] + + +def _duplicate_issues(files: list[MaskFile]) -> list[ValidationIssue]: + seen: dict[tuple[str, str], MaskFile] = {} + duplicates = [] + for f in files: + key = (f.image, f.mask_type.casefold()) + if key in seen: + duplicates += [seen[key].path, f.path] + seen[key] = f + if not duplicates: + return [] + return [ + error( + "duplicate_mask", + "Some images have two masks of the same type (for example a file in the masks " + "folder and a label image of the same name, or two files with the same name and " + "different extensions).", + sorted(set(duplicates)), + "Keep one mask per image and type.", + ) + ] + + +def _header_issues( + files: list[MaskFile], images: Mapping[str, TiffMetadata], folder: Path | None, root: Path +) -> list[ValidationIssue]: + """Read each mask's header (shape, pixel type) and check it against its image.""" + problems: dict[str, list[str]] = {} + for index, mask in enumerate(files): + try: + shape, dtype = _mask_header(mask, folder, root) + except ( + OSError, + ValueError, + tifffile.TiffFileError, + NgffError, + ZarrUnavailableError, + ) as exc: + logger.info("Mask %s could not be opened: %s", mask.path, exc) + problems.setdefault("unreadable_mask", []).append(mask.path) + continue + files[index] = MaskFile(mask.image, mask.mask_type, mask.path, mask.source, shape, dtype) + if len(shape) != 2: + problems.setdefault("mask_not_2d", []).append(mask.path) + elif shape != _image_yx(images[mask.image]): + problems.setdefault("mask_size_mismatch", []).append(mask.path) + if dtype.kind not in "biu": + problems.setdefault("mask_pixel_type", []).append(mask.path) + return [ + error(code, _TEXTS[code][0], paths, _TEXTS[code][1]) + for code, paths in sorted(problems.items()) + ] + + +_TEXTS = { + "unreadable_mask": ( + "Some masks could not be opened.", + "Check that they are TIFF files (or OME-Zarr label images) with one image each and " + "not damaged; the technical details are in the log.", + ), + "mask_not_2d": ( + "Some masks are not 2D: a mask must be one label image per position.", + "Save each mask as a single 2D label image.", + ), + "mask_size_mismatch": ( + "Some masks do not have the size (Y and X) of their image.", + "Check that each mask was made from its image, at full resolution.", + ), + "mask_pixel_type": ( + "Some masks do not hold whole numbers: mask pixels must be integers (labels) or " + "true/false.", + "Save the masks as integer label images.", + ), +} + + +def _mask_header(mask: MaskFile, folder: Path | None, root: Path) -> tuple: + """Return a mask's shape without its length-1 axes, and its pixel type; no pixels read.""" + if mask.source == FOLDER: + with tifffile.TiffFile(Path(folder) / mask.path) as tif: + if len(tif.series) != 1: + raise ValueError(f"{len(tif.series)} image series") + shape, dtype = tif.series[0].shape, tif.series[0].dtype + else: + image, name = mask.path.rsplit("/labels/", 1) + array = open_ngff_labels(root / image, name).array + shape, dtype = array.shape, np.dtype(array.dtype) + return tuple(int(n) for n in shape if n != 1), dtype.newbyteorder("=") + + +def _image_yx(meta: TiffMetadata) -> tuple[int, int] | None: + """Return the image's Y and X sizes from its header, or None when they are not known.""" + if not meta.readable or "Y" not in meta.axes or "X" not in meta.axes: + return None + return (meta.shape[meta.axes.index("Y")], meta.shape[meta.axes.index("X")]) + + +def _missing_issues( + types: tuple[str, ...], files: list[MaskFile], images: Mapping[str, TiffMetadata] +) -> list[ValidationIssue]: + have = {(f.image, f.mask_type) for f in files} + issues = [] + for mask_type in types: + missing = sorted(image for image in images if (image, mask_type) not in have) + if missing: + issues.append( + warning( + "mask_missing", + f"{plural(len(missing), 'image')} {_have(missing)} no {mask_type} mask; " + f"their {mask_type} objects cannot be measured.", + missing, + f"Add the {mask_type} masks of these images, or exclude the images.", + ) + ) + return issues + + +def _are(items: list) -> str: + return "is" if len(items) == 1 else "are" + + +def _match(items: list) -> str: + return "matches" if len(items) == 1 else "match" + + +def _have(items: list) -> str: + return "has" if len(items) == 1 else "have" diff --git a/src/morphobricks_preprocessing/report.py b/src/morphobricks_preprocessing/report.py index 1c342e83..f47286d3 100644 --- a/src/morphobricks_preprocessing/report.py +++ b/src/morphobricks_preprocessing/report.py @@ -1,4 +1,4 @@ -"""Plain-language summaries of a preflight, a run and the group statistics. +"""Plain-language summaries of a preflight, a run, the group statistics and the masks. The command-line interface prints these texts and the notebook shows them, so both tell the user the same thing. Nothing here reads images or writes files. @@ -18,6 +18,8 @@ from .dataset import BalanceTable from .design import ExperimentDesign from .discovery import DiscoveryResult +from .mask_scan import MaskReport, PairRelation +from .masks import MaskInventory from .run import Preflight, RunResult from .validation import Severity, ValidationIssue, plural @@ -328,3 +330,96 @@ def _number(value: float | None) -> str: if value is None: return "–" return f"{value:.6g}" + + +_MASK_KINDS = { + "instances": "one object per label", + "grouped": "each label in several pieces (for example organelles grouped by cell)", + "one_value": "one value per image (a binary mask or a single object)", + "empty": "no objects", +} +_RELATIONS = { + "inside": "inside", + "grouped_by": "grouped by the labels of", + "same_objects": "the same objects as", +} + + +def mask_summary(masks: MaskInventory) -> str: + """Summarise the masks found: each mask type and how many images have it.""" + if not masks.types: + return "Masks: none found." + lines = [ + f"Masks: {plural(len(masks.types), 'mask type')} for {plural(len(masks.images), 'image')}:" + ] + for name in masks.types: + found = [f for f in masks.files if f.mask_type == name] + sources = " and ".join(sorted({f.source for f in found})) + lines.append(f" {name}: {plural(len(found), 'image')} ({sources})") + return "\n".join(lines) + + +def mask_report_summary(report: MaskReport) -> str: + """Describe the mask types, their relationships and the hierarchy in plain words. + + Every statement is a suggestion from the masks' pixels; the text ends by asking the + user to name the masks and check the relationships. + """ + lines = [f"Mask types ({plural(report.images, 'image')} scanned):"] + for name, summary in sorted(report.types.items()): + lines.append( + f" {name}: {_MASK_KINDS[summary.kind]}; {plural(summary.objects, 'object')}, " + f"median {summary.median_pieces_per_label:g} piece(s) per label" + ) + relationships = [pair for pair in report.relationships if pair.relation != "same_objects"] + if relationships: + lines += ["", "Relationships found:"] + lines += [_relationship_line(pair) for pair in relationships] + if report.same_objects: + lines += ["", "Masks that look like two segmentations of the same objects:"] + lines += [f" {first} and {second}" for first, second in report.same_objects] + lines += ["", "Hierarchy (direct parents only):"] + lines += _hierarchy_lines(report.hierarchy) + lines += [ + "", + "These are suggestions from the masks' pixels: name each mask, then check the " + "relationships before using them.", + ] + return "\n".join(lines) + + +def _relationship_line(pair: PairRelation) -> str: + """One relationship with its evidence, e.g. how many parents have 0, 1, 2, 3+ children.""" + e = pair.evidence + text = ( + f" {pair.child} {_RELATIONS[pair.relation]} {pair.parent}" + f"{'' if pair.direct else ' (through another mask)'}: " + f"{e['inside']} of {plural(e['children'], 'object')} inside; " + f"{pair.parent} with 0/1/2/3+ {pair.child}: {e['parents_with_0_children']}/" + f"{e['parents_with_1_child']}/{e['parents_with_2_children']}/" + f"{e['parents_with_3_or_more_children']}" + ) + extra = [] + if e["outside_parents"]: + extra.append(f"{e['outside_parents']} outside every {pair.parent}") + if e["split_between_parents"]: + extra.append(f"{e['split_between_parents']} split between {pair.parent}") + return text + (f"; {', '.join(extra)}" if extra else "") + + +def _hierarchy_lines(hierarchy) -> list[str]: + """Return the mask types as an indented tree; a type with two parents shows twice.""" + children: dict[str, list[str]] = {name: [] for name in hierarchy} + for name, parents in hierarchy.items(): + for parent in parents: + children[parent].append(name) + lines: list[str] = [] + + def show(name: str, depth: int) -> None: + lines.append(f"{' ' * (depth + 1)}{name}") + for child in sorted(children[name]): + show(child, depth + 1) + + for name in sorted(n for n, parents in hierarchy.items() if not parents): + show(name, 0) + return lines diff --git a/src/morphobricks_preprocessing/run.py b/src/morphobricks_preprocessing/run.py index 1d5ada94..15b42f11 100644 --- a/src/morphobricks_preprocessing/run.py +++ b/src/morphobricks_preprocessing/run.py @@ -59,6 +59,7 @@ read_manifest, write_manifest, ) +from .masks import MaskInventory, find_masks from .metadata import MetadataReport, TiffMetadata, inspect_images from .planes import ( DEFAULT_MEMORY_LIMIT, @@ -104,6 +105,9 @@ class Preflight: dataset Images per condition and biological repeat, and the warnings about their balance. + masks + The masks of the images (contract §13), found from their headers, when + masks were asked for; None otherwise. """ discovery: DiscoveryResult @@ -112,6 +116,7 @@ class Preflight: exclusion_issues: tuple[ValidationIssue, ...] = () suggested_mapping: AxisMapping | None = None dataset: DatasetReport = field(default_factory=DatasetReport) + masks: MaskInventory | None = None @property def issues(self) -> tuple[ValidationIssue, ...]: @@ -121,6 +126,7 @@ def issues(self) -> tuple[ValidationIssue, ...]: *self.exclusion_issues, *self.metadata.issues, *self.dataset.issues, + *(self.masks.issues if self.masks is not None else ()), ] return tuple(sorted(merged, key=lambda issue: not issue.is_error)) @@ -258,6 +264,8 @@ def preflight( output_folder: str | os.PathLike | None = None, *, exclude: Iterable[str] = (), + masks_folder: str | os.PathLike | None = None, + ome_zarr_labels: bool = False, ) -> Preflight: """Pass 1: check the folder and the image headers without reading pixels. @@ -270,6 +278,14 @@ def preflight( exclude Relative paths of images to leave out (``/`` or backslash separators), for example damaged files. + masks_folder + The masks folder (contract §13): one subfolder per mask type, each + mirroring the experiment's folders. The masks are paired with the images + and checked from their headers; the label images of OME-Zarr images are + included. + ome_zarr_labels + Use the label images of OME-Zarr images as masks, also without a masks + folder. Returns ------- @@ -281,6 +297,9 @@ def preflight( excluded = tuple(sorted({_normalise(path) for path in exclude})) images = [path for path in discovery.images if path not in excluded] metadata = inspect_images(root, images) + masks = None + if masks_folder is not None or ome_zarr_labels: + masks = find_masks(masks_folder, discovery, metadata.files) return Preflight( discovery=discovery, metadata=metadata, @@ -288,6 +307,7 @@ def preflight( exclusion_issues=tuple(_exclusion_issues(discovery, excluded)), suggested_mapping=prefill_experiment_mapping(metadata), dataset=check_dataset(discovery, excluded=excluded), + masks=masks, ) @@ -305,6 +325,8 @@ def run_experiment( convert_tiffs: bool = False, conversion: ConversionSettings | None = None, design: ExperimentDesign | None = None, + masks_folder: str | os.PathLike | None = None, + ome_zarr_labels: bool = False, ) -> RunResult: """Pass 2: write or update the histogram table and the run manifest. @@ -350,6 +372,10 @@ def run_experiment( with the table, and its SHA-256 recorded in the manifest. Without a design nothing is asked, the manifest records none, and an existing ``experiment_design.csv`` is left as it is. + masks_folder, ome_zarr_labels + As :func:`preflight`. The masks are recorded in the manifest (``masks``, + contract §9); an error about them stops the run. Their pixels are not + read: the mask step (:func:`.mask_scan.scan_masks`) does that. Returns ------- @@ -358,7 +384,9 @@ def run_experiment( """ root = Path(experiment_root).resolve() output = default_output_folder(root) if output_folder is None else Path(output_folder).resolve() - pre = preflight(root, output, exclude=exclude) + pre = preflight( + root, output, exclude=exclude, masks_folder=masks_folder, ome_zarr_labels=ome_zarr_labels + ) review = review_axes(pre.metadata, mapping, missing_as_singleton=missing_as_singleton) dataset = pre.dataset design_problems: list[ValidationIssue] = [] @@ -373,6 +401,7 @@ def run_experiment( *review.combined_issues(pre.metadata), *design_problems, *dataset.issues, + *(pre.masks.issues if pre.masks is not None else ()), ] dtype = pre.metadata.dtype bins = None @@ -557,6 +586,7 @@ def stop(**counts) -> RunResult: versions=package_versions(), written_at=datetime.now(UTC).isoformat(timespec="seconds"), experiment_design_sha256=design_sha256, + masks=None if pre.masks is None else pre.masks.to_manifest(), ) write_manifest(manifest, staged / MANIFEST_FILENAME) logger.info( diff --git a/tests/test_mask_scan.py b/tests/test_mask_scan.py new file mode 100644 index 00000000..177b08ae --- /dev/null +++ b/tests/test_mask_scan.py @@ -0,0 +1,247 @@ +"""The mask step: kinds of masks, relationships with their evidence, hierarchy (contract §13).""" + +import json +from collections import deque + +import numpy as np + +from morphobricks_preprocessing import ( + pair_counts, + piece_counts, + preflight, + scan_masks, + write_mask_report, +) +from tiff_tree import write_data + +SHAPE = (64, 64) +# Cells in four quadrants; cell 1 has two nuclei, cell 3 none. +CELLS = {1: (2, 30, 2, 30), 2: (2, 30, 34, 62), 3: (34, 62, 2, 30), 4: (34, 62, 34, 62)} +NUCLEI = {10: (5, 11, 5, 11), 11: (18, 24, 18, 24), 12: (10, 16, 40, 46), 13: (40, 46, 40, 46)} + + +def boxes(spec): + mask = np.zeros(SHAPE, dtype=np.uint16) + for label, (r0, r1, c0, c1) in spec.items(): + mask[r0:r1, c0:c1] = label + return mask + + +def hierarchy_masks(): + """Organoid > Cells > Nuclei > Nucleoli; NucleiAlt = the nuclei segmented again; Vesicles + = single pixels carrying their cell's label; Left and Top = two overlapping regions.""" + nucleoli = {100 + n: (r0 + 2, r0 + 4, c0 + 2, c0 + 4) for n, (r0, _, c0, _) in NUCLEI.items()} + alt = {200 + n: (r0 - 1, r1 + 1, c0 - 1, c1 + 1) for n, (r0, r1, c0, c1) in NUCLEI.items()} + vesicles = np.zeros(SHAPE, dtype=np.uint16) + for label, (r0, _, c0, _) in CELLS.items(): + for dr, dc in ((1, 1), (1, 20), (20, 1), (20, 20)): + vesicles[r0 + dr, c0 + dc] = label + organoid = np.zeros(SHAPE, dtype=np.uint8) + organoid[1:63, 1:63] = 1 + return { + "Organoid": organoid, + "Cells": boxes(CELLS), + "Nuclei": boxes(NUCLEI), + "NucleiAlt": boxes(alt), + "Nucleoli": boxes(nucleoli), + "Vesicles": vesicles, + } + + +def write_experiment(tmp_path, masks, images=("ctrl/a.tif", "drug/a.tif")): + root, folder = tmp_path / "exp", tmp_path / "masks" + for rel in images: + write_data(root / rel, np.zeros(SHAPE, dtype=np.uint8)) + for name, mask in masks.items(): + write_data(folder / name / rel, mask) + return root, folder + + +def scan(tmp_path, masks): + root, folder = write_experiment(tmp_path, masks) + inventory = preflight(root, masks_folder=folder).masks + assert inventory.ok, inventory.issues + return scan_masks(inventory, root) + + +def by_pair(report): + return {(p.child, p.parent): p for p in report.pairs} + + +# Pieces ------------------------------------------------------------------------------- + + +def reference_pieces(mask): + """Pieces per label by breadth-first search (independent of the brick's code).""" + seen = np.zeros(mask.shape, dtype=bool) + pieces = {} + for start in zip(*np.nonzero(mask), strict=True): + if seen[start]: + continue + label = mask[start] + pieces[int(label)] = pieces.get(int(label), 0) + 1 + queue = deque([start]) + seen[start] = True + while queue: + r, c = queue.popleft() + for nr, nc in ((r + 1, c), (r - 1, c), (r, c + 1), (r, c - 1)): + inside = 0 <= nr < mask.shape[0] and 0 <= nc < mask.shape[1] + if inside and not seen[nr, nc] and mask[nr, nc] == label: + seen[nr, nc] = True + queue.append((nr, nc)) + return pieces + + +def test_pieces_are_4_connected_and_counted_per_label(): + rng = np.random.default_rng(0) + for _ in range(50): + mask = rng.integers(0, 4, size=tuple(rng.integers(1, 20, 2))).astype(np.int32) + labels, counts = piece_counts(mask) + assert dict(zip(labels.tolist(), counts.tolist(), strict=True)) == reference_pieces(mask) + diagonal = np.eye(3, dtype=np.uint8) # touching only at corners: three pieces + assert piece_counts(diagonal)[1].tolist() == [3] + assert piece_counts(np.zeros((3, 3), np.uint8))[0].size == 0 + assert piece_counts(np.full((2, 2), 7, np.uint8))[1].tolist() == [1] + + +# Pair counts -------------------------------------------------------------------------- + + +def test_several_children_per_parent_outside_and_split_children_are_counted(): + cells = boxes({1: (0, 10, 0, 10), 2: (0, 10, 12, 22), 3: (12, 22, 0, 10)}) + nuclei = boxes( + { + 5: (1, 3, 1, 3), # in cell 1 + 6: (5, 7, 5, 7), # in cell 1: two nuclei in one cell + 7: (2, 4, 14, 16), # in cell 2 + 8: (30, 32, 30, 32), # outside every cell + 9: (4, 6, 9, 13), # split between cells 1 and 2 + } + ) + counts = pair_counts(nuclei, cells) + assert counts["children"] == 5 and counts["parents"] == 3 + assert counts["inside"] == 3 and counts["inside_parent_larger"] == 3 + assert counts["outside_parents"] == 1 + assert counts["split_between_parents"] == 1 + assert ( + counts["parents_with_0_children"], + counts["parents_with_1_child"], + counts["parents_with_2_children"], + counts["parents_with_3_or_more_children"], + ) == (1, 1, 1, 0) + assert counts["children_matched_one_to_one"] == 0 + + +def test_pair_counts_of_an_empty_or_disjoint_child(): + cells = boxes({1: (0, 10, 0, 10)}) + assert pair_counts(np.zeros(SHAPE, np.uint16), cells)["parents_with_0_children"] == 1 + far = boxes({3: (40, 42, 40, 42)}) + counts = pair_counts(far, cells) + assert (counts["outside_parents"], counts["parents_with_0_children"]) == (1, 1) + + +# Whole experiment --------------------------------------------------------------------- + + +def test_kinds_of_masks(tmp_path): + report = scan(tmp_path, hierarchy_masks()) + kinds = {name: summary.kind for name, summary in report.types.items()} + assert kinds == { + "Cells": "instances", + "Nuclei": "instances", + "NucleiAlt": "instances", + "Nucleoli": "instances", + "Organoid": "one_value", + "Vesicles": "grouped", + } + assert report.images == 2 + assert report.types["Cells"].objects == 8 # 4 per image + assert report.types["Vesicles"].median_pieces_per_label == 4 + + +def test_relationships_and_their_evidence(tmp_path): + pairs = by_pair(scan(tmp_path, hierarchy_masks())) + nuclei = pairs[("Nuclei", "Cells")] + assert nuclei.relation == "inside" and nuclei.direct + # Over two images: cell 1 has two nuclei, cells 2 and 4 one, cell 3 none. + e = nuclei.evidence + assert (e["children"], e["inside"]) == (8, 8) + assert ( + e["parents_with_0_children"], + e["parents_with_1_child"], + e["parents_with_2_children"], + e["parents_with_3_or_more_children"], + ) == (2, 4, 2, 0) + assert pairs[("Vesicles", "Cells")].relation == "grouped_by" + assert pairs[("Cells", "Organoid")].relation == "inside" + assert pairs[("Nuclei", "NucleiAlt")].relation == "same_objects" + assert pairs[("NucleiAlt", "Nuclei")].relation == "same_objects" + assert pairs[("Cells", "Nuclei")].relation == "none" # parents are never inside children + + +def test_nested_hierarchy_keeps_direct_parents_only(tmp_path): + report = scan(tmp_path, hierarchy_masks()) + pairs = by_pair(report) + # Nucleoli are inside cells too, but through their nuclei. + assert pairs[("Nucleoli", "Cells")].relation == "inside" + assert not pairs[("Nucleoli", "Cells")].direct + assert dict(report.hierarchy) == { + "Cells": ("Organoid",), + "Nuclei": ("Cells",), + "NucleiAlt": ("Cells",), + "Nucleoli": ("Nuclei", "NucleiAlt"), + "Organoid": (), + "Vesicles": ("Cells",), + } + assert report.same_objects == (("Nuclei", "NucleiAlt"),) + + +def test_a_mask_may_have_two_unrelated_direct_parents(tmp_path): + left = np.zeros(SHAPE, np.uint8) + left[:, :40] = 1 + top = np.zeros(SHAPE, np.uint8) + top[:40, :] = 1 + dots = boxes({1: (5, 7, 5, 7), 2: (20, 22, 30, 32), 3: (30, 32, 10, 12)}) + report = scan(tmp_path, {"Left": left, "Top": top, "Dots": dots}) + assert report.hierarchy["Dots"] == ("Left", "Top") + assert by_pair(report)[("Left", "Top")].relation == "none" # same size: not nested + + +def test_label_numbers_shared_with_the_parent_do_not_make_a_mask_grouped(tmp_path): + cells = boxes(CELLS) + nuclei = boxes( + {label: (r0 + 5, r0 + 10, c0 + 5, c0 + 10) for label, (r0, _, c0, _) in CELLS.items()} + ) + pairs = by_pair(scan(tmp_path, {"Cells": cells, "Nuclei": nuclei})) + assert pairs[("Nuclei", "Cells")].relation == "inside" + assert pairs[("Nuclei", "Cells")].evidence["child_pixels_with_parent_label"] == 2 * 4 * 25 + + +def test_report_file(tmp_path): + report = scan(tmp_path, hierarchy_masks()) + path = write_mask_report(report, tmp_path / "out" / "mask_relationships.json") + raw = path.read_bytes() + assert b"\r\n" not in raw + data = json.loads(raw) + assert data["report_version"] == "1" + assert data["thresholds"]["inside_object_share"] == 0.9 + assert data["hierarchy"]["Nucleoli"] == ["Nuclei", "NucleiAlt"] + assert data["same_objects"] == [["Nuclei", "NucleiAlt"]] + assert {r["relation"] for r in data["relationships"]} == { + "inside", + "grouped_by", + "same_objects", + } + assert len(data["pairs"]) == 6 * 5 + assert "morphobricks-preprocessing" in data["versions"] + second = write_mask_report(report, tmp_path / "again.json") + assert json.loads(second.read_bytes())["pairs"] == data["pairs"] + + +def test_progress_is_reported_per_image(tmp_path): + root, folder = write_experiment(tmp_path, {"Cells": boxes(CELLS)}) + calls = [] + scan_masks( + preflight(root, masks_folder=folder).masks, root, progress=lambda *a: calls.append(a) + ) + assert calls == [(0, 2, "ctrl/a.tif"), (1, 2, "drug/a.tif"), (2, 2, "")] diff --git a/tests/test_masks.py b/tests/test_masks.py new file mode 100644 index 00000000..2ae5caa0 --- /dev/null +++ b/tests/test_masks.py @@ -0,0 +1,239 @@ +"""Masks: pairing every image with its masks and checking them from their headers (contract §13).""" + +import json + +import numpy as np +import pytest +import tifffile + +from morphobricks_preprocessing import ( + prefill_mapping, + preflight, + read_manifest, + run_experiment, +) +from morphobricks_preprocessing.cli import main +from tiff_tree import write_data, write_file + +IMAGES = ["ctrl/rep1/a.tif", "ctrl/rep1/b.tif", "drug/rep1/a.tif"] + + +def label_image(seed, shape=(16, 20)): + mask = np.zeros(shape, dtype=np.uint16) + mask[2:8, 3:9] = 1 + seed + mask[10:14, 10:18] = 2 + seed + return mask + + +@pytest.fixture +def experiment(tmp_path): + root = tmp_path / "exp" + for seed, rel in enumerate(IMAGES): + write_data(root / rel, np.full((2, 16, 20), seed, dtype=np.uint8), "CYX") + return root + + +@pytest.fixture +def masks(tmp_path): + folder = tmp_path / "masks" + for kind in ("Cells", "Nuclei"): + for seed, rel in enumerate(IMAGES): + write_data(folder / kind / rel, label_image(seed)) + return folder + + +def inventory(experiment, masks, **kwargs): + pre = preflight(experiment, masks_folder=masks, **kwargs) + return pre.masks + + +def codes(issues): + return [issue.code for issue in issues] + + +@pytest.fixture +def no_pixel_reads(monkeypatch): + def fail(*args, **kwargs): + raise AssertionError("pixel data was read") + + for cls in (tifffile.TiffFile, tifffile.TiffPage, tifffile.TiffFrame, tifffile.TiffPageSeries): + monkeypatch.setattr(cls, "asarray", fail) + + +# Pairing ------------------------------------------------------------------------------ + + +def test_every_image_gets_its_masks_by_type(experiment, masks, no_pixel_reads): + found = inventory(experiment, masks) + assert found.ok and found.issues == () + assert found.types == ("Cells", "Nuclei") + assert found.images == tuple(IMAGES) + by_type = found.for_image("ctrl/rep1/b.tif") + assert sorted(by_type) == ["Cells", "Nuclei"] + assert by_type["Nuclei"].path == "Nuclei/ctrl/rep1/b.tif" + assert (by_type["Nuclei"].shape, by_type["Nuclei"].dtype) == ((16, 20), np.dtype("uint16")) + assert found.for_image("unknown.tif") == {} + + +def test_any_number_of_mask_types(experiment, masks): + for kind in ("Nucleoli", "Vesicles", "Organoid", "Background"): + for seed, rel in enumerate(IMAGES): + write_data(masks / kind / rel, label_image(seed)) + found = inventory(experiment, masks) + assert found.types == ("Background", "Cells", "Nuclei", "Nucleoli", "Organoid", "Vesicles") + assert len(found.files) == 6 * len(IMAGES) + + +def test_masks_match_images_by_position_whatever_their_extension(tmp_path): + root = tmp_path / "exp" + write_data(root / "img_01.ome.tif", np.zeros((16, 20), dtype=np.uint8)) + write_data(tmp_path / "masks" / "Cells" / "IMG_01.TIFF", label_image(0)) + found = inventory(root, tmp_path / "masks") + assert [(f.image, f.path) for f in found.files] == [("img_01.ome.tif", "Cells/IMG_01.TIFF")] + + +def test_missing_and_unpaired_masks_are_warnings(experiment, masks): + (masks / "Nuclei" / "drug" / "rep1" / "a.tif").unlink() + write_data(masks / "Nuclei" / "drug" / "rep2" / "c.tif", label_image(0)) + write_file(masks / "notes.txt") + found = inventory(experiment, masks) + assert found.ok + assert codes(found.issues) == ["mask_missing", "mask_outside_type_folder", "mask_without_image"] + missing, loose, unpaired = found.issues + assert missing.paths == ("drug/rep1/a.tif",) and "no Nuclei mask" in missing.message + assert loose.paths == ("notes.txt",) + assert unpaired.paths == ("Nuclei/drug/rep2/c.tif",) + + +def test_masks_of_excluded_images_are_not_used(experiment, masks): + found = inventory(experiment, masks, exclude=["drug/rep1/a.tif"]) + assert "drug/rep1/a.tif" not in found.images + assert codes(found.issues) == ["mask_without_image"] + + +# Errors ------------------------------------------------------------------------------- + + +def test_the_masks_folder_must_exist_and_be_outside_the_experiment(experiment, tmp_path): + assert codes(inventory(experiment, tmp_path / "nowhere").issues) == ["masks_folder_not_found"] + write_data(experiment / "masks" / "Cells" / "ctrl" / "rep1" / "a.tif", label_image(0)) + inside = inventory(experiment, experiment / "masks") + assert "masks_inside_experiment" in codes(inside.issues) + + +@pytest.mark.parametrize( + ("data", "code"), + [ + (np.zeros((16, 21), dtype=np.uint16), "mask_size_mismatch"), + (np.zeros((2, 16, 20), dtype=np.uint16), "mask_not_2d"), + (np.zeros((16, 20), dtype=np.float32), "mask_pixel_type"), + ], +) +def test_masks_are_checked_from_their_headers(experiment, masks, data, code, no_pixel_reads): + write_data(masks / "Cells" / "ctrl" / "rep1" / "a.tif", data) + found = inventory(experiment, masks) + assert codes(found.issues) == [code] + assert found.issues[0].paths == ("Cells/ctrl/rep1/a.tif",) + assert not found.ok + + +def test_a_length_one_axis_and_boolean_masks_are_accepted(experiment, masks): + write_data(masks / "Cells" / "ctrl" / "rep1" / "a.tif", np.zeros((1, 16, 20), dtype=np.uint8)) + write_data(masks / "Nuclei" / "ctrl" / "rep1" / "a.tif", np.zeros((16, 20), dtype=bool)) + assert inventory(experiment, masks).ok + + +def test_unreadable_masks_are_errors(experiment, masks): + write_file(masks / "Cells" / "ctrl" / "rep1" / "a.tif", b"not a tiff") + assert codes(inventory(experiment, masks).issues) == ["unreadable_mask"] + + +def test_mask_names_must_differ_by_more_than_case(experiment, masks): + write_data(masks / "cells" / "ctrl" / "rep1" / "a.tif", label_image(0)) + found = inventory(experiment, masks) + assert "mask_names_differ_in_case" in codes(found.issues) + + +def test_one_mask_per_image_and_type(experiment, masks): + write_data(masks / "Cells" / "ctrl" / "rep1" / "a.tiff", label_image(0)) + found = inventory(experiment, masks) + duplicate = next(i for i in found.issues if i.code == "duplicate_mask") + assert duplicate.paths == ("Cells/ctrl/rep1/a.tif", "Cells/ctrl/rep1/a.tiff") + + +# OME-Zarr label images ---------------------------------------------------------------- + + +def test_label_images_of_ome_zarr_images_are_masks(tmp_path): + pytest.importorskip("zarr") + from zarr_tree import write_labels, write_ngff + + path = write_ngff(tmp_path / "exp" / "ctrl" / "a.ome.zarr", np.zeros((16, 20), np.uint8), "yx") + write_labels(path, {"cells": label_image(0), "nuclei": label_image(1)}, "yx") + found = preflight(tmp_path / "exp", ome_zarr_labels=True).masks + assert found.ok + assert found.types == ("cells", "nuclei") + assert [f.path for f in found.files] == [ + "ctrl/a.ome.zarr/labels/cells", + "ctrl/a.ome.zarr/labels/nuclei", + ] + assert preflight(tmp_path / "exp").masks is None # only when masks are asked for + + +# Run, manifest and command line ------------------------------------------------------- + + +def test_the_run_records_the_masks_in_the_manifest(experiment, masks, tmp_path): + out = tmp_path / "out" + result = run_experiment( + experiment, prefill_mapping("CYX"), output_folder=out, masks_folder=masks + ) + assert result.written + recorded = read_manifest(out / "run_manifest.json").masks + assert recorded["types"] == ["Cells", "Nuclei"] + assert recorded["files"][0] == { + "image": "ctrl/rep1/a.tif", + "type": "Cells", + "path": "Cells/ctrl/rep1/a.tif", + "source": "folder", + } + # Without masks the manifest has no masks key, as before. + run_experiment(experiment, prefill_mapping("CYX"), output_folder=out) + data = json.loads((out / "run_manifest.json").read_text(encoding="utf-8")) + assert "masks" not in data and read_manifest(out / "run_manifest.json").masks is None + + +def test_a_mask_error_stops_the_run(experiment, masks, tmp_path): + write_data(masks / "Cells" / "ctrl" / "rep1" / "a.tif", np.zeros((3, 3), dtype=np.uint8)) + result = run_experiment( + experiment, prefill_mapping("CYX"), output_folder=tmp_path / "out", masks_folder=masks + ) + assert not result.written + assert "mask_size_mismatch" in codes(result.issues) + + +def test_check_lists_the_masks(experiment, masks, capsys): + assert main(["check", str(experiment), "--masks", str(masks)]) == 0 + out = capsys.readouterr().out + assert "Masks: 2 mask types for 3 images:" in out + assert " Nuclei: 3 images (folder)" in out + + +def test_masks_command_writes_the_report(experiment, masks, tmp_path, capsys): + out = tmp_path / "out" + assert main(["masks", str(experiment), "--masks", str(masks), "--output", str(out)]) == 0 + printed = capsys.readouterr().out + assert "[3/3] drug/rep1/a.tif" in printed + assert "Hierarchy (direct parents only):" in printed + report = json.loads((out / "mask_relationships.json").read_text(encoding="utf-8")) + assert report["report_version"] == "1" + assert sorted(report["masks"]) == ["Cells", "Nuclei"] + assert report["same_objects"] == [["Cells", "Nuclei"]] + + +def test_masks_command_stops_on_mask_errors_or_without_masks(experiment, masks, capsys): + assert main(["masks", str(experiment)]) == 1 + assert "Masks: none found." in capsys.readouterr().out + write_data(masks / "Cells" / "ctrl" / "rep1" / "a.tif", np.zeros((3, 3), dtype=np.uint8)) + assert main(["masks", str(experiment), "--masks", str(masks)]) == 1 + assert "do not have the size" in capsys.readouterr().out From d905baccdf02d933bc1a6a7b8f62adbe63305b6c Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 15:30:23 +0000 Subject: [PATCH 2/5] Skip the mask-name case test on case-insensitive file systems Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01U6vJfhGWQ7pTpsJhw48Zf7 --- tests/test_masks.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_masks.py b/tests/test_masks.py index 2ae5caa0..86b27db1 100644 --- a/tests/test_masks.py +++ b/tests/test_masks.py @@ -149,6 +149,10 @@ def test_unreadable_masks_are_errors(experiment, masks): def test_mask_names_must_differ_by_more_than_case(experiment, masks): + try: + (masks / "cells").mkdir() + except FileExistsError: + pytest.skip("the file system does not distinguish upper and lower case") write_data(masks / "cells" / "ctrl" / "rep1" / "a.tif", label_image(0)) found = inventory(experiment, masks) assert "mask_names_differ_in_case" in codes(found.issues) From 8177b8731b69eabd4ce876fac9d9d30507d9276a Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 15:56:19 +0000 Subject: [PATCH 3/5] Let users change the mask-relationship thresholds (owner's answer, Q29) The seven thresholds become defaults (DEFAULT_THRESHOLDS): scan_masks and pair_counts take thresholds=, the masks command one option each (--inside-object-share 0.8, ...). mask_thresholds checks them before any mask is read (shares and IoU in (0, 1], median pieces >= 1), and the values used are written in the report and printed. The owner confirmed 2D masks at the image's YX, integer/bool pixels and the masks folder outside the experiment; plan Q29 records it, with the v0.2.0 tag proposal and the next steps (notebook step, saved answers). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01U6vJfhGWQ7pTpsJhw48Zf7 --- AGENTS.md | 5 +- docs/contract.md | 12 ++- docs/plan.md | 17 +++- docs/user_guide.md | 14 ++- src/morphobricks_preprocessing/__init__.py | 4 + src/morphobricks_preprocessing/cli.py | 38 +++++++- src/morphobricks_preprocessing/mask_scan.py | 98 ++++++++++++++++----- src/morphobricks_preprocessing/report.py | 9 +- tests/test_mask_scan.py | 56 ++++++++++++ tests/test_masks.py | 12 +++ 10 files changed, 231 insertions(+), 34 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 519c9fab..819e2f1d 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -95,7 +95,8 @@ src/morphobricks_preprocessing/ checks (2D, image YX, integer/bool), MaskInventory (to_manifest), read_mask mask_scan.py the optional mask step (Q29): piece_counts (scipy.ndimage), pair_counts, scan_masks -> MaskReport (kinds, relations, direct-parent hierarchy, evidence), - write_mask_report (mask_relationships.json); fixed THRESHOLDS + write_mask_report (mask_relationships.json); DEFAULT_THRESHOLDS, + mask_thresholds (checked changes) zarr_input.py OME-Zarr input: open_ngff (OME-NGFF 0.4/0.5 metadata -> level-0 array, axis letters and physical sizes from the scale and axis units), ngff_label_names/open_ngff_labels (label images, for later bricks), @@ -414,7 +415,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | Dataset consistency (Q21, P2) | Step 1 "Dataset consistency" panel: balance table (images and technical-repeat folders per condition x biological repeat), `unbalanced_images` (< 0.5 or > 2 x median), `repeat_missing_conditions`, `different_pixel_sizes`, `missing_pixel_size` (OME-Zarr scale included), `different_plane_counts`; `different_sizes` narrowed to width/height; consistency issues shown only in that panel, not in Messages. Warnings only; outputs unchanged | | Experiment design (Q22, P3) | Step 1 after the labels: pairing dropdowns, one yellow "Yes, they are comparable" question per repeat group (single whole-experiment question without biological-repeat folders), control checkboxes; `experiment_design.csv` written after every answer; reruns load it and fold the questions behind a Change button; step 3 gated. CLI `run` reads it or `--design`, else writes a template and exits 1. `run_experiment(design=...)` optional; manifest `experiment_design_sha256` (optional key, version 1). Two folders of one condition in a group block that group | | Position identity (Q20, D02) | `position_id` = image name without extension, in every CSV row and manifest file entry; FOV key = labels + `position_id`; duplicates in one leaf folder are the error `duplicate_position_id`. CSV columns and schema 2 unchanged; manifest version 1 with an optional key | -| Masks (Q29) | Optional masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, paired by leaf folder + `position_id`; OME-Zarr label images are masks too. Header-only errors (`masks_folder_not_found`, `masks_inside_experiment`, `mask_names_differ_in_case`, `duplicate_mask`, `unreadable_mask`, `mask_not_2d`, `mask_size_mismatch`, `mask_pixel_type`) and warnings (`mask_missing`, `mask_without_image`, `mask_outside_type_folder`) in `check`/`run`; manifest key `masks` (optional, version 1). Optional mask step (`masks` command, `scan_masks`): kinds, relations with evidence, direct-parent hierarchy, `mask_relationships.json` (report version 1); fixed thresholds; suggestions only, the user names and confirms. SciPy core dependency. Notebook step later. To confirm: 2D masks, integer/bool pixels, fixed thresholds | +| Masks (Q29) | Optional masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, paired by leaf folder + `position_id`; OME-Zarr label images are masks too. Header-only errors (`masks_folder_not_found`, `masks_inside_experiment`, `mask_names_differ_in_case`, `duplicate_mask`, `unreadable_mask`, `mask_not_2d`, `mask_size_mismatch`, `mask_pixel_type`) and warnings (`mask_missing`, `mask_without_image`, `mask_outside_type_folder`) in `check`/`run`; manifest key `masks` (optional, version 1). Optional mask step (`masks` command, `scan_masks`): kinds, relations with evidence, direct-parent hierarchy, `mask_relationships.json` (report version 1); default thresholds the user may change (`thresholds=` keyword, one CLI option each, checked by `mask_thresholds`, written in the report); suggestions only, the user names and confirms. SciPy core dependency. Confirmed by the owner: 2D masks at the image's YX, integer/bool pixels, masks folder outside the experiment. Next: notebook step and a file of the user's answers; 3D/time-lapse masks later; tag v0.2.0 proposed after merge | | Notebook style | Follow the lab's NucleiSky and LabChronicle notebooks (see §7) | | Widget widths | **99%**, never 100% | | Paths and filenames in widget HTML and notebook text | **`…` / `**…**`, never `` or inline code** | diff --git a/docs/contract.md b/docs/contract.md index 0d4286aa..f23b3ea3 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -283,7 +283,7 @@ the user names the masks and confirms or corrects the relationships. For a child type and a parent type, the counts are summed over the images with both masks. Each child object's parent is the parent object sharing most of its pixels. -| Relation of child to parent | Rule (fixed thresholds, written in the report) | +| Relation of child to parent | Rule (default thresholds; the values used are written in the report) | |---|---| | `same_objects` | At least 90% of the children and 90% of the parents are matched one to one at an intersection over union of at least 0.5: two segmentations of the same objects. | | `grouped_by` | The child is `grouped`, at least 90% of its objects are inside as below, and at least 99% of its pixels carry the label of the parent object under them. | @@ -299,6 +299,14 @@ Each child object's parent is the parent object sharing most of its pixels. | `child_pixels`, `child_pixels_with_parent_label` | Child pixels; those whose label equals the parent label under them. | | `children_matched_one_to_one`, `parents_matched_one_to_one` | Objects with a partner at an intersection over union of at least 0.5. | +**Thresholds.** The numbers above are defaults (`DEFAULT_THRESHOLDS`), and each may be +changed: `inside_pixel_share` (0.5), `inside_object_share` (0.9), `parent_larger_share` +(0.9), `same_objects_iou` (0.5), `same_objects_share` (0.9), `grouped_label_share` (0.99) +and `grouped_median_pieces` (2). Shares and the intersection over union must be above 0 +and at most 1, the median number of pieces at least 1; an unknown name or another value +raises `ValueError` before any mask is read (`mask_thresholds`). The command line takes +one option per threshold (`--inside-object-share 0.8`, ...). + **Hierarchy.** A relation `inside` or `grouped_by` is **direct** unless another type lies between: nucleoli inside nuclei inside cells have the direct parent nuclei, and the relation nucleoli inside cells is kept in the report but not in the hierarchy. A type may @@ -310,7 +318,7 @@ The report is UTF-8 JSON (two-space indent, LF), written atomically to the outpu | Key | Meaning | |---|---| | `report_version` | `1`. | -| `thresholds` | The thresholds above. | +| `thresholds` | The thresholds used, all seven. | | `images` | Images scanned (every image with at least one mask). | | `masks` | Per mask type: `kind`, `images`, `objects` (labels over all images) and `median_pieces_per_label`. | | `hierarchy` | Per mask type, its direct parents (sorted; empty for none). | diff --git a/docs/plan.md b/docs/plan.md index 52d6aef9..3e230001 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -713,7 +713,16 @@ Raised during review; to be decided before the milestones that need them: line and documentation first; the notebook step follows in a later change. SciPy becomes a core dependency (the owner's agreement: "can use scipy if good reasons"): counting the connected pieces of every label of one 2048 x 2048 HPMEC cell mask took - about 5 s with a NumPy union-find, and 0.3 to 0.9 s per HPMEC mask with `scipy.ndimage`. Proposed defaults, to confirm in - review: masks are 2D with the Y and X sizes of their image (3D masks later), their - pixels integers or booleans, and the thresholds fixed (written in every report), not - settings.* + about 5 s with a NumPy union-find, and 0.3 to 0.9 s per HPMEC mask with + `scipy.ndimage`. Confirmed by the owner (27 October 2026): masks are 2D with the Y and + X sizes of their image (3D and time-lapse masks are a later follow-up), their pixels + integers or booleans, and the masks folder outside the experiment. The thresholds are + defaults the user may change, `scan_masks(..., thresholds=...)` in Python and one + option per threshold on the `masks` command (`--inside-object-share 0.8`, ...); + `mask_thresholds` checks them (shares and the intersection over union above 0 and at + most 1, the median number of pieces at least 1) before any mask is read, and the values + used are written in every report. Proposed to the brick's maintainer: tag v0.2.0 once + this is merged, so the other bricks can pin it (the follow-ups P4 and P5 then go to + v0.3.0). Next, in their own changes: a notebook step for the masks (name each mask, + say what it contains, confirm or correct the relationships) and a file saving those + answers for the other bricks, as `experiment_design.csv` does for the repeats.* diff --git a/docs/user_guide.md b/docs/user_guide.md index 48e41aea..8b46c112 100644 --- a/docs/user_guide.md +++ b/docs/user_guide.md @@ -357,6 +357,12 @@ pixels, not decisions: you name each kind of mask and check the relationships in that use them. The report is written to `mask_relationships.json` in the output folder (§6). The command takes `--output`, `--exclude` and `--log` like the others, and exits with status 1 when an error blocks it or when no masks are found. + +The rules behind each suggestion use thresholds you can change, for example +`--inside-object-share 0.8` to call a kind of mask "inside" another when 80% of its objects +(not 90%) lie inside, or `--same-objects-iou 0.7` to ask for a closer match before two +kinds of mask count as the same objects. `morphobricks-preprocessing masks --help` lists +all seven with their defaults; the values used are printed and written in the report. `python -m morphobricks_preprocessing` works the same way. ## 6. Output files @@ -532,9 +538,11 @@ report.hierarchy # {"Nuclei": ("Cells",), "Cells": (), ...}: direct parents write_mask_report(report, "/data/exp1_MorphoBricks_Preprocessing/mask_relationships.json") ``` -`preflight(..., ome_zarr_labels=True)` uses the label images of OME-Zarr images without a -masks folder; `read_mask` reads one mask's pixels; `pair_counts` and `piece_counts` give the -counts behind the report for one pair of masks or one mask. +`scan_masks(..., thresholds={"inside_object_share": 0.8})` changes a threshold (the others +keep their defaults, `DEFAULT_THRESHOLDS`). `preflight(..., ome_zarr_labels=True)` uses the +label images of OME-Zarr images without a masks folder; `read_mask` reads one mask's pixels; +`pair_counts` and `piece_counts` give the counts behind the report for one pair of masks or +one mask. Reading the table, group statistics and normalisation parameters come from the normalization brick (`morphobricks_normalization`, installed with this package): diff --git a/src/morphobricks_preprocessing/__init__.py b/src/morphobricks_preprocessing/__init__.py index 071e7514..fe518a33 100644 --- a/src/morphobricks_preprocessing/__init__.py +++ b/src/morphobricks_preprocessing/__init__.py @@ -67,10 +67,12 @@ write_manifest, ) from .mask_scan import ( + DEFAULT_THRESHOLDS, MASK_REPORT_FILENAME, MaskReport, MaskTypeSummary, PairRelation, + mask_thresholds, pair_counts, piece_counts, scan_masks, @@ -148,6 +150,7 @@ def __dir__() -> list[str]: __all__ = [ "CONSISTENCY_CODES", "DEFAULT_MEMORY_LIMIT", + "DEFAULT_THRESHOLDS", "DESIGN_COLUMNS", "DESIGN_FILENAME", "DESIGN_SCHEMA_VERSION", @@ -239,6 +242,7 @@ def __dir__() -> list[str]: "iter_xyc_planes", "load_design", "make_plane_id", + "mask_thresholds", "min_max", "ngff_label_names", "ome_zarr_relative_path", diff --git a/src/morphobricks_preprocessing/cli.py b/src/morphobricks_preprocessing/cli.py index c810d94a..49f33538 100644 --- a/src/morphobricks_preprocessing/cli.py +++ b/src/morphobricks_preprocessing/cli.py @@ -34,7 +34,13 @@ write_experiment_design, ) from .discovery import check_output_folder, default_output_folder -from .mask_scan import MASK_REPORT_FILENAME, scan_masks, write_mask_report +from .mask_scan import ( + DEFAULT_THRESHOLDS, + MASK_REPORT_FILENAME, + mask_thresholds, + scan_masks, + write_mask_report, +) from .planes import DEFAULT_MEMORY_LIMIT, format_bytes from .report import ( balance_summary, @@ -201,9 +207,31 @@ def build_parser() -> argparse.ArgumentParser: ) _common(masks) masks.add_argument("--log", metavar="FILE", help="write technical details to this file") + limits = masks.add_argument_group( + "relationship thresholds", "change the defaults (contract §13.1); shares are 0 to 1" + ) + for name, text in _THRESHOLD_HELP.items(): + limits.add_argument( + "--" + name.replace("_", "-"), + dest=name, + type=float, + metavar="X", + help=f"{text} (default {DEFAULT_THRESHOLDS[name]:g})", + ) return parser +_THRESHOLD_HELP = { + "inside_pixel_share": "share of a child's pixels in one parent for it to be inside", + "inside_object_share": "share of the children inside for a mask to be inside another", + "parent_larger_share": "share of the inside children whose parent is larger", + "same_objects_iou": "intersection over union for two objects to match one to one", + "same_objects_share": "share of the objects of both masks matched for the same objects", + "grouped_label_share": "share of a grouped mask's pixels carrying the parent's label", + "grouped_median_pieces": "median pieces per label above which a mask is grouped", +} + + def main(argv: Sequence[str] | None = None) -> int: """Run the command line; return the exit status.""" _safe_output() @@ -278,6 +306,12 @@ def _check(args: argparse.Namespace) -> int: def _masks(args: argparse.Namespace) -> int: + changes = {name: getattr(args, name) for name in _THRESHOLD_HELP} + try: + thresholds = mask_thresholds({k: v for k, v in changes.items() if v is not None}) + except ValueError as exc: + print(exc) + return 1 pre = preflight( args.experiment, args.output, @@ -301,7 +335,7 @@ def progress(done: int, total: int, path: str) -> None: print(f"[{done + 1}/{total}] {path}", flush=True) print() - report = scan_masks(pre.masks, pre.discovery.root, progress=progress) + report = scan_masks(pre.masks, pre.discovery.root, progress=progress, thresholds=thresholds) output = Path(args.output) if args.output else default_output_folder(pre.discovery.root) path = write_mask_report(report, output / MASK_REPORT_FILENAME) print() diff --git a/src/morphobricks_preprocessing/mask_scan.py b/src/morphobricks_preprocessing/mask_scan.py index 9196c19d..4434216f 100644 --- a/src/morphobricks_preprocessing/mask_scan.py +++ b/src/morphobricks_preprocessing/mask_scan.py @@ -26,6 +26,10 @@ nuclei). A type may have several direct parents (cells inside an organoid and inside a region of the slide). +The numbers above are the default thresholds (:data:`DEFAULT_THRESHOLDS`); each can be changed +(``scan_masks(..., thresholds=...)``, or the ``masks`` command's options), and the values used +are written in the report. + Nothing is decided for the user: the report is a suggestion with its evidence; the user names the masks and confirms or corrects the relationships. Pieces are counted with 4-connectivity (``scipy.ndimage.label``). @@ -56,7 +60,7 @@ MASK_REPORT_VERSION = "1" """Version of the report's format.""" -THRESHOLDS = MappingProxyType( +DEFAULT_THRESHOLDS = MappingProxyType( { "inside_pixel_share": 0.5, "inside_object_share": 0.9, @@ -64,10 +68,13 @@ "same_objects_iou": 0.5, "same_objects_share": 0.9, "grouped_label_share": 0.99, - "grouped_median_pieces": 2, + "grouped_median_pieces": 2.0, } ) -"""Fixed thresholds of the relationships (contract §13); written in every report.""" +"""Default thresholds of the kinds and relationships (contract §13.1); each may be changed +(:func:`mask_thresholds`), and the values used are written in every report.""" + +_SHARES = frozenset(DEFAULT_THRESHOLDS) - {"grouped_median_pieces"} _COUNTS = ( "children", @@ -147,11 +154,14 @@ class MaskReport: One summary per mask type, by name. pairs Every ordered pair of types seen together in at least one image, sorted. + thresholds + The thresholds the kinds and relations were decided with. """ images: int types: Mapping[str, MaskTypeSummary] = field(default_factory=lambda: MappingProxyType({})) pairs: tuple[PairRelation, ...] = () + thresholds: Mapping[str, float] = field(default_factory=lambda: DEFAULT_THRESHOLDS) @property def relationships(self) -> tuple[PairRelation, ...]: @@ -181,7 +191,7 @@ def to_dict(self) -> dict: """Return the report as JSON-compatible data (contract §13).""" return { "report_version": MASK_REPORT_VERSION, - "thresholds": dict(THRESHOLDS), + "thresholds": dict(self.thresholds), "images": self.images, "masks": { name: { @@ -200,6 +210,34 @@ def to_dict(self) -> dict: } +def mask_thresholds(changes: Mapping[str, float] | None = None) -> Mapping[str, float]: + """Return the default thresholds with some replaced, after checking them. + + Parameters + ---------- + changes + New values by name (:data:`DEFAULT_THRESHOLDS`); the others keep their default. + + Raises + ------ + ValueError + For an unknown name, a share or intersection over union outside (0, 1], or a median + number of pieces below 1. + """ + thresholds = dict(DEFAULT_THRESHOLDS) + for name, value in (changes or {}).items(): + if name not in thresholds: + known = ", ".join(DEFAULT_THRESHOLDS) + raise ValueError(f"Unknown mask threshold {name!r}; the thresholds are {known}.") + value = float(value) + if name in _SHARES and not 0 < value <= 1: + raise ValueError(f"The mask threshold {name} must be above 0 and at most 1.") + if name not in _SHARES and not value >= 1: + raise ValueError(f"The mask threshold {name} must be at least 1.") + thresholds[name] = value + return MappingProxyType(thresholds) + + def piece_counts(mask: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """Return the labels of a 2D mask and the number of 4-connected pieces of each. @@ -223,11 +261,15 @@ def piece_counts(mask: np.ndarray) -> tuple[np.ndarray, np.ndarray]: return values, counts -def pair_counts(child: np.ndarray, parent: np.ndarray) -> dict[str, int]: +def pair_counts( + child: np.ndarray, parent: np.ndarray, *, thresholds: Mapping[str, float] | None = None +) -> dict[str, int]: """Return the counts relating the objects of one mask (child) to another (parent). - Both are 2D label images of one position; label 0 is background. + Both are 2D label images of one position; label 0 is background. ``thresholds`` changes + the default ``inside_pixel_share`` and ``same_objects_iou`` (:func:`mask_thresholds`). """ + limits = mask_thresholds(thresholds) counts = dict.fromkeys(_COUNTS, 0) on = child != 0 child_labels, child_of = np.unique(child[on], return_inverse=True) @@ -252,13 +294,13 @@ def pair_counts(child: np.ndarray, parent: np.ndarray) -> dict[str, int]: counts["parents_with_0_children"] = slots return counts union = child_sizes[kid] + parent_sizes[host] - shared - matched = shared / union >= THRESHOLDS["same_objects_iou"] + matched = shared / union >= limits["same_objects_iou"] counts["children_matched_one_to_one"] = int(np.unique(kid[matched]).size) counts["parents_matched_one_to_one"] = int(np.unique(host[matched]).size) best = np.lexsort((-shared, kid)) lead = np.r_[True, kid[best][1:] != kid[best][:-1]] best_kid, best_host, best_shared = kid[best][lead], host[best][lead], shared[best][lead] - inside = best_shared / child_sizes[best_kid] >= THRESHOLDS["inside_pixel_share"] + inside = best_shared / child_sizes[best_kid] >= limits["inside_pixel_share"] counts["inside"] = int(inside.sum()) counts["inside_parent_larger"] = int( np.count_nonzero(parent_sizes[best_host[inside]] > child_sizes[best_kid[inside]]) @@ -278,6 +320,7 @@ def scan_masks( experiment_root: str | os.PathLike, *, progress: Callable[[int, int, str], None] | None = None, + thresholds: Mapping[str, float] | None = None, ) -> MaskReport: """Read every image's masks and describe the mask types and their relationships. @@ -289,11 +332,19 @@ def scan_masks( The experiment folder (for OME-Zarr label images). progress Called as ``progress(done, total, image)`` before each image is read. + thresholds + Thresholds to change from their defaults, by name (:func:`mask_thresholds`). Returns ------- MaskReport + + Raises + ------ + ValueError + For an invalid threshold, before any mask is read. """ + limits = mask_thresholds(thresholds) root = Path(experiment_root).resolve() images = inventory.images pieces: dict[str, list[np.ndarray]] = {name: [] for name in inventory.types} @@ -313,17 +364,18 @@ def scan_masks( for child in masks: for parent in masks: if child != parent: - found = pair_counts(masks[child], masks[parent]) + found = pair_counts(masks[child], masks[parent], thresholds=limits) total = totals.setdefault((child, parent), dict.fromkeys(_COUNTS, 0)) for key, value in found.items(): total[key] += value if progress is not None: progress(len(images), len(images), "") - types = {name: _summary(name, pieces[name], single[name]) for name in inventory.types} + types = {name: _summary(name, pieces[name], single[name], limits) for name in inventory.types} return MaskReport( images=len(images), types=MappingProxyType(types), - pairs=_relations(types, totals), + pairs=_relations(types, totals, limits), + thresholds=limits, ) @@ -343,7 +395,9 @@ def write_mask_report(report: MaskReport, path: str | os.PathLike) -> Path: return target -def _summary(name: str, pieces: list[np.ndarray], single: list[bool]) -> MaskTypeSummary: +def _summary( + name: str, pieces: list[np.ndarray], single: list[bool], limits: Mapping[str, float] +) -> MaskTypeSummary: """Return a mask type's summary from its pieces per label in every image.""" every = np.concatenate(pieces) if pieces else np.zeros(0, dtype=np.int64) median = float(np.median(every)) if every.size else 0.0 @@ -351,42 +405,46 @@ def _summary(name: str, pieces: list[np.ndarray], single: list[bool]) -> MaskTyp kind = "empty" elif all(single): kind = "one_value" - elif median > THRESHOLDS["grouped_median_pieces"]: + elif median > limits["grouped_median_pieces"]: kind = "grouped" else: kind = "instances" return MaskTypeSummary(name, kind, len(pieces), int(every.size), median) -def _relation(child: MaskTypeSummary, counts: Mapping[str, int]) -> str: +def _relation( + child: MaskTypeSummary, counts: Mapping[str, int], limits: Mapping[str, float] +) -> str: """Return the relation of a child type to a parent type from their counts.""" children, parents = counts["children"], counts["parents"] if not children or not parents: return "none" - share = THRESHOLDS["same_objects_share"] + share = limits["same_objects_share"] if ( counts["children_matched_one_to_one"] >= share * children and counts["parents_matched_one_to_one"] >= share * parents ): return "same_objects" - inside = counts["inside"] >= THRESHOLDS["inside_object_share"] * children + inside = counts["inside"] >= limits["inside_object_share"] * children if ( child.kind == "grouped" and inside and counts["child_pixels_with_parent_label"] - >= THRESHOLDS["grouped_label_share"] * counts["child_pixels"] + >= limits["grouped_label_share"] * counts["child_pixels"] ): return "grouped_by" - larger = counts["inside_parent_larger"] >= THRESHOLDS["parent_larger_share"] * counts["inside"] + larger = counts["inside_parent_larger"] >= limits["parent_larger_share"] * counts["inside"] return "inside" if inside and larger else "none" def _relations( - types: Mapping[str, MaskTypeSummary], totals: Mapping[tuple[str, str], Mapping[str, int]] + types: Mapping[str, MaskTypeSummary], + totals: Mapping[tuple[str, str], Mapping[str, int]], + limits: Mapping[str, float], ) -> tuple[PairRelation, ...]: """Return every pair's relation, marking the direct ones (module docstring).""" found = { - (child, parent): _relation(types[child], counts) + (child, parent): _relation(types[child], counts, limits) for (child, parent), counts in totals.items() } nested = {pair for pair, relation in found.items() if relation in ("inside", "grouped_by")} diff --git a/src/morphobricks_preprocessing/report.py b/src/morphobricks_preprocessing/report.py index f47286d3..3f9c754d 100644 --- a/src/morphobricks_preprocessing/report.py +++ b/src/morphobricks_preprocessing/report.py @@ -18,7 +18,7 @@ from .dataset import BalanceTable from .design import ExperimentDesign from .discovery import DiscoveryResult -from .mask_scan import MaskReport, PairRelation +from .mask_scan import DEFAULT_THRESHOLDS, MaskReport, PairRelation from .masks import MaskInventory from .run import Preflight, RunResult from .validation import Severity, ValidationIssue, plural @@ -380,6 +380,13 @@ def mask_report_summary(report: MaskReport) -> str: lines += [f" {first} and {second}" for first, second in report.same_objects] lines += ["", "Hierarchy (direct parents only):"] lines += _hierarchy_lines(report.hierarchy) + changed = [ + f"{name} {value:g}" + for name, value in report.thresholds.items() + if value != DEFAULT_THRESHOLDS[name] + ] + used = ", ".join(changed) + "; the others the defaults" if changed else "the defaults" + lines += ["", f"Thresholds: {used}"] lines += [ "", "These are suggestions from the masks' pixels: name each mask, then check the " diff --git a/tests/test_mask_scan.py b/tests/test_mask_scan.py index 177b08ae..a9cf078b 100644 --- a/tests/test_mask_scan.py +++ b/tests/test_mask_scan.py @@ -4,8 +4,11 @@ from collections import deque import numpy as np +import pytest from morphobricks_preprocessing import ( + DEFAULT_THRESHOLDS, + mask_thresholds, pair_counts, piece_counts, preflight, @@ -245,3 +248,56 @@ def test_progress_is_reported_per_image(tmp_path): preflight(root, masks_folder=folder).masks, root, progress=lambda *a: calls.append(a) ) assert calls == [(0, 2, "ctrl/a.tif"), (1, 2, "drug/a.tif"), (2, 2, "")] + + +# Thresholds --------------------------------------------------------------------------- + + +def test_thresholds_can_be_changed_and_are_checked(): + assert dict(mask_thresholds()) == dict(DEFAULT_THRESHOLDS) + changed = mask_thresholds({"inside_object_share": 0.75, "grouped_median_pieces": 3}) + assert changed["inside_object_share"] == 0.75 and changed["grouped_median_pieces"] == 3 + assert changed["same_objects_iou"] == DEFAULT_THRESHOLDS["same_objects_iou"] + for bad in ( + {"inside_share": 0.5}, + {"same_objects_iou": 0}, + {"inside_object_share": 1.2}, + {"same_objects_share": float("nan")}, + {"grouped_median_pieces": 0.5}, + ): + with pytest.raises(ValueError): + mask_thresholds(bad) + + +def test_changed_thresholds_change_the_relations_and_are_reported(tmp_path): + # Half of the dots lie inside the region: not inside at 90%, inside at 50%. + region = np.zeros(SHAPE, np.uint8) + region[:, :32] = 1 + dots = boxes({1: (5, 7, 5, 7), 2: (20, 22, 10, 12), 3: (5, 7, 40, 42), 4: (20, 22, 50, 52)}) + root, folder = write_experiment(tmp_path, {"Region": region, "Dots": dots}) + inventory = preflight(root, masks_folder=folder).masks + assert by_pair(scan_masks(inventory, root))[("Dots", "Region")].relation == "none" + report = scan_masks(inventory, root, thresholds={"inside_object_share": 0.5}) + assert by_pair(report)[("Dots", "Region")].relation == "inside" + assert report.to_dict()["thresholds"]["inside_object_share"] == 0.5 + # A stricter one-to-one match makes the two nuclei segmentations different objects. + pairs = by_pair(scan(tmp_path / "strict", hierarchy_masks())) + assert pairs[("Nuclei", "NucleiAlt")].relation == "same_objects" + root, folder = write_experiment(tmp_path / "iou", hierarchy_masks()) + strict = scan_masks( + preflight(root, masks_folder=folder).masks, root, thresholds={"same_objects_iou": 0.9} + ) + assert by_pair(strict)[("Nuclei", "NucleiAlt")].relation != "same_objects" + + +def test_invalid_thresholds_stop_before_any_mask_is_read(tmp_path): + root, folder = write_experiment(tmp_path, {"Cells": boxes(CELLS)}) + calls = [] + with pytest.raises(ValueError, match="at most 1"): + scan_masks( + preflight(root, masks_folder=folder).masks, + root, + progress=lambda *a: calls.append(a), + thresholds={"inside_pixel_share": 2}, + ) + assert calls == [] diff --git a/tests/test_masks.py b/tests/test_masks.py index 86b27db1..97430e42 100644 --- a/tests/test_masks.py +++ b/tests/test_masks.py @@ -241,3 +241,15 @@ def test_masks_command_stops_on_mask_errors_or_without_masks(experiment, masks, write_data(masks / "Cells" / "ctrl" / "rep1" / "a.tif", np.zeros((3, 3), dtype=np.uint8)) assert main(["masks", str(experiment), "--masks", str(masks)]) == 1 assert "do not have the size" in capsys.readouterr().out + + +def test_masks_command_takes_thresholds(experiment, masks, tmp_path, capsys): + out = tmp_path / "out" + args = ["masks", str(experiment), "--masks", str(masks), "--output", str(out)] + assert main([*args, "--same-objects-iou", "0.8", "--inside-object-share", "0.7"]) == 0 + printed = capsys.readouterr().out + assert "Thresholds: inside_object_share 0.7, same_objects_iou 0.8; the others" in printed + report = json.loads((out / "mask_relationships.json").read_text(encoding="utf-8")) + assert report["thresholds"]["same_objects_iou"] == 0.8 + assert main([*args, "--grouped-label-share", "1.5"]) == 1 + assert "must be above 0 and at most 1" in capsys.readouterr().out From 07a63bc369acbd5452d20b9fd686f5d126a55bcd Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 27 Oct 2026 17:20:45 +0000 Subject: [PATCH 4/5] Masks may also be named after their image: _.tif (owner) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit MorphoFeature's end-to-end trial on MorphoRelate's JUMP-Target-2 conversion found masks written as _.tif beside the images' folders. The masks folder may now mirror the experiment with the type after the image name, as well as one folder per type; both layouts may be mixed. A file that could name two images (one image name starting another's) or fit both layouts is the error ambiguous_mask, never guessed. On the four JUMP smoke fields both layouts give the same report. Contract §8 and §13, user guide §2 and §8, plan Q29, AGENTS. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01U6vJfhGWQ7pTpsJhw48Zf7 --- AGENTS.md | 5 +- docs/contract.md | 8 +- docs/plan.md | 8 +- docs/user_guide.md | 13 +++- src/morphobricks_preprocessing/masks.py | 97 ++++++++++++++++++------- tests/test_mask_scan.py | 16 ++++ tests/test_masks.py | 51 +++++++++++++ 7 files changed, 161 insertions(+), 37 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 819e2f1d..80eb3c90 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -91,7 +91,8 @@ src/morphobricks_preprocessing/ report.py plain-language summaries shared by CLI and notebook, group statistics tables cli.py, __main__ `morphobricks-preprocessing check|run|masks` (also `python -m morphobricks_preprocessing`) masks.py masks (Q29): find_masks pairs every image with its masks (masks folder with one - folder per type mirroring the experiment, OME-Zarr label images), header-only + folder per type mirroring the experiment, or _.tif in the + experiment's folders, mixed; OME-Zarr label images), header-only checks (2D, image YX, integer/bool), MaskInventory (to_manifest), read_mask mask_scan.py the optional mask step (Q29): piece_counts (scipy.ndimage), pair_counts, scan_masks -> MaskReport (kinds, relations, direct-parent hierarchy, evidence), @@ -415,7 +416,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | Dataset consistency (Q21, P2) | Step 1 "Dataset consistency" panel: balance table (images and technical-repeat folders per condition x biological repeat), `unbalanced_images` (< 0.5 or > 2 x median), `repeat_missing_conditions`, `different_pixel_sizes`, `missing_pixel_size` (OME-Zarr scale included), `different_plane_counts`; `different_sizes` narrowed to width/height; consistency issues shown only in that panel, not in Messages. Warnings only; outputs unchanged | | Experiment design (Q22, P3) | Step 1 after the labels: pairing dropdowns, one yellow "Yes, they are comparable" question per repeat group (single whole-experiment question without biological-repeat folders), control checkboxes; `experiment_design.csv` written after every answer; reruns load it and fold the questions behind a Change button; step 3 gated. CLI `run` reads it or `--design`, else writes a template and exits 1. `run_experiment(design=...)` optional; manifest `experiment_design_sha256` (optional key, version 1). Two folders of one condition in a group block that group | | Position identity (Q20, D02) | `position_id` = image name without extension, in every CSV row and manifest file entry; FOV key = labels + `position_id`; duplicates in one leaf folder are the error `duplicate_position_id`. CSV columns and schema 2 unchanged; manifest version 1 with an optional key | -| Masks (Q29) | Optional masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, paired by leaf folder + `position_id`; OME-Zarr label images are masks too. Header-only errors (`masks_folder_not_found`, `masks_inside_experiment`, `mask_names_differ_in_case`, `duplicate_mask`, `unreadable_mask`, `mask_not_2d`, `mask_size_mismatch`, `mask_pixel_type`) and warnings (`mask_missing`, `mask_without_image`, `mask_outside_type_folder`) in `check`/`run`; manifest key `masks` (optional, version 1). Optional mask step (`masks` command, `scan_masks`): kinds, relations with evidence, direct-parent hierarchy, `mask_relationships.json` (report version 1); default thresholds the user may change (`thresholds=` keyword, one CLI option each, checked by `mask_thresholds`, written in the report); suggestions only, the user names and confirms. SciPy core dependency. Confirmed by the owner: 2D masks at the image's YX, integer/bool pixels, masks folder outside the experiment. Next: notebook step and a file of the user's answers; 3D/time-lapse masks later; tag v0.2.0 proposed after merge | +| Masks (Q29) | Optional masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, or the experiment's folders with `_.tif` (both may be mixed; `ambiguous_mask` when a name could fit two images or both layouts), paired by leaf folder + `position_id`; OME-Zarr label images are masks too. Header-only errors (`masks_folder_not_found`, `masks_inside_experiment`, `mask_names_differ_in_case`, `duplicate_mask`, `ambiguous_mask`, `unreadable_mask`, `mask_not_2d`, `mask_size_mismatch`, `mask_pixel_type`) and warnings (`mask_missing`, `mask_without_image`, `mask_outside_type_folder`) in `check`/`run`; manifest key `masks` (optional, version 1). Optional mask step (`masks` command, `scan_masks`): kinds, relations with evidence, direct-parent hierarchy, `mask_relationships.json` (report version 1); default thresholds the user may change (`thresholds=` keyword, one CLI option each, checked by `mask_thresholds`, written in the report); suggestions only, the user names and confirms. SciPy core dependency. Confirmed by the owner: 2D masks at the image's YX, integer/bool pixels, masks folder outside the experiment. Next: notebook step and a file of the user's answers; 3D/time-lapse masks later; tag v0.2.0 proposed after merge | | Notebook style | Follow the lab's NucleiSky and LabChronicle notebooks (see §7) | | Widget widths | **99%**, never 100% | | Paths and filenames in widget HTML and notebook text | **`…` / `**…**`, never `` or inline code** | diff --git a/docs/contract.md b/docs/contract.md index f23b3ea3..1f5291d1 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -107,7 +107,7 @@ method is an offset and a scale, `(pixels − offset) / scale`. | Severity | Situations | Effect | |---|---|---| -| Error | Experiment design not confirmed (`design_not_confirmed`) or a repeat group holding two folders of one condition (`repeat_group_conflict`), when a design is used (§12); mixed dtypes; images with different numbers of channels (`different_channel_counts`, §2); unreadable TIFF (header or pixel data) or OME-Zarr folder in a leaf folder that is not excluded; OME-Zarr images, or converting TIFF images to OME-Zarr, without the `zarr` extra; OME-Zarr images whose declared axes do not fit the experiment's roles; every image excluded; more than one image series in a file; unsupported sample format or axis; missing/duplicate X or Y, duplicate roles, unassigned non-singleton axis, inconsistent mapping; plane that cannot be read within the memory limit, not even in blocks (for example a plane stored compressed in one piece); folder named `single` in another case, or beside other folders (§1); sibling folders differing only in case; images in one leaf folder with the same position identifier (§1); folder below the technical-repeat level; output folder inside the experiment; no usable planes; when masks are given, a masks folder missing or inside the experiment, mask names differing only in case, two masks of one type for an image, or a mask that is unreadable, not 2D, of another size than its image or not of integers or booleans (§13) | Processing is blocked | +| Error | Experiment design not confirmed (`design_not_confirmed`) or a repeat group holding two folders of one condition (`repeat_group_conflict`), when a design is used (§12); mixed dtypes; images with different numbers of channels (`different_channel_counts`, §2); unreadable TIFF (header or pixel data) or OME-Zarr folder in a leaf folder that is not excluded; OME-Zarr images, or converting TIFF images to OME-Zarr, without the `zarr` extra; OME-Zarr images whose declared axes do not fit the experiment's roles; every image excluded; more than one image series in a file; unsupported sample format or axis; missing/duplicate X or Y, duplicate roles, unassigned non-singleton axis, inconsistent mapping; plane that cannot be read within the memory limit, not even in blocks (for example a plane stored compressed in one piece); folder named `single` in another case, or beside other folders (§1); sibling folders differing only in case; images in one leaf folder with the same position identifier (§1); folder below the technical-repeat level; output folder inside the experiment; no usable planes; when masks are given, a masks folder missing or inside the experiment, mask names differing only in case, two masks of one type for an image, a mask whose name could belong to two images or types, or a mask that is unreadable, not 2D, of another size than its image or not of integers or booleans (§13) | Processing is blocked | | Warning | Files that are not images; empty folders; images outside leaf folders; images excluded by the user, or excluded paths that match no image; channel names differing between files; non-finite pixels; different image widths or heights; different numbers of depth slices, time points or other planes; physical pixel sizes that differ or are known for only some images; unbalanced numbers of images; biological repeats missing from some conditions (§11); planes without finite pixels while others remain; planes over the memory limit that are read in blocks; images without a mask of a type other images have, masks matching no analysed image, files outside the mask-type folders (§13) | Processing continues; the item is reported | Every message states what happened in plain language, lists the affected paths, and @@ -258,12 +258,12 @@ masks/ masks folder, outside the experiment folder | Rule | Detail | |---|---| -| Layout | One subfolder per mask type, each mirroring the experiment's folders. Any number of mask types. The label images of an OME-Zarr image (its OME-NGFF `labels` group, with the `zarr` extra) are masks too, of the type named by the label image; they are used whenever masks are asked for, with or without a masks folder. | +| Layout | Either one subfolder per mask type, each mirroring the experiment's folders (`//.../.tif`), or the experiment's folders with the type after the image name (`/.../_.tif`, as segmentation tools often write it, e.g. MorphoRelate's JUMP-Target-2 conversion); the two may be mixed. Any number of mask types. The label images of an OME-Zarr image (its OME-NGFF `labels` group, with the `zarr` extra) are masks too, of the type named by the label image; they are used whenever masks are asked for, with or without a masks folder. | | Pairing | A mask belongs to the image in the same leaf folder with the same `position_id` (§1), ignoring upper/lower case; its file is a TIFF (`.tif`, `.tiff`, optionally `.ome`, any case). Excluded images have no masks. | -| Names | A mask type is named by its folder or label image. Two names differing only in case are the error `mask_names_differ_in_case`; one image with two masks of one type (two extensions, or a file and a label image) is the error `duplicate_mask`. | +| Names | A mask type is named by its folder, by the suffix after the image name, or by its label image. A file that could belong to two images (an image name that starts another's, `a_b_cells.tif` beside images `a` and `a_b`) or fit both layouts is the error `ambiguous_mask`. Two names differing only in case are the error `mask_names_differ_in_case`; one image with two masks of one type (two extensions, or a file and a label image) is the error `duplicate_mask`. | | Header checks | From the headers only, no pixels read (`check`, `run`, `preflight`): a mask is 2D once its length-1 axes are removed (`mask_not_2d`), has the Y and X sizes of its image (`mask_size_mismatch`), and integer or boolean pixels (`mask_pixel_type`); a mask that cannot be opened, or holds several series, is `unreadable_mask`. All four are errors. | | Folder | The masks folder must exist (`masks_folder_not_found`) and be outside the experiment folder (`masks_inside_experiment`), where its folders would be read as labels. Both errors. | -| Warnings | `mask_missing` (images without a mask of a type that other images have: their objects of that type cannot be measured), `mask_without_image` (files matching no analysed image), `mask_outside_type_folder` (files directly in the masks folder). | +| Warnings | `mask_missing` (images without a mask of a type that other images have: their objects of that type cannot be measured), `mask_without_image` (files matching no analysed image), `mask_outside_type_folder` (files directly in the masks folder that are not TIFF masks). | | Run | `run_experiment(..., masks_folder=...)` stops on a mask error and records the masks in the manifest's `masks` key (§9). Without masks, the run and its outputs are unchanged. | ### 13.1 The mask step (`mask_relationships.json`, report version 1) diff --git a/docs/plan.md b/docs/plan.md index 3e230001..ba2e9dc3 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -725,4 +725,10 @@ Raised during review; to be decided before the milestones that need them: this is merged, so the other bricks can pin it (the follow-ups P4 and P5 then go to v0.3.0). Next, in their own changes: a notebook step for the masks (name each mask, say what it contains, confirm or correct the relationships) and a file saving those - answers for the other bricks, as `experiment_design.csv` does for the repeats.* + answers for the other bricks, as `experiment_design.csv` does for the repeats. + Also decided (owner, 27 October 2026, after MorphoFeature's end-to-end trial on + MorphoRelate's JUMP-Target-2 conversion, which writes `_.tif` beside the + images' folders): the masks folder may instead mirror the experiment with the type after + the image name, `/.../_.tif`, and the two layouts may be mixed. A + file that could name two images or fit both layouts is the error `ambiguous_mask`, never + guessed. On the four JUMP smoke fields, both layouts give the same report.* diff --git a/docs/user_guide.md b/docs/user_guide.md index 8b46c112..7d111e64 100644 --- a/docs/user_guide.md +++ b/docs/user_guide.md @@ -109,11 +109,17 @@ masks/ control/bio_rep1/tech_rep1/img_02.tif ``` +- Or keep the masks in the same folders as the experiment, each named after its image + followed by `_` and its kind (`img_01_cells.tif`, `img_01_nuclei.tif`), as many + segmentation tools write them; the kind is the part after the image name. Both ways may be + mixed in one masks folder. - A mask is a 2D TIFF label image with the width and height of its image: background 0, each object its own whole number (or a true/false mask). The extension may differ from the image's (`img_01.tiff` for `img_01.ome.tif`). -- You may have any number of kinds of mask. Their folder names are their names: they must - differ by more than upper/lower case. +- You may have any number of kinds of mask. Their folder names (or the part after the image + name) are their names: they must differ by more than upper/lower case. If one image's name + starts another's (`a.tif` and `a_b.tif`), a mask named `a_b_cells.tif` could belong to + either and is refused: rename the masks or the images. - The label images stored inside an OME-Zarr image (its `labels` group) are used as masks too, named after the label image. - An image without a mask of one kind is reported; its objects of that kind cannot be @@ -491,7 +497,8 @@ continues. | `mask_size_mismatch` | Error | A mask does not have the width and height of its image. Check that it was made from that image at full resolution. | | `mask_pixel_type` | Error | A mask's pixels are not whole numbers (for example decimals). Save it as an integer label image. | | `mask_missing` | Warning | Some images have no mask of a kind that other images have; their objects of that kind cannot be measured. Add the masks, or exclude the images. | -| `mask_without_image`, `mask_outside_type_folder` | Warning | Files in the masks folder that match no analysed image (check the folder and file names), or that are not in a folder of their kind. They are not used. | +| `ambiguous_mask` | Error | A mask's name could belong to two images (one image name starts another's) or fit both layouts. Rename the masks or the images so each names one image and one kind. | +| `mask_without_image`, `mask_outside_type_folder` | Warning | Files in the masks folder that match no analysed image (check the folder and file names), or other files directly in the masks folder. They are not used. | | `zarr_not_installed` | Error | OME-Zarr images, or converting TIFF images, need the `zarr` extra. Install it, exclude the OME-Zarr images, or read the TIFF images directly. | ## 9. Python interface diff --git a/src/morphobricks_preprocessing/masks.py b/src/morphobricks_preprocessing/masks.py index 75fc50e6..375b093c 100644 --- a/src/morphobricks_preprocessing/masks.py +++ b/src/morphobricks_preprocessing/masks.py @@ -6,11 +6,16 @@ /////.tif A mask belongs to the image with the same position (same leaf folder, same ``position_id``, -contract §1); its file may be a TIFF (``.tif``, ``.tiff``, any case, with an optional -``.ome``). The label images of an OME-Zarr image (its OME-NGFF ``labels`` group) are masks too, -of the type named by the label image, when the ``zarr`` extra is installed. +contract §1). The type may also follow the image name instead of being a folder:: -Any number of mask types is accepted; each is named by its folder (or label image), and two + ////_.tif + +and the two layouts may be mixed. Its file may be a TIFF (``.tif``, ``.tiff``, any case, with +an optional ``.ome``). The label images of an OME-Zarr image (its OME-NGFF ``labels`` group) +are masks too, of the type named by the label image, when the ``zarr`` extra is installed. + +Any number of mask types is accepted; each is named by its folder, by the suffix after the +image name, or by its label image, and two names may not differ only in case. Only headers are read here, never pixels: a mask must be 2D (length-1 axes aside) with the image's Y and X sizes, and its pixels integers or booleans. Reading the pixels and relating the masks to each other is :mod:`.mask_scan`. @@ -218,39 +223,77 @@ def _folder_issues(folder: Path, experiment_root: Path) -> list[ValidationIssue] def _folder_masks( folder: Path, discovery: DiscoveryResult, images: Mapping[str, TiffMetadata] ) -> tuple[list[MaskFile], list[ValidationIssue]]: - """Pair the TIFF files of every mask-type folder with the images.""" - wanted = {} + """Pair the TIFF files of the masks folder with the images, in either layout. + + A file is a mask of an image when its path is ``//`` (one + folder per mask type) or ``/_`` (the type after the + image name, as segmentation tools often write it). Each file is read the way that matches + an analysed image; a file matching two images or both layouts is refused, never guessed. + """ + by_key, by_leaf = {}, {} for leaf in discovery.leaves: + prefix = "" if leaf.relative_path == "." else leaf.relative_path + "/" for image in leaf.images: if image in images: - prefix = "" if leaf.relative_path == "." else leaf.relative_path + "/" - wanted[(prefix + position_id(image.rsplit("/", 1)[-1])).casefold()] = image - files, unpaired, loose = [], [], [] - for entry in sorted(folder.iterdir()): - if is_system_name(entry.name): + stem = position_id(image.rsplit("/", 1)[-1]) + by_key[(prefix + stem).casefold()] = image + by_leaf.setdefault(prefix.casefold(), []).append((stem, image)) + files, unpaired, loose, ambiguous = [], [], [], [] + for path in sorted(folder.rglob("*")): + relative = path.relative_to(folder) + if any(is_system_name(part) for part in relative.parts) or not path.is_file(): continue - if not entry.is_dir(): - loose.append(entry.name) + mask_path = relative.as_posix() + if not is_tiff_name(path.name): + (loose if len(relative.parts) == 1 else unpaired).append(mask_path) continue - for path in sorted(entry.rglob("*")): - relative = path.relative_to(entry) - if any(is_system_name(part) for part in relative.parts) or not path.is_file(): - continue - key = relative.with_name(position_id(path.name)).as_posix().casefold() - mask_path = path.relative_to(folder).as_posix() - if is_tiff_name(path.name) and key in wanted: - files.append(MaskFile(wanted[key], entry.name, mask_path, FOLDER)) - else: - unpaired.append(mask_path) + found = _suffix_matches(relative, by_leaf) + if len(relative.parts) > 1: + key = Path(*relative.parts[1:]).with_name(position_id(path.name)).as_posix() + if key.casefold() in by_key: + found.append((by_key[key.casefold()], relative.parts[0])) + if len(found) == 1: + files.append(MaskFile(found[0][0], found[0][1], mask_path, FOLDER)) + else: + (ambiguous if found else unpaired).append(mask_path) + issues = _layout_issues(loose, unpaired) + if ambiguous: + issues.append( + error( + "ambiguous_mask", + f"{plural(len(ambiguous), 'mask')} could belong to more than one image or mask " + "type (an image name that starts another image's name, or a path that fits both " + "layouts), so they cannot be paired.", + ambiguous, + "Rename the masks (or the images) so that each mask names one image and one type.", + ) + ) + return files, issues + + +def _suffix_matches(relative: Path, by_leaf: Mapping[str, list]) -> list[tuple[str, str]]: + """Return ``(image, mask type)`` for every image whose name, then ``_``, starts the file's.""" + prefix = "" if len(relative.parts) == 1 else relative.parent.as_posix() + "/" + stem = position_id(relative.name) + found = [] + for image_stem, image in by_leaf.get(prefix.casefold(), []): + head = image_stem + "_" + if len(stem) > len(head) and stem.casefold().startswith(head.casefold()): + found.append((image, stem[len(head) :])) + return found + + +def _layout_issues(loose: list[str], unpaired: list[str]) -> list[ValidationIssue]: issues = [] if loose: issues.append( warning( "mask_outside_type_folder", - f"{plural(len(loose), 'file')} in the masks folder {_are(loose)} not in a " - "mask-type folder and will not be used.", + f"{plural(len(loose), 'file')} in the masks folder {_are(loose)} not masks " + "and will not be used.", loose, - "Put each mask in the folder of its type, e.g. Nuclei//.tif.", + "Keep only masks there: Nuclei//.../.tif, or " + "/.../_Nuclei.tif.", ) ) if unpaired: @@ -264,7 +307,7 @@ def _folder_masks( "masks of excluded images.", ) ) - return files, issues + return issues def _label_masks(root: Path, images: Mapping[str, TiffMetadata]) -> list[MaskFile]: diff --git a/tests/test_mask_scan.py b/tests/test_mask_scan.py index a9cf078b..85253550 100644 --- a/tests/test_mask_scan.py +++ b/tests/test_mask_scan.py @@ -301,3 +301,19 @@ def test_invalid_thresholds_stop_before_any_mask_is_read(tmp_path): thresholds={"inside_pixel_share": 2}, ) assert calls == [] + + +def test_both_layouts_give_the_same_report(tmp_path): + masks = hierarchy_masks() + by_folder = scan(tmp_path / "folders", masks) + root = tmp_path / "suffix" / "exp" + folder = tmp_path / "suffix" / "masks" + for rel in ("ctrl/a.tif", "drug/a.tif"): + write_data(root / rel, np.zeros(SHAPE, dtype=np.uint8)) + for name, mask in masks.items(): + write_data(folder / rel.replace("a.tif", f"a_{name}.tif"), mask) + inventory = preflight(root, masks_folder=folder).masks + assert inventory.ok, inventory.issues + by_suffix = scan_masks(inventory, root) + assert by_suffix.to_dict()["pairs"] == by_folder.to_dict()["pairs"] + assert dict(by_suffix.hierarchy) == dict(by_folder.hierarchy) diff --git a/tests/test_masks.py b/tests/test_masks.py index 97430e42..005c8af1 100644 --- a/tests/test_masks.py +++ b/tests/test_masks.py @@ -111,6 +111,57 @@ def test_masks_of_excluded_images_are_not_used(experiment, masks): assert codes(found.issues) == ["mask_without_image"] +# The type after the image name ------------------------------------------------------- + + +def test_masks_named_image_then_type_are_paired(experiment, tmp_path, no_pixel_reads): + folder = tmp_path / "masks" + for seed, rel in enumerate(IMAGES): + leaf, name = rel.rsplit("/", 1) + stem = name.rsplit(".", 1)[0] + write_data(folder / leaf / f"{stem}_cells.tif", label_image(seed)) + write_data(folder / leaf / f"{stem}_nuclei.tif", label_image(seed)) + found = inventory(experiment, folder) + assert found.ok and found.issues == () + assert found.types == ("cells", "nuclei") + by_type = found.for_image("ctrl/rep1/b.tif") + assert by_type["nuclei"].path == "ctrl/rep1/b_nuclei.tif" + + +def test_the_two_layouts_may_be_mixed(experiment, tmp_path): + folder = tmp_path / "masks" + for seed, rel in enumerate(IMAGES): + write_data(folder / "Cells" / rel, label_image(seed)) + leaf, name = rel.rsplit("/", 1) + write_data(folder / leaf / f"{name.rsplit('.', 1)[0]}_Nuclei.tif", label_image(seed)) + found = inventory(experiment, folder) + assert found.ok and found.types == ("Cells", "Nuclei") + assert len(found.files) == 2 * len(IMAGES) + + +def test_images_in_the_experiment_folder_itself(tmp_path): + root = tmp_path / "exp" + write_data(root / "img_01.tif", np.zeros((16, 20), dtype=np.uint8)) + write_data(tmp_path / "masks" / "img_01_cells.tiff", label_image(0)) + write_data(tmp_path / "masks" / "img_02_cells.tif", label_image(0)) + found = inventory(root, tmp_path / "masks") + assert [(f.image, f.mask_type) for f in found.files] == [("img_01.tif", "cells")] + assert codes(found.issues) == ["mask_without_image"] + + +def test_a_mask_that_could_name_two_images_is_refused(tmp_path): + root = tmp_path / "exp" + for name in ("a.tif", "a_b.tif"): + write_data(root / "ctrl" / name, np.zeros((16, 20), dtype=np.uint8)) + write_data(tmp_path / "masks" / "ctrl" / "a_b_cells.tif", label_image(0)) # a + b_cells? + write_data(tmp_path / "masks" / "ctrl" / "a_nuclei.tif", label_image(0)) # a only + found = inventory(root, tmp_path / "masks") + ambiguous = next(issue for issue in found.issues if issue.code == "ambiguous_mask") + assert ambiguous.paths == ("ctrl/a_b_cells.tif",) + assert not found.ok + assert [(f.image, f.mask_type) for f in found.files] == [("ctrl/a.tif", "nuclei")] + + # Errors ------------------------------------------------------------------------------- From 127cd1e43b6cf972c84c805ec97f27afb88316cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Hidalgo?= <47150974+IvanHCenalmor@users.noreply.github.com> Date: Tue, 27 Oct 2026 23:05:06 +0200 Subject: [PATCH 5/5] Renumber the masks decision to plan question 32 PR #27 (folder step without the experiment design) takes question 31, and goes in first. Co-Authored-By: Claude Opus 5.5 --- AGENTS.md | 10 +++++----- docs/plan.md | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 95a12156..ca299568 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -96,11 +96,11 @@ src/morphobricks_preprocessing/ run.py preflight(), run_experiment() (incremental), sample_planes() for previews report.py plain-language summaries shared by CLI and notebook, group statistics tables cli.py, __main__ `morphobricks-preprocessing check|run|masks` (also `python -m morphobricks_preprocessing`) - masks.py masks (Q31): find_masks pairs every image with its masks (masks folder with one + masks.py masks (Q32): find_masks pairs every image with its masks (masks folder with one folder per type mirroring the experiment, or _.tif in the experiment's folders, mixed; OME-Zarr label images), header-only checks (2D, image YX, integer/bool), MaskInventory (to_manifest), read_mask - mask_scan.py the optional mask step (Q31): piece_counts (scipy.ndimage), pair_counts, + mask_scan.py the optional mask step (Q32): piece_counts (scipy.ndimage), pair_counts, scan_masks -> MaskReport (kinds, relations, direct-parent hierarchy, evidence), write_mask_report (mask_relationships.json); DEFAULT_THRESHOLDS, mask_thresholds (checked changes) @@ -429,7 +429,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | Topic | Decision | |---|---| | Binning | Per-plane min–max, same bin count per experiment, pooled counts for groups (binning study) | -| Core dependencies | `numpy`, `tifffile`, `scipy` (the mask step's piece counts, Q31) and `morphobricks-normalization` (git dependency pinned to its release tag v0.2.0; `[tool.hatch.metadata] allow-direct-references`). `ipywidgets` and `matplotlib` in the `notebook` extra; `zarr` in the `zarr` extra; `pytest` and `ruff` in `dev` | +| Core dependencies | `numpy`, `tifffile`, `scipy` (the mask step's piece counts, Q32) and `morphobricks-normalization` (git dependency pinned to its release tag v0.2.0; `[tool.hatch.metadata] allow-direct-references`). `ipywidgets` and `matplotlib` in the `notebook` extra; `zarr` in the `zarr` extra; `pytest` and `ruff` in `dev` | | OME-Zarr (Q12) | Input may be TIFF and/or OME-Zarr (folders ending in `.zarr`), mixed. OME-Zarr axes come from their metadata and are matched by role; level 0 of a single image only, several images per store is `multiple_series`. TIFF choice once per experiment (`convert_tiffs`): read directly (default) or convert to `/ome_zarr` (OME-NGFF 0.4) and compute rows from the copies; switching recomputes TIFF rows and removes unused copies. Reads OME-NGFF 0.4 and 0.5 (0.5 needs zarr 3). In the notebook: its own step-3 card "(Optional) Image conversion" after Settings (Q13): radio buttons (read directly = recommended), a disk-space warning (amber, red when the estimate exceeds free space, never blocking), and conversion settings (chunk size, compression + level, pyramid levels) shown only when converting. CLI `--convert-tiff` with `--chunk-size`, `--compression`, `--compression-level`, `--pyramid-levels`. No overlap setting: chunks never overlap | | System files (Q5) | Ignored silently; `single` check and sibling uniqueness are case-insensitive | | Too-deep folders (Q7) | Error | @@ -451,7 +451,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | Image type (Q29, P4) | Step 2 card after the axis roles: Fluorescence or RGB, suggestion preselected, confirmed with **Yes, this is the right image type**; RGB suggested for TIFF photometric RGB with three samples on the C axis; RGB needs three integer channels (`check_rgb`), else disabled with the reason and `rgb_not_possible` in a run. Manifest keys `image_type`, `channels_from_rgb_samples` (optional, version 1); type changes recompute nothing. CLI `--image-type`, default the printed suggestion | | Intensity checks (Q30, P5) | Warnings from the table alone with the brick's `image_summaries`/`intensity_checks`, shown automatically after step 3 and in step 4 for any table; limits under a collapsed **Check settings** panel, recorded in the manifest (`intensity_checks`); per-channel box plots with flagged images and cumulative distributions per biological repeat, notebook only; the CLI prints the warnings (`--outlier-limit`, `--saturation-limit`, `--min-images`) | | Position identity (Q20, D02) | `position_id` = image name without extension, in every CSV row and manifest file entry; FOV key = labels + `position_id`; duplicates in one leaf folder are the error `duplicate_position_id`. CSV columns and schema 2 unchanged; manifest version 1 with an optional key | -| Masks (Q31) | Optional masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, or the experiment's folders with `_.tif` (both may be mixed; `ambiguous_mask` when a name could fit two images or both layouts), paired by leaf folder + `position_id`; OME-Zarr label images are masks too. Header-only errors (`masks_folder_not_found`, `masks_inside_experiment`, `mask_names_differ_in_case`, `duplicate_mask`, `ambiguous_mask`, `unreadable_mask`, `mask_not_2d`, `mask_size_mismatch`, `mask_pixel_type`) and warnings (`mask_missing`, `mask_without_image`, `mask_outside_type_folder`) in `check`/`run`; manifest key `masks` (optional, version 1). Optional mask step (`masks` command, `scan_masks`): kinds, relations with evidence, direct-parent hierarchy, `mask_relationships.json` (report version 1); default thresholds the user may change (`thresholds=` keyword, one CLI option each, checked by `mask_thresholds`, written in the report); suggestions only, the user names and confirms. SciPy core dependency. Confirmed by the owner: 2D masks at the image's YX, integer/bool pixels, masks folder outside the experiment. Next: notebook step and a file of the user's answers; 3D/time-lapse masks later; tag v0.2.0 proposed after merge | +| Masks (Q32) | Optional masks folder outside the experiment, one subfolder per mask type mirroring the experiment's folders, or the experiment's folders with `_.tif` (both may be mixed; `ambiguous_mask` when a name could fit two images or both layouts), paired by leaf folder + `position_id`; OME-Zarr label images are masks too. Header-only errors (`masks_folder_not_found`, `masks_inside_experiment`, `mask_names_differ_in_case`, `duplicate_mask`, `ambiguous_mask`, `unreadable_mask`, `mask_not_2d`, `mask_size_mismatch`, `mask_pixel_type`) and warnings (`mask_missing`, `mask_without_image`, `mask_outside_type_folder`) in `check`/`run`; manifest key `masks` (optional, version 1). Optional mask step (`masks` command, `scan_masks`): kinds, relations with evidence, direct-parent hierarchy, `mask_relationships.json` (report version 1); default thresholds the user may change (`thresholds=` keyword, one CLI option each, checked by `mask_thresholds`, written in the report); suggestions only, the user names and confirms. SciPy core dependency. Confirmed by the owner: 2D masks at the image's YX, integer/bool pixels, masks folder outside the experiment. Next: notebook step and a file of the user's answers; 3D/time-lapse masks later; tag v0.2.0 proposed after merge | | Notebook style | Follow the lab's NucleiSky and LabChronicle notebooks (see §7) | | Widget widths | **99%**, never 100% | | Paths and filenames in widget HTML and notebook text | **`…` / `**…**`, never `` or inline code** | @@ -742,7 +742,7 @@ jupyter lab notebooks/MorphoBricks_Preprocessing.ipynb P5. - **Since v0.1.0:** follow-ups P4 (the image type, question 29) and P5 (intensity checks and plots, question 30) of the normalization plan. The owner tags v0.2.0 next. -- **Masks (question 31, in review):** pairing every image with its masks, header checks, +- **Masks (question 32, in review):** pairing every image with its masks, header checks, the manifest's `masks` key, and the mask step with `mask_relationships.json`; library, command line and documentation (the notebook step later). - **Not yet verified:** diff --git a/docs/plan.md b/docs/plan.md index 7e1a4c50..26dd2884 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -748,7 +748,7 @@ Raised during review; to be decided before the milestones that need them: kind with the action once per kind, takes `--outlier-limit`, `--saturation-limit` and `--min-images`, and draws no plots. Comparisons stay within a condition, as the brick does. `plane_histograms.csv` and `experiment_design.csv` do not change.* -31. Can this brick find each image's masks and tell how they relate, so that a user who is +32. Can this brick find each image's masks and tell how they relate, so that a user who is not a programmer only has to name the masks (MorphoFeature decision D78)? Bricks that measure objects (MorphoObjects, MorphoRelate, CellPatcher, MorphoFeature's workflows) need every image paired with its masks and the parent-child relationships between mask