diff --git a/AGENTS.md b/AGENTS.md index 8763c917..27b3ed5f 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -160,7 +160,8 @@ src/morphobricks_preprocessing/ report.py plain-language summaries shared by CLI and notebook, group statistics tables mask_text.py the masks' texts (mask_summary, mask_report_summary) and the notebook's tables of the mask step (mask_type_rows, mask_relationship_rows, - mask_hierarchy_lines); moved from report.py in PP-10d + mask_hierarchy_lines, with the mask_hierarchy_cycle warning, #105); + moved from report.py in PP-10d discovery_report.py the run report (discovery_report.json and .txt, design v0.4 §4): status, manifest SHA-256, what was found, reused, recomputed and replaced cli.py, __main__ `morphobricks-preprocessing check|run|masks|calibrate|design|convert` (also `python -m morphobricks_preprocessing`); @@ -230,7 +231,7 @@ src/morphobricks_preprocessing/ notebook_axis_step.py step 2: role dropdowns, axis check, missing-axes card (_AxisStep) notebook_axis_preview.py step 2: preview, confirmation, image type (_AxisPreview) notebook_run.py step 3: settings, what will happen, the run (_RunStep) - notebook_results.py step 4: pooled statistics per group (_ResultsStep) + notebook_results.py step 4: pooled statistics per group and channel (_ResultsStep) notebook_plots.py the PNG plots and IntensityPanel (the intensity checks card) notebook_render.py HTML helpers: issues, balance, labels, role effect, captions, tiles notebook_options.py the steps' choices and fixed settings (ROLE_OPTIONS, ...) @@ -444,8 +445,12 @@ functions the tests use. histogram equals the version-2 row's, whole and in blocks (`tests/test_schema_v2_equivalence.py` with the committed fixture in `tests/data`, written by `tests/data/make_v2_fixture.py` with the pre-P3b code). -- **Step 4** pools every channel of a group together, as with version 2 - (`report.channel_histograms`), so its table and plot are unchanged. +- **Step 4** summarises each channel of a group on its own (issue #104): one + `GroupSummary` per group and channel (`channel`, `xyc_plane_count`, + `channel_histogram_count`), a `channel` column in the table and one plot panel per + channel. Pooling every channel together (as with version 2) only with + `group_summaries(..., pool_channels=True)`, labelled `all channels (pooled)`; + `GroupSummary.plane_count` is a deprecated alias of `channel_histogram_count`. - **`position_id`** is filled in every row (D02, plan §12 question 20). Earlier tables with it blank are filled on the next run without reading pixels. - **Reading:** labels are always read back as text (`NA`, `001`, `TRUE` stay as typed). diff --git a/CHANGELOG.md b/CHANGELOG.md index b5ea7082..feb8fa51 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,37 @@ The design is `docs/restart/preprocessing-design-v0.4.md` in MorphoFeature. The are in [docs/plan.md](docs/plan.md) §12. Upgrading from v0.3: see the [migration guide](docs/migration-0.3-to-0.4.md). +## Unreleased + +- **Group statistics per channel (#104).** `group_summaries` and the notebook's step 4 no + longer pool distinct channels together: each group has one summary per channel, and the + intensity curves are drawn in one panel per channel. Single-channel tables give the same + numbers as before. + - `GroupSummary` gains `channel` (the channel index, or `ALL_CHANNELS`/None for an + aggregate), `xyc_plane_count` (XYC planes) and `channel_histogram_count` (histograms + pooled). `plane_count`, which counted each channel of a plane as a plane, remains as a + deprecated property equal to `channel_histogram_count` and warns + (`DeprecationWarning`). The positional constructor changed: build summaries by + keyword, or use `group_summaries`. + - `group_table` adds the columns `channel`, `XYC planes` and `channel histograms` (in + place of `planes`). + - The earlier all-channel aggregate is available on request, + `group_summaries(..., pool_channels=True)`, and is labelled `all channels (pooled)`. + - A group whose planes have different numbers of channels is refused (`ValueError`). + - **Migration:** no file changes. `plane_histograms.csv` (schema 4), the manifest and + the other outputs are unchanged; step 4 writes no file. Code reading + `GroupSummary.plane_count` should read `xyc_plane_count` or `channel_histogram_count` + and select rows by `channel`; to reproduce a pre-#104 table, pass + `pool_channels=True`. No other MorphoBricks repository uses these names. + +- **Cyclic mask hierarchies are named (#105).** When the direct-parent relations of the mask + step form a cycle (possible because pairs of types may be seen together in different + images), the hierarchy text no longer prints an empty or truncated tree: it starts with + the warning `mask_hierarchy_cycle`, naming every type and link of the cycle, and shows + the group as one `[cycle: …]` entry; every type is listed, and no direction is guessed. + New `MaskReport.hierarchy_cycles` and `mask_scan.hierarchy_cycles`. Acyclic hierarchies + print as before; `mask_relationships.json` is unchanged (report version 1). + ## v0.4.0 (5 November 2026) The final v0.4.0, with the small changes since v0.4.0rc5 (#147) and the release step: the diff --git a/docs/contract.md b/docs/contract.md index f0cc9299..edd998e1 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -116,6 +116,14 @@ on the next run without reading pixels. Defined by the [normalization brick's contract][normalization contract] §1.2.3: groups pool their planes' (representative, count) pairs and never average per-plane statistics. +This brick's optional group statistics (`group_summaries`, the notebook's step 4) pool +**one channel at a time** (issue #104): one summary per group and channel, with +`xyc_plane_count` (XYC planes) and `channel_histogram_count` (histograms pooled) reported +apart. Channels are pooled together only when asked (`pool_channels=True`), and the +summary then says so (`channel` None, shown as `all channels (pooled)`). A group whose +planes have different numbers of channels is refused (`ValueError`). These statistics are +shown only; no file is written, so no file schema depends on them. + ## 7. Normalization parameters Defined by the [normalization brick's contract][normalization contract] §1.2.4: every @@ -538,6 +546,15 @@ relation nucleoli inside cells is kept in the report but not in the hierarchy. A 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. +**Cycles (issue #105).** Relations come from pairwise evidence summed over the images with +both masks, so types seen together in different images can form a cycle (A inside B, B +inside C, C inside A). No parent-child hierarchy can then be inferred, and none is guessed: +`MaskReport.hierarchy_cycles` names each group of types in a cycle, and the hierarchy text +(`mask_report_summary`, `mask_hierarchy_lines`, the notebook) starts with the named warning +`mask_hierarchy_cycle` listing every type and link of the cycle, then shows the group as one +`[cycle: A, B, C]` entry; every type is listed. `mask_relationships.json` is unchanged +(report version 1): its `hierarchy` keeps every direct link, cycles included. + The report is UTF-8 JSON (two-space indent, LF), written atomically to the output folder: | Key | Meaning | diff --git a/docs/plan.md b/docs/plan.md index ef653865..bb463ede 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -1000,7 +1000,13 @@ Raised during review; to be decided before the milestones that need them: schema version 1 or 2 is recomputed completely on the next run; the manifest stays version 1, records schema version 3, and a file's `plane_count` counts its XYC planes. Step 4 reads the new layout and, as before, pools every channel of a group together, - so its table (whose "planes" still counts channel planes) and plot do not change. The + so its table (whose "planes" still counts channel planes) and plot do not change. + *Amended (6 November 2026, issue #104): step 4 no longer pools channels. Distinct + channels are different measurements (two fluorescence markers, or red, green and + blue), so each group is summarised and plotted channel by channel; the table has a + `channel` column, counts `XYC planes` and `channel histograms` apart, and an + all-channel aggregate is only made on explicit request (`pool_channels=True`) and + labelled `all channels (pooled)`. No file changes: step 4 writes nothing.* The warning `no_finite_pixels` counts the planes with a channel without finite pixels. The deprecated re-exports of the moved names (question 19) keep pointing at the brick's current objects, now with its Milestone 2b signatures; names it added later are not diff --git a/docs/user_guide.md b/docs/user_guide.md index 202bca9b..e1e872c8 100644 --- a/docs/user_guide.md +++ b/docs/user_guide.md @@ -543,11 +543,14 @@ results or comparing another grouping. - A spinning wheel shows that work is in progress while the table is read, statistics are calculated and the plot is drawn. The button says **Calculating…** and stays disabled until the calculation finishes or an error is shown; you can then try again. - - The table lists pooled statistics per group: minimum, maximum, mean, standard - deviation, median and percentiles. Every channel of a group is pooled together, as - before; its **planes** column counts each channel of a plane as one plane. - - The plot compares the groups' intensity distributions: curves further right are - brighter. + - The table lists pooled statistics per group **and channel**: minimum, maximum, mean, + standard deviation, median and percentiles. Each channel (**C0**, **C1**, …) has its + own row, since different fluorescence markers (or the red, green and blue of an RGB + image) are different measurements and are never mixed. **XYC planes** counts the + planes (an image at one time point) and **channel histograms** the histograms pooled + in that row. + - The plot compares the groups' intensity distributions, one panel per channel: curves + further right are brighter. - **Intensity checks** shows the same warnings, settings and plots as step 3 for the chosen table, so tables of other experiments can be checked too. It needs only the table, not the images. @@ -1102,6 +1105,7 @@ continues. | `channel_without_finite_pixels` | Warning (intensity check) | An image's channel has no finite pixels (every value NaN or infinite), so its intensities were not checked and it cannot be normalised from its own pixels. Exclude the image, or fix the channel in the source data and run again. | | `channel_names_differ` | Error | Channel names differ between images at the same channel number, so different stains would be pooled as one channel. Spaces at the start or end of a name are ignored (and noted in the manifest); otherwise names must match exactly, and names that differ only in inner spaces are quoted and reported as such. Check the channel order; export the images again with the same channel order and names, or exclude the images whose channels differ. | | `non_finite_pixels` | Warning | NaN or infinite pixels are left out of the histograms and counted. | +| `mask_hierarchy_cycle` | Warning (mask step) | The direct parents of some mask types form a cycle (for example A inside B in some images, B inside C in others and C inside A), so no parent-child hierarchy can be inferred for them. The hierarchy lists them together as `[cycle: …]` with every link. Check those masks; nothing is guessed. | | `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. | @@ -1230,6 +1234,7 @@ 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 +report.hierarchy_cycles # (): types whose direct parents form a cycle, if any write_mask_report(report, "/data/exp1_MorphoBricks_Preprocessing/mask_relationships.json") ``` diff --git a/src/morphobricks_preprocessing/mask_scan.py b/src/morphobricks_preprocessing/mask_scan.py index f5ef2174..b59c0ed0 100644 --- a/src/morphobricks_preprocessing/mask_scan.py +++ b/src/morphobricks_preprocessing/mask_scan.py @@ -49,7 +49,7 @@ import json import logging import os -from collections.abc import Callable, Mapping +from collections.abc import Callable, Mapping, Sequence from dataclasses import dataclass, field from pathlib import Path from types import MappingProxyType @@ -195,6 +195,16 @@ def hierarchy(self) -> Mapping[str, tuple[str, ...]]: parents[pair.child].append(pair.parent) return MappingProxyType({name: tuple(sorted(found)) for name, found in parents.items()}) + @property + def hierarchy_cycles(self) -> tuple[tuple[str, ...], ...]: + """The groups of mask types whose direct parents form a cycle (issue #105). + + Pairwise evidence from different images can make A inside B, B inside C and C + inside A: no parent-child hierarchy can then be inferred for them. Each group is + sorted, and the groups too; empty when the hierarchy has no cycle. + """ + return hierarchy_cycles(self.hierarchy) + @property def same_objects(self) -> tuple[tuple[str, str], ...]: """Pairs of types that segment the same objects, each pair once.""" @@ -228,6 +238,34 @@ def to_dict(self) -> dict: } +def hierarchy_cycles(hierarchy: Mapping[str, Sequence[str]]) -> tuple[tuple[str, ...], ...]: + """Return the groups of names whose parents form a cycle, in ``{child: parents}``. + + A group is a strongly connected set of two or more names, or one name that is its own + parent. Each group is sorted, and the groups too; empty for an acyclic hierarchy. + """ + names = {*hierarchy, *(parent for parents in hierarchy.values() for parent in parents)} + reach = {name: _ancestors(hierarchy, name) for name in names} + groups = { + tuple(sorted(other for other in found if name in reach[other])) + for name, found in reach.items() + if name in found + } + return tuple(sorted(groups)) + + +def _ancestors(hierarchy: Mapping[str, Sequence[str]], name: str) -> set[str]: + """Return the names reached from ``name`` through one or more parent links.""" + found: set[str] = set() + todo = list(hierarchy.get(name, ())) + while todo: + parent = todo.pop() + if parent not in found: + found.add(parent) + todo.extend(hierarchy.get(parent, ())) + return found + + def mask_thresholds(changes: Mapping[str, float] | None = None) -> Mapping[str, float]: """Return the default thresholds with some replaced, after checking them. diff --git a/src/morphobricks_preprocessing/mask_text.py b/src/morphobricks_preprocessing/mask_text.py index e849d081..a1a5c71c 100644 --- a/src/morphobricks_preprocessing/mask_text.py +++ b/src/morphobricks_preprocessing/mask_text.py @@ -10,7 +10,12 @@ from collections.abc import Mapping, Sequence -from .mask_scan import DEFAULT_THRESHOLDS, MaskReport, PairRelation +from .mask_scan import ( + DEFAULT_THRESHOLDS, + MaskReport, + PairRelation, + hierarchy_cycles, +) from .masks import MaskInventory from .validation import plural @@ -185,18 +190,47 @@ def _relationship_line(pair: PairRelation) -> str: def _hierarchy_lines(hierarchy: Mapping[str, Sequence[str]]) -> 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} + """Return the mask types as an indented tree; a type with two parents shows twice. + + Every type is listed. Types whose direct parents form a cycle (issue #105) have no + parent-child hierarchy: a named warning ``mask_hierarchy_cycle`` comes first, with + every type and link of the cycle, and the tree shows them together as one + ``[cycle: ...]`` entry (the strongly connected components), never guessing a direction. + """ + cycles = hierarchy_cycles(hierarchy) + group = {name: (name,) for name in hierarchy} + group.update({parent: (parent,) for parents in hierarchy.values() for parent in parents}) + group.update({name: cycle for cycle in cycles for name in cycle}) + children: dict[tuple[str, ...], set[tuple[str, ...]]] = {g: set() for g in group.values()} + has_parent: set[tuple[str, ...]] = set() 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]): + if group[parent] != group[name]: + children[group[parent]].add(group[name]) + has_parent.add(group[name]) + lines = [_cycle_warning(cycle, hierarchy) for cycle in cycles] + + def show(node: tuple[str, ...], depth: int) -> None: + label = f"[cycle: {', '.join(node)}]" if node in cycles else node[0] + lines.append(f"{' ' * (depth + 1)}{label}") + for child in sorted(children[node]): show(child, depth + 1) - for name in sorted(n for n, parents in hierarchy.items() if not parents): - show(name, 0) + for node in sorted(set(group.values()) - has_parent): + show(node, 0) return lines + + +def _cycle_warning(cycle: tuple[str, ...], hierarchy: Mapping[str, Sequence[str]]) -> str: + """Name a cycle's types and links, and say no hierarchy can be inferred for them.""" + links = [ + f"{name} in {parent}" + for name in cycle + for parent in hierarchy.get(name, ()) + if parent in cycle + ] + return ( + f" Warning (mask_hierarchy_cycle): the direct parents of {', '.join(cycle)} form a " + f"cycle ({'; '.join(links)}), so no parent-child hierarchy can be inferred for " + "them; they are shown together below. Check these masks." + ) diff --git a/src/morphobricks_preprocessing/notebook_plots.py b/src/morphobricks_preprocessing/notebook_plots.py index f2917a04..152f1a56 100644 --- a/src/morphobricks_preprocessing/notebook_plots.py +++ b/src/morphobricks_preprocessing/notebook_plots.py @@ -40,7 +40,6 @@ from .notebook_render import _label from .notebook_values import display_range from .planes import Plane -from .report import channel_histograms from .validation import plural @@ -433,33 +432,45 @@ def plane_png(plane: Plane, size: float = 4.2) -> bytes: def distribution_png(groups: dict[tuple[str, ...], list]) -> bytes | None: - """Plot the cumulative intensity distribution of each group (up to 12 groups).""" - curves = [] - for key, members in list(groups.items())[:_MAX_PLOTTED_GROUPS]: - try: - pooled = PooledHistogram(channel_histograms(members)) - except NoValidPixelsError: - continue - cumulative = cumulative_share(pooled) - label = " / ".join(key[-2:]) if len(key) > 1 else key[0] - curves.append((label, pooled, cumulative)) - if not curves: + """Plot the cumulative intensity distribution of each group, one panel per channel. + + Channels are never pooled together (issue #104): each panel holds one channel's + curves, one per group (up to 12 groups), each with its own intensity axis. + """ + shown = list(groups.items())[:_MAX_PLOTTED_GROUPS] + channels = max((r.channel_count for _, members in shown for r in members), default=0) + panels = [] + for channel in range(channels): + curves = [] + for key, members in shown: + try: + pooled = PooledHistogram([r.channel_histogram(channel) for r in members]) + except NoValidPixelsError: + continue + label = " / ".join(key[-2:]) if len(key) > 1 else key[0] + curves.append((label, pooled, cumulative_share(pooled))) + if curves: + panels.append((channel, curves)) + if not panels: return None - low = min(p.percentile(0.1) for _, p, _ in curves) - high = max(p.percentile(99.9) for _, p, _ in curves) - figure = Figure(figsize=(7.2, 3.4)) - axes = figure.add_subplot() - for label, pooled, cumulative in curves: - axes.step(pooled.values, cumulative, where="post", label=label, linewidth=1.4) - if high > low: - axes.set_xlim(low, high) - axes.set_ylim(0, 1) - axes.set_xlabel("Intensity") - axes.set_ylabel("Share of pixels at or below") - axes.grid(alpha=0.3) - axes.legend(fontsize=8, loc="lower right", frameon=False) - if len(groups) > _MAX_PLOTTED_GROUPS: - axes.set_title(f"First {_MAX_PLOTTED_GROUPS} of {len(groups)} groups", fontsize=9) + figure = Figure(figsize=(7.2, 3.4 * len(panels))) + for index, (channel, curves) in enumerate(panels): + axes = figure.add_subplot(len(panels), 1, index + 1) + low = min(p.percentile(0.1) for _, p, _ in curves) + high = max(p.percentile(99.9) for _, p, _ in curves) + for label, pooled, cumulative in curves: + axes.step(pooled.values, cumulative, where="post", label=label, linewidth=1.4) + if high > low: + axes.set_xlim(low, high) + axes.set_ylim(0, 1) + axes.set_xlabel(f"Intensity, channel C{channel}") + axes.set_ylabel("Share of pixels at or below") + axes.grid(alpha=0.3) + axes.legend(fontsize=8, loc="lower right", frameon=False) + title = f"Channel C{channel}" + if len(groups) > _MAX_PLOTTED_GROUPS: + title += f": first {_MAX_PLOTTED_GROUPS} of {len(groups)} groups" + axes.set_title(title, fontsize=9) figure.tight_layout() buffer = io.BytesIO() figure.savefig(buffer, format="png", dpi=100) diff --git a/src/morphobricks_preprocessing/notebook_results.py b/src/morphobricks_preprocessing/notebook_results.py index 3443a542..bc1b3706 100644 --- a/src/morphobricks_preprocessing/notebook_results.py +++ b/src/morphobricks_preprocessing/notebook_results.py @@ -217,16 +217,18 @@ def _calculate_statistics(self, path: Path, percentiles: Sequence[float]) -> Non first = records[0].channel_histogram(0) if records else None self.stats_tiles.value = ui.tiles( [ - ("🧪", f"{len(summaries):,}", f"{LEVEL_NAMES[level]} groups"), + ("🧪", f"{len(groups):,}", f"{LEVEL_NAMES[level]} groups"), ("🖼️", f"{len(images):,}", "images"), - ("🧩", f"{len(records):,}", "planes"), + ("🧩", f"{len(records):,}", "XYC planes"), + ("🎨", f"{records[0].channel_count:,}" if records else "-", "channels"), ("🔢", first.dtype.name if first else "-", "pixel type"), ("📊", f"{first.bin_count:,}" if first else "-", "bins per plane"), ] ) rows = group_table(summaries, level) depth = len(summaries[0].labels) if summaries else 1 - self.statistics.value = ui.table(rows, numeric_from=depth) + # The labels and the channel column are text; the counts and statistics numbers. + self.statistics.value = ui.table(rows, numeric_from=depth + 1) ui.set_visible(self.statistics_card, True) self.stats_status.value = ui.working("Drawing the intensity distributions…") png = distribution_png(groups) @@ -236,7 +238,7 @@ def _calculate_statistics(self, path: Path, percentiles: Sequence[float]) -> Non self.table_checks.show(records) self.stats_status.value = ui.status( "ok", - f"Statistics of {plural(len(summaries), LEVEL_NAMES[level])} calculated.", + f"Statistics of {plural(len(groups), LEVEL_NAMES[level])} calculated.", [f"From the table {path}"], ) self.stats_next.value = ui.next_steps( diff --git a/src/morphobricks_preprocessing/report.py b/src/morphobricks_preprocessing/report.py index 108b67f2..4c68b2ca 100644 --- a/src/morphobricks_preprocessing/report.py +++ b/src/morphobricks_preprocessing/report.py @@ -7,6 +7,7 @@ from __future__ import annotations +import warnings from collections import Counter from collections.abc import Iterable, Sequence from dataclasses import dataclass @@ -317,27 +318,39 @@ def _stopped_summary(result: RunResult, design_template: bool) -> str: return "\n".join(lines) +ALL_CHANNELS = None +"""The ``channel`` of a deliberate all-channel aggregate (``pool_channels=True``).""" + + @dataclass(frozen=True) class GroupSummary: - """Pooled statistics of one group of planes. + """Pooled statistics of one channel of one group of XYC planes (issue #104). Attributes ---------- labels The group's labels, from condition down to the chosen level. - plane_count - Planes in the group, each channel of an XYC plane counted as one plane, - including planes without finite pixels. + channel + The channel summarised (0 for C0), or :data:`ALL_CHANNELS` (None) for a + deliberate aggregate of every channel, asked for with ``pool_channels=True``. + xyc_plane_count + XYC planes (rows of the table: an image at one Z and T) in the group, including + planes without finite pixels. + channel_histogram_count + Channel histograms pooled: equal to ``xyc_plane_count`` for one channel, and the + planes times their channels for an all-channel aggregate. pixel_count Finite pixels pooled. minimum, maximum, mean, std, median - Pooled statistics; None if the group has no finite pixels. + Pooled statistics; None if the group has no finite pixels in that channel. percentiles ``{q: value}`` for the requested percentiles. """ labels: tuple[str, ...] - plane_count: int + channel: int | None + xyc_plane_count: int + channel_histogram_count: int pixel_count: int minimum: float | None = None maximum: float | None = None @@ -346,13 +359,34 @@ class GroupSummary: median: float | None = None percentiles: dict[float, float] | None = None + @property + def plane_count(self) -> int: + """Deprecated name of ``channel_histogram_count`` (before issue #104). + + It counted each channel of an XYC plane as one plane; use ``xyc_plane_count`` for + planes and ``channel_histogram_count`` for histograms. + """ + warnings.warn( + "GroupSummary.plane_count is deprecated: use xyc_plane_count (XYC planes) or " + "channel_histogram_count (channel histograms pooled).", + DeprecationWarning, + stacklevel=2, + ) + return self.channel_histogram_count + def group_summaries( records: Sequence[PlaneRecord], level: str, percentiles: Sequence[float] = (1.0, 99.8), + *, + pool_channels: bool = False, ) -> list[GroupSummary]: - """Pool the planes of every group at ``level`` and summarise them. + """Summarise every channel of every group at ``level``, each channel on its own. + + Distinct channels are different measurements (two fluorescence markers, or the red, + green and blue of an RGB image), so by default each group has one summary per + channel, channels in order, and nothing is pooled across channels (issue #104). Parameters ---------- @@ -363,49 +397,96 @@ def group_summaries( ``"image"``. percentiles Percentiles to report, from 0 to 100. + pool_channels + True pools every channel of a group into one summary whose ``channel`` is + :data:`ALL_CHANNELS`, the behaviour before issue #104. Only on explicit request: + the table labels it as an aggregate of all channels. + + Raises + ------ + ValueError + If the rows of one group have different numbers of channels: which channel is + which cannot be known, so nothing is summarised. """ summaries = [] for key, members in group_records(list(records), level).items(): - histograms = channel_histograms(members) - try: - pooled = PooledHistogram(histograms) - except NoValidPixelsError: - summaries.append(GroupSummary(key, len(histograms), 0)) - continue - summaries.append( - GroupSummary( - labels=key, - plane_count=len(histograms), - pixel_count=pooled.pixel_count, - minimum=pooled.minimum, - maximum=pooled.maximum, - mean=pooled.mean, - std=pooled.std, - median=pooled.median, - percentiles={q: pooled.percentile(q) for q in percentiles}, + counts = {record.channel_count for record in members} + if len(counts) > 1: + raise ValueError( + f"The planes of group {' / '.join(key)} have different numbers of channels " + f"({', '.join(map(str, sorted(counts)))}); their statistics cannot be " + "matched channel by channel." ) - ) + if pool_channels: + selections = [(ALL_CHANNELS, channel_histograms(members))] + else: + selections = [ + (c, [record.channel_histogram(c) for record in members]) + for c in range(counts.pop()) + ] + summaries += [ + _summary(key, channel, len(members), histograms, percentiles) + for channel, histograms in selections + ] return summaries +def _summary( + labels: tuple[str, ...], + channel: int | None, + planes: int, + histograms: list[PlaneHistogram], + percentiles: Sequence[float], +) -> GroupSummary: + """Pool ``histograms`` (one channel of a group, or all of them) and summarise them.""" + try: + pooled = PooledHistogram(histograms) + except NoValidPixelsError: + return GroupSummary(labels, channel, planes, len(histograms), 0) + return GroupSummary( + labels=labels, + channel=channel, + xyc_plane_count=planes, + channel_histogram_count=len(histograms), + pixel_count=pooled.pixel_count, + minimum=pooled.minimum, + maximum=pooled.maximum, + mean=pooled.mean, + std=pooled.std, + median=pooled.median, + percentiles={q: pooled.percentile(q) for q in percentiles}, + ) + + def channel_histograms(records: Iterable[PlaneRecord]) -> list[PlaneHistogram]: """Return the histogram of every channel of every row, rows first. - The group statistics pool every channel of a group together, as they did - when the table had one row per channel (schema version 2), so the table and - plot of step 4 are unchanged. + Only for a deliberate all-channel aggregate (``pool_channels=True``): the default + statistics keep each channel apart (issue #104). """ return [histogram for record in records for histogram in record.histograms] +def channel_label(channel: int | None) -> str: + """Return how the tables name a channel: ``C0``, or ``all channels (pooled)``.""" + return "all channels (pooled)" if channel is ALL_CHANNELS else f"C{channel}" + + def group_table(summaries: Sequence[GroupSummary], level: str) -> list[list[str]]: - """Turn the summaries into rows of text, with a header row, for printing or HTML.""" + """Turn the summaries into rows of text, with a header row, for printing or HTML. + + One row per group and channel: the labels, the ``channel`` (``C0``, or ``all + channels (pooled)`` for a deliberate aggregate), the ``XYC planes`` and the ``channel + histograms`` pooled, then the statistics. + """ label_columns = [LEVEL_NAMES[name] for name in (*LEVELS, "image")] depth = (*LEVELS, "image").index(level) + 1 qs = next((list(s.percentiles) for s in summaries if s.percentiles), []) header = [ *label_columns[:depth], - "planes", + "channel", + "XYC planes", + "channel histograms", "pixels", "min", "max", @@ -418,7 +499,15 @@ def group_table(summaries: Sequence[GroupSummary], level: str) -> list[list[str] for s in summaries: values = [s.minimum, s.maximum, s.mean, s.std, s.median] values += [None if not s.percentiles else s.percentiles[q] for q in qs] - rows.append([*s.labels, f"{s.plane_count:,}", f"{s.pixel_count:,}", *map(_number, values)]) + counts = (s.xyc_plane_count, s.channel_histogram_count, s.pixel_count) + rows.append( + [ + *s.labels, + channel_label(s.channel), + *(f"{n:,}" for n in counts), + *map(_number, values), + ] + ) return rows diff --git a/tests/test_mask_hierarchy_cycles.py b/tests/test_mask_hierarchy_cycles.py new file mode 100644 index 00000000..e9ea70c8 --- /dev/null +++ b/tests/test_mask_hierarchy_cycles.py @@ -0,0 +1,112 @@ +"""Cyclic mask hierarchies are named and shown whole, never as an empty tree (issue #105).""" + +from types import MappingProxyType + +import numpy as np + +from morphobricks_preprocessing import preflight, scan_masks +from morphobricks_preprocessing.mask_scan import MaskReport, MaskTypeSummary, PairRelation +from morphobricks_preprocessing.mask_text import mask_hierarchy_lines, mask_report_summary +from tiff_tree import write_data + +EVIDENCE = { + "children": 1, + "inside": 1, + "outside_parents": 0, + "split_between_parents": 0, + "ambiguous_parent": 0, + "parents_with_0_children": 0, + "parents_with_1_child": 1, + "parents_with_2_children": 0, + "parents_with_3_or_more_children": 0, +} + + +def report(*links, types=()): + """A MaskReport whose direct relations are the ``(child, parent)`` links given.""" + names = sorted({*types, *(name for link in links for name in link)}) + summaries = {name: MaskTypeSummary(name, "instances", 1, 1, 1.0) for name in names} + pairs = tuple(PairRelation(c, p, "inside", True, MappingProxyType(EVIDENCE)) for c, p in links) + return MaskReport(images=1, types=MappingProxyType(summaries), pairs=pairs) + + +def warnings(lines): + return [line for line in lines if "mask_hierarchy_cycle" in line] + + +def test_two_masks_inside_each_other_are_a_named_cycle(): + r = report(("A", "B"), ("B", "A")) + assert r.hierarchy_cycles == (("A", "B"),) + lines = mask_hierarchy_lines(r) + (warning,) = warnings(lines) + assert "A, B form a cycle (A in B; B in A)" in warning + assert "no parent-child hierarchy can be inferred" in warning + assert lines[1:] == [" [cycle: A, B]"] + + +def test_a_three_mask_cycle_is_not_an_empty_tree(): + # Before #105 no type was without a parent, so the tree printed nothing at all. + r = report(("A", "B"), ("B", "C"), ("C", "A"), types=("Free",)) + assert r.hierarchy_cycles == (("A", "B", "C"),) + lines = mask_hierarchy_lines(r) + assert "(A in B; B in C; C in A)" in warnings(lines)[0] + assert lines[1:] == [" [cycle: A, B, C]", " Free"] + text = mask_report_summary(r) + tree = text.split("Hierarchy (direct parents only):\n")[1] + assert "Warning (mask_hierarchy_cycle)" in tree and "[cycle: A, B, C]" in tree + + +def test_a_cycle_below_a_root_is_shown_without_endless_recursion(): + # Cells is a root; Nuclei and Spots are inside each other and both inside Cells. + r = report(("Nuclei", "Cells"), ("Nuclei", "Spots"), ("Spots", "Nuclei"), ("Dots", "Spots")) + assert r.hierarchy_cycles == (("Nuclei", "Spots"),) + lines = mask_hierarchy_lines(r) + assert len(warnings(lines)) == 1 + assert lines[1:] == [" Cells", " [cycle: Nuclei, Spots]", " Dots"] + + +def test_a_mask_that_is_its_own_parent_is_a_cycle(): + r = report(("A", "A"), ("B", "A")) + assert r.hierarchy_cycles == (("A",),) + lines = mask_hierarchy_lines(r) + assert "(A in A)" in warnings(lines)[0] + assert lines[1:] == [" [cycle: A]", " B"] + + +def test_a_valid_hierarchy_is_unchanged(): + r = report(("Cells", "Organoid"), ("Nuclei", "Cells"), types=("Other",)) + assert r.hierarchy_cycles == () + assert mask_hierarchy_lines(r) == [ + " Organoid", + " Cells", + " Nuclei", + " Other", + ] + + +def test_a_type_with_two_parents_still_shows_under_each(): + r = report(("Dots", "Left"), ("Dots", "Top")) + assert r.hierarchy_cycles == () + assert mask_hierarchy_lines(r) == [" Left", " Dots", " Top", " Dots"] + + +def test_sparse_pairwise_evidence_from_a_scan_makes_a_named_cycle(tmp_path): + # Each pair of types is seen together in one image only: A in B, B in C, C in A. + def box(size): + mask = np.zeros((32, 32), np.uint16) + start = (32 - size) // 2 + mask[start : start + size, start : start + size] = 1 + return mask + + root, folder = tmp_path / "exp", tmp_path / "masks" + for image, (small, big) in {"one": ("A", "B"), "two": ("B", "C"), "three": ("C", "A")}.items(): + write_data(root / image / "x.tif", np.zeros((32, 32), np.uint8)) + write_data(folder / small / image / "x.tif", box(6)) + write_data(folder / big / image / "x.tif", box(24)) + inventory = preflight(root, masks_folder=folder).masks + assert inventory.ok, inventory.issues + scanned = scan_masks(inventory, root) + assert dict(scanned.hierarchy) == {"A": ("B",), "B": ("C",), "C": ("A",)} + assert scanned.hierarchy_cycles == (("A", "B", "C"),) + lines = mask_hierarchy_lines(scanned) + assert len(warnings(lines)) == 1 and lines[1:] == [" [cycle: A, B, C]"] diff --git a/tests/test_notebook_debt.py b/tests/test_notebook_debt.py index 20a4d26d..28f9f185 100644 --- a/tests/test_notebook_debt.py +++ b/tests/test_notebook_debt.py @@ -76,3 +76,30 @@ def test_the_unused_pyramid_options_are_removed(): assert not hasattr(notebook_options, "PYRAMID_OPTIONS") assert "PYRAMID_OPTIONS" not in notebook_ui.__all__ + + +def test_distribution_png_draws_one_panel_per_channel(monkeypatch): + # Issue #104: the curves are drawn per channel, never from every channel pooled. + from morphobricks_normalization.histogram import compute_plane_histogram + from morphobricks_normalization.records import PlaneRecord + + from morphobricks_preprocessing import notebook_plots + + channels = (np.full((8, 8), 10, np.uint16), np.full((8, 8), 1000, np.uint16)) + record = PlaneRecord( + "c/r/t/a.tif", + "c", + "r", + "t", + histograms=tuple(compute_plane_histogram(p, 16) for p in channels), + ) + used = [] + real = notebook_plots.cumulative_share + + def spy(pooled): + used.append(pooled) + return real(pooled) + + monkeypatch.setattr(notebook_plots, "cumulative_share", spy) + assert notebook_plots.distribution_png({("c",): [record]}) is not None + assert [(p.minimum, p.maximum, p.pixel_count) for p in used] == [(10, 10, 64), (1000, 1000, 64)] diff --git a/tests/test_notebook_ui.py b/tests/test_notebook_ui.py index 08d9290f..ea72aa5b 100644 --- a/tests/test_notebook_ui.py +++ b/tests/test_notebook_ui.py @@ -1248,8 +1248,8 @@ def test_full_run_and_statistics(experiment, tmp_path): assert "The table from step 3 is selected" in nb.stats_status.value assert not ui.is_visible(nb.statistics_card) nb.statistics_button.click() - assert ui.is_visible(nb.statistics_card) and "ctrl" in nb.statistics.value - assert "Statistics of 2 conditions calculated" in nb.stats_status.value + assert ui.is_visible(nb.statistics_card) and nb.statistics.value.count("ctrl") == 2 + assert "Statistics of 2 conditions calculated" in nb.stats_status.value # 2 channels: #104 assert "plane_histograms.csv" in nb.stats_status.value assert ui.is_visible(nb.distribution_card) assert bytes(nb.distribution.value).startswith(b"\x89PNG") diff --git a/tests/test_report.py b/tests/test_report.py index 73634d84..bdca0445 100644 --- a/tests/test_report.py +++ b/tests/test_report.py @@ -1,8 +1,11 @@ +from dataclasses import replace + import numpy as np import pytest from morphobricks_normalization.csv_io import read_plane_histograms +from morphobricks_normalization.histogram import compute_plane_histogram from morphobricks_normalization.pooling import PooledHistogram -from morphobricks_normalization.records import group_records +from morphobricks_normalization.records import PlaneRecord, group_records from morphobricks_preprocessing import ( discover_experiment, @@ -12,6 +15,7 @@ run_experiment, ) from morphobricks_preprocessing.report import ( + ALL_CHANNELS, format_issues, format_table, group_table, @@ -123,6 +127,7 @@ def test_run_summary_when_blocked(experiment): def test_group_summaries_pool_the_planes(experiment): + # One channel: one summary per group, the same numbers as before issue #104. result = run_experiment(experiment, prefill_mapping("TYX")) records = read_plane_histograms(result.csv_path) summaries = group_summaries(records, "biological_repeat", percentiles=(5, 95)) @@ -131,24 +136,124 @@ def test_group_summaries_pool_the_planes(experiment): summaries, group_records(records, "biological_repeat").values(), strict=True ): pooled = PooledHistogram(members) - assert summary.plane_count == len(members) + assert summary.channel == 0 + assert summary.xyc_plane_count == summary.channel_histogram_count == len(members) assert summary.pixel_count == pooled.pixel_count assert summary.mean == pooled.mean and summary.median == pooled.median assert summary.percentiles == {5: pooled.percentile(5), 95: pooled.percentile(95)} + # One channel pooled "across channels" is the same summary, labelled as an aggregate. + pooled = group_summaries(records, "biological_repeat", (5, 95), pool_channels=True) + assert [replace(s, channel=0) for s in pooled] == summaries def test_group_table_and_text(experiment): result = run_experiment(experiment, prefill_mapping("TYX")) records = read_plane_histograms(result.csv_path) rows = group_table(group_summaries(records, "condition"), "condition") - assert rows[0] == ["condition", "planes", "pixels", "min", "max", "mean", "std", "median"] + [ + assert rows[0] == [ + "condition", + "channel", + "XYC planes", + "channel histograms", + "pixels", + "min", + "max", + "mean", + "std", + "median", "p1", "p99.8", ] assert [row[0] for row in rows[1:]] == ["ctrl", "drug"] - assert rows[1][1:3] == ["6", "864"] + assert rows[1][1:5] == ["C0", "6", "6", "864"] text = format_table(rows) - assert text.splitlines()[0].startswith("condition planes") + assert text.splitlines()[0].startswith("condition channel XYC planes") + + +def channel_records(levels, condition="c", repeat="r", planes=2, size=8, seed=0): + """Records of ``planes`` XYC planes whose channel ``c`` holds values around ``levels[c]``.""" + rng = np.random.default_rng(seed) + records = [] + for i in range(planes): + histograms = tuple( + compute_plane_histogram( + (level + rng.integers(0, 3, (size, size))).astype(np.uint16), 256 + ) + for level in levels + ) + records.append( + PlaneRecord(f"{condition}/{repeat}/x.tif", condition, repeat, "single", histograms, t=i) + ) + return records + + +@pytest.mark.parametrize( + "levels", [(10, 1000), (10, 1000, 30000), (5, 60000, 200, 4000)], ids=["2", "3", "4"] +) +def test_group_summaries_keep_fluorescence_channels_apart(levels): + # Issue #104: strongly separated channels are summarised one by one, never mixed. + records = channel_records(levels) + channel_records(levels, condition="d", seed=1) + summaries = group_summaries(records, "condition", percentiles=(50,)) + assert [(s.labels, s.channel) for s in summaries] == [ + ((condition,), c) for condition in ("c", "d") for c in range(len(levels)) + ] + for s in summaries: + level = levels[s.channel] + assert level <= s.minimum and s.maximum <= level + 2 + assert level <= s.mean <= level + 2 + # Two XYC planes, one histogram of this channel each: not planes x channels. + assert s.xyc_plane_count == s.channel_histogram_count == 2 + assert s.pixel_count == 2 * 64 + rows = group_table(summaries, "condition") + assert [row[1] for row in rows[1:]] == [f"C{c}" for c in range(len(levels))] * 2 + assert {row[2] for row in rows[1:]} == {"2"} + + +def test_an_all_channel_aggregate_is_explicit_and_labelled(): + records = channel_records((10, 1000)) + (summary,) = group_summaries(records, "condition", pool_channels=True) + assert summary.channel is ALL_CHANNELS + assert (summary.xyc_plane_count, summary.channel_histogram_count) == (2, 4) + assert 10 < summary.mean < 1000 # the mixed value the default never reports + row = group_table([summary], "condition")[1] + assert row[1:4] == ["all channels (pooled)", "2", "4"] + + +def test_a_channel_without_finite_pixels_is_reported_alone(): + finite = compute_plane_histogram(np.full((8, 8), 5.0, np.float32), 16) + empty = compute_plane_histogram(np.full((8, 8), np.nan, np.float32), 16) + records = [PlaneRecord("c/r/x.tif", "c", "r", "single", (finite, empty), t=t) for t in range(2)] + first, second = group_summaries(records, "condition") + assert first.channel == 0 and first.mean == 5.0 and first.pixel_count == 128 + assert second.channel == 1 and second.pixel_count == 0 and second.mean is None + assert second.xyc_plane_count == 2 + assert group_table([first, second], "condition")[2][1:] == ["C1", "2", "2", "0"] + ["–"] * 7 + + +def test_rgb_channels_are_summarised_one_by_one(tmp_path): + root = tmp_path / "exp" + rgb = np.stack([np.full((8, 8), v, np.uint8) for v in (20, 120, 240)]) + write_data(root / "a/x.tif", rgb, axes="CYX") + result = run_experiment(root, prefill_mapping("CYX"), image_type="rgb") + assert result.written + summaries = group_summaries(read_plane_histograms(result.csv_path), "condition") + assert [(s.channel, s.mean, s.xyc_plane_count) for s in summaries] == [ + (0, 20.0, 1), + (1, 120.0, 1), + (2, 240.0, 1), + ] + + +def test_a_group_with_different_channel_counts_is_refused(): + records = channel_records((10, 1000)) + channel_records((10,), seed=1) + with pytest.raises(ValueError, match="different numbers of channels"): + group_summaries(records, "condition") + + +def test_plane_count_is_a_deprecated_alias_of_the_histogram_count(): + (summary,) = group_summaries(channel_records((10, 1000)), "condition", pool_channels=True) + with pytest.warns(DeprecationWarning, match="xyc_plane_count"): + assert summary.plane_count == 4 def test_group_without_finite_pixels(tmp_path): @@ -159,4 +264,4 @@ def test_group_without_finite_pixels(tmp_path): summaries = group_summaries(read_plane_histograms(result.csv_path), "condition") empty = summaries[1] assert empty.labels == ("b",) and empty.pixel_count == 0 and empty.mean is None - assert group_table(summaries, "condition")[2][3:] == ["–"] * 7 + assert group_table(summaries, "condition")[2][5:] == ["–"] * 7 diff --git a/tests/test_schema_v2_equivalence.py b/tests/test_schema_v2_equivalence.py index 0de2f044..e6a285c8 100644 --- a/tests/test_schema_v2_equivalence.py +++ b/tests/test_schema_v2_equivalence.py @@ -14,6 +14,7 @@ import csv import json import shutil +from dataclasses import replace import numpy as np import pytest @@ -161,10 +162,33 @@ def test_step_4_statistics_are_those_of_version_2(tmp_path, pixels, dtype, level root = tmp_path / "exp" bins, axes = write_experiment(root, dtype, pixels) result = run_experiment(root, prefill_mapping(axes), bin_count=bins) - now = group_summaries(read_plane_histograms(result.csv_path), level, (1.0, 99.8)) - before = group_summaries(v2_records(dtype), level, (1.0, 99.8)) + records = read_plane_histograms(result.csv_path) + # Issue #104: each channel on its own, from the version-2 rows of that channel. + now = group_summaries(records, level, (1.0, 99.8)) + before = group_summaries(v2_planes(dtype), level, (1.0, 99.8)) assert now == before assert group_table(now, level) == group_table(before, level) + # The all-channel aggregate, asked for, has the numbers step 4 showed before #104. + pooled = group_summaries(records, level, (1.0, 99.8), pool_channels=True) + legacy = group_summaries(v2_records(dtype), level, (1.0, 99.8), pool_channels=True) + statistics = ("pixel_count", "minimum", "maximum", "mean", "std", "median", "percentiles") + for p, old in zip(pooled, legacy, strict=True): + assert p.labels == old.labels + assert p.channel_histogram_count == old.xyc_plane_count # a v2 row was one channel + assert [getattr(p, f) for f in statistics] == [getattr(old, f) for f in statistics] + + +def v2_planes(dtype): + """The version-2 rows gathered into XYC planes, channels in order (``c``).""" + planes = {} + for record, row in zip(v2_records(dtype), read_rows(DATA / f"v2_{dtype}.csv"), strict=True): + key = (record.source_relative_path, row["z"], row["t"]) + planes.setdefault(key, {})[int(row["c"])] = record + return [ + replace(first, histograms=tuple(channels[c].histograms[0] for c in sorted(channels))) + for channels in planes.values() + for first in [channels[min(channels)]] + ] def test_a_version_2_table_is_recomputed_on_the_next_run(tmp_path, pixels):