diff --git a/AGENTS.md b/AGENTS.md index a740bfd5..9d02a324 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -459,7 +459,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | 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`) | | Histogram range (Q34, P6) | Step 3 **Settings**: the brick's `HistogramSettings` (`minmax` default, or `robust` with adjustable quantiles); recorded in the manifest's `settings.histogram`, any change recomputes every row. Every plane binned with `histogram_with_range_check` (memory estimate with `range_check=True`, through `planes._Binning`); `narrow_histogram_range` warnings from `range_check_issues`, one per channel, with the narrow planes recorded in the manifest's `narrow_histogram_ranges` so kept rows still report them. CSV schema 4 (underflow and overflow columns, zero with `minmax`); schema 3 tables recomputed. CLI `--histogram-range`, `--robust-quantiles`. Normalization brick pinned to v0.3.1 | | 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 (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 | +| 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; tag v0.2.0 proposed after merge. Plane identity (issue #36): masks may have Z and/or T named in their own metadata (never guessed: unnamed axes are `mask_not_2d`); the masks of an image share their planes (`mask_planes_differ`) and `run` checks them against the image's confirmed roles (`mask_planes_mismatch`, `mask_plane_issues`), so a 2D mask is refused for a Z-stack or time-lapse; manifest mask entries gain `axes` and `shape` (optional keys, version 1); `read_mask(z=, t=)` reads one plane; `scan_masks` scans plane by plane | | 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/README.md b/README.md index b10e952c..403e212f 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,10 @@ Pre-alpha. Implemented and tested: 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`) + (`mask_scan`, command `masks`). A mask is 2D, or has depth (Z) and/or time (T) planes + named in its own metadata, one per image plane; a 2D mask of a Z-stack or time-lapse is + refused, projections included (`mask_planes_mismatch`), and the manifest records each + mask's `axes` and `shape` - A step-by-step notebook ([notebooks/MorphoBricks_Preprocessing.ipynb](notebooks/MorphoBricks_Preprocessing.ipynb)) and a command line (`morphobricks-preprocessing check` / `run` / `masks`) diff --git a/docs/contract.md b/docs/contract.md index 5f9b6630..a461f70a 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -118,7 +118,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); image type RGB for images that are not three channels of an integer type (`rgb_not_possible`, §13); mixed dtypes; images with different numbers of channels (`different_channel_counts`, §2) or different channel names at the same channel number (`channel_names_differ`, §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 (§15) | 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); image type RGB for images that are not three channels of an integer type (`rgb_not_possible`, §13); mixed dtypes; images with different numbers of channels (`different_channel_counts`, §2) or different channel names at the same channel number (`channel_names_differ`, §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, has axes other than Y, X and named Z/T (`mask_not_2d`), is of another size than its image, is not of integers or booleans, has other depth slices or time points than the image's other masks or, in `run`, than the image (§15) | Processing is blocked | | Warning | Intensity checks of the written table (`image_intensity_outlier`, `repeat_intensity_shift`, `saturated_image`, `no_intensity_spread`, §14), reported apart from the run's issues; files that are not images; empty folders; images outside leaf folders; images excluded by the user, or excluded paths that match no image; 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; channels whose central intensity range covers too few histogram bins (`narrow_histogram_range`, §4); 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 (§15) | Processing continues; the item is reported | Every message states what happened in plain language, lists the affected paths, and @@ -161,7 +161,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), `channel_names_trimmed` (true when leading or trailing whitespace was removed from the image's channel names as they were read, §2; `channel_names` holds the trimmed names; false when absent). 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 (§15), 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. | +| `masks` | The masks paired with the images (§15), 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`, `axes` and `shape`: the mask's axes and lengths without its length-1 axes, `YX` for a 2D mask, with `Z` and/or `T` for a mask with depth slices or time points, e.g. `TYX` and `[50, 512, 512]`; §15). Absent when the run had no masks. Optional key: the manifest stays version 1, and the table does not depend on the masks. `axes` and `shape` are absent in manifests written before them, whose masks are all 2D; a later brick must read `axes` before using a mask and never apply a 2D mask to several planes. | | `image_type` | `fluorescence` or `rgb`, the image type the user confirmed (§13); null when the run was given none, or absent in manifests written before it. Optional key: the manifest stays version 1. Changing it recomputes nothing. | | `channels_from_rgb_samples` | `true` when the C role is the RGB samples axis of images stored as RGB colour (§13), `false` otherwise; absent in manifests written before it. Optional key. | | `intensity_checks` | The limits of the intensity checks run on the table (§14): `IntensityCheckSettings.to_dict()` of the normalization brick (`library_version`, `outlier_limit`, `saturation_limit`, `min_images`); absent in manifests written before it. Optional key. Changing the limits recomputes nothing. | @@ -311,7 +311,8 @@ masks/ masks folder, outside the experiment 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, 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. | +| Header checks | From the headers only, no pixels read (`check`, `run`, `preflight`): once its length-1 axes are removed, a mask has Y and X and at most a depth axis Z and a time axis T, each named so in its own metadata (TIFF axes as `tifffile` reads them, e.g. an ImageJ hyperstack or OME-TIFF; OME-NGFF axes for a label image); any other or unnamed axis, such as a plain stack that could be depth or time, is `mask_not_2d` and is never guessed. A mask has the Y and X sizes of its image (`mask_size_mismatch`) and integer or boolean pixels (`mask_pixel_type`); the masks of one image have the same depth slices and time points (`mask_planes_differ`); a mask that cannot be opened, or holds several series, is `unreadable_mask`. All five are errors. | +| Planes | A mask belongs to the image's planes (issue #36): its Z and T lengths (1 for an axis it lacks) must equal the image's under the confirmed axis roles, checked by `run` (`mask_planes_mismatch`, error; `mask_plane_issues`). So a 2D mask belongs only to an image with one depth slice and one time point, and is refused for a Z-stack or a time-lapse, even as a projection (a decision for the owner if wanted later); a mask's plane at index `i` along Z or T belongs to the image's rows with `z = i` or `t = i`. `read_mask(..., z=, t=)` reads one plane and requires `z` and `t` for a mask with those axes. `check` and the `masks` command have no confirmed roles and do not compare with the image. | | 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 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. | @@ -330,7 +331,12 @@ the user names the masks and confirms or corrects the relationships. | `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. +A mask with depth slices or time points is scanned one plane at a time, each plane like an +image of its own: objects are never matched across planes, `objects` counts a label once per +plane it is in, and `images` still counts images. + +For a child type and a parent type, the counts are summed over the images (and planes) with +both masks. Each child object's parent is the parent object sharing most of its pixels. | Relation of child to parent | Rule (default thresholds; the values used are written in the report) | diff --git a/docs/plan.md b/docs/plan.md index 3a18267f..36e76a0d 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -841,6 +841,15 @@ Raised during review; to be decided before the milestones that need them: 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.* + *Follow-up (October 2026, data-flow audit, issue #36): masks get a plane identity. A + mask may have a depth axis Z and/or a time axis T, named in its own metadata (never + guessed: an unnamed stack axis is `mask_not_2d`); the masks of one image share their + planes (`mask_planes_differ`), and `run` checks them against the image's confirmed + roles (`mask_planes_mismatch`). A 2D mask is therefore refused for a Z-stack or a + time-lapse; whether a 2D mask may stand for a projection of such an image is left to + the owner, so it is refused for now. The manifest's mask entries gain `axes` and + `shape` (optional keys, manifest version 1); `read_mask(z=, t=)` reads one plane and + `scan_masks` scans plane by plane.* 33. What does this brick's second tag, v0.2.0, hold? *Decided (27 October 2026): v0.2.0 holds everything merged after v0.1.0: the image type (follow-up P4, question 29), the intensity checks and plots (follow-up P5, question diff --git a/docs/user_guide.md b/docs/user_guide.md index a4849e65..4bb6ef55 100644 --- a/docs/user_guide.md +++ b/docs/user_guide.md @@ -113,9 +113,13 @@ masks/ 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, +- A mask is a 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`). +- A 2D mask belongs to an image with one plane. For a Z-stack or a time-lapse, the mask needs + one plane per depth slice or time point, with its axes named in its metadata (an ImageJ + hyperstack or OME-TIFF with axes `ZYX` or `TYX`); a 2D mask for such an image, or a plain + stack whose axis could be depth or time, is refused by `run`/`check`. - 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 @@ -478,7 +482,8 @@ The default output folder is next to the experiment folder: relationship between every two kinds with its counts, and the hierarchy. The rules and thresholds are in [contract.md](contract.md) §15. -When masks are given, `run_manifest.json` also lists every image's masks (`masks`). +When masks are given, `run_manifest.json` also lists every image's masks (`masks`), with +each mask's axes and size (for example `TYX` and 50 x 512 x 512 for a mask per time point). 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 @@ -565,7 +570,9 @@ continues. | `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_not_2d` | Error | A mask has dimensions other than Y and X that its metadata does not name as depth (Z) or time (T), for example a plain stack. Save it as a 2D label image, or with its axes named (an ImageJ hyperstack or OME-TIFF with axes TYX or ZYX). | +| `mask_planes_differ` | Error | The masks of one image have different numbers of depth slices or time points. Make every mask of an image with the image's planes. | +| `mask_planes_mismatch` | Error (`run`) | A mask does not have the depth slices and time points of its image under the confirmed axis roles: a 2D mask fits only an image with one plane (not a Z-stack or time-lapse), and a mask with Z or T needs as many slices or time points as its image. Save one mask plane per image plane, with the axes named, or check the axis roles. | | `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. | diff --git a/src/morphobricks_preprocessing/__init__.py b/src/morphobricks_preprocessing/__init__.py index 41e52fe4..621091a9 100644 --- a/src/morphobricks_preprocessing/__init__.py +++ b/src/morphobricks_preprocessing/__init__.py @@ -96,7 +96,7 @@ scan_masks, write_mask_report, ) -from .masks import MaskFile, MaskInventory, find_masks, read_mask +from .masks import MaskFile, MaskInventory, find_masks, mask_plane_issues, read_mask from .metadata import ( ImageMetadata, MetadataReport, @@ -271,6 +271,7 @@ def __dir__() -> list[str]: "iter_xyc_planes", "load_design", "make_plane_id", + "mask_plane_issues", "mask_thresholds", "min_max", "ngff_label_names", diff --git a/src/morphobricks_preprocessing/mask_scan.py b/src/morphobricks_preprocessing/mask_scan.py index 3a64d55f..1a286839 100644 --- a/src/morphobricks_preprocessing/mask_scan.py +++ b/src/morphobricks_preprocessing/mask_scan.py @@ -33,6 +33,10 @@ 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``). + +Masks with depth slices or time points (``Z``/``T`` in :attr:`.masks.MaskFile.axes`) are +scanned one plane at a time, each plane counted like an image of its own: objects are never +matched across planes, and the masks of one image must have the same planes. """ from __future__ import annotations @@ -107,7 +111,8 @@ class MaskTypeSummary: images Images with a mask of this type. objects - Distinct labels, summed over the images. + Distinct labels, summed over the images (and over the planes of masks with depth + slices or time points: a label present in 5 time points counts 5 times). median_pieces_per_label Median number of 4-connected pieces per label over all objects; 0 when empty. """ @@ -342,35 +347,51 @@ def scan_masks( Raises ------ ValueError - For an invalid threshold, before any mask is read. + For an invalid threshold, before any mask is read; for an image whose masks have + different depth slices or time points; for a mask that cannot be read as the planes + its header gave (:func:`.masks.read_mask`). """ limits = mask_thresholds(thresholds) 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} + seen: dict[str, int] = dict.fromkeys(inventory.types, 0) 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], thresholds=limits) - total = totals.setdefault((child, parent), dict.fromkeys(_COUNTS, 0)) - for key, value in found.items(): - total[key] += value + found_masks = inventory.for_image(image) + planes = {tuple(mask.planes.items()) for mask in found_masks.values()} + if len(planes) != 1: + raise ValueError( + f"the masks of {image} have different depth slices or time points " + "(mask_planes_differ)" + ) + counts_of = dict(next(iter(planes))) + for name in found_masks: + seen[name] += 1 + # A mask with depth slices or time points is scanned one plane at a time, as if + # each plane were an image: objects are never matched across planes. + for z in range(counts_of["Z"]): + for t in range(counts_of["T"]): + masks = { + name: read_mask( + mask, + inventory.folder, + root, + z=z if "Z" in mask.axes else None, + t=t if "T" in mask.axes else None, + ) + for name, mask in found_masks.items() + } + _scan_plane(masks, pieces, single, totals, limits) if progress is not None: progress(len(images), len(images), "") - types = {name: _summary(name, pieces[name], single[name], limits) for name in inventory.types} + types = { + name: _summary(name, seen[name], pieces[name], single[name], limits) + for name in inventory.types + } return MaskReport( images=len(images), types=MappingProxyType(types), @@ -395,10 +416,35 @@ def write_mask_report(report: MaskReport, path: str | os.PathLike) -> Path: return target +def _scan_plane( + masks: Mapping[str, np.ndarray], + pieces: dict[str, list[np.ndarray]], + single: dict[str, list[bool]], + totals: dict[tuple[str, str], dict[str, int]], + limits: Mapping[str, float], +) -> None: + """Add the pieces of every mask of one plane, and the counts of every pair, to the totals.""" + 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], thresholds=limits) + total = totals.setdefault((child, parent), dict.fromkeys(_COUNTS, 0)) + for key, value in found.items(): + total[key] += value + + def _summary( - name: str, pieces: list[np.ndarray], single: list[bool], limits: Mapping[str, float] + name: str, + images: int, + 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.""" + """Return a mask type's summary from its pieces per label in every image plane.""" 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: @@ -409,7 +455,7 @@ def _summary( kind = "grouped" else: kind = "instances" - return MaskTypeSummary(name, kind, len(pieces), int(every.size), median) + return MaskTypeSummary(name, kind, images, int(every.size), median) def _relation( diff --git a/src/morphobricks_preprocessing/masks.py b/src/morphobricks_preprocessing/masks.py index 2b45abb7..c87be4fb 100644 --- a/src/morphobricks_preprocessing/masks.py +++ b/src/morphobricks_preprocessing/masks.py @@ -16,14 +16,19 @@ 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`. +names may not differ only in case. Only headers are read here, never pixels: a mask has the +image's Y and X sizes and integer or boolean pixels. Its other dimensions (length-1 axes +aside) must be labelled depth (Z) or time (T) in its own metadata, never guessed, and every +mask of one image has the same depth slices and time points. Whether those are the image's +planes needs the confirmed axis roles: :func:`mask_plane_issues`, run by +:func:`.run.run_experiment` (issue #36). Reading the pixels and relating the masks to each +other is :mod:`.mask_scan`. """ from __future__ import annotations import logging +import math import os from collections.abc import Iterable, Mapping from dataclasses import dataclass @@ -33,6 +38,7 @@ import numpy as np import tifffile +from .axes import AxisMapping 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 @@ -46,6 +52,9 @@ OME_ZARR = "ome-zarr" """Source of a mask kept as a label image inside an OME-Zarr image.""" +PLANE_AXES = ("Z", "T") +"""The dimensions a mask may have besides Y and X: depth slices and time points.""" + @dataclass(frozen=True) class MaskFile: @@ -63,9 +72,13 @@ class MaskFile: source :data:`FOLDER` or :data:`OME_ZARR`. shape - Shape of the mask without its length-1 axes, ``(Y, X)`` when it is valid. + Shape of the mask without its length-1 axes, ``(Y, X)`` for a 2D mask. dtype Pixel type, or None when the header could not be read. + axes + One letter per entry of ``shape``, from the mask's own metadata: ``"YX"`` for a 2D + mask, with ``Z`` and/or ``T`` for a mask with several depth slices or time points + (e.g. ``"TYX"``). Empty when the header could not be read. """ image: str @@ -74,6 +87,15 @@ class MaskFile: source: str shape: tuple[int, ...] = () dtype: np.dtype | None = None + axes: str = "" + + @property + def planes(self) -> dict[str, int]: + """Number of depth slices (``"Z"``) and time points (``"T"``); 1 for an axis it lacks.""" + return { + axis: self.shape[self.axes.index(axis)] if axis in self.axes else 1 + for axis in PLANE_AXES + } @dataclass(frozen=True) @@ -117,7 +139,14 @@ def to_manifest(self) -> dict: "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} + { + "image": f.image, + "type": f.mask_type, + "path": f.path, + "source": f.source, + "axes": f.axes, + "shape": list(f.shape), + } for f in self.files ], } @@ -178,26 +207,97 @@ def find_masks( ) -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). +def read_mask( + mask: MaskFile, + masks_folder: Path | None, + experiment_root: Path, + *, + z: int | None = None, + t: int | None = None, +) -> np.ndarray: + """Read one plane of a mask as a 2D array. + + Parameters + ---------- + mask + The mask, as found by :func:`find_masks`. + masks_folder, experiment_root + Where its file is (:attr:`MaskInventory.folder` for a TIFF mask, the experiment for + a label image). + z, t + The depth slice and time point to read, zero-based, the image's ``z`` and ``t``. Each + is required when the mask has that axis (``mask.axes``), and may only be None or 0 + otherwise: a mask's planes are never guessed (issue #36). Raises ------ ValueError - If the mask is not 2D once its length-1 axes are removed. + If the mask's dimensions are not the ones its header gave (``mask.axes``), or ``z`` + or ``t`` is missing or outside the mask. """ + wanted = {"Z": z, "T": t} + for axis, index in wanted.items(): + length = mask.planes[axis] if mask.axes else 1 + if index is None and length > 1: + raise ValueError( + f"the mask {mask.path} has {length} planes along {axis}; say which one to read" + ) + if index is not None and not 0 <= index < length: + raise ValueError(f"the mask {mask.path} has no plane {axis}={index}") if mask.source == FOLDER: - pixels = tifffile.imread(Path(masks_folder) / mask.path) + pixels, axes = _read_tiff_plane(Path(masks_folder) / mask.path, wanted) 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})") + labels = open_ngff_labels(image, name) + selection = tuple( + (wanted[a] or 0) if a in wanted and n > 1 else (0 if n == 1 else slice(None)) + for a, n in zip(labels.axes, labels.array.shape, strict=True) + ) + pixels = np.asarray(labels.array[selection]) + axes = "".join(a for a, s in zip(labels.axes, selection, strict=True) if s == slice(None)) + if axes != "YX" or (mask.axes and pixels.shape != tuple(_yx(mask.axes, mask.shape))): + raise ValueError(f"the mask {mask.path} is not the 2D plane its header gave ({axes})") return pixels +def _read_tiff_plane(path: Path, wanted: Mapping[str, int | None]) -> tuple[np.ndarray, str]: + """Read the plane at ``wanted`` Z and T of a TIFF mask; return it and its axes. + + One page is read when the planes are stored one per page (the usual layout); otherwise + the whole mask is read and the plane taken from it. + """ + with tifffile.TiffFile(path) as tif: + series = tif.series[0] + axes, shape = series.axes, tuple(series.shape) + index = tuple( + (wanted.get(a) or 0) if n > 1 and a in wanted else (0 if n == 1 else slice(None)) + for a, n in zip(axes, shape, strict=True) + ) + kept = "".join(a for a, i in zip(axes, index, strict=True) if i == slice(None)) + pages = len(series.pages) + split = next( + (k for k in range(len(shape) + 1) if math.prod(shape[:k]) == pages), + None, + ) + if ( + split is not None + and all(isinstance(i, int) for i in index[:split]) + and all(n == 1 or a in "YX" for a, n in zip(axes[split:], shape[split:], strict=True)) + and pages > 1 + ): + key = int(np.ravel_multi_index(index[:split], shape[:split])) + pixels = tif.asarray(key=key, series=0).reshape(shape[split:]) + pixels = pixels[index[split:]] + else: + pixels = series.asarray()[index] + return pixels, kept + + +def _yx(axes: str, shape: tuple[int, ...]) -> tuple[int, int]: + return shape[axes.index("Y")], shape[axes.index("X")] + + def _folder_issues(folder: Path, experiment_root: Path) -> list[ValidationIssue]: if not folder.is_dir(): return [ @@ -368,11 +468,14 @@ def _duplicate_issues(files: list[MaskFile]) -> list[ValidationIssue]: 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.""" + """Read each mask's header (axes, shape, pixel type) and check it against its image. + + Also check that the masks of one image have the same depth slices and time points. + """ problems: dict[str, list[str]] = {} for index, mask in enumerate(files): try: - shape, dtype = _mask_header(mask, folder, root) + axes, shape, dtype = _mask_header(mask, folder, root) except ( OSError, ValueError, @@ -383,13 +486,23 @@ def _header_issues( 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: + if not _known_axes(axes): problems.setdefault("mask_not_2d", []).append(mask.path) - elif shape != _image_yx(images[mask.image]): + continue + files[index] = MaskFile( + mask.image, mask.mask_type, mask.path, mask.source, shape, dtype, axes + ) + if _yx(axes, 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) + by_image: dict[str, list[MaskFile]] = {} + for mask in files: + if mask.axes: + by_image.setdefault(mask.image, []).append(mask) + for masks in by_image.values(): + if len({tuple(m.planes.items()) for m in masks}) > 1: + problems.setdefault("mask_planes_differ", []).extend(m.path for m in masks) return [ error(code, _TEXTS[code][0], paths, _TEXTS[code][1]) for code, paths in sorted(problems.items()) @@ -403,8 +516,15 @@ def _header_issues( "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.", + "Some masks have dimensions other than Y and X that their metadata does not name as " + "depth (Z) or time (T), so their planes cannot be matched to the image's.", + "Save each mask as a 2D label image, or with its depth and time axes named in its " + "metadata (for example an ImageJ hyperstack or OME-TIFF with axes TYX or ZYX).", + ), + "mask_planes_differ": ( + "The masks of some images have different numbers of depth slices or time points, so " + "they cannot all belong to the same planes of their image.", + "Make every mask of an image with the image's depth slices and time points.", ), "mask_size_mismatch": ( "Some masks do not have the size (Y and X) of their image.", @@ -419,17 +539,91 @@ def _header_issues( 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.""" + """Return a mask's axes and shape without its length-1 axes, and its pixel type. + + The axes are the mask's own metadata (``tifffile``'s series axes, or the OME-NGFF axes of + the label image); no pixels are 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 + axes, shape, dtype = tif.series[0].axes, 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("=") + labels = open_ngff_labels(root / image, name) + axes, shape, dtype = labels.axes, labels.array.shape, np.dtype(labels.array.dtype) + kept = [(a, int(n)) for a, n in zip(axes, shape, strict=True) if n != 1] + return ( + "".join(a for a, _ in kept), + tuple(n for _, n in kept), + dtype.newbyteorder("="), + ) + + +def _known_axes(axes: str) -> bool: + """Return True if the axes are Y and X, each once, plus at most one Z and one T.""" + return ( + axes.count("Y") == 1 + and axes.count("X") == 1 + and all(axes.count(a) <= 1 for a in PLANE_AXES) + and set(axes) <= {"Y", "X", *PLANE_AXES} + ) + + +def mask_plane_issues( + inventory: MaskInventory, + images: Iterable[TiffMetadata], + file_mappings: Mapping[str, AxisMapping], +) -> list[ValidationIssue]: + """Check that every mask has the depth slices and time points of its image (issue #36). + + An image's depth slices and time points come from its confirmed axis roles; a mask's + from its own metadata (:attr:`MaskFile.planes`). They must be equal: a 2D mask belongs + to an image with one plane, never to every plane of a Z-stack or time-lapse. + + Parameters + ---------- + inventory + The masks found by :func:`find_masks`. + images + Header information of the analysed images. + file_mappings + Each image's own axis mapping (:attr:`.axes.AxisReview.file_mappings`); images + without one are left to the axis review's errors. + + Returns + ------- + list of ValidationIssue + The error ``mask_planes_mismatch`` naming the masks that do not fit, or nothing. + """ + headers = {meta.relative_path: meta for meta in images} + mismatched = [] + for mask in inventory.files: + mapping = file_mappings.get(mask.image) + if not mask.axes or mapping is None or mask.image not in headers: + continue + shape = headers[mask.image].shape + planes = {} + for axis in PLANE_AXES: + index = mapping.index_of(axis) + planes[axis] = 1 if index is None else int(shape[index]) + if planes != mask.planes: + mismatched.append(mask.path) + if not mismatched: + return [] + return [ + error( + "mask_planes_mismatch", + f"{plural(len(mismatched), 'mask')} {_do(mismatched)} not have the depth slices " + f"and time points of {_their(mismatched)} image: a 2D mask fits only an image with " + "one plane, and a mask with depth (Z) or time (T) needs the same number of slices " + "or time points as its image.", + mismatched, + "Save one mask plane for every depth slice and time point of the image, with the " + "axes named in the mask's metadata, or check the image's axis roles.", + ) + ] def _image_yx(meta: TiffMetadata) -> tuple[int, int] | None: @@ -459,6 +653,14 @@ def _missing_issues( return issues +def _do(items: list) -> str: + return "does" if len(items) == 1 else "do" + + +def _their(items: list) -> str: + return "its" if len(items) == 1 else "their" + + def _are(items: list) -> str: return "is" if len(items) == 1 else "are" diff --git a/src/morphobricks_preprocessing/run.py b/src/morphobricks_preprocessing/run.py index cc413170..256311ca 100644 --- a/src/morphobricks_preprocessing/run.py +++ b/src/morphobricks_preprocessing/run.py @@ -69,7 +69,7 @@ read_manifest, write_manifest, ) -from .masks import MaskInventory, find_masks +from .masks import MaskInventory, find_masks, mask_plane_issues from .metadata import MetadataReport, TiffMetadata, inspect_images from .planes import ( DEFAULT_MEMORY_LIMIT, @@ -396,8 +396,10 @@ def run_experiment( ``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. + contract §9); an error about them stops the run, including a mask whose + depth slices and time points are not its image's under the confirmed + roles (``mask_planes_mismatch``). Their pixels are not read: the mask + step (:func:`.mask_scan.scan_masks`) does that. image_type The image type the user confirmed, ``"fluorescence"`` or ``"rgb"`` (:mod:`.image_type`), recorded in the manifest. ``"rgb"`` for images that @@ -451,6 +453,11 @@ def run_experiment( *design_problems, *dataset.issues, *(pre.masks.issues if pre.masks is not None else ()), + *( + mask_plane_issues(pre.masks, pre.metadata.files, review.file_mappings) + if pre.masks is not None + else () + ), ] dtype = pre.metadata.dtype bins = None diff --git a/tests/test_mask_scan.py b/tests/test_mask_scan.py index 06d0a8ef..3c206f86 100644 --- a/tests/test_mask_scan.py +++ b/tests/test_mask_scan.py @@ -317,3 +317,36 @@ def test_both_layouts_give_the_same_report(tmp_path): 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) + + +# Masks with time points (issue #36) --------------------------------------------------- + + +def test_masks_with_time_points_are_scanned_plane_by_plane(tmp_path): + frames = hierarchy_masks() + # Three time points: the same masks, then the nuclei moved, then no nuclei at all. + moved = np.roll(frames["Nuclei"], 30, axis=1) + series = { + "Cells": np.stack([frames["Cells"]] * 3), + "Nuclei": np.stack([frames["Nuclei"], moved, np.zeros_like(moved)]), + } + root, folder = tmp_path / "exp", tmp_path / "masks" + write_data(root / "ctrl" / "a.tif", np.zeros((3, *SHAPE), dtype=np.uint8), "TYX") + for name, data in series.items(): + write_data(folder / name / "ctrl" / "a.tif", data, "TYX") + inventory = preflight(root, masks_folder=folder).masks + assert inventory.ok, inventory.issues + report = scan_masks(inventory, root) + # The same as three single-plane images, one per time point. + flat_root, flat_folder = tmp_path / "flat", tmp_path / "flat_masks" + for t in range(3): + write_data(flat_root / "ctrl" / f"a{t}.tif", np.zeros(SHAPE, dtype=np.uint8)) + for name, data in series.items(): + write_data(flat_folder / name / "ctrl" / f"a{t}.tif", data[t]) + flat = scan_masks(preflight(flat_root, masks_folder=flat_folder).masks, flat_root) + assert report.images == 1 and flat.images == 3 + assert report.types["Nuclei"].objects == flat.types["Nuclei"].objects == 8 + assert report.types["Cells"].images == 1 + assert [(p.child, p.parent, dict(p.evidence)) for p in report.pairs] == [ + (p.child, p.parent, dict(p.evidence)) for p in flat.pairs + ] diff --git a/tests/test_masks.py b/tests/test_masks.py index 508ea129..8ff44f7c 100644 --- a/tests/test_masks.py +++ b/tests/test_masks.py @@ -251,6 +251,8 @@ def test_the_run_records_the_masks_in_the_manifest(experiment, masks, tmp_path): "type": "Cells", "path": "Cells/ctrl/rep1/a.tif", "source": "folder", + "axes": "YX", + "shape": [16, 20], } # Without masks the manifest has no masks key, as before. run_experiment(experiment, prefill_mapping("CYX"), output_folder=out) @@ -304,3 +306,132 @@ def test_masks_command_takes_thresholds(experiment, masks, tmp_path, capsys): 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 + + +# Planes: depth slices and time points (issue #36) ------------------------------------- + + +def time_lapse(tmp_path, n_t=3, image_axes="TYX"): + """An experiment of one time-lapse image (T, Y, X) and a masks folder.""" + root = tmp_path / "exp" + write_data(root / "ctrl" / "a.tif", np.zeros((n_t, 16, 20), dtype=np.uint8), image_axes) + return root, tmp_path / "masks" + + +def run_with_masks(root, masks, roles, tmp_path): + return run_experiment( + root, prefill_mapping(roles), output_folder=tmp_path / "out", masks_folder=masks + ) + + +def test_a_2d_mask_of_a_time_lapse_or_z_stack_is_refused(tmp_path): + for roles in ("TYX", "ZYX"): + root, masks = time_lapse(tmp_path / roles, image_axes=roles) + write_data(masks / "Cells" / "ctrl" / "a.tif", label_image(0)) + assert preflight(root, masks_folder=masks).masks.ok # the roles are not confirmed yet + result = run_with_masks(root, masks, roles, tmp_path / roles) + assert not result.written + mismatch = next(i for i in result.issues if i.code == "mask_planes_mismatch") + assert mismatch.paths == ("Cells/ctrl/a.tif",) + assert "a 2D mask fits only an image with one plane" in mismatch.message + + +def test_a_mask_with_the_images_time_points_is_accepted_and_recorded(tmp_path, no_pixel_reads): + root, masks = time_lapse(tmp_path) + data = np.stack([label_image(t) for t in range(3)]) + write_data(masks / "Cells" / "ctrl" / "a.tif", data, "TYX") + found = preflight(root, masks_folder=masks).masks + assert found.ok + cells = found.for_image("ctrl/a.tif")["Cells"] + assert (cells.axes, cells.shape, cells.planes) == ("TYX", (3, 16, 20), {"Z": 1, "T": 3}) + assert found.to_manifest()["files"][0]["axes"] == "TYX" + assert found.to_manifest()["files"][0]["shape"] == [3, 16, 20] + + +def test_the_run_records_each_masks_planes(tmp_path): + root, masks = time_lapse(tmp_path) + write_data(masks / "Cells" / "ctrl" / "a.tif", np.stack([label_image(0)] * 3), "TYX") + result = run_with_masks(root, masks, "TYX", tmp_path) + assert result.written, result.issues + entry = read_manifest(tmp_path / "out" / "run_manifest.json").masks["files"][0] + assert (entry["axes"], entry["shape"]) == ("TYX", [3, 16, 20]) + + +def test_a_mask_with_other_planes_than_its_image_is_refused(tmp_path): + root, masks = time_lapse(tmp_path) + write_data(masks / "Cells" / "ctrl" / "a.tif", np.stack([label_image(0)] * 2), "TYX") + write_data(masks / "Nuclei" / "ctrl" / "a.tif", np.stack([label_image(0)] * 3), "ZYX") + result = run_with_masks(root, masks, "TYX", tmp_path) + assert not result.written + found = {i.code: i.paths for i in result.issues} + assert found["mask_planes_mismatch"] == ("Cells/ctrl/a.tif", "Nuclei/ctrl/a.tif") + # The two masks of one image do not even agree with each other: found from the headers. + assert found["mask_planes_differ"] == ("Cells/ctrl/a.tif", "Nuclei/ctrl/a.tif") + + +def test_a_mask_whose_extra_axis_is_not_named_is_refused(tmp_path): + root, masks = time_lapse(tmp_path) + # No axes in the metadata: tifffile calls the stack Q, which could be depth or time. + write_data(masks / "Cells" / "ctrl" / "a.tif", np.stack([label_image(0)] * 3)) + found = preflight(root, masks_folder=masks).masks + assert codes(found.issues) == ["mask_not_2d"] + assert "does not name as depth (Z) or time (T)" in found.issues[0].message + + +def test_read_mask_reads_one_plane_and_never_guesses_which(tmp_path): + from morphobricks_preprocessing import read_mask + + root, masks = time_lapse(tmp_path) + data = np.stack([label_image(t) for t in range(3)]) + write_data(masks / "Cells" / "ctrl" / "a.tif", data, "TYX") + found = preflight(root, masks_folder=masks).masks + cells = found.for_image("ctrl/a.tif")["Cells"] + for t in range(3): + np.testing.assert_array_equal(read_mask(cells, found.folder, root, t=t), data[t]) + with pytest.raises(ValueError, match="3 planes along T"): + read_mask(cells, found.folder, root) + with pytest.raises(ValueError, match="no plane T=3"): + read_mask(cells, found.folder, root, t=3) + with pytest.raises(ValueError, match="no plane Z=1"): + read_mask(cells, found.folder, root, z=1, t=0) + + +def test_read_mask_of_a_volume_stored_in_one_page(tmp_path): + from morphobricks_preprocessing import read_mask + + root = tmp_path / "exp" + write_data(root / "ctrl" / "a.tif", np.zeros((3, 16, 20), dtype=np.uint8), "ZYX") + data = np.stack([label_image(z) for z in range(3)]) + write_data( + tmp_path / "masks" / "Cells" / "ctrl" / "a.tif", + data, + "ZYX", + volumetric=True, + tile=(3, 16, 16), + ) + found = preflight(root, masks_folder=tmp_path / "masks").masks + cells = found.for_image("ctrl/a.tif")["Cells"] + assert cells.planes == {"Z": 3, "T": 1} + np.testing.assert_array_equal(read_mask(cells, found.folder, root, z=2), data[2]) + + +def test_label_images_with_time_points(tmp_path): + pytest.importorskip("zarr") + from morphobricks_preprocessing import read_mask + from zarr_tree import write_labels, write_ngff + + path = write_ngff(tmp_path / "exp" / "a.ome.zarr", np.zeros((3, 16, 20), np.uint8), "tyx") + data = np.stack([label_image(t) for t in range(3)]) + write_labels(path, {"cells": data}, "tyx") + found = preflight(tmp_path / "exp", ome_zarr_labels=True).masks + assert found.ok + cells = found.for_image("a.ome.zarr")["cells"] + assert (cells.axes, cells.planes) == ("TYX", {"Z": 1, "T": 3}) + np.testing.assert_array_equal(read_mask(cells, None, tmp_path / "exp", t=1), data[1]) + result = run_experiment( + tmp_path / "exp", + prefill_mapping("TYX"), + output_folder=tmp_path / "out", + ome_zarr_labels=True, + ) + assert result.written, result.issues