From df5b7720847714f4f90bf21eb08e2169981799c0 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 07:16:14 +0000 Subject: [PATCH 01/12] Add RGB-to-scalar stain presets (H&E first) New stain module with rgb_to_scalar and a read-only registry of four P0 presets (he_eosin, he_hematoxylin, luminance_inverted, lab_a). 8-bit sRGB, axes YXC or CYX, refusals instead of repairs, optional histology extra (scikit-image, imported lazily). Contract section 3.4 and plan N67 record the scope change; rgb.py and the file formats are untouched. Refs #94 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- AGENTS.md | 5 +- CHANGELOG.md | 8 + docs/contract.md | 44 +++++ docs/plan.md | 14 +- pyproject.toml | 2 + src/morphobricks_normalization/stain.py | 235 +++++++++++++++++++++++ tests/test_stain.py | 236 ++++++++++++++++++++++++ 7 files changed, 542 insertions(+), 2 deletions(-) create mode 100644 src/morphobricks_normalization/stain.py create mode 100644 tests/test_stain.py diff --git a/AGENTS.md b/AGENTS.md index da806db..6f4eb9e 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -56,6 +56,8 @@ src/morphobricks_normalization/ scopes.py group_key, group_name, scope_groups: groups at each scope, repeat groups rgb.py StarDist's percentile normalization: rgb_percentile, apply_rgb (CSBDeep's float32 arithmetic), check_rgb + stain.py RGB-to-scalar stain presets: rgb_to_scalar, PRESETS (H&E first; optional + scikit-image extra `histology`; contract §3.4, N67) fitting.py fit, by channel ID parameters.py Parameters (apply, select, to_columns, from_columns with Normalization_Offset_ and Normalization_Scale_) @@ -113,7 +115,8 @@ range check. 7. **Parameters as numbers.** Parameters are numeric columns, never JSON; only the method's columns are written; settings are stored once per run (N21, N23). `from_columns(to_columns(p))` must equal `p` exactly. -8. **Dependencies.** Required: `numpy` only. No GPU, deep-learning, widget or file-format +8. **Dependencies.** Required: `numpy` only; the optional extra `histology` (scikit-image) is + imported lazily by `stain.py` alone (N67). No GPU, deep-learning, widget or file-format dependency beyond what the statistics file needs. ## 3. How we work diff --git a/CHANGELOG.md b/CHANGELOG.md index 22047e3..557fc10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,14 @@ Release notes and migration notes of `morphobricks-normalization`. Numbers are p issues of [CellMigrationLab/morphobricks-normalization](https://github.com/CellMigrationLab/morphobricks-normalization). The decisions are in [docs/plan.md](docs/plan.md) §12. +## Unreleased + +- **RGB-to-scalar stain presets** (issue #94, plan N67, contract §3.4): new `stain` module with + `rgb_to_scalar(pixels, *, axes, preset)` and the preset registry `PRESETS` (`he_eosin`, + `he_hematoxylin`, `luminance_inverted`, `lab_a`). 8-bit sRGB, `axes` `"YXC"` or `"CYX"`. + scikit-image is the new optional extra `histology`. `rgb.py`, the parameter columns and + `plane_histograms.csv` are unchanged. + ## v0.5.0 (2026-11-11) Breaking release. Stage 3 S3.2 (MorphoBricks-Spec roadmap; ADR 0003, plan N61–N66, PR #90): diff --git a/docs/contract.md b/docs/contract.md index b27d4d2..ab92fce 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -488,6 +488,50 @@ y_c = (x_c − offset) / scale for c = R, G, B --- +### 3.4 RGB-to-scalar stain presets (`stain.py`, issue #94, plan N67) + +A method family of its own, separate from §3.3 (which is unchanged): `rgb_to_scalar(pixels, *, +axes, preset)` returns `ScalarResult(image, exclusion, provenance)` for consumers that need one +channel from a stained image. normalization owns the conversion; consumers call it and never +reimplement it. No file reading, tiling, segmentation or WSI handling is done here. + +| `preset` (version 1) | `image` | `exclusion` | +|---|---|---| +| `he_eosin` | `separate_stains(rgb, hed_from_rgb)`, component 1 (eosin) | component 0 (haematoxylin) | +| `he_hematoxylin` | the same call, component 0 | `None` | +| `luminance_inverted` | `1 - rgb2gray(rgb)` (baseline) | `None` | +| `lab_a` | `rgb2lab(rgb)[..., 1]` (CIELAB a*, signed, not clipped; comparator) | `None` | + +Every preset has `stain` (`"he"` or `None`), `method`, `matrix`, `target_component`, +`exclusion_component`, `polarity = "bright"` (higher means more stain), `fixed_range = None` and +`version = 1`. `fixed_range` stays `None` until a stain-OD range is calibrated from a reference +set and recorded; it is not hard-coded from one image. No trichrome presets yet. + +**Inputs (refused, never repaired).** `axes` is exactly `"YXC"` or `"CYX"`, with no default; +the array is 3-D with exactly three samples (R, G, B) on the colour axis (an alpha channel, or +any other count, is refused with a message, never dropped); the dtype is `uint8` (8-bit sRGB +only; other dtypes are refused, 16-bit and ICC handling are later); an unknown preset is +refused with the known ones listed. The source array is never modified and the result is +deterministic; `YXC` and `CYX` of the same image give bitwise identical outputs. `image` and +`exclusion` are float32 `(Y, X)` and finite (`uint8` cannot hold NaN or infinity). + +**scikit-image behaviour accepted in P0.** `separate_stains` floors every pixel value at 1e-6, +treats 255 as white, divides by `-log(1e-6)` and sets negative component values to 0. The +`clipped_fraction` of the provenance is the fraction of pixels whose target component is +exactly 0 after the call (a pixel clipped from a negative value and one that was 0 cannot be told +apart); it is `None` for `luminance_inverted` and `lab_a`, which do not clip. A self-computed +optical-density path is later work. + +**Provenance** (read-only mapping): `preset`, `preset_version`, `stain`, `method`, `matrix`, +`target_component`, `exclusion_component`, `polarity`, `input_dtype`, `axes`, `colour_space` +(`"srgb"`), `skimage_version`, `clipped_fraction`. + +**Dependency.** scikit-image is the optional extra `histology` +(`pip install morphobricks-normalization[histology]`), imported only inside the conversion; +without it `rgb_to_scalar` raises `ImportError` naming the extra. `import +morphobricks_normalization` needs NumPy only. The module is not re-exported from the package +root; import it as `morphobricks_normalization.stain`. + ## 4. Scopes, groups and weights ### 4.1 Scopes diff --git a/docs/plan.md b/docs/plan.md index 1618a6c..a12db6e 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -218,7 +218,9 @@ calling brick. ## 9. Dependencies `numpy` only. No GPU or deep-learning dependency; StarDist and CSBDeep are not dependencies, -their formula is implemented here. +their formula is implemented here. The one exception is the optional extra `histology` +(`scikit-image>=0.22`) for the RGB-to-scalar stain presets of `stain.py` (N67); the core +imports and works without it. ## 10. Milestones @@ -534,6 +536,16 @@ Added with the project owner on 11 November 2026, for Stage 3 (S3.2, pull reques (`ChannelScopeIncompleteError`, the Spec's `A_SCOPE_INCOMPLETE`); the user may select the planes that have it with `partial_channels="planes_with_channel"`, stored with the settings and reported per fit (the notebook asks explicitly) (contract §1.5, §8.3). +61. **N67** (13 November 2026, issue #94) RGB-to-scalar stain presets are now in scope, H&E + first, and normalization owns the conversion; consumers (FiberFeatures #210) call it and + never reimplement it. This reverses the withdrawal of stain methods (N5, N17-N18, N24-N28) + for this method family only: new module `stain.py`, `rgb_to_scalar` with a code-defined + registry of frozen, versioned presets (P0: `he_eosin`, `he_hematoxylin`, + `luminance_inverted`, `lab_a`; contract §3.4). scikit-image is the optional extra + `histology` (N5 stays for the core). scikit-image's clipping in `separate_stains` is + accepted in P0 and reported as `clipped_fraction`; `fixed_range` stays `None` until + calibrated. The StarDist percentile normalization (`rgb.py`, §3.1) is unchanged. Later + (not in P0): trichrome and custom stain vectors (P1), RGB colour normalization (P2). Implementation choices of S3.2 where the Spec leaves the detail open (open questions 1-7 of pull request #90), each taken on the side the Spec states or the safer one, and listed in that diff --git a/pyproject.toml b/pyproject.toml index 0b1e568..6b321b4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,6 +24,7 @@ dependencies = [ ] [project.optional-dependencies] +histology = ["scikit-image>=0.22"] dev = [ "pytest>=7.4", "ruff>=0.6", @@ -122,3 +123,4 @@ baseline = "tools/code_rules_baseline.json" "morphobricks_normalization.bin_reachability" = "pure" "morphobricks_normalization.tail_moments" = "pure" "morphobricks_normalization.sums" = "pure" +"morphobricks_normalization.stain" = "pure" diff --git a/src/morphobricks_normalization/stain.py b/src/morphobricks_normalization/stain.py new file mode 100644 index 0000000..1d3e151 --- /dev/null +++ b/src/morphobricks_normalization/stain.py @@ -0,0 +1,235 @@ +"""RGB-to-scalar conversion with stain presets (contract §3.4, plan N67, issue #94). + +:func:`rgb_to_scalar` turns an 8-bit sRGB image into one float32 map, chosen by a named, +versioned preset (:data:`PRESETS`), for consumers that need a single channel from a stained +image (FiberFeatures first). It is a method family of its own: the RGB joint-percentile +normalization of ``rgb.py`` is unchanged and unrelated. + +Presets (P0): ``he_eosin``, ``he_hematoxylin`` (Ruifrok-Johnston H&E vectors through +``skimage.color.separate_stains`` and ``hed_from_rgb``), ``luminance_inverted`` (baseline) and +``lab_a`` (CIELAB a*, comparator). A new stain is a new preset record, not a new code path. + +Input contract. ``axes`` is exactly ``"YXC"`` or ``"CYX"`` and is never guessed; the array is +3-D with exactly three samples on the colour axis (R, G, B in that order; an alpha channel is +refused, never dropped); the dtype is ``uint8`` (sRGB; other dtypes are refused, 16-bit and ICC +handling are later work). Nothing is reordered, rescaled or repaired silently, and the input +array is never modified. ``uint8`` cannot hold NaN or infinity, so the output is finite. + +Behaviour inherited from scikit-image, accepted in P0 and recorded in the provenance: +``separate_stains`` floors every pixel value at 1e-6, treats 255 as white, divides by +``-log(1e-6)`` and **sets negative component values to 0**. The provenance reports +``clipped_fraction``: the fraction of the pixels whose target component is exactly 0 after the +call (a pixel that was negative or truly 0 cannot be told apart). ``lab_a`` is signed and not +clipped, and ``luminance_inverted`` is not clipped either; both report ``None``. + +``fixed_range`` of a preset is ``None`` for now: the stain-OD range a consumer needs for fixed +intensity scaling must be calibrated from a reference set and recorded, not hard-coded from one +image (issue #94). + +scikit-image is the optional ``histology`` extra, imported only inside the conversion; without +it :func:`rgb_to_scalar` raises :class:`ImportError`, and importing this module needs NumPy only. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from dataclasses import dataclass +from types import MappingProxyType, ModuleType +from typing import Any + +import numpy as np + +AXES = ("YXC", "CYX") +SAMPLE_COUNT = 3 +COLOUR_SPACE = "srgb" +INSTALL_HINT = "pip install morphobricks-normalization[histology]" + + +@dataclass(frozen=True) +class StainPreset: + """A frozen, versioned recipe for one scalar map from an RGB image. + + Attributes + ---------- + id + Preset identifier, the key in :data:`PRESETS`. + version + Version of the recipe; any change of numbers is a new version. + stain + Stain declared in the MorphoBricks Spec (``"he"``), or ``None`` for a colour projection. + method + ``fixed_deconvolution``, ``luminance`` or ``cielab``. + matrix + Name of the scikit-image stain matrix (``hed_from_rgb``), or ``None``. + target_component + Index of the component (deconvolution) or CIELAB channel returned, or ``None``. + exclusion_component + Component returned as the optional exclusion map (haematoxylin), or ``None``. + polarity + ``"bright"``: a higher value means more stain. + fixed_range + Stain range for fixed intensity scaling; ``None`` until calibrated from a reference set. + """ + + id: str + version: int = 1 + stain: str | None = None + method: str = "fixed_deconvolution" + matrix: str | None = None + target_component: int | None = None + exclusion_component: int | None = None + polarity: str = "bright" + fixed_range: tuple[float, float] | None = None + + +@dataclass(frozen=True) +class ScalarResult: + """The scalar map of :func:`rgb_to_scalar`. + + Attributes + ---------- + image + float32 ``(Y, X)``, finite. + exclusion + float32 ``(Y, X)`` (haematoxylin for ``he_eosin``), or ``None`` if the preset has none. + provenance + Read-only mapping that records how the map was made. + """ + + image: np.ndarray + exclusion: np.ndarray | None + provenance: Mapping[str, Any] + + +PRESETS: Mapping[str, StainPreset] = MappingProxyType( + { + preset.id: preset + for preset in ( + StainPreset( + "he_eosin", + stain="he", + matrix="hed_from_rgb", + target_component=1, + exclusion_component=0, + ), + StainPreset("he_hematoxylin", stain="he", matrix="hed_from_rgb", target_component=0), + StainPreset("luminance_inverted", method="luminance"), + StainPreset("lab_a", method="cielab", target_component=1), + ) + } +) +"""The four P0 presets, read-only. No trichrome presets yet (issue #94, P1).""" + + +def _to_yxc(pixels: np.ndarray, axes: str) -> np.ndarray: + """Return a C-contiguous ``(Y, X, 3)`` copy after checking the input contract.""" + if not isinstance(pixels, np.ndarray): + raise TypeError(f"pixels must be a NumPy array; got {type(pixels).__name__}.") + if not isinstance(axes, str) or axes not in AXES: + raise ValueError(f"axes must be exactly 'YXC' or 'CYX', never guessed; got {axes!r}.") + if pixels.ndim != 3: + raise ValueError(f"Expected a 3-D array with axes {axes}; got {pixels.ndim}-D.") + colour_axis = axes.index("C") + if pixels.shape[colour_axis] != SAMPLE_COUNT: + extra = ( + " An alpha channel is refused, not dropped." if pixels.shape[colour_axis] == 4 else "" + ) + raise ValueError( + f"Expected exactly 3 samples (R, G, B) on the colour axis of {axes}; " + f"got {pixels.shape[colour_axis]}.{extra}" + ) + if pixels.dtype != np.uint8: + raise TypeError( + f"Only 8-bit sRGB (uint8) is supported; got {pixels.dtype}. " + "Nothing is rescaled or converted silently." + ) + return np.ascontiguousarray(np.moveaxis(pixels, colour_axis, 2)) + + +def _skimage_color() -> tuple[ModuleType, ModuleType]: + """Import ``skimage`` and ``skimage.color`` lazily; name the extra when missing.""" + try: + import skimage + from skimage import color + except ImportError as error: + raise ImportError(f"RGB-to-scalar presets need scikit-image: {INSTALL_HINT}") from error + return skimage, color + + +def _convert(rgb: np.ndarray, preset: StainPreset) -> tuple[np.ndarray, np.ndarray | None, str]: + """Return ``(image, exclusion, skimage version)`` as float32, from a ``(Y, X, 3)`` uint8.""" + skimage, color = _skimage_color() + exclusion = None + if preset.method == "fixed_deconvolution": + matrix = getattr(color, preset.matrix) + components = color.separate_stains(rgb, matrix) + image = components[..., preset.target_component] + if preset.exclusion_component is not None: + exclusion = components[..., preset.exclusion_component] + elif preset.method == "luminance": + image = 1.0 - color.rgb2gray(rgb) + elif preset.method == "cielab": + image = color.rgb2lab(rgb)[..., preset.target_component] + else: + raise ValueError(f"Preset {preset.id!r} has an unknown method {preset.method!r}.") + as_float32 = np.ascontiguousarray(image, dtype=np.float32) + excluded = None if exclusion is None else np.ascontiguousarray(exclusion, dtype=np.float32) + return as_float32, excluded, skimage.__version__ + + +def rgb_to_scalar(pixels: np.ndarray, *, axes: str, preset: str) -> ScalarResult: + """Convert an 8-bit sRGB image to one float32 map with a stain preset. + + Parameters + ---------- + pixels + ``uint8`` array, 3-D, with exactly three samples (R, G, B) on the colour axis. Never + modified. + axes + ``"YXC"`` or ``"CYX"``; required, never guessed. + preset + A key of :data:`PRESETS`. + + Returns + ------- + ScalarResult + ``image`` float32 ``(Y, X)``; ``exclusion`` (``he_eosin`` only) or ``None``; + ``provenance`` with preset id and version, stain, method, matrix, target and exclusion + component, polarity, input dtype, axes, colour space, scikit-image version and + ``clipped_fraction`` (fraction of target-component pixels exactly 0, ``None`` where the + preset does not clip). + + Raises + ------ + ValueError + Unknown preset, axes other than ``"YXC"``/``"CYX"``, not 3-D, or not three samples + (an alpha channel included). + TypeError + ``pixels`` is not a NumPy array or its dtype is not ``uint8``. + ImportError + scikit-image is not installed (``pip install morphobricks-normalization[histology]``). + """ + if preset not in PRESETS: + raise ValueError(f"Unknown preset {preset!r}; known presets: {', '.join(PRESETS)}.") + recipe = PRESETS[preset] + rgb = _to_yxc(pixels, axes) + image, exclusion, skimage_version = _convert(rgb, recipe) + clipped = float(np.mean(image == 0)) if recipe.method == "fixed_deconvolution" else None + provenance = MappingProxyType( + { + "preset": recipe.id, + "preset_version": recipe.version, + "stain": recipe.stain, + "method": recipe.method, + "matrix": recipe.matrix, + "target_component": recipe.target_component, + "exclusion_component": recipe.exclusion_component, + "polarity": recipe.polarity, + "input_dtype": str(pixels.dtype), + "axes": axes, + "colour_space": COLOUR_SPACE, + "skimage_version": skimage_version, + "clipped_fraction": clipped, + } + ) + return ScalarResult(image=image, exclusion=exclusion, provenance=provenance) diff --git a/tests/test_stain.py b/tests/test_stain.py new file mode 100644 index 0000000..6c369a1 --- /dev/null +++ b/tests/test_stain.py @@ -0,0 +1,236 @@ +import subprocess +import sys +import textwrap + +import numpy as np +import pytest + +from morphobricks_normalization.stain import PRESETS, ScalarResult, StainPreset, rgb_to_scalar + +EOSIN = (230, 130, 190) +HEMATOXYLIN = (80, 60, 150) + + +def random_image(seed=0, shape=(7, 9)): + return np.random.default_rng(seed).integers(0, 256, (*shape, 3), dtype=np.uint8) + + +def patches(left, right): + """A 4x8 YXC image: left colour in the first four columns, right colour in the rest.""" + image = np.empty((4, 8, 3), dtype=np.uint8) + image[:, :4] = left + image[:, 4:] = right + return image + + +def test_the_module_imports_and_fails_clearly_without_skimage(): + code = textwrap.dedent( + """ + import sys + sys.modules["skimage"] = None + sys.modules["skimage.color"] = None + import numpy as np + import morphobricks_normalization + from morphobricks_normalization.stain import rgb_to_scalar + try: + rgb_to_scalar(np.zeros((2, 2, 3), np.uint8), axes="YXC", preset="he_eosin") + except ImportError as error: + print(error) + else: + raise SystemExit("no ImportError") + """ + ) + result = subprocess.run( + [sys.executable, "-c", code], capture_output=True, text=True, check=False + ) + assert result.returncode == 0, result.stderr + assert "pip install morphobricks-normalization[histology]" in result.stdout + + +def test_the_registry_has_exactly_the_four_p0_presets(): + assert list(PRESETS) == ["he_eosin", "he_hematoxylin", "luminance_inverted", "lab_a"] + assert all( + p.version == 1 and p.polarity == "bright" and p.fixed_range is None + for p in PRESETS.values() + ) + + +def test_presets_are_immutable(): + with pytest.raises(TypeError): + PRESETS["new"] = StainPreset("new") # type: ignore[index] + with pytest.raises(AttributeError): + PRESETS["he_eosin"].version = 2 # type: ignore[misc] + + +@pytest.mark.parametrize("preset", list(PRESETS)) +def test_yxc_and_cyx_give_bitwise_identical_outputs(preset): + pytest.importorskip("skimage") + yxc = random_image() + cyx = np.moveaxis(yxc, 2, 0) + a = rgb_to_scalar(yxc, axes="YXC", preset=preset) + b = rgb_to_scalar(cyx, axes="CYX", preset=preset) + assert a.image.tobytes() == b.image.tobytes() + if a.exclusion is not None: + assert a.exclusion.tobytes() == b.exclusion.tobytes() + assert a.provenance["clipped_fraction"] == b.provenance["clipped_fraction"] + + +@pytest.mark.parametrize("axes", ["yxc", "XYC", "YXCZ", "", "CXY", None, ("Y", "X", "C")]) +def test_wrong_axes_are_refused(axes): + with pytest.raises(ValueError, match="axes"): + rgb_to_scalar(random_image(), axes=axes, preset="he_eosin") + + +def test_axes_are_required(): + with pytest.raises(TypeError): + rgb_to_scalar(random_image(), preset="he_eosin") # type: ignore[call-arg] + + +@pytest.mark.parametrize( + ("array", "axes", "message"), + [ + (np.zeros((5, 5), np.uint8), "YXC", "3-D"), + (np.zeros((5, 5, 4), np.uint8), "YXC", "alpha"), + (np.zeros((4, 5, 5), np.uint8), "CYX", "alpha"), + (np.zeros((5, 5, 2), np.uint8), "YXC", "exactly 3"), + (np.zeros((5, 5, 3), np.uint8), "CYX", "exactly 3"), + (np.zeros((2, 5, 5, 3), np.uint8), "YXC", "3-D"), + ], +) +def test_wrong_shapes_are_refused(array, axes, message): + with pytest.raises(ValueError, match=message): + rgb_to_scalar(array, axes=axes, preset="he_eosin") + + +@pytest.mark.parametrize("dtype", ["float32", "float64", "uint16", "int16", "bool"]) +def test_other_dtypes_are_refused(dtype): + with pytest.raises(TypeError, match="uint8"): + rgb_to_scalar(np.zeros((4, 4, 3), dtype), axes="YXC", preset="he_eosin") + + +def test_a_non_array_is_refused(): + with pytest.raises(TypeError, match="NumPy array"): + rgb_to_scalar([[[0, 0, 0]]], axes="YXC", preset="he_eosin") # type: ignore[arg-type] + + +def test_unknown_preset_is_refused_listing_the_known_ones(): + with pytest.raises(ValueError, match="he_eosin.*he_hematoxylin.*luminance_inverted.*lab_a"): + rgb_to_scalar(random_image(), axes="YXC", preset="trichrome_blue_bpx") + + +@pytest.mark.parametrize("preset", list(PRESETS)) +def test_input_is_unchanged_and_output_is_float32_yx_finite(preset): + pytest.importorskip("skimage") + image = random_image(shape=(6, 11)) + before = image.copy() + result = rgb_to_scalar(image, axes="YXC", preset=preset) + assert isinstance(result, ScalarResult) + assert np.array_equal(image, before) and image.dtype == np.uint8 + assert result.image.dtype == np.float32 and result.image.shape == (6, 11) + assert np.isfinite(result.image).all() + if result.exclusion is not None: + assert result.exclusion.dtype == np.float32 and result.exclusion.shape == (6, 11) + assert np.isfinite(result.exclusion).all() + assert not np.shares_memory(result.exclusion, image) + assert not np.shares_memory(result.image, image) + + +def test_cyx_input_is_unchanged(): + pytest.importorskip("skimage") + cyx = np.moveaxis(random_image(), 2, 0).copy() + before = cyx.copy() + rgb_to_scalar(cyx, axes="CYX", preset="he_eosin") + assert np.array_equal(cyx, before) + + +def test_eosin_and_hematoxylin_presets_separate_the_two_stains(): + pytest.importorskip("skimage") + image = patches(EOSIN, HEMATOXYLIN) # eosin on the left, haematoxylin on the right + eosin = rgb_to_scalar(image, axes="YXC", preset="he_eosin").image + hema = rgb_to_scalar(image, axes="YXC", preset="he_hematoxylin").image + assert eosin[:, :4].mean() > eosin[:, 4:].mean() + assert hema[:, 4:].mean() > hema[:, :4].mean() + + +def test_equal_luminance_colours_are_separated_by_lab_a_not_by_luminance(): + pytest.importorskip("skimage") + # rgb2gray = 0.2125 R + 0.7154 G + 0.0721 B: both colours give 121.3 / 255 within 0.001. + red, green = (200, 100, 100), (100, 130, 100) + image = patches(red, green) + luminance = rgb_to_scalar(image, axes="YXC", preset="luminance_inverted").image + lab_a = rgb_to_scalar(image, axes="YXC", preset="lab_a").image + # 0.005 is about one uint8 grey step (1/255 = 0.0039) of the weighted sum. + assert luminance.max() - luminance.min() < 0.005 + assert lab_a[:, :4].mean() - lab_a[:, 4:].mean() > 40 # CIELAB a* units, red vs green + + +def test_lab_a_is_signed_and_not_clipped(): + pytest.importorskip("skimage") + result = rgb_to_scalar(patches(EOSIN, (60, 160, 60)), axes="YXC", preset="lab_a") + assert result.image.min() < 0 < result.image.max() + assert result.provenance["clipped_fraction"] is None + + +def test_he_eosin_exclusion_is_the_hematoxylin_image_and_others_have_none(): + pytest.importorskip("skimage") + image = random_image() + eosin = rgb_to_scalar(image, axes="YXC", preset="he_eosin") + hema = rgb_to_scalar(image, axes="YXC", preset="he_hematoxylin") + assert eosin.exclusion is not None + assert eosin.exclusion.tobytes() == hema.image.tobytes() + for preset in ("he_hematoxylin", "luminance_inverted", "lab_a"): + assert rgb_to_scalar(image, axes="YXC", preset=preset).exclusion is None + + +def test_provenance_contents_and_clipped_fraction(): + skimage = pytest.importorskip("skimage") + image = patches(EOSIN, HEMATOXYLIN) + result = rgb_to_scalar(image, axes="YXC", preset="he_eosin") + expected = { + "preset": "he_eosin", + "preset_version": 1, + "stain": "he", + "method": "fixed_deconvolution", + "matrix": "hed_from_rgb", + "target_component": 1, + "exclusion_component": 0, + "polarity": "bright", + "input_dtype": "uint8", + "axes": "YXC", + "colour_space": "srgb", + "skimage_version": skimage.__version__, + "clipped_fraction": float(np.mean(result.image == 0)), + } + assert dict(result.provenance) == expected + assert 0.0 <= result.provenance["clipped_fraction"] <= 1.0 + lab = rgb_to_scalar(image, axes="YXC", preset="lab_a").provenance + assert lab["stain"] is None and lab["matrix"] is None and lab["method"] == "cielab" + assert ( + rgb_to_scalar(image, axes="YXC", preset="luminance_inverted").provenance["method"] + == "luminance" + ) + + +def test_the_clipped_fraction_counts_zeros_of_the_target_component(): + pytest.importorskip("skimage") + # Left half pure white: every component is 0 there; right half the haematoxylin colour. + image = patches((255, 255, 255), HEMATOXYLIN) + result = rgb_to_scalar(image, axes="YXC", preset="he_hematoxylin") + assert result.provenance["clipped_fraction"] == pytest.approx(0.5) + + +def test_provenance_is_immutable(): + pytest.importorskip("skimage") + provenance = rgb_to_scalar(random_image(), axes="YXC", preset="he_eosin").provenance + with pytest.raises(TypeError): + provenance["preset"] = "other" # type: ignore[index] + + +@pytest.mark.parametrize("preset", list(PRESETS)) +def test_the_conversion_is_deterministic(preset): + pytest.importorskip("skimage") + image = random_image(3) + a = rgb_to_scalar(image, axes="YXC", preset=preset) + b = rgb_to_scalar(image, axes="YXC", preset=preset) + assert a.image.tobytes() == b.image.tobytes() + assert dict(a.provenance) == dict(b.provenance) From b152472898e17f2086e944563800e320b527547f Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 07:19:03 +0000 Subject: [PATCH 02/12] Run the stain tests in CI and copy the input explicitly The main CI matrix installs only .[dev], so tests/test_stain.py was skipped there; a dedicated histology job installs the extra and runs it. _to_yxc returned a view for contiguous YXC input (np.ascontiguousarray does not copy); it now always copies, as its docstring says. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- .github/workflows/ci.yml | 12 ++++++++++++ src/morphobricks_normalization/stain.py | 7 +++++-- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 790011e..25808bc 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -58,3 +58,15 @@ jobs: python-version: "3.12" - run: python -m pip install -e ".[dev]" csbdeep - run: python -m pytest tests/test_rgb.py -rs + + # The RGB-to-scalar stain presets (tests/test_stain.py) need the optional extra + # `histology` (scikit-image); the main matrix runs without it, so these tests run here. + histology: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + - run: python -m pip install -e ".[dev,histology]" + - run: python -m pytest tests/test_stain.py -rs diff --git a/src/morphobricks_normalization/stain.py b/src/morphobricks_normalization/stain.py index 1d3e151..72077f2 100644 --- a/src/morphobricks_normalization/stain.py +++ b/src/morphobricks_normalization/stain.py @@ -122,7 +122,10 @@ class ScalarResult: def _to_yxc(pixels: np.ndarray, axes: str) -> np.ndarray: - """Return a C-contiguous ``(Y, X, 3)`` copy after checking the input contract.""" + """Return a new C-contiguous ``(Y, X, 3)`` array after checking the input contract. + + Always a copy, never a view of ``pixels``, so no later step can modify the caller's array. + """ if not isinstance(pixels, np.ndarray): raise TypeError(f"pixels must be a NumPy array; got {type(pixels).__name__}.") if not isinstance(axes, str) or axes not in AXES: @@ -143,7 +146,7 @@ def _to_yxc(pixels: np.ndarray, axes: str) -> np.ndarray: f"Only 8-bit sRGB (uint8) is supported; got {pixels.dtype}. " "Nothing is rescaled or converted silently." ) - return np.ascontiguousarray(np.moveaxis(pixels, colour_axis, 2)) + return np.array(np.moveaxis(pixels, colour_axis, 2), order="C", copy=True) def _skimage_color() -> tuple[ModuleType, ModuleType]: From beadf6238b5a745c3b5852f9894a7e201d2fc191 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 08:14:27 +0000 Subject: [PATCH 03/12] Add the preset catalogue: to_scalar, custom_od, passthrough and hints to_scalar(pixels, *, axes, preset, white, stain_vectors, polarizer) is the general entry point; rgb_to_scalar delegates to it and the four P0 presets give bitwise identical outputs. New presets: fluorescence and shg (passthrough), psr_brightfield, he_custom and trichrome_custom (custom_od, optical density against a declared white in NumPy), psr_polarized (rgb2gray, declared polarizer), trichrome_blue (bpx_from_rgb) and trichrome_green (minus CIELAB a*). Each preset carries a modality and consumer hints. Helpers white_from_region and stain_vectors_from_regions derive the white and the stain vectors from labelled regions. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- pyproject.toml | 2 + .../optical_density.py | 293 ++++++++++ .../scalar_inputs.py | 94 ++++ src/morphobricks_normalization/stain.py | 514 ++++++++++++++---- tests/stain_simulation.py | 55 ++ tests/test_optical_density.py | 227 ++++++++ tests/test_stain.py | 65 ++- tests/test_stain_catalogue.py | 487 +++++++++++++++++ 8 files changed, 1624 insertions(+), 113 deletions(-) create mode 100644 src/morphobricks_normalization/optical_density.py create mode 100644 src/morphobricks_normalization/scalar_inputs.py create mode 100644 tests/stain_simulation.py create mode 100644 tests/test_optical_density.py create mode 100644 tests/test_stain_catalogue.py diff --git a/pyproject.toml b/pyproject.toml index 6b321b4..10dff82 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -124,3 +124,5 @@ baseline = "tools/code_rules_baseline.json" "morphobricks_normalization.tail_moments" = "pure" "morphobricks_normalization.sums" = "pure" "morphobricks_normalization.stain" = "pure" +"morphobricks_normalization.optical_density" = "pure" +"morphobricks_normalization.scalar_inputs" = "pure" diff --git a/src/morphobricks_normalization/optical_density.py b/src/morphobricks_normalization/optical_density.py new file mode 100644 index 0000000..4ff5b18 --- /dev/null +++ b/src/morphobricks_normalization/optical_density.py @@ -0,0 +1,293 @@ +"""Optical density and stain unmixing against a declared white (contract §3.4, plan N68). + +The in-house ``custom_od`` path of :mod:`morphobricks_normalization.stain`, in NumPy only (no +scikit-image), so that slide-specific stain vectors (the TWOMBLI approach) work without the +``histology`` extra. Every function takes a ``(Y, X, 3)`` RGB array already checked by the +caller and returns new arrays; nothing is modified in place. + +Beer-Lambert law per channel ``c``: ``I_c = white_c * 10 ** (-OD_c)`` and ``OD = C @ M``, where +the rows of ``M`` are the unit optical-density vectors of the stains and ``C`` holds one +concentration per stain. Hence: + +- ``OD_c = -log10(max(I_c, 1) / white_c)``. A pixel value of 0 has no finite optical density; + it is floored at 1 (the smallest non-zero 8-bit value) and the fraction of pixels with a + floored channel is reported, never hidden. +- ``C = OD @ inv(M)``. Negative concentrations are kept, not clipped: they mean the pixel does + not fit the declared stains (noise, a third colour, a white that is too dark), and the caller + reports their fraction. Pixels brighter than the white have a negative optical density in + that channel and are reported too. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from dataclasses import dataclass +from types import MappingProxyType +from typing import Any + +import numpy as np + +CHANNELS = 3 +WHITE_MAX = 255.0 +"""Upper limit of a declared white: the largest 8-bit value.""" +FLOOR = 1.0 +"""Pixel values below 1 (that is, 0) are raised to 1 before the logarithm.""" +RESIDUAL = "residual" +"""Name of the third stain vector, declared or completed as a cross product.""" +MIN_VECTOR_NORM = 1e-6 +"""A stain vector shorter than this has no direction; optical densities are of order 0.1-1.""" +MIN_ANGLE_DEGREES = 5.0 +"""Two stain vectors closer than this (in either direction) are refused as collinear. + +For two unit vectors at angle ``theta`` completed by their cross product the condition number +of ``M`` is ``cot(theta / 2)``: 5 degrees gives 22.9, so this rule is the binding one for two +vectors and :data:`MAX_CONDITION` the one for three nearly coplanar vectors. +""" +MAX_CONDITION = 50.0 +"""Largest accepted 2-norm condition number of the unit stain matrix. + +The relative error of the concentrations is at most ``cond(M)`` times the relative error of +the optical densities. One 8-bit grey level at a mid-range value (``I`` about 80, ``OD`` about +0.5) is an optical-density error of about 1 %, so 50 bounds the worst-case amplification at +about 50 %: beyond it the unmixed maps are dominated by quantisation noise. Ruifrok-Johnston +H&E-DAB has 3.8 and scikit-image's Methyl Blue/Ponceau Fuchsin matrix 2.0. +""" +MIN_REGION_OD = 0.01 +"""Smallest optical-density norm of a region's mean for it to define a stain vector. + +0.01 is about 2 % absorption, five grey levels below a white of 240: a region below it is +background, not stain, and its direction would be noise. +""" + + +@dataclass(frozen=True) +class StainMatrix: + """Unit stain vectors as rows, in the order target, exclusion, residual. + + Attributes + ---------- + names + The three stain names, row by row. + rows + The unit vectors, as nested tuples (read-only). + condition + 2-norm condition number of the matrix. + completed + Whether the residual row is the normalised cross product of the first two. + """ + + names: tuple[str, str, str] + rows: tuple[tuple[float, float, float], ...] + condition: float + completed: bool + + def array(self) -> np.ndarray: + """Return the matrix as a new float64 ``(3, 3)`` array.""" + return np.array(self.rows, dtype=np.float64) + + +def check_white(white: object) -> np.ndarray: + """Return a declared white as a float64 ``(3,)`` array after checking it. + + Parameters + ---------- + white + Background RGB of the 8-bit image: three finite numbers in ``(0, 255]``. + + Returns + ------- + numpy.ndarray + float64 ``(3,)``. + + Raises + ------ + ValueError + Not three finite numbers in ``(0, 255]``. + """ + values = np.asarray(white) + if values.dtype.kind not in "iuf" or values.shape != (CHANNELS,): + raise ValueError( + "white must be three numbers (the background R, G, B of the 8-bit image); " + f"got {white!r}." + ) + values = values.astype(np.float64) + if not np.isfinite(values).all() or (values <= 0).any() or (values > WHITE_MAX).any(): + raise ValueError( + f"Each value of white must be finite and in (0, 255]; got {tuple(values.tolist())}. " + "Measure it on background, for example with white_from_region." + ) + return values + + +def _unit_vector(name: str, values: object) -> np.ndarray: + """Return the stain vector ``name`` normalised to unit length, or refuse it.""" + vector = np.asarray(values) + if vector.dtype.kind not in "iuf" or vector.shape != (CHANNELS,): + raise ValueError( + f"Stain vector {name!r} must be three numbers (R, G, B optical densities); " + f"got {values!r}." + ) + vector = vector.astype(np.float64) + if not np.isfinite(vector).all(): + raise ValueError(f"Stain vector {name!r} is not finite: {tuple(vector.tolist())}.") + norm = float(np.linalg.norm(vector)) + if norm < MIN_VECTOR_NORM: + raise ValueError( + f"Stain vector {name!r} has (almost) zero length ({norm:.3g}); it has no direction." + ) + return vector / norm + + +def _check_names(vectors: Mapping[str, Any], target: str, exclusion: str) -> None: + """Refuse a mapping without the target and exclusion stains or with other names.""" + if not isinstance(vectors, Mapping): + raise TypeError(f"stain_vectors must be a mapping name -> 3 values; got {vectors!r}.") + allowed = {target, exclusion, RESIDUAL} + missing = [name for name in (target, exclusion) if name not in vectors] + extra = sorted(str(name) for name in vectors if name not in allowed) + if missing or extra: + raise ValueError( + f"stain_vectors must name {target!r} and {exclusion!r}, and optionally " + f"{RESIDUAL!r}; missing {missing}, not allowed {extra}." + ) + + +def _angle_degrees(a: np.ndarray, b: np.ndarray) -> float: + """Return the angle between two unit vectors, ignoring their sign (0 to 90 degrees).""" + return float(np.degrees(np.arccos(min(abs(float(a @ b)), 1.0)))) + + +def stain_matrix(vectors: Mapping[str, Any], target: str, exclusion: str) -> StainMatrix: + """Build and check the unit stain matrix of a ``custom_od`` preset. + + Parameters + ---------- + vectors + Mapping stain name -> three optical-density components (R, G, B), any positive scale. + Exactly ``target`` and ``exclusion``, plus optionally ``"residual"``; without it the + third row is the normalised cross product of the first two. + target, exclusion + The preset's stain names; their rows come first and second. + + Returns + ------- + StainMatrix + + Raises + ------ + ValueError + Missing or extra names; a vector that is not three finite numbers or has zero length; + two vectors within :data:`MIN_ANGLE_DEGREES`; a condition number above + :data:`MAX_CONDITION`. + """ + _check_names(vectors, target, exclusion) + names = (target, exclusion, RESIDUAL) + rows = [_unit_vector(name, vectors[name]) for name in names[:2]] + completed = RESIDUAL not in vectors + if not completed: + rows.append(_unit_vector(RESIDUAL, vectors[RESIDUAL])) + for i in range(len(rows)): + for j in range(i + 1, len(rows)): + angle = _angle_degrees(rows[i], rows[j]) + if angle < MIN_ANGLE_DEGREES: + raise ValueError( + f"Stain vectors {names[i]!r} and {names[j]!r} are {angle:.2f} degrees apart " + f"(minimum {MIN_ANGLE_DEGREES}); collinear stains cannot be separated." + ) + if completed: + cross = np.cross(rows[0], rows[1]) + rows.append(cross / np.linalg.norm(cross)) + matrix = np.array(rows) + condition = float(np.linalg.cond(matrix)) + if not np.isfinite(condition) or condition > MAX_CONDITION: + raise ValueError( + f"The stain matrix is ill-conditioned (condition number {condition:.3g}, maximum " + f"{MAX_CONDITION}); the vectors are too close to a plane to be separated." + ) + as_tuples = tuple(tuple(float(v) for v in row) for row in matrix) + return StainMatrix(names, as_tuples, condition, completed) + + +def optical_density(rgb: np.ndarray, white: np.ndarray) -> tuple[np.ndarray, float]: + """Return ``(OD, floored_fraction)`` of a ``(Y, X, 3)`` image against ``white``. + + ``OD = -log10(max(I, 1) / white)`` per channel, float64. ``floored_fraction`` is the + fraction of pixels with at least one channel below 1 (raised to 1). + """ + values = rgb.astype(np.float64) + floored = values < FLOOR + od = -np.log10(np.maximum(values, FLOOR) / white) + return od, float(np.mean(floored.any(axis=-1))) + + +def above_white_fraction(rgb: np.ndarray, white: np.ndarray) -> float: + """Return the fraction of pixels brighter than ``white`` in at least one channel.""" + return float(np.mean((rgb.astype(np.float64) > white).any(axis=-1))) + + +def concentrations(od: np.ndarray, matrix: StainMatrix) -> np.ndarray: + """Return ``OD @ inv(M)``: one float64 concentration per stain, signed, never clipped.""" + return od @ np.linalg.inv(matrix.array()) + + +def _check_region(region: object, shape: tuple[int, int], name: str) -> np.ndarray: + """Return a boolean ``(Y, X)`` region, or refuse a wrong type, shape or empty region.""" + if not isinstance(region, np.ndarray) or region.dtype != np.bool_: + raise TypeError(f"Region {name!r} must be a boolean NumPy array (Y, X).") + if region.shape != shape: + raise ValueError(f"Region {name!r} has shape {region.shape}; the image is {shape} (Y, X).") + if not region.any(): + raise ValueError(f"Region {name!r} is empty; select at least one pixel.") + return region + + +def region_white(rgb: np.ndarray, region: object) -> tuple[float, float, float]: + """Return the median R, G, B of a ``(Y, X, 3)`` image over a background region. + + Raises + ------ + ValueError + Empty region, wrong shape, or a median of 0 or 255 (no light, or saturated: the + true white is then unknown). + TypeError + The region is not a boolean array. + """ + mask = _check_region(region, rgb.shape[:2], "background") + median = np.median(rgb[mask].astype(np.float64), axis=0) + if (median <= 0).any() or (median >= WHITE_MAX).any(): + raise ValueError( + f"The background median {tuple(median.tolist())} has a channel at 0 or 255 " + "(dark or saturated); the white cannot be measured from it. Lower the exposure " + "or choose another background region." + ) + return tuple(float(v) for v in median) + + +def region_vectors( + rgb: np.ndarray, white: np.ndarray, regions: Mapping[str, Any] +) -> Mapping[str, tuple[float, float, float]]: + """Return the unit mean optical-density vector of each labelled region. + + Raises + ------ + ValueError + No regions, an empty or wrongly shaped region, or a region whose mean optical density + is below :data:`MIN_REGION_OD` (background, not stain). + TypeError + ``regions`` is not a mapping, or a region is not a boolean array. + """ + if not isinstance(regions, Mapping) or not regions: + raise TypeError("regions must be a non-empty mapping stain name -> boolean (Y, X) array.") + od, _ = optical_density(rgb, white) + vectors = {} + for name, region in regions.items(): + mask = _check_region(region, rgb.shape[:2], name) + mean = od[mask].mean(axis=0) + norm = float(np.linalg.norm(mean)) + if norm < MIN_REGION_OD: + raise ValueError( + f"Region {name!r} has a mean optical density of {norm:.4f} (minimum " + f"{MIN_REGION_OD}); it looks like background, not stain." + ) + vectors[name] = tuple(float(v) for v in mean / norm) + return MappingProxyType(vectors) diff --git a/src/morphobricks_normalization/scalar_inputs.py b/src/morphobricks_normalization/scalar_inputs.py new file mode 100644 index 0000000..004a62d --- /dev/null +++ b/src/morphobricks_normalization/scalar_inputs.py @@ -0,0 +1,94 @@ +"""Input contracts of the image-to-scalar presets (contract §3.4, plan N67, N68). + +:func:`rgb_yxc` checks an 8-bit sRGB image with declared axes and returns a ``(Y, X, 3)`` +copy; :func:`plane_float32` checks a single ``"YX"`` plane for the ``passthrough`` presets and +returns a float32 copy. Nothing is reordered, rescaled or repaired silently, and the caller's +array is never returned or modified. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + +AXES = ("YXC", "CYX") +PLANE_AXES = "YX" +SAMPLE_COUNT = 3 + + +def rgb_yxc(pixels: np.ndarray, axes: str) -> np.ndarray: + """Return a new C-contiguous ``(Y, X, 3)`` array after checking the RGB input contract. + + Always a copy, never a view of ``pixels``, so no later step can modify the caller's array. + """ + if not isinstance(pixels, np.ndarray): + raise TypeError(f"pixels must be a NumPy array; got {type(pixels).__name__}.") + if not isinstance(axes, str) or axes not in AXES: + raise ValueError(f"axes must be exactly 'YXC' or 'CYX', never guessed; got {axes!r}.") + if pixels.ndim != 3: + raise ValueError(f"Expected a 3-D array with axes {axes}; got {pixels.ndim}-D.") + colour_axis = axes.index("C") + if pixels.shape[colour_axis] != SAMPLE_COUNT: + extra = ( + " An alpha channel is refused, not dropped." if pixels.shape[colour_axis] == 4 else "" + ) + raise ValueError( + f"Expected exactly 3 samples (R, G, B) on the colour axis of {axes}; " + f"got {pixels.shape[colour_axis]}.{extra}" + ) + if pixels.dtype != np.uint8: + raise TypeError( + f"Only 8-bit sRGB (uint8) is supported; got {pixels.dtype}. " + "Nothing is rescaled or converted silently." + ) + return np.array(np.moveaxis(pixels, colour_axis, 2), order="C", copy=True) + + +def plane_float32(pixels: np.ndarray, axes: str) -> tuple[np.ndarray, dict[str, Any]]: + """Return a float32 copy of a finite 2-D ``"YX"`` plane and whether the cast was exact. + + Raises + ------ + ValueError + axes other than ``"YX"``, not 2-D, NaN or infinity, or values beyond float32. + TypeError + Not a NumPy array, or a dtype that is not integer or float (bool, complex). + """ + if not isinstance(pixels, np.ndarray): + raise TypeError(f"pixels must be a NumPy array; got {type(pixels).__name__}.") + if axes != PLANE_AXES: + raise ValueError(f"axes must be exactly 'YX' for a single plane; got {axes!r}.") + if pixels.ndim != 2: + raise ValueError(f"Expected a 2-D array with axes YX; got {pixels.ndim}-D.") + if pixels.dtype.kind not in "iuf": + raise TypeError( + f"Expected a real numeric dtype (integer or float); got {pixels.dtype}. " + "bool and complex planes are refused." + ) + if pixels.dtype.kind == "f" and not np.isfinite(pixels).all(): + raise ValueError("The plane has NaN or infinite values; they are refused, not repaired.") + with np.errstate(over="ignore"): # an overflow becomes infinity, refused just below + image = pixels.astype(np.float32) # always a new array + if not np.isfinite(image).all(): + raise ValueError( + "The plane has values beyond the float32 range (about 3.4e38); they are refused." + ) + return image, {"float32_exact": float32_exact(image, pixels)} + + +def float32_exact(image: np.ndarray, pixels: np.ndarray) -> bool: + """Return whether the float32 ``image`` holds exactly the values of ``pixels``. + + Every float32 is exact in float64, so float64 compares floats and integers up to 32 bits + exactly. 64-bit integers above 2**53 are not exact in float64; they are compared as + integers after checking that the float32 values fit the integer type. + """ + back = image.astype(np.float64) + if pixels.dtype.kind == "f" or pixels.dtype.itemsize < 8: + return bool(np.array_equal(back, pixels.astype(np.float64))) + info = np.iinfo(pixels.dtype) + # float(info.max) rounds up to 2**63 or 2**64, the first value that does not fit. + if (back < info.min).any() or (back >= float(info.max)).any(): + return False + return bool(np.array_equal(back.astype(pixels.dtype), pixels)) diff --git a/src/morphobricks_normalization/stain.py b/src/morphobricks_normalization/stain.py index 72077f2..d94836e 100644 --- a/src/morphobricks_normalization/stain.py +++ b/src/morphobricks_normalization/stain.py @@ -1,53 +1,83 @@ -"""RGB-to-scalar conversion with stain presets (contract §3.4, plan N67, issue #94). - -:func:`rgb_to_scalar` turns an 8-bit sRGB image into one float32 map, chosen by a named, -versioned preset (:data:`PRESETS`), for consumers that need a single channel from a stained -image (FiberFeatures first). It is a method family of its own: the RGB joint-percentile -normalization of ``rgb.py`` is unchanged and unrelated. - -Presets (P0): ``he_eosin``, ``he_hematoxylin`` (Ruifrok-Johnston H&E vectors through -``skimage.color.separate_stains`` and ``hed_from_rgb``), ``luminance_inverted`` (baseline) and -``lab_a`` (CIELAB a*, comparator). A new stain is a new preset record, not a new code path. - -Input contract. ``axes`` is exactly ``"YXC"`` or ``"CYX"`` and is never guessed; the array is -3-D with exactly three samples on the colour axis (R, G, B in that order; an alpha channel is -refused, never dropped); the dtype is ``uint8`` (sRGB; other dtypes are refused, 16-bit and ICC -handling are later work). Nothing is reordered, rescaled or repaired silently, and the input -array is never modified. ``uint8`` cannot hold NaN or infinity, so the output is finite. - -Behaviour inherited from scikit-image, accepted in P0 and recorded in the provenance: -``separate_stains`` floors every pixel value at 1e-6, treats 255 as white, divides by -``-log(1e-6)`` and **sets negative component values to 0**. The provenance reports -``clipped_fraction``: the fraction of the pixels whose target component is exactly 0 after the -call (a pixel that was negative or truly 0 cannot be told apart). ``lab_a`` is signed and not -clipped, and ``luminance_inverted`` is not clipped either; both report ``None``. - -``fixed_range`` of a preset is ``None`` for now: the stain-OD range a consumer needs for fixed -intensity scaling must be calibrated from a reference set and recorded, not hard-coded from one -image (issue #94). - -scikit-image is the optional ``histology`` extra, imported only inside the conversion; without -it :func:`rgb_to_scalar` raises :class:`ImportError`, and importing this module needs NumPy only. +"""Image-to-scalar conversion with stain and modality presets (contract §3.4, N67, N68). + +:func:`to_scalar` turns one declared imaging unit into one float32 ``(Y, X)`` map, chosen by a +named, versioned preset (:data:`PRESETS`), for consumers that need a single channel +(FiberFeatures first). It is a method family of its own: the RGB joint-percentile +normalization of ``rgb.py`` is unchanged and unrelated. :func:`rgb_to_scalar` is the P0 entry +point, kept unchanged; it calls :func:`to_scalar` without the optional arguments. + +Methods (a new stain is a new preset record, not a new code path): + +- ``passthrough`` (``fluorescence``, ``shg``): a 2-D ``"YX"`` plane of any real numeric dtype, + returned as a float32 copy. Values that are not finite, or not finite in float32, are refused. + The cast rounds integers above 2**24 and float64 values to the nearest float32; the + provenance reports ``float32_exact``. +- ``fixed_deconvolution`` (``he_eosin``, ``he_hematoxylin``, ``trichrome_blue``): + ``skimage.color.separate_stains`` with a fixed scikit-image matrix. It floors every pixel at + 1e-6, treats 255 as white, divides by ``-log(1e-6)`` and **sets negative component values to + 0**; ``clipped_fraction`` is the fraction of pixels whose target component is exactly 0 + after the call (a pixel clipped from a negative value and one that was 0 cannot be told + apart). +- ``luminance`` (``luminance_inverted``: ``1 - rgb2gray``; ``psr_polarized``: ``rgb2gray``, not + inverted, since birefringent fibres are bright on a dark field) and ``cielab`` (``lab_a``: + a*; ``trichrome_green``: minus a*), signed or bounded, not clipped. +- ``custom_od`` (``psr_brightfield``, ``he_custom``, ``trichrome_custom``): optical density + against a declared ``white`` and user-supplied ``stain_vectors``, in NumPy + (:mod:`morphobricks_normalization.optical_density`), negative concentrations kept. + +RGB input contract. ``axes`` is exactly ``"YXC"`` or ``"CYX"`` and is never guessed; the array +is 3-D with exactly three samples on the colour axis (R, G, B in that order; an alpha channel is +refused, never dropped); the dtype is ``uint8`` (sRGB; 16-bit and ICC handling are later work). +Nothing is reordered, rescaled or repaired silently, and the input array is never modified. + +Arguments a preset does not use (``white`` for ``lab_a``, ``polarizer`` for ``he_eosin``) are +refused, never ignored. ``psr_polarized`` needs ``polarizer`` (``"circular"`` or ``"linear"``): +a linear polarizer darkens fibres at some angles to its axis, so the provenance hints record +``orientation_valid = False`` for it and consumers must refuse orientation features then. + +Each preset carries consumer ``hints`` (``polarity``, ``scaling``, ``mask``, ``measures``): +recommendations for consumers such as FiberFeatures (#214), recorded in the provenance and +never applied here. ``fixed_range`` is ``None`` until calibrated from a reference set. + +scikit-image is the optional ``histology`` extra, imported only by the presets that call it; +``passthrough`` and ``custom_od`` need NumPy only, and importing this module needs NumPy only. """ from __future__ import annotations from collections.abc import Mapping -from dataclasses import dataclass +from dataclasses import dataclass, field from types import MappingProxyType, ModuleType from typing import Any import numpy as np -AXES = ("YXC", "CYX") -SAMPLE_COUNT = 3 +from morphobricks_normalization import optical_density as od +from morphobricks_normalization import scalar_inputs as inputs + +AXES = inputs.AXES +PLANE_AXES = inputs.PLANE_AXES +SAMPLE_COUNT = inputs.SAMPLE_COUNT COLOUR_SPACE = "srgb" INSTALL_HINT = "pip install morphobricks-normalization[histology]" +MODALITIES = ("fluorescence", "shg", "brightfield", "polarized") +SCALING_HINTS = ("per_image", "fixed") +MASK_HINTS = ( + "exclude_nuclei_rbc_empty", + "exclude_nuclei", + "exclude_counterstain", + "background", + "none", + "optional", +) +POLARIZERS = ("circular", "linear") +OPTIONAL_ARGUMENTS = ("white", "stain_vectors", "polarizer") +EMPTY: Mapping[str, Any] = MappingProxyType({}) @dataclass(frozen=True) class StainPreset: - """A frozen, versioned recipe for one scalar map from an RGB image. + """A frozen, versioned recipe for one scalar map from one imaging unit. Attributes ---------- @@ -56,19 +86,32 @@ class StainPreset: version Version of the recipe; any change of numbers is a new version. stain - Stain declared in the MorphoBricks Spec (``"he"``), or ``None`` for a colour projection. + Stain as proposed for the MorphoBricks Spec (``"he"``, ``"picrosirius_red"``, + ``"trichrome_blue"``, ``"trichrome_green"``), or ``None``. method - ``fixed_deconvolution``, ``luminance`` or ``cielab``. + ``passthrough``, ``fixed_deconvolution``, ``luminance``, ``cielab`` or ``custom_od``. matrix - Name of the scikit-image stain matrix (``hed_from_rgb``), or ``None``. + Name of the scikit-image stain matrix (``hed_from_rgb``, ``bpx_from_rgb``), or ``None``. target_component Index of the component (deconvolution) or CIELAB channel returned, or ``None``. exclusion_component Component returned as the optional exclusion map (haematoxylin), or ``None``. polarity - ``"bright"``: a higher value means more stain. + ``"bright"``: a higher value means more of the measured structure. fixed_range - Stain range for fixed intensity scaling; ``None`` until calibrated from a reference set. + Range for fixed intensity scaling; ``None`` until calibrated from a reference set. + modality + One of :data:`MODALITIES`. + hints + Read-only consumer hints: ``polarity``, ``scaling`` (:data:`SCALING_HINTS`), ``mask`` + (:data:`MASK_HINTS`) and ``measures`` (short text). Recommendations, not applied here. + needs + The optional arguments of :func:`to_scalar` this preset requires; the others are + refused. + inverted + ``luminance``: ``1 - gray`` instead of ``gray``; ``cielab``: minus the channel. + target_stain, exclusion_stain + ``custom_od``: names of the target and exclusion vectors in ``stain_vectors``. """ id: str @@ -80,18 +123,25 @@ class StainPreset: exclusion_component: int | None = None polarity: str = "bright" fixed_range: tuple[float, float] | None = None + modality: str = "brightfield" + hints: Mapping[str, str] = field(default=EMPTY, hash=False) + needs: frozenset[str] = frozenset() + inverted: bool = False + target_stain: str | None = None + exclusion_stain: str | None = None @dataclass(frozen=True) class ScalarResult: - """The scalar map of :func:`rgb_to_scalar`. + """The scalar map of :func:`to_scalar`. Attributes ---------- image float32 ``(Y, X)``, finite. exclusion - float32 ``(Y, X)`` (haematoxylin for ``he_eosin``), or ``None`` if the preset has none. + float32 ``(Y, X)`` (haematoxylin for ``he_eosin``, the exclusion stain for + ``custom_od``), or ``None`` if the preset has none. provenance Read-only mapping that records how the map was made. """ @@ -101,6 +151,16 @@ class ScalarResult: provenance: Mapping[str, Any] +def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: + """Return read-only consumer hints; every preset is ``bright``.""" + return MappingProxyType( + {"polarity": "bright", "scaling": scaling, "mask": mask, "measures": measures} + ) + + +_HE_MASK = "exclude_nuclei_rbc_empty" +_CUSTOM = frozenset({"white", "stain_vectors"}) + PRESETS: Mapping[str, StainPreset] = MappingProxyType( { preset.id: preset @@ -111,42 +171,100 @@ class ScalarResult: matrix="hed_from_rgb", target_component=1, exclusion_component=0, + hints=_hints("fixed", _HE_MASK, "eosin fibrillar texture, not collagen-specific"), + ), + StainPreset( + "he_hematoxylin", + stain="he", + matrix="hed_from_rgb", + target_component=0, + hints=_hints("fixed", "none", "haematoxylin, mainly nuclei"), + ), + StainPreset( + "luminance_inverted", + method="luminance", + inverted=True, + hints=_hints("fixed", _HE_MASK, "inverted brightness (baseline), nuclei-dominated"), + ), + StainPreset( + "lab_a", + method="cielab", + target_component=1, + hints=_hints("fixed", _HE_MASK, "CIELAB a* redness (comparator), low contrast"), + ), + StainPreset( + "fluorescence", + method="passthrough", + modality="fluorescence", + hints=_hints("per_image", "optional", "the labelled structure"), + ), + StainPreset( + "shg", + method="passthrough", + modality="shg", + hints=_hints("fixed", "none", "fibrillar collagen (mainly type I)"), + ), + StainPreset( + "psr_brightfield", + stain="picrosirius_red", + method="custom_od", + needs=_CUSTOM, + target_stain="psr", + exclusion_stain="counterstain", + hints=_hints("fixed", "exclude_counterstain", "collagen"), + ), + StainPreset( + "psr_polarized", + stain="picrosirius_red", + method="luminance", + modality="polarized", + needs=frozenset({"polarizer"}), + hints=_hints("fixed", "background", "birefringent packed collagen"), + ), + StainPreset( + "trichrome_blue", + stain="trichrome_blue", + matrix="bpx_from_rgb", + target_component=0, + hints=_hints("fixed", "exclude_nuclei", "collagen (methyl blue)"), + ), + StainPreset( + "trichrome_green", + stain="trichrome_green", + method="cielab", + target_component=1, + inverted=True, + hints=_hints("fixed", "exclude_nuclei", "collagen (green, minus CIELAB a*)"), + ), + StainPreset( + "he_custom", + stain="he", + method="custom_od", + needs=_CUSTOM, + target_stain="eosin", + exclusion_stain="hematoxylin", + hints=_hints("fixed", _HE_MASK, "eosin fibrillar texture, not collagen-specific"), + ), + StainPreset( + "trichrome_custom", + method="custom_od", + needs=_CUSTOM, + target_stain="collagen", + exclusion_stain="nuclei", + hints=_hints("fixed", "exclude_nuclei", "collagen"), ), - StainPreset("he_hematoxylin", stain="he", matrix="hed_from_rgb", target_component=0), - StainPreset("luminance_inverted", method="luminance"), - StainPreset("lab_a", method="cielab", target_component=1), ) } ) -"""The four P0 presets, read-only. No trichrome presets yet (issue #94, P1).""" - - -def _to_yxc(pixels: np.ndarray, axes: str) -> np.ndarray: - """Return a new C-contiguous ``(Y, X, 3)`` array after checking the input contract. +"""The preset catalogue (issue #94, owner decision of 13 November 2026), read-only.""" - Always a copy, never a view of ``pixels``, so no later step can modify the caller's array. - """ - if not isinstance(pixels, np.ndarray): - raise TypeError(f"pixels must be a NumPy array; got {type(pixels).__name__}.") - if not isinstance(axes, str) or axes not in AXES: - raise ValueError(f"axes must be exactly 'YXC' or 'CYX', never guessed; got {axes!r}.") - if pixels.ndim != 3: - raise ValueError(f"Expected a 3-D array with axes {axes}; got {pixels.ndim}-D.") - colour_axis = axes.index("C") - if pixels.shape[colour_axis] != SAMPLE_COUNT: - extra = ( - " An alpha channel is refused, not dropped." if pixels.shape[colour_axis] == 4 else "" - ) - raise ValueError( - f"Expected exactly 3 samples (R, G, B) on the colour axis of {axes}; " - f"got {pixels.shape[colour_axis]}.{extra}" - ) - if pixels.dtype != np.uint8: - raise TypeError( - f"Only 8-bit sRGB (uint8) is supported; got {pixels.dtype}. " - "Nothing is rescaled or converted silently." - ) - return np.array(np.moveaxis(pixels, colour_axis, 2), order="C", copy=True) +PRESETS_BY_MODALITY: Mapping[str, tuple[str, ...]] = MappingProxyType( + { + modality: tuple(key for key, preset in PRESETS.items() if preset.modality == modality) + for modality in MODALITIES + } +) +"""Preset ids by modality, read-only, in registry order.""" def _skimage_color() -> tuple[ModuleType, ModuleType]: @@ -170,9 +288,11 @@ def _convert(rgb: np.ndarray, preset: StainPreset) -> tuple[np.ndarray, np.ndarr if preset.exclusion_component is not None: exclusion = components[..., preset.exclusion_component] elif preset.method == "luminance": - image = 1.0 - color.rgb2gray(rgb) + gray = color.rgb2gray(rgb) + image = 1.0 - gray if preset.inverted else gray elif preset.method == "cielab": - image = color.rgb2lab(rgb)[..., preset.target_component] + channel = color.rgb2lab(rgb)[..., preset.target_component] + image = -channel if preset.inverted else channel else: raise ValueError(f"Preset {preset.id!r} has an unknown method {preset.method!r}.") as_float32 = np.ascontiguousarray(image, dtype=np.float32) @@ -180,8 +300,152 @@ def _convert(rgb: np.ndarray, preset: StainPreset) -> tuple[np.ndarray, np.ndarr return as_float32, excluded, skimage.__version__ +def _custom_od( + rgb: np.ndarray, preset: StainPreset, white: object, stain_vectors: object +) -> tuple[np.ndarray, np.ndarray, dict[str, Any]]: + """Unmix ``rgb`` with the declared white and vectors; return image, exclusion, provenance.""" + white_values = od.check_white(white) + matrix = od.stain_matrix(stain_vectors, preset.target_stain, preset.exclusion_stain) + density, floored = od.optical_density(rgb, white_values) + amounts = od.concentrations(density, matrix) + image = np.ascontiguousarray(amounts[..., 0], dtype=np.float32) + exclusion = np.ascontiguousarray(amounts[..., 1], dtype=np.float32) + extra = { + "white": tuple(white_values.tolist()), + "stain_names": matrix.names, + "stain_matrix": matrix.rows, + "condition_number": matrix.condition, + "residual_completed": matrix.completed, + "od_definition": "-log10(max(I, 1) / white)", + "floored_fraction": floored, + "negative_fraction": float(np.mean(amounts[..., 0] < 0)), + "above_white_fraction": od.above_white_fraction(rgb, white_values), + } + return image, exclusion, extra + + +def _check_arguments(recipe: StainPreset, given: Mapping[str, Any]) -> None: + """Refuse a missing required argument, an unused one, or an unknown polarizer.""" + for name in OPTIONAL_ARGUMENTS: + if name in recipe.needs and given[name] is None: + raise ValueError(f"Preset {recipe.id!r} needs {name}; it has no default.") + if name not in recipe.needs and given[name] is not None: + raise ValueError( + f"Preset {recipe.id!r} does not use {name}; it is refused rather than ignored." + ) + polarizer = given["polarizer"] + if "polarizer" in recipe.needs and polarizer not in POLARIZERS: + raise ValueError(f"polarizer must be 'circular' or 'linear'; got {polarizer!r}.") + + +def _provenance( + recipe: StainPreset, pixels: np.ndarray, axes: str, extra: Mapping[str, Any] +) -> Mapping[str, Any]: + """Return the read-only provenance of one conversion.""" + hints = dict(recipe.hints) + if "polarizer" in recipe.needs: + # A linear polarizer darkens fibres at some angles to its axis (orientation bias). + hints["orientation_valid"] = extra["polarizer"] == "circular" + record = { + "preset": recipe.id, + "preset_version": recipe.version, + "stain": recipe.stain, + "modality": recipe.modality, + "method": recipe.method, + "matrix": recipe.matrix, + "target_component": recipe.target_component, + "exclusion_component": recipe.exclusion_component, + "polarity": recipe.polarity, + "hints": MappingProxyType(hints), + "input_dtype": str(pixels.dtype), + "axes": axes, + "colour_space": None if recipe.method == "passthrough" else COLOUR_SPACE, + "skimage_version": None, + "clipped_fraction": None, + } + record.update(extra) + return MappingProxyType(record) + + +def to_scalar( + pixels: np.ndarray, + *, + axes: str, + preset: str, + white: object = None, + stain_vectors: Mapping[str, Any] | None = None, + polarizer: str | None = None, +) -> ScalarResult: + """Convert one imaging unit to one float32 map with a preset of :data:`PRESETS`. + + Parameters + ---------- + pixels + ``passthrough`` presets: a 2-D real numeric plane. Others: a ``uint8`` 3-D array with + exactly three samples (R, G, B) on the colour axis. Never modified. + axes + ``"YX"`` for ``passthrough``; ``"YXC"`` or ``"CYX"`` otherwise. Required, never guessed. + preset + A key of :data:`PRESETS`. + white + ``custom_od`` only: background R, G, B of the 8-bit image, three finite values in + ``(0, 255]`` (see :func:`white_from_region`). + stain_vectors + ``custom_od`` only: mapping stain name -> three R, G, B optical-density components; the + preset's target and exclusion names, and optionally ``"residual"`` (otherwise the cross + product of the first two). See :func:`stain_vectors_from_regions`. + polarizer + ``psr_polarized`` only: ``"circular"`` or ``"linear"``. + + Returns + ------- + ScalarResult + ``image`` float32 ``(Y, X)``; ``exclusion`` or ``None``; read-only ``provenance`` + (contract §3.4): preset id, version, stain, modality, method, matrix, components, + polarity, hints, input dtype, axes, colour space, scikit-image version (``None`` when + not used), ``clipped_fraction``, and the method's own entries (``custom_od``: white, + stain names, matrix, condition number, floored, negative and above-white fractions; + ``psr_polarized``: polarizer; ``passthrough``: ``float32_exact``). + + Raises + ------ + ValueError + Unknown preset; wrong axes, shape or values; a missing, unused or invalid optional + argument; stain vectors that are missing, extra, collinear or ill-conditioned. + TypeError + ``pixels`` is not a NumPy array or has a refused dtype. + ImportError + The preset needs scikit-image and it is not installed. + """ + if preset not in PRESETS: + raise ValueError(f"Unknown preset {preset!r}; known presets: {', '.join(PRESETS)}.") + recipe = PRESETS[preset] + given = {"white": white, "stain_vectors": stain_vectors, "polarizer": polarizer} + _check_arguments(recipe, given) + exclusion = None + if recipe.method == "passthrough": + image, extra = inputs.plane_float32(pixels, axes) + elif recipe.method == "custom_od": + image, exclusion, extra = _custom_od( + inputs.rgb_yxc(pixels, axes), recipe, white, stain_vectors + ) + else: + image, exclusion, version = _convert(inputs.rgb_yxc(pixels, axes), recipe) + clipped = float(np.mean(image == 0)) if recipe.method == "fixed_deconvolution" else None + extra = {"skimage_version": version, "clipped_fraction": clipped} + if "polarizer" in recipe.needs: + extra["polarizer"] = polarizer + provenance = _provenance(recipe, pixels, axes, extra) + return ScalarResult(image=image, exclusion=exclusion, provenance=provenance) + + def rgb_to_scalar(pixels: np.ndarray, *, axes: str, preset: str) -> ScalarResult: - """Convert an 8-bit sRGB image to one float32 map with a stain preset. + """Convert an 8-bit sRGB image to one float32 map with a preset that needs no argument. + + The P0 entry point (N67): :func:`to_scalar` without ``white``, ``stain_vectors`` or + ``polarizer``, so a preset that needs one of them is refused with a message naming it. + Outputs of ``he_eosin``, ``he_hematoxylin``, ``luminance_inverted`` and ``lab_a`` are + unchanged from N67. Parameters ---------- @@ -196,43 +460,81 @@ def rgb_to_scalar(pixels: np.ndarray, *, axes: str, preset: str) -> ScalarResult Returns ------- ScalarResult - ``image`` float32 ``(Y, X)``; ``exclusion`` (``he_eosin`` only) or ``None``; - ``provenance`` with preset id and version, stain, method, matrix, target and exclusion - component, polarity, input dtype, axes, colour space, scikit-image version and - ``clipped_fraction`` (fraction of target-component pixels exactly 0, ``None`` where the - preset does not clip). + As :func:`to_scalar`. Raises ------ ValueError - Unknown preset, axes other than ``"YXC"``/``"CYX"``, not 3-D, or not three samples - (an alpha channel included). + Unknown preset, a preset that needs an argument, axes other than ``"YXC"``/``"CYX"``, + not 3-D, or not three samples (an alpha channel included). TypeError ``pixels`` is not a NumPy array or its dtype is not ``uint8``. ImportError scikit-image is not installed (``pip install morphobricks-normalization[histology]``). """ - if preset not in PRESETS: - raise ValueError(f"Unknown preset {preset!r}; known presets: {', '.join(PRESETS)}.") - recipe = PRESETS[preset] - rgb = _to_yxc(pixels, axes) - image, exclusion, skimage_version = _convert(rgb, recipe) - clipped = float(np.mean(image == 0)) if recipe.method == "fixed_deconvolution" else None - provenance = MappingProxyType( - { - "preset": recipe.id, - "preset_version": recipe.version, - "stain": recipe.stain, - "method": recipe.method, - "matrix": recipe.matrix, - "target_component": recipe.target_component, - "exclusion_component": recipe.exclusion_component, - "polarity": recipe.polarity, - "input_dtype": str(pixels.dtype), - "axes": axes, - "colour_space": COLOUR_SPACE, - "skimage_version": skimage_version, - "clipped_fraction": clipped, - } - ) - return ScalarResult(image=image, exclusion=exclusion, provenance=provenance) + return to_scalar(pixels, axes=axes, preset=preset) + + +def white_from_region( + pixels: np.ndarray, *, axes: str, region: np.ndarray +) -> tuple[float, float, float]: + """Return the median R, G, B of a boolean background region: the declared ``white``. + + Parameters + ---------- + pixels + ``uint8`` RGB array as for :func:`to_scalar`. Never modified. + axes + ``"YXC"`` or ``"CYX"``. + region + Boolean ``(Y, X)`` array selecting background (no tissue). + + Returns + ------- + tuple of float + The white, valid for ``to_scalar(..., white=...)``. + + Raises + ------ + ValueError + Empty region, shape mismatch, or a median channel at 0 or 255 (dark or saturated). + TypeError + Wrong pixel dtype, or the region is not boolean. + """ + return od.region_white(inputs.rgb_yxc(pixels, axes), region) + + +def stain_vectors_from_regions( + pixels: np.ndarray, *, axes: str, white: object, regions: Mapping[str, np.ndarray] +) -> Mapping[str, tuple[float, float, float]]: + """Return one unit optical-density vector per labelled region (the TWOMBLI approach). + + Each vector is the mean of ``-log10(max(I, 1) / white)`` over the region, normalised to + unit length: label the matrix stain and the nuclear stain; the background gives ``white`` + (:func:`white_from_region`). + + Parameters + ---------- + pixels + ``uint8`` RGB array as for :func:`to_scalar`. Never modified. + axes + ``"YXC"`` or ``"CYX"``. + white + Declared background R, G, B, three finite values in ``(0, 255]``. + regions + Mapping stain name -> boolean ``(Y, X)`` array of pixels of that stain only. + + Returns + ------- + Mapping + Read-only, stain name -> unit vector, valid as ``stain_vectors``. + + Raises + ------ + ValueError + Invalid white; an empty or wrongly shaped region; a region whose mean optical density + is about 0 (background, not stain). + TypeError + Wrong pixel dtype, ``regions`` not a non-empty mapping, or a region not boolean. + """ + return od.region_vectors(inputs.rgb_yxc(pixels, axes), od.check_white(white), regions) diff --git a/tests/stain_simulation.py b/tests/stain_simulation.py new file mode 100644 index 0000000..5ab7d35 --- /dev/null +++ b/tests/stain_simulation.py @@ -0,0 +1,55 @@ +"""Synthetic stained images from the Beer-Lambert law, with their exact error bounds. + +``I = white * 10 ** (-C @ M)`` per pixel, with the rows of ``M`` unit optical-density vectors. +Rounding to uint8 moves each value by at most 0.5, so the optical density of a rounded value +``r`` (at least 1) is off by at most ``log10(r / (r - 0.5))`` per channel (the larger of the +two one-sided bounds, ``log10((r + 0.5) / r)`` being smaller). Concentrations are +``OD @ inv(M)``, so stain ``j`` is off by at most ``sum_c bound_c * |inv(M)[c, j]|``. +""" + +import numpy as np + +WHITE = (240.0, 235.0, 245.0) +# Unit-free optical-density directions (R, G, B absorbance), normalised where used. +HEMATOXYLIN = (0.65, 0.70, 0.29) # Ruifrok-Johnston +EOSIN = (0.07, 0.99, 0.11) # Ruifrok-Johnston +DAB = (0.27, 0.57, 0.78) # Ruifrok-Johnston +PSR = (0.10, 0.80, 0.55) # picrosirius red: absorbs green and blue, passes red +COUNTERSTAIN = HEMATOXYLIN + + +def unit(vector): + vector = np.asarray(vector, dtype=np.float64) + return vector / np.linalg.norm(vector) + + +def unit_matrix(*vectors): + """Rows: the unit vectors; with two, the third is their normalised cross product.""" + rows = [unit(v) for v in vectors] + if len(rows) == 2: + rows.append(unit(np.cross(rows[0], rows[1]))) + return np.array(rows) + + +def concentrations(seed, shape, stains, high=0.6): + """Non-negative concentrations, uniform in [0, high], one per stain.""" + return np.random.default_rng(seed).uniform(0.0, high, (*shape, stains)) + + +def beer_lambert(amounts, matrix, white=WHITE): + """Float intensities of concentrations ``amounts`` (..., k) with the first k rows of M.""" + od = amounts @ matrix[: amounts.shape[-1]] + return np.asarray(white) * 10.0 ** (-od) + + +def to_uint8(intensities): + assert intensities.min() >= 1.5 and intensities.max() <= 254.5, "no floor, no saturation" + return np.rint(intensities).astype(np.uint8) + + +def concentration_bound(rounded, matrix): + """Per-pixel, per-stain bound of the concentration error caused by uint8 rounding.""" + r = rounded.astype(np.float64) + od_bound = np.log10(r / (r - 0.5)) + # float32 output and float64 arithmetic: a few float32 ulps of values below about 2. + return od_bound @ np.abs(np.linalg.inv(matrix)) + 1e-6 diff --git a/tests/test_optical_density.py b/tests/test_optical_density.py new file mode 100644 index 0000000..8efb366 --- /dev/null +++ b/tests/test_optical_density.py @@ -0,0 +1,227 @@ +"""Optical density and stain unmixing (NumPy only): analytic checks of the Beer-Lambert maths.""" + +import math + +import numpy as np +import pytest + +from morphobricks_normalization import optical_density as od +from stain_simulation import ( + COUNTERSTAIN, + DAB, + EOSIN, + HEMATOXYLIN, + PSR, + WHITE, + beer_lambert, + concentration_bound, + concentrations, + to_uint8, + unit, + unit_matrix, +) + +WHITE_ARRAY = np.array(WHITE) + + +def three(): + return {"eosin": EOSIN, "hematoxylin": HEMATOXYLIN, "residual": DAB} + + +def test_float_round_trip_is_exact_with_three_vectors(): + matrix = od.stain_matrix(three(), "eosin", "hematoxylin") + amounts = concentrations(0, (8, 9), 3) + image = beer_lambert(amounts, matrix.array()) + density, floored = od.optical_density(image, WHITE_ARRAY) + recovered = od.concentrations(density, matrix) + assert floored == 0.0 + np.testing.assert_allclose(recovered, amounts, rtol=0, atol=1e-12) + + +def test_float_round_trip_is_exact_with_two_vectors_and_the_cross_product(): + matrix = od.stain_matrix({"psr": PSR, "counterstain": COUNTERSTAIN}, "psr", "counterstain") + assert matrix.completed and matrix.names == ("psr", "counterstain", "residual") + expected = unit_matrix(PSR, COUNTERSTAIN) + np.testing.assert_allclose(matrix.array(), expected, rtol=0, atol=1e-15) + amounts = concentrations(1, (8, 9), 2) + density, _ = od.optical_density(beer_lambert(amounts, expected), WHITE_ARRAY) + recovered = od.concentrations(density, matrix) + np.testing.assert_allclose(recovered[..., :2], amounts, rtol=0, atol=1e-12) + np.testing.assert_allclose(recovered[..., 2], 0.0, rtol=0, atol=1e-12) + + +@pytest.mark.parametrize("vectors", [three(), {"eosin": EOSIN, "hematoxylin": HEMATOXYLIN}]) +def test_uint8_round_trip_is_within_the_rounding_bound(vectors): + matrix = od.stain_matrix(vectors, "eosin", "hematoxylin") + m = matrix.array() + amounts = concentrations(2, (20, 20), len(vectors)) + rounded = to_uint8(beer_lambert(amounts, m)) + density, _ = od.optical_density(rounded, WHITE_ARRAY) + error = np.abs(od.concentrations(density, matrix)[..., : len(vectors)] - amounts) + bound = concentration_bound(rounded, m)[..., : len(vectors)] + assert (error <= bound).all() + assert error.max() > 1e-4 # the quantisation is real, so the bound is exercised + + +def test_optical_density_floors_zero_and_reports_it(): + image = np.full((2, 5, 3), 120, np.uint8) + image[0, 0] = (0, 50, 50) + image[1, 4] = (0, 0, 0) + density, floored = od.optical_density(image, WHITE_ARRAY) + assert floored == pytest.approx(2 / 10) + assert density[0, 0, 0] == pytest.approx(-math.log10(1 / WHITE[0])) + assert np.isfinite(density).all() + + +def test_above_white_fraction_counts_any_brighter_channel(): + image = np.full((2, 4, 3), 100, np.uint8) + image[0, 0] = (241, 100, 100) # red above 240 + image[0, 1] = (240, 235, 245) # equal to white: not above + image[1, 3] = (255, 255, 255) + assert od.above_white_fraction(image, WHITE_ARRAY) == pytest.approx(2 / 8) + density, _ = od.optical_density(image, WHITE_ARRAY) + assert density[0, 0, 0] < 0 and np.all(density[0, 1] == 0) + + +def test_two_vector_condition_number_is_cot_half_angle(): + for degrees in (10.0, 30.0, 60.0, 90.0): + theta = math.radians(degrees) + a, b = (1.0, 0.0, 0.0), (math.cos(theta), math.sin(theta), 0.0) + matrix = od.stain_matrix({"t": a, "e": b}, "t", "e") + assert matrix.condition == pytest.approx(1 / math.tan(theta / 2), rel=1e-12) + # The 5-degree rule keeps two-vector matrices well below MAX_CONDITION. + assert 1 / math.tan(math.radians(od.MIN_ANGLE_DEGREES) / 2) < od.MAX_CONDITION + + +def test_vectors_are_normalised_and_scale_free(): + a = od.stain_matrix({"t": PSR, "e": COUNTERSTAIN}, "t", "e") + b = od.stain_matrix({"t": np.multiply(PSR, 7), "e": np.multiply(COUNTERSTAIN, 0.1)}, "t", "e") + np.testing.assert_allclose(a.array(), b.array(), rtol=0, atol=1e-15) + np.testing.assert_allclose(np.linalg.norm(a.array(), axis=1), 1.0, rtol=0, atol=1e-15) + + +@pytest.mark.parametrize( + ("vectors", "message"), + [ + ({"t": PSR}, "missing \\['e'\\]"), + ({"e": PSR}, "missing \\['t'\\]"), + ({"t": PSR, "e": COUNTERSTAIN, "dab": DAB}, "not allowed \\['dab'\\]"), + ({"t": PSR, "e": (0.0, 0.0, 0.0)}, "zero length"), + ({"t": PSR, "e": (np.nan, 0.5, 0.5)}, "not finite"), + ({"t": PSR, "e": (np.inf, 0.5, 0.5)}, "not finite"), + ({"t": PSR, "e": (0.5, 0.5)}, "three numbers"), + ({"t": PSR, "e": ("a", "b", "c")}, "three numbers"), + ({"t": PSR, "e": np.multiply(PSR, 2.0)}, "collinear"), + ({"t": PSR, "e": np.multiply(PSR, -1.0)}, "collinear"), + ({"t": (1.0, 0.0, 0.0), "e": (1.0, 0.07, 0.0)}, "4.00 degrees"), + ({"t": PSR, "e": COUNTERSTAIN, "residual": PSR}, "collinear"), + ], +) +def test_invalid_stain_vectors_are_refused(vectors, message): + with pytest.raises(ValueError, match=message): + od.stain_matrix(vectors, "t", "e") + + +def test_an_angle_just_above_the_limit_is_accepted(): + theta = math.radians(od.MIN_ANGLE_DEGREES + 0.01) + od.stain_matrix({"t": (1, 0, 0), "e": (math.cos(theta), math.sin(theta), 0)}, "t", "e") + + +def test_three_nearly_coplanar_vectors_are_refused_as_ill_conditioned(): + third = unit((1.0, 1.0, 0.02)) # 45 degrees from both, but almost in their plane + rows = np.array([(1, 0, 0), (0, 1, 0), third]) + assert np.linalg.cond(rows) > od.MAX_CONDITION + with pytest.raises(ValueError, match="ill-conditioned"): + od.stain_matrix({"t": rows[0], "e": rows[1], "residual": third}, "t", "e") + + +def test_a_non_mapping_is_refused(): + with pytest.raises(TypeError, match="mapping"): + od.stain_matrix([PSR, COUNTERSTAIN], "t", "e") + + +@pytest.mark.parametrize( + "white", + [(0, 200, 200), (256, 200, 200), (-1, 200, 200), (np.nan, 200, 200), (np.inf, 1, 1)], +) +def test_white_out_of_range_is_refused(white): + with pytest.raises(ValueError, match=r"\(0, 255\]"): + od.check_white(white) + + +@pytest.mark.parametrize( + "white", [(200, 200), (200, 200, 200, 255), 240, "white", (True, True, True), [[1, 2, 3]]] +) +def test_white_of_wrong_shape_or_type_is_refused(white): + with pytest.raises(ValueError, match="three numbers"): + od.check_white(white) + + +def test_white_accepts_255_and_fractions(): + assert od.check_white([255, 0.5, 170]).tolist() == [255.0, 0.5, 170.0] + + +def test_region_white_is_the_background_median(): + rng = np.random.default_rng(3) + image = rng.integers(20, 120, (10, 12, 3)).astype(np.uint8) + background = np.zeros((10, 12), bool) + background[:4] = True + image[background] = rng.integers(160, 181, (int(background.sum()), 3)) + expected = tuple(np.median(image[background].astype(float), axis=0)) + assert od.region_white(image, background) == expected + + +@pytest.mark.parametrize("value", [0, 255]) +def test_region_white_refuses_dark_or_saturated_backgrounds(value): + image = np.full((4, 4, 3), 200, np.uint8) + image[..., 1] = value + with pytest.raises(ValueError, match="dark or saturated"): + od.region_white(image, np.ones((4, 4), bool)) + + +@pytest.mark.parametrize( + ("region", "error", "message"), + [ + (np.zeros((4, 4), bool), ValueError, "empty"), + (np.ones((4, 5), bool), ValueError, "shape"), + (np.ones((4, 4), np.uint8), TypeError, "boolean"), + ([[True]], TypeError, "boolean"), + ], +) +def test_region_white_refuses_bad_regions(region, error, message): + with pytest.raises(error, match=message): + od.region_white(np.full((4, 4, 3), 200, np.uint8), region) + + +def test_region_vectors_recover_the_generating_vectors(): + m = unit_matrix(PSR, COUNTERSTAIN) + amounts = np.zeros((20, 20, 2)) + amounts[:, :10, 0] = concentrations(4, (20, 10), 1, high=0.8)[..., 0] + 0.2 + amounts[:, 10:, 1] = concentrations(5, (20, 10), 1, high=0.8)[..., 0] + 0.2 + image = to_uint8(beer_lambert(amounts, m)) + left = np.zeros((20, 20), bool) + left[:, :10] = True + vectors = od.region_vectors(image, WHITE_ARRAY, {"psr": left, "counterstain": ~left}) + for name, row in (("psr", m[0]), ("counterstain", m[1])): + angle = math.degrees(math.acos(min(1.0, float(np.dot(vectors[name], row))))) + assert angle < 0.2, (name, angle) + assert np.linalg.norm(vectors[name]) == pytest.approx(1.0) + with pytest.raises(TypeError): + vectors["psr"] = (1.0, 0.0, 0.0) + + +def test_region_vectors_refuse_background_like_and_empty_regions(): + image = np.empty((4, 4, 3), np.uint8) + image[...] = (239, 234, 244) # one grey level below white: OD about 0.002 per channel + region = np.ones((4, 4), bool) + with pytest.raises(ValueError, match="background"): + od.region_vectors(image, WHITE_ARRAY, {"psr": region}) + with pytest.raises(ValueError, match="empty"): + od.region_vectors(image, WHITE_ARRAY, {"psr": np.zeros((4, 4), bool)}) + with pytest.raises(TypeError, match="non-empty mapping"): + od.region_vectors(image, WHITE_ARRAY, {}) + + +def test_eosin_dab_hematoxylin_are_well_conditioned(): + matrix = od.stain_matrix(three(), "eosin", "hematoxylin") + assert matrix.condition < 10 diff --git a/tests/test_stain.py b/tests/test_stain.py index 6c369a1..8857935 100644 --- a/tests/test_stain.py +++ b/tests/test_stain.py @@ -7,6 +7,7 @@ from morphobricks_normalization.stain import PRESETS, ScalarResult, StainPreset, rgb_to_scalar +P0 = ["he_eosin", "he_hematoxylin", "luminance_inverted", "lab_a"] EOSIN = (230, 130, 190) HEMATOXYLIN = (80, 60, 150) @@ -47,14 +48,55 @@ def test_the_module_imports_and_fails_clearly_without_skimage(): assert "pip install morphobricks-normalization[histology]" in result.stdout -def test_the_registry_has_exactly_the_four_p0_presets(): - assert list(PRESETS) == ["he_eosin", "he_hematoxylin", "luminance_inverted", "lab_a"] +def test_the_four_p0_presets_are_unchanged(): + assert list(PRESETS)[:4] == P0 + numbers = { + key: (p.version, p.stain, p.method, p.matrix, p.target_component, p.exclusion_component) + for key, p in PRESETS.items() + if key in P0 + } + assert numbers == { + "he_eosin": (1, "he", "fixed_deconvolution", "hed_from_rgb", 1, 0), + "he_hematoxylin": (1, "he", "fixed_deconvolution", "hed_from_rgb", 0, None), + "luminance_inverted": (1, None, "luminance", None, None, None), + "lab_a": (1, None, "cielab", None, 1, None), + } assert all( - p.version == 1 and p.polarity == "bright" and p.fixed_range is None - for p in PRESETS.values() + p.polarity == "bright" and p.fixed_range is None and not p.needs + for key, p in PRESETS.items() + if key in P0 ) +@pytest.mark.parametrize("preset", P0) +def test_p0_outputs_equal_the_direct_scikit_image_calls(preset): + # Regression for N68: the outputs of b152472 were these calls cast to float32. + color = pytest.importorskip("skimage.color") + rgb = random_image(5, shape=(13, 17)) + hed = color.separate_stains(rgb, color.hed_from_rgb) + expected = { + "he_eosin": (hed[..., 1], hed[..., 0]), + "he_hematoxylin": (hed[..., 0], None), + "luminance_inverted": (1.0 - color.rgb2gray(rgb), None), + "lab_a": (color.rgb2lab(rgb)[..., 1], None), + }[preset] + result = rgb_to_scalar(rgb, axes="YXC", preset=preset) + assert result.image.tobytes() == expected[0].astype(np.float32).tobytes() + if expected[1] is None: + assert result.exclusion is None + else: + assert result.exclusion.tobytes() == expected[1].astype(np.float32).tobytes() + + +@pytest.mark.parametrize( + ("preset", "needed"), + [("psr_brightfield", "white"), ("he_custom", "white"), ("psr_polarized", "polarizer")], +) +def test_rgb_to_scalar_refuses_presets_that_need_an_argument(preset, needed): + with pytest.raises(ValueError, match=f"needs {needed}"): + rgb_to_scalar(random_image(), axes="YXC", preset=preset) + + def test_presets_are_immutable(): with pytest.raises(TypeError): PRESETS["new"] = StainPreset("new") # type: ignore[index] @@ -62,7 +104,7 @@ def test_presets_are_immutable(): PRESETS["he_eosin"].version = 2 # type: ignore[misc] -@pytest.mark.parametrize("preset", list(PRESETS)) +@pytest.mark.parametrize("preset", P0) def test_yxc_and_cyx_give_bitwise_identical_outputs(preset): pytest.importorskip("skimage") yxc = random_image() @@ -118,7 +160,7 @@ def test_unknown_preset_is_refused_listing_the_known_ones(): rgb_to_scalar(random_image(), axes="YXC", preset="trichrome_blue_bpx") -@pytest.mark.parametrize("preset", list(PRESETS)) +@pytest.mark.parametrize("preset", P0) def test_input_is_unchanged_and_output_is_float32_yx_finite(preset): pytest.importorskip("skimage") image = random_image(shape=(6, 11)) @@ -190,11 +232,18 @@ def test_provenance_contents_and_clipped_fraction(): "preset": "he_eosin", "preset_version": 1, "stain": "he", + "modality": "brightfield", "method": "fixed_deconvolution", "matrix": "hed_from_rgb", "target_component": 1, "exclusion_component": 0, "polarity": "bright", + "hints": { + "polarity": "bright", + "scaling": "fixed", + "mask": "exclude_nuclei_rbc_empty", + "measures": "eosin fibrillar texture, not collagen-specific", + }, "input_dtype": "uint8", "axes": "YXC", "colour_space": "srgb", @@ -224,9 +273,11 @@ def test_provenance_is_immutable(): provenance = rgb_to_scalar(random_image(), axes="YXC", preset="he_eosin").provenance with pytest.raises(TypeError): provenance["preset"] = "other" # type: ignore[index] + with pytest.raises(TypeError): + provenance["hints"]["scaling"] = "per_image" # type: ignore[index] -@pytest.mark.parametrize("preset", list(PRESETS)) +@pytest.mark.parametrize("preset", P0) def test_the_conversion_is_deterministic(preset): pytest.importorskip("skimage") image = random_image(3) diff --git a/tests/test_stain_catalogue.py b/tests/test_stain_catalogue.py new file mode 100644 index 0000000..c84f6fa --- /dev/null +++ b/tests/test_stain_catalogue.py @@ -0,0 +1,487 @@ +"""The preset catalogue of N68: passthrough, custom_od, polarized PSR, trichrome, hints.""" + +import math +import subprocess +import sys +import textwrap + +import numpy as np +import pytest + +from morphobricks_normalization.stain import ( + MASK_HINTS, + MODALITIES, + PRESETS, + PRESETS_BY_MODALITY, + SCALING_HINTS, + stain_vectors_from_regions, + to_scalar, + white_from_region, +) +from stain_simulation import ( + COUNTERSTAIN, + DAB, + EOSIN, + HEMATOXYLIN, + PSR, + WHITE, + beer_lambert, + concentration_bound, + concentrations, + to_uint8, + unit_matrix, +) + +CUSTOM = { + "psr_brightfield": ("psr", "counterstain", PSR, COUNTERSTAIN), + "he_custom": ("eosin", "hematoxylin", EOSIN, HEMATOXYLIN), + "trichrome_custom": ("collagen", "nuclei", (0.80, 0.59, 0.11), HEMATOXYLIN), +} +SKIMAGE_RGB = ["he_eosin", "he_hematoxylin", "luminance_inverted", "lab_a"] +SKIMAGE_RGB += ["psr_polarized", "trichrome_blue", "trichrome_green"] +RGB = SKIMAGE_RGB + list(CUSTOM) + + +def arguments(preset): + """The optional arguments a preset needs, with valid values.""" + if preset in CUSTOM: + target, exclusion, a, b = CUSTOM[preset] + return {"white": WHITE, "stain_vectors": {target: a, exclusion: b}} + if preset == "psr_polarized": + return {"polarizer": "circular"} + return {} + + +def random_rgb(seed=0, shape=(7, 9)): + return np.random.default_rng(seed).integers(0, 256, (*shape, 3), dtype=np.uint8) + + +def needs_skimage(preset): + if preset in SKIMAGE_RGB: + pytest.importorskip("skimage") + + +def test_the_registry_holds_the_catalogue_in_order(): + assert list(PRESETS) == [ + "he_eosin", + "he_hematoxylin", + "luminance_inverted", + "lab_a", + "fluorescence", + "shg", + "psr_brightfield", + "psr_polarized", + "trichrome_blue", + "trichrome_green", + "he_custom", + "trichrome_custom", + ] + assert set(RGB) | {"fluorescence", "shg"} == set(PRESETS) + with pytest.raises(TypeError): + PRESETS_BY_MODALITY["shg"] = () # type: ignore[index] + + +def test_presets_by_modality(): + assert dict(PRESETS_BY_MODALITY) == { + "fluorescence": ("fluorescence",), + "shg": ("shg",), + "brightfield": ( + "he_eosin", + "he_hematoxylin", + "luminance_inverted", + "lab_a", + "psr_brightfield", + "trichrome_blue", + "trichrome_green", + "he_custom", + "trichrome_custom", + ), + "polarized": ("psr_polarized",), + } + + +@pytest.mark.parametrize("preset", list(PRESETS)) +def test_every_preset_has_valid_read_only_hints(preset): + recipe = PRESETS[preset] + assert recipe.modality in MODALITIES and recipe.version == 1 and recipe.fixed_range is None + assert set(recipe.hints) == {"polarity", "scaling", "mask", "measures"} + assert recipe.hints["polarity"] == recipe.polarity == "bright" + assert recipe.hints["scaling"] in SCALING_HINTS and recipe.hints["mask"] in MASK_HINTS + assert recipe.hints["measures"] + with pytest.raises(TypeError): + recipe.hints["scaling"] = "per_image" # type: ignore[index] + + +def test_specific_hints_of_the_owner_table(): + hints = {key: (p.hints["scaling"], p.hints["mask"]) for key, p in PRESETS.items()} + assert hints["fluorescence"] == ("per_image", "optional") + assert hints["shg"] == ("fixed", "none") + assert hints["he_eosin"] == ("fixed", "exclude_nuclei_rbc_empty") + assert hints["psr_brightfield"] == ("fixed", "exclude_counterstain") + assert hints["psr_polarized"] == ("fixed", "background") + assert hints["trichrome_blue"] == hints["trichrome_green"] == ("fixed", "exclude_nuclei") + assert "not collagen-specific" in PRESETS["he_eosin"].hints["measures"] + + +# --- passthrough ------------------------------------------------------------------------ + + +@pytest.mark.parametrize("preset", ["fluorescence", "shg"]) +@pytest.mark.parametrize("dtype", ["uint8", "uint16", "int16", "int32", "float16", "float64"]) +def test_passthrough_returns_the_same_values_as_float32_copy(preset, dtype): + plane = (np.arange(30).reshape(5, 6) * 3 - 20).astype(dtype) + if dtype.startswith("uint"): + plane = np.arange(30, dtype=dtype).reshape(5, 6) + before = plane.copy() + result = to_scalar(plane, axes="YX", preset=preset) + assert result.image.dtype == np.float32 and result.image.shape == plane.shape + assert np.array_equal(result.image.astype(np.float64), plane.astype(np.float64)) + assert not np.shares_memory(result.image, plane) + assert np.array_equal(plane, before) + assert result.exclusion is None + provenance = result.provenance + assert provenance["method"] == "passthrough" and provenance["float32_exact"] is True + assert provenance["colour_space"] is None and provenance["skimage_version"] is None + assert provenance["axes"] == "YX" and provenance["input_dtype"] == dtype + + +def test_passthrough_float32_input_is_copied_not_aliased(): + plane = np.linspace(-1, 1, 12, dtype=np.float32).reshape(3, 4) + result = to_scalar(plane, axes="YX", preset="fluorescence") + assert result.image.tobytes() == plane.tobytes() + assert not np.shares_memory(result.image, plane) + + +@pytest.mark.parametrize( + ("plane", "exact"), + [ + (np.array([[2**24, 2**24 + 1]], np.int32), False), + (np.array([[2**24, 2**25]], np.int64), True), + (np.array([[2**63 - 1, 0]], np.int64), False), + (np.array([[2**64 - 1, 0]], np.uint64), False), + (np.array([[2**60, -(2**62)]], np.int64), True), + (np.array([[0.1, 0.5]], np.float64), False), + (np.array([[0.25, 0.5]], np.float64), True), + ], +) +def test_passthrough_reports_whether_the_float32_cast_is_exact(plane, exact): + result = to_scalar(plane, axes="YX", preset="shg") + assert result.provenance["float32_exact"] is exact + + +@pytest.mark.parametrize( + ("plane", "axes", "error", "message"), + [ + (np.zeros((3, 4)), "YXC", ValueError, "'YX'"), + (np.zeros((3, 4)), "XY", ValueError, "'YX'"), + (np.zeros((3, 4, 3)), "YX", ValueError, "2-D"), + (np.zeros((3, 4), bool), "YX", TypeError, "bool"), + (np.zeros((3, 4), complex), "YX", TypeError, "complex"), + (np.array([[np.nan, 1.0]]), "YX", ValueError, "NaN"), + (np.array([[np.inf, 1.0]]), "YX", ValueError, "NaN or infinite"), + (np.array([[1e39, 1.0]]), "YX", ValueError, "float32 range"), + ([[1.0, 2.0]], "YX", TypeError, "NumPy array"), + ], +) +def test_passthrough_refusals(plane, axes, error, message): + with pytest.raises(error, match=message): + to_scalar(plane, axes=axes, preset="fluorescence") + + +# --- arguments ---------------------------------------------------------------------------- + +UNUSED = [ + ("lab_a", "white"), + ("he_eosin", "polarizer"), + ("he_eosin", "stain_vectors"), + ("fluorescence", "white"), + ("shg", "polarizer"), + ("psr_polarized", "white"), + ("psr_brightfield", "polarizer"), + ("trichrome_blue", "stain_vectors"), +] + + +@pytest.mark.parametrize(("preset", "name"), UNUSED) +def test_unused_arguments_are_refused_not_ignored(preset, name): + given = arguments(preset) + given[name] = {"white": WHITE, "stain_vectors": {"a": PSR}, "polarizer": "circular"}[name] + pixels = np.zeros((3, 4)) if preset in ("fluorescence", "shg") else random_rgb() + axes = "YX" if preset in ("fluorescence", "shg") else "YXC" + with pytest.raises(ValueError, match=f"does not use {name}"): + to_scalar(pixels, axes=axes, preset=preset, **given) + + +@pytest.mark.parametrize("preset", list(CUSTOM)) +@pytest.mark.parametrize("missing", ["white", "stain_vectors"]) +def test_custom_presets_need_white_and_vectors(preset, missing): + given = arguments(preset) + del given[missing] + with pytest.raises(ValueError, match=f"needs {missing}"): + to_scalar(random_rgb(), axes="YXC", preset=preset, **given) + + +@pytest.mark.parametrize("polarizer", [None, "Circular", "none", "", 1]) +def test_psr_polarized_needs_a_declared_polarizer(polarizer): + with pytest.raises(ValueError, match="polarizer"): + to_scalar(random_rgb(), axes="YXC", preset="psr_polarized", polarizer=polarizer) + + +@pytest.mark.parametrize("preset", list(CUSTOM)) +def test_custom_presets_refuse_missing_and_extra_names(preset): + target, exclusion, a, b = CUSTOM[preset] + with pytest.raises(ValueError, match=f"missing \\['{exclusion}'\\]"): + to_scalar(random_rgb(), axes="YXC", preset=preset, white=WHITE, stain_vectors={target: a}) + with pytest.raises(ValueError, match="not allowed \\['dab'\\]"): + vectors = {target: a, exclusion: b, "dab": DAB} + to_scalar(random_rgb(), axes="YXC", preset=preset, white=WHITE, stain_vectors=vectors) + + +@pytest.mark.parametrize("white", [(240, 235), (0, 235, 245), (240, 235, 300), (np.nan, 1, 1)]) +def test_custom_presets_refuse_an_invalid_white(white): + vectors = arguments("psr_brightfield")["stain_vectors"] + with pytest.raises(ValueError, match="white"): + to_scalar( + random_rgb(), axes="YXC", preset="psr_brightfield", white=white, stain_vectors=vectors + ) + + +def test_custom_presets_refuse_collinear_and_ill_conditioned_vectors(): + with pytest.raises(ValueError, match="collinear"): + to_scalar( + random_rgb(), + axes="YXC", + preset="psr_brightfield", + white=WHITE, + stain_vectors={"psr": PSR, "counterstain": np.multiply(PSR, 3)}, + ) + third = np.add(PSR, COUNTERSTAIN) / 2 + (0.003, -0.003, 0.0) + with pytest.raises(ValueError, match="ill-conditioned"): + to_scalar( + random_rgb(), + axes="YXC", + preset="psr_brightfield", + white=WHITE, + stain_vectors={"psr": PSR, "counterstain": COUNTERSTAIN, "residual": third}, + ) + + +# --- custom_od ---------------------------------------------------------------------------- + + +@pytest.mark.parametrize("preset", list(CUSTOM)) +@pytest.mark.parametrize("explicit_third", [False, True]) +def test_custom_od_beer_lambert_round_trip(preset, explicit_third): + target, exclusion, a, b = CUSTOM[preset] + vectors = {target: a, exclusion: b} + m = unit_matrix(a, b) + stains = 2 + if explicit_third: + vectors["residual"] = DAB + m = unit_matrix(a, b, DAB) + stains = 3 + amounts = concentrations(6, (24, 24), stains) + rgb = to_uint8(beer_lambert(amounts, m)) + cyx = np.moveaxis(rgb, 2, 0).copy() + result = to_scalar(cyx, axes="CYX", preset=preset, white=WHITE, stain_vectors=vectors) + bound = concentration_bound(rgb, m) + assert (np.abs(result.image - amounts[..., 0]) <= bound[..., 0]).all() + assert (np.abs(result.exclusion - amounts[..., 1]) <= bound[..., 1]).all() + provenance = result.provenance + assert provenance["residual_completed"] is (not explicit_third) + np.testing.assert_allclose(provenance["stain_matrix"], m, rtol=0, atol=1e-15) + assert provenance["condition_number"] == pytest.approx(np.linalg.cond(m), rel=1e-12) + assert provenance["stain_names"] == (target, exclusion, "residual") + assert provenance["white"] == WHITE and provenance["skimage_version"] is None + assert provenance["floored_fraction"] == provenance["above_white_fraction"] == 0.0 + + +def test_custom_od_keeps_negative_concentrations_and_reports_the_fractions(): + m = unit_matrix(PSR, COUNTERSTAIN) + image = np.empty((2, 5, 3), np.uint8) + image[...] = to_uint8(beer_lambert(np.array([0.5, 0.2]), m)) # ordinary stained pixel + # Target -0.15 with exclusion 0.6: outside the declared stains, kept negative. + image[0, :2] = to_uint8(beer_lambert(np.array([-0.15, 0.6]), m)) + image[1, 0] = (250, 200, 200) # brighter than the white in red + image[1, 1] = (0, 90, 90) # floored red + image[1, 2] = (255, 255, 0) # above white in R and G, floored in B + vectors = {"psr": PSR, "counterstain": COUNTERSTAIN} + result = to_scalar( + image, axes="YXC", preset="psr_brightfield", white=WHITE, stain_vectors=vectors + ) + provenance = result.provenance + assert result.image[0, 0] == pytest.approx(-0.15, abs=0.01) + assert provenance["negative_fraction"] == pytest.approx(np.mean(result.image < 0)) + assert (result.image[0, :2] < 0).all() and provenance["negative_fraction"] >= 0.2 + assert provenance["above_white_fraction"] == pytest.approx(2 / 10) + assert provenance["floored_fraction"] == pytest.approx(2 / 10) + assert provenance["clipped_fraction"] is None + assert np.isfinite(result.image).all() and np.isfinite(result.exclusion).all() + + +def test_custom_od_uses_the_declared_white_not_255(): + # A background of 170 (as in Zenodo 4550300) is stain-free against white 170, not 255. + image = np.full((3, 3, 3), 170, np.uint8) + vectors = {"eosin": EOSIN, "hematoxylin": HEMATOXYLIN} + zero = to_scalar( + image, axes="YXC", preset="he_custom", white=(170, 170, 170), stain_vectors=vectors + ) + assert np.array_equal(zero.image, np.zeros((3, 3), np.float32)) + false = to_scalar( + image, axes="YXC", preset="he_custom", white=(255, 255, 255), stain_vectors=vectors + ) + assert (false.exclusion > 0.3).all() # false haematoxylin density in a stain-free area + + +# --- scikit-image presets ----------------------------------------------------------------- + + +def test_psr_polarized_is_rgb2gray_not_inverted(): + color = pytest.importorskip("skimage.color") + rgb = random_rgb(7) + for polarizer, valid in (("circular", True), ("linear", False)): + result = to_scalar(rgb, axes="YXC", preset="psr_polarized", polarizer=polarizer) + assert result.image.tobytes() == color.rgb2gray(rgb).astype(np.float32).tobytes() + assert result.provenance["polarizer"] == polarizer + assert result.provenance["hints"]["orientation_valid"] is valid + assert result.provenance["method"] == "luminance" and result.exclusion is None + bright = to_scalar( + np.full((2, 2, 3), 250, np.uint8), axes="YXC", preset="psr_polarized", polarizer="circular" + ) + assert (bright.image > 0.9).all() # a bright fibre stays bright + + +def test_trichrome_blue_is_methyl_blue_of_bpx(): + color = pytest.importorskip("skimage.color") + rgb = random_rgb(8) + result = to_scalar(rgb, axes="YXC", preset="trichrome_blue") + expected = color.separate_stains(rgb, color.bpx_from_rgb)[..., 0].astype(np.float32) + assert result.image.tobytes() == expected.tobytes() + assert result.exclusion is None + assert result.provenance["clipped_fraction"] == float(np.mean(expected == 0)) + assert result.provenance["matrix"] == "bpx_from_rgb" + blue = to_scalar( + np.full((2, 2, 3), (40, 90, 200), np.uint8), axes="YXC", preset="trichrome_blue" + ).image + red = to_scalar( + np.full((2, 2, 3), (200, 60, 90), np.uint8), axes="YXC", preset="trichrome_blue" + ).image + assert blue.mean() > red.mean() + + +def test_trichrome_green_is_minus_lab_a_signed(): + color = pytest.importorskip("skimage.color") + rgb = random_rgb(9) + result = to_scalar(rgb, axes="YXC", preset="trichrome_green") + expected = (-color.rgb2lab(rgb)[..., 1]).astype(np.float32) + assert result.image.tobytes() == expected.tobytes() + assert result.image.min() < 0 < result.image.max() + assert result.provenance["clipped_fraction"] is None + + +# --- shared contract ---------------------------------------------------------------------- + + +@pytest.mark.parametrize("preset", RGB) +def test_yxc_and_cyx_are_bitwise_identical_and_input_is_unchanged(preset): + needs_skimage(preset) + yxc = random_rgb(10) + cyx = np.moveaxis(yxc, 2, 0).copy() + before = (yxc.copy(), cyx.copy()) + a = to_scalar(yxc, axes="YXC", preset=preset, **arguments(preset)) + b = to_scalar(cyx, axes="CYX", preset=preset, **arguments(preset)) + assert a.image.tobytes() == b.image.tobytes() + assert (a.exclusion is None) == (b.exclusion is None) + if a.exclusion is not None: + assert a.exclusion.tobytes() == b.exclusion.tobytes() + assert np.array_equal(yxc, before[0]) and np.array_equal(cyx, before[1]) + assert a.image.dtype == np.float32 and a.image.shape == (7, 9) + assert np.isfinite(a.image).all() + assert not np.shares_memory(a.image, yxc) + differing = {"axes"} + assert {k for k in a.provenance if a.provenance[k] != b.provenance[k]} == differing + + +@pytest.mark.parametrize("preset", list(PRESETS)) +def test_provenance_is_read_only_and_records_the_hints(preset): + needs_skimage(preset) + if preset in ("fluorescence", "shg"): + result = to_scalar(np.ones((2, 3), np.uint16), axes="YX", preset=preset) + else: + result = to_scalar(random_rgb(), axes="YXC", preset=preset, **arguments(preset)) + provenance = result.provenance + assert provenance["preset"] == preset and provenance["modality"] == PRESETS[preset].modality + for key, value in PRESETS[preset].hints.items(): + assert provenance["hints"][key] == value + with pytest.raises(TypeError): + provenance["preset"] = "x" # type: ignore[index] + with pytest.raises(TypeError): + provenance["hints"]["mask"] = "none" # type: ignore[index] + + +# --- helpers ------------------------------------------------------------------------------ + + +def test_helpers_give_a_white_and_vectors_that_unmix_the_image(): + m = unit_matrix(PSR, COUNTERSTAIN) + amounts = np.zeros((30, 30, 2)) + amounts[10:20, :, 0] = 0.4 + concentrations(11, (10, 30), 1, high=0.5)[..., 0] + amounts[20:, :, 1] = 0.4 + concentrations(12, (10, 30), 1, high=0.5)[..., 0] + rgb = to_uint8(beer_lambert(amounts, m)) + cyx = np.moveaxis(rgb, 2, 0).copy() + rows = np.arange(30)[:, None] * np.ones((1, 30), int) + background, psr, nuclei = rows < 10, (rows >= 10) & (rows < 20), rows >= 20 + white = white_from_region(cyx, axes="CYX", region=background) + assert white == tuple(float(round(v)) for v in WHITE) + vectors = stain_vectors_from_regions( + cyx, axes="CYX", white=white, regions={"psr": psr, "counterstain": nuclei} + ) + for name, row in (("psr", m[0]), ("counterstain", m[1])): + angle = math.degrees(math.acos(min(1.0, float(np.dot(vectors[name], row))))) + assert angle < 0.5, (name, angle) + result = to_scalar( + cyx, axes="CYX", preset="psr_brightfield", white=white, stain_vectors=vectors + ) + assert np.abs(result.image[background]).max() < 0.01 + assert result.image[psr].mean() > 0.5 and np.abs(result.image[nuclei]).mean() < 0.02 + + +def test_helpers_check_their_inputs(): + rgb = random_rgb() + with pytest.raises(ValueError, match="axes"): + white_from_region(rgb, axes="YX", region=np.ones((7, 9), bool)) + with pytest.raises(ValueError, match="white"): + stain_vectors_from_regions( + rgb, axes="YXC", white=(0, 1, 1), regions={"psr": np.ones((7, 9), bool)} + ) + + +def test_numpy_only_presets_work_without_scikit_image(): + code = textwrap.dedent( + """ + import sys + sys.modules["skimage"] = None + sys.modules["skimage.color"] = None + import numpy as np + from morphobricks_normalization.stain import to_scalar + plane = to_scalar(np.arange(6.0).reshape(2, 3), axes="YX", preset="fluorescence") + rgb = np.full((2, 2, 3), 120, np.uint8) + vectors = {"psr": (0.1, 0.8, 0.55), "counterstain": (0.65, 0.7, 0.29)} + od = to_scalar(rgb, axes="YXC", preset="psr_brightfield", white=(240, 235, 245), + stain_vectors=vectors) + assert "skimage" not in [m for m in sys.modules if sys.modules[m] is not None] + for preset, extra in (("trichrome_blue", {}), ("psr_polarized", {"polarizer": "linear"})): + try: + to_scalar(rgb, axes="YXC", preset=preset, **extra) + except ImportError as error: + assert "[histology]" in str(error) + else: + raise SystemExit("no ImportError for " + preset) + print(plane.image.dtype, od.image.dtype) + """ + ) + result = subprocess.run( + [sys.executable, "-c", code], capture_output=True, text=True, check=False + ) + assert result.returncode == 0, result.stderr + assert result.stdout.split() == ["float32", "float32"] From 30bab3dfc80cd6126bb7181bc8273230a9a1433c Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 08:16:27 +0000 Subject: [PATCH 04/12] =?UTF-8?q?Document=20the=20preset=20catalogue=20(co?= =?UTF-8?q?ntract=20=C2=A73.4,=20plan=20N68)=20and=20run=20its=20tests=20i?= =?UTF-8?q?n=20CI?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Contract §3.4 gets the preset table with modality, method, target and exclusion, hints, refusals, the optical-density definitions and limits. Plan N68 records the owner's catalogue decision; CHANGELOG and the AGENTS.md repository map list the new modules. The histology CI job now runs the catalogue and optical-density tests too. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- .github/workflows/ci.yml | 4 +- AGENTS.md | 9 ++- CHANGELOG.md | 10 +++ docs/contract.md | 146 ++++++++++++++++++++++++++++----------- docs/plan.md | 17 +++++ 5 files changed, 143 insertions(+), 43 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 25808bc..f3ef0cf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -59,7 +59,7 @@ jobs: - run: python -m pip install -e ".[dev]" csbdeep - run: python -m pytest tests/test_rgb.py -rs - # The RGB-to-scalar stain presets (tests/test_stain.py) need the optional extra + # The stain presets (tests/test_stain*.py, test_optical_density.py) need the optional extra # `histology` (scikit-image); the main matrix runs without it, so these tests run here. histology: runs-on: ubuntu-latest @@ -69,4 +69,4 @@ jobs: with: python-version: "3.12" - run: python -m pip install -e ".[dev,histology]" - - run: python -m pytest tests/test_stain.py -rs + - run: python -m pytest tests/test_stain.py tests/test_stain_catalogue.py tests/test_optical_density.py -rs diff --git a/AGENTS.md b/AGENTS.md index 6f4eb9e..3a7c9fc 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -56,8 +56,13 @@ src/morphobricks_normalization/ scopes.py group_key, group_name, scope_groups: groups at each scope, repeat groups rgb.py StarDist's percentile normalization: rgb_percentile, apply_rgb (CSBDeep's float32 arithmetic), check_rgb - stain.py RGB-to-scalar stain presets: rgb_to_scalar, PRESETS (H&E first; optional - scikit-image extra `histology`; contract §3.4, N67) + stain.py image-to-scalar presets: to_scalar, rgb_to_scalar, PRESETS (catalogue with + modality and consumer hints), white_from_region, stain_vectors_from_regions + (optional scikit-image extra `histology`; contract §3.4, N67, N68) + optical_density.py optical density against a declared white, stain matrices (normalised, + cross-product completion, collinearity and conditioning checks), unmixing + (NumPy only; the `custom_od` presets) + scalar_inputs.py input checks of the presets: RGB with declared axes, single YX planes fitting.py fit, by channel ID parameters.py Parameters (apply, select, to_columns, from_columns with Normalization_Offset_ and Normalization_Scale_) diff --git a/CHANGELOG.md b/CHANGELOG.md index 557fc10..a223636 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,16 @@ The decisions are in [docs/plan.md](docs/plan.md) §12. `he_hematoxylin`, `luminance_inverted`, `lab_a`). 8-bit sRGB, `axes` `"YXC"` or `"CYX"`. scikit-image is the new optional extra `histology`. `rgb.py`, the parameter columns and `plane_histograms.csv` are unchanged. +- **Preset catalogue** (issue #94, plan N68, contract §3.4): general entry point + `to_scalar(pixels, *, axes, preset, white=None, stain_vectors=None, polarizer=None)`; + `rgb_to_scalar` delegates to it and the four P0 presets give bitwise identical outputs. New + presets `fluorescence`, `shg` (passthrough), `psr_brightfield`, `he_custom`, + `trichrome_custom` (`custom_od`: optical density against a declared white, user stain vectors, + NumPy only), `psr_polarized` (needs `polarizer`), `trichrome_blue`, `trichrome_green`. + `StainPreset` gains `modality`, read-only consumer `hints`, `needs`, `inverted`, + `target_stain` and `exclusion_stain`; the provenance gains `modality` and `hints`. New + helpers `white_from_region`, `stain_vectors_from_regions`; new pure modules `optical_density` + and `scalar_inputs`; `PRESETS_BY_MODALITY`. ## v0.5.0 (2026-11-11) diff --git a/docs/contract.md b/docs/contract.md index ab92fce..6dfdbdf 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -488,49 +488,117 @@ y_c = (x_c − offset) / scale for c = R, G, B --- -### 3.4 RGB-to-scalar stain presets (`stain.py`, issue #94, plan N67) - -A method family of its own, separate from §3.3 (which is unchanged): `rgb_to_scalar(pixels, *, -axes, preset)` returns `ScalarResult(image, exclusion, provenance)` for consumers that need one -channel from a stained image. normalization owns the conversion; consumers call it and never -reimplement it. No file reading, tiling, segmentation or WSI handling is done here. - -| `preset` (version 1) | `image` | `exclusion` | -|---|---|---| -| `he_eosin` | `separate_stains(rgb, hed_from_rgb)`, component 1 (eosin) | component 0 (haematoxylin) | -| `he_hematoxylin` | the same call, component 0 | `None` | -| `luminance_inverted` | `1 - rgb2gray(rgb)` (baseline) | `None` | -| `lab_a` | `rgb2lab(rgb)[..., 1]` (CIELAB a*, signed, not clipped; comparator) | `None` | - -Every preset has `stain` (`"he"` or `None`), `method`, `matrix`, `target_component`, -`exclusion_component`, `polarity = "bright"` (higher means more stain), `fixed_range = None` and -`version = 1`. `fixed_range` stays `None` until a stain-OD range is calibrated from a reference -set and recorded; it is not hard-coded from one image. No trichrome presets yet. - -**Inputs (refused, never repaired).** `axes` is exactly `"YXC"` or `"CYX"`, with no default; +### 3.4 Image-to-scalar presets (`stain.py`, issue #94, plan N67, N68) + +A method family of its own, separate from §3.3 (which is unchanged): `to_scalar(pixels, *, +axes, preset, white=None, stain_vectors=None, polarizer=None)` returns `ScalarResult(image, +exclusion, provenance)`, one float32 `(Y, X)` map of one declared imaging unit, for consumers +that need one channel (FiberFeatures first). normalization owns the conversion; consumers call +it and never reimplement it. `rgb_to_scalar(pixels, *, axes, preset)` (N67) is `to_scalar` +without the optional arguments. No file reading, tiling, segmentation or WSI handling is done +here. Modules: `stain.py` (API, registry), `optical_density.py` (OD maths, NumPy), +`scalar_inputs.py` (input checks). + +**Catalogue** (every preset version 1, `polarity = "bright"`, `fixed_range = None`): + +| `preset` | modality | method | `image` | `exclusion` | needs | scaling | mask | measures | +|---|---|---|---|---|---|---|---|---| +| `he_eosin` | brightfield | `fixed_deconvolution` | `separate_stains(rgb, hed_from_rgb)[..., 1]` | component 0 (haematoxylin) | — | fixed | `exclude_nuclei_rbc_empty` | eosin fibrillar texture, not collagen-specific | +| `he_hematoxylin` | brightfield | `fixed_deconvolution` | the same call, component 0 | `None` | — | fixed | `none` | haematoxylin, mainly nuclei | +| `luminance_inverted` | brightfield | `luminance` | `1 - rgb2gray(rgb)` (baseline) | `None` | — | fixed | `exclude_nuclei_rbc_empty` | inverted brightness, nuclei-dominated | +| `lab_a` | brightfield | `cielab` | `rgb2lab(rgb)[..., 1]` (a*, signed; comparator) | `None` | — | fixed | `exclude_nuclei_rbc_empty` | a* redness, low contrast | +| `fluorescence` | fluorescence | `passthrough` | the `"YX"` plane as float32 | `None` | — | per_image | `optional` | the labelled structure | +| `shg` | shg | `passthrough` | the `"YX"` plane as float32 | `None` | — | fixed | `none` | fibrillar collagen (mainly type I) | +| `psr_brightfield` | brightfield | `custom_od` | concentration of `"psr"` | concentration of `"counterstain"` | white, stain_vectors | fixed | `exclude_counterstain` | collagen | +| `psr_polarized` | polarized | `luminance` | `rgb2gray(rgb)`, **not inverted** | `None` | polarizer | fixed | `background` | birefringent packed collagen | +| `trichrome_blue` | brightfield | `fixed_deconvolution` | `separate_stains(rgb, bpx_from_rgb)[..., 0]` (Methyl Blue) | `None` | — | fixed | `exclude_nuclei` | collagen (methyl blue) | +| `trichrome_green` | brightfield | `cielab` | `-rgb2lab(rgb)[..., 1]` (minus a*, signed) | `None` | — | fixed | `exclude_nuclei` | collagen (green) | +| `he_custom` | brightfield | `custom_od` | concentration of `"eosin"` | concentration of `"hematoxylin"` | white, stain_vectors | fixed | `exclude_nuclei_rbc_empty` | eosin fibrillar texture, not collagen-specific | +| `trichrome_custom` | brightfield | `custom_od` | concentration of `"collagen"` | concentration of `"nuclei"` | white, stain_vectors | fixed | `exclude_nuclei` | collagen | + +`stain` is `"he"`, `"picrosirius_red"`, `"trichrome_blue"`, `"trichrome_green"` (the values +proposed in MorphoBricks-Spec #36) or `None` (colour projections, passthrough, +`trichrome_custom`). `PRESETS_BY_MODALITY` lists the ids by modality. `trichrome_blue` uses +scikit-image's Methyl Blue/Ponceau Fuchsin matrix, not a universal Masson profile; for a +slide-specific profile use `trichrome_custom`. + +**Hints** (`preset.hints`, read-only, recorded in the provenance): `polarity` (`"bright"`), +`scaling` (`"per_image"` or `"fixed"`), `mask` (`exclude_nuclei_rbc_empty`, `exclude_nuclei`, +`exclude_counterstain`, `background`, `none`, `optional`) and `measures` (short text). They +are recommendations for consumers (FiberFeatures #214), **never applied here**. + +**Arguments.** A preset's `needs` lists the optional arguments it requires: each one missing is +refused, and each argument it does not use (`white` for `lab_a`, `polarizer` for `he_eosin`) is +refused, never ignored. `polarizer` is `"circular"` or `"linear"`; the provenance hints of +`psr_polarized` record `orientation_valid = (polarizer == "circular")`: a linear polarizer +darkens fibres at some angles to its axis, so consumers must refuse orientation features then. +Hue is not used (later, as a separate feature). + +**RGB inputs (refused, never repaired).** `axes` is exactly `"YXC"` or `"CYX"`, with no default; the array is 3-D with exactly three samples (R, G, B) on the colour axis (an alpha channel, or any other count, is refused with a message, never dropped); the dtype is `uint8` (8-bit sRGB -only; other dtypes are refused, 16-bit and ICC handling are later); an unknown preset is -refused with the known ones listed. The source array is never modified and the result is -deterministic; `YXC` and `CYX` of the same image give bitwise identical outputs. `image` and -`exclusion` are float32 `(Y, X)` and finite (`uint8` cannot hold NaN or infinity). - -**scikit-image behaviour accepted in P0.** `separate_stains` floors every pixel value at 1e-6, -treats 255 as white, divides by `-log(1e-6)` and sets negative component values to 0. The -`clipped_fraction` of the provenance is the fraction of pixels whose target component is -exactly 0 after the call (a pixel clipped from a negative value and one that was 0 cannot be told -apart); it is `None` for `luminance_inverted` and `lab_a`, which do not clip. A self-computed -optical-density path is later work. - -**Provenance** (read-only mapping): `preset`, `preset_version`, `stain`, `method`, `matrix`, -`target_component`, `exclusion_component`, `polarity`, `input_dtype`, `axes`, `colour_space` -(`"srgb"`), `skimage_version`, `clipped_fraction`. +only; 16-bit and ICC handling are later); an unknown preset is refused with the known ones +listed. The source array is never modified and the result is deterministic; `YXC` and `CYX` of +the same image give bitwise identical outputs. `image` and `exclusion` are float32 `(Y, X)` and +finite. + +**Passthrough inputs.** `axes` is exactly `"YX"`; the array is 2-D of a real numeric dtype +(integer or float; bool and complex refused); NaN and infinity are refused, and so are values +that overflow float32 (above about 3.4e38). The result is a new float32 array (never a view of +the input) with the same values up to the float32 cast, which rounds integers above 2**24 and +float64 values; the provenance reports `float32_exact`. + +**scikit-image behaviour accepted (`fixed_deconvolution`).** `separate_stains` floors every pixel +value at 1e-6, treats 255 as white, divides by `-log(1e-6)` and sets negative component values +to 0. The `clipped_fraction` of the provenance is the fraction of pixels whose target component +is exactly 0 after the call (a pixel clipped from a negative value and one that was 0 cannot be +told apart); it is `None` for every other method. + +**`custom_od`** (NumPy, no scikit-image; `optical_density.py`): + +- `white`: the declared background R, G, B of the 8-bit image, three finite values in + `(0, 255]` (`white_from_region` measures it). It is never assumed to be 255: a background + of 170 against 255 would be read as stain. +- Optical density per channel: `OD = -log10(max(I, 1) / white)`. A value of 0 has no finite + OD and is floored at 1; `floored_fraction` is the fraction of pixels with a floored channel. + `above_white_fraction` is the fraction of pixels brighter than `white` in at least one + channel (negative OD there). +- `stain_vectors`: mapping stain name -> three R, G, B OD components, any positive scale; exactly + the preset's target and exclusion names, plus optionally `"residual"`. Missing or other names + are refused. Each vector is normalised to unit length (non-finite values, a wrong length or + a norm below 1e-6 are refused). Without `"residual"` the third row is the normalised cross + product of the first two (scikit-image's completion of `rgb_from_hed`). +- Rows of `M`: target, exclusion, residual. Two vectors closer than 5 degrees (either + direction) are refused as collinear; for two vectors completed by the cross product + `cond(M) = cot(angle / 2)`, so 5 degrees gives 22.9. A 2-norm condition number above 50 is + refused as ill-conditioned (three nearly coplanar vectors): concentration errors are at most + `cond(M)` times the relative OD error, about 1 % per grey level at mid range. +- Concentrations `C = OD @ inv(M)`; `image` is the target, `exclusion` the exclusion stain. + **Negative concentrations are kept, not clipped**; `negative_fraction` is the fraction of + pixels with a negative target. +- `stain_vectors_from_regions(pixels, *, axes, white, regions)`: for each boolean region, the + mean OD over the region normalised to unit length (the TWOMBLI approach: the user labels the + matrix stain and the nuclear stain; the background gives `white`). An empty or wrongly + shaped region, or a mean OD norm below 0.01 (background-like), is refused. +- `white_from_region(pixels, *, axes, region)`: the per-channel median over a boolean + background region; an empty or wrongly shaped region, or a median of 0 or 255 (dark or + saturated, so the true white is unknown), is refused. + +**Provenance** (read-only mapping): `preset`, `preset_version`, `stain`, `modality`, `method`, +`matrix`, `target_component`, `exclusion_component`, `polarity`, `hints` (read-only), +`input_dtype`, `axes`, `colour_space` (`"srgb"`, `None` for passthrough), `skimage_version` +(`None` when scikit-image is not called), `clipped_fraction`; plus, for `custom_od`: `white`, +`stain_names`, `stain_matrix` (unit rows), `condition_number`, `residual_completed`, +`od_definition`, `floored_fraction`, `negative_fraction`, `above_white_fraction`; for +`psr_polarized`: `polarizer` (and `hints["orientation_valid"]`); for passthrough: +`float32_exact`. **Dependency.** scikit-image is the optional extra `histology` -(`pip install morphobricks-normalization[histology]`), imported only inside the conversion; -without it `rgb_to_scalar` raises `ImportError` naming the extra. `import -morphobricks_normalization` needs NumPy only. The module is not re-exported from the package -root; import it as `morphobricks_normalization.stain`. +(`pip install morphobricks-normalization[histology]`), imported only by the presets that call +it (`fixed_deconvolution`, `luminance`, `cielab`); without it they raise `ImportError` naming the +extra, while `passthrough` and `custom_od` work with NumPy only. `import +morphobricks_normalization` needs NumPy only. The modules are not re-exported from the package +root; import `morphobricks_normalization.stain`. ## 4. Scopes, groups and weights diff --git a/docs/plan.md b/docs/plan.md index a12db6e..370031d 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -546,6 +546,23 @@ Added with the project owner on 11 November 2026, for Stage 3 (S3.2, pull reques accepted in P0 and reported as `clipped_fraction`; `fixed_range` stays `None` until calibrated. The StarDist percentile normalization (`rgb.py`, §3.1) is unchanged. Later (not in P0): trichrome and custom stain vectors (P1), RGB colour normalization (P2). +62. **N68** (13 November 2026, issue #94, owner's preset catalogue) Presets are image + preprocessing steps owned here: each turns one declared imaging unit into one float32 + `(Y, X)` map with provenance and carries a modality and **consumer hints** (`polarity`, + `scaling` per image or fixed, `mask`, `measures`) that are recommendations for consumers + (FiberFeatures #214), never applied here. New general entry point `to_scalar` (`rgb_to_scalar` + delegates, P0 outputs bitwise unchanged). Catalogue: `fluorescence` and `shg` (passthrough + of a `"YX"` plane as float32), `psr_brightfield`, `he_custom`, `trichrome_custom` + (`custom_od`), `psr_polarized`, `trichrome_blue` (`bpx_from_rgb`), `trichrome_green` (minus + CIELAB a*). `custom_od` is implemented here in NumPy, not via scikit-image: optical density + against a **declared white** (never assumed 255; measured once with `white_from_region` and + recorded), user-supplied stain vectors (two, the third completed by a cross product, or + three), collinear (< 5 degrees) and ill-conditioned (condition number > 50) matrices refused, + negative concentrations kept and reported; `stain_vectors_from_regions` derives the vectors + from labelled regions (TWOMBLI). `psr_polarized` is luminance, not inverted, and requires a + declared polarizer; only `"circular"` gives `orientation_valid` (linear polarizers darken + fibres at some angles). Arguments a preset does not use are refused. Colour normalization + stays later, with a frozen reference only (contract §3.4). Implementation choices of S3.2 where the Spec leaves the detail open (open questions 1-7 of pull request #90), each taken on the side the Spec states or the safer one, and listed in that From 5040a67d191c31d597462cd85ec8eeda830abb36 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 08:19:45 +0000 Subject: [PATCH 05/12] Use default_factory for the preset hints (Python 3.11) Python 3.11's dataclasses refuse a MappingProxyType default as mutable (3.13 accepts it), so test collection failed on the 3.11 CI jobs. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- src/morphobricks_normalization/stain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/morphobricks_normalization/stain.py b/src/morphobricks_normalization/stain.py index d94836e..66ca8e3 100644 --- a/src/morphobricks_normalization/stain.py +++ b/src/morphobricks_normalization/stain.py @@ -124,7 +124,7 @@ class StainPreset: polarity: str = "bright" fixed_range: tuple[float, float] | None = None modality: str = "brightfield" - hints: Mapping[str, str] = field(default=EMPTY, hash=False) + hints: Mapping[str, str] = field(default_factory=lambda: EMPTY, hash=False) needs: frozenset[str] = frozenset() inverted: bool = False target_stain: str | None = None From 3a7af3e941a36e1c15301d7cd755f396ca52aea4 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 08:22:24 +0000 Subject: [PATCH 06/12] Refuse masked and empty images, rename zero_fraction, clarify polarity Review of b152472 (ChatGPT, verified against this head): - np.ma.MaskedArray passed the ndarray check and np.array dropped its mask; masked arrays are now refused for RGB and plane presets. - An image with Y or X of 0 reached np.mean and gave a NaN provenance; it is now refused. - clipped_fraction counted every zero component, not only clipped ones; renamed zero_fraction (unreleased API). - polarity 'bright' is a signal direction, not a stain amount. Tests: masked and empty refusals, fixed reference values (white, black, neutral grey). The pixelwise scikit-image parity test of the P0 presets was already present since 5040a67's parent. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- docs/contract.md | 11 ++-- docs/plan.md | 6 ++- .../scalar_inputs.py | 13 +++++ src/morphobricks_normalization/stain.py | 12 +++-- tests/test_stain.py | 14 +++--- tests/test_stain_catalogue.py | 50 +++++++++++++++++-- 6 files changed, 86 insertions(+), 20 deletions(-) diff --git a/docs/contract.md b/docs/contract.md index 6dfdbdf..f561c1b 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -522,7 +522,9 @@ proposed in MorphoBricks-Spec #36) or `None` (colour projections, passthrough, scikit-image's Methyl Blue/Ponceau Fuchsin matrix, not a universal Masson profile; for a slide-specific profile use `trichrome_custom`. -**Hints** (`preset.hints`, read-only, recorded in the provenance): `polarity` (`"bright"`), +**Hints** (`preset.hints`, read-only, recorded in the provenance): `polarity` (`"bright"`: the +structure of interest has higher values than its surroundings; a signal direction, not a stain +amount), `scaling` (`"per_image"` or `"fixed"`), `mask` (`exclude_nuclei_rbc_empty`, `exclude_nuclei`, `exclude_counterstain`, `background`, `none`, `optional`) and `measures` (short text). They are recommendations for consumers (FiberFeatures #214), **never applied here**. @@ -538,7 +540,8 @@ Hue is not used (later, as a separate feature). the array is 3-D with exactly three samples (R, G, B) on the colour axis (an alpha channel, or any other count, is refused with a message, never dropped); the dtype is `uint8` (8-bit sRGB only; 16-bit and ICC handling are later); an unknown preset is refused with the known ones -listed. The source array is never modified and the result is deterministic; `YXC` and `CYX` of +listed. For every preset, RGB or plane, a masked array (`np.ma.MaskedArray`, whose mask would be +dropped) and an image with no pixels (Y or X of 0) are refused. The source array is never modified and the result is deterministic; `YXC` and `CYX` of the same image give bitwise identical outputs. `image` and `exclusion` are float32 `(Y, X)` and finite. @@ -550,7 +553,7 @@ float64 values; the provenance reports `float32_exact`. **scikit-image behaviour accepted (`fixed_deconvolution`).** `separate_stains` floors every pixel value at 1e-6, treats 255 as white, divides by `-log(1e-6)` and sets negative component values -to 0. The `clipped_fraction` of the provenance is the fraction of pixels whose target component +to 0. The `zero_fraction` of the provenance is the fraction of pixels whose target component is exactly 0 after the call (a pixel clipped from a negative value and one that was 0 cannot be told apart); it is `None` for every other method. @@ -587,7 +590,7 @@ told apart); it is `None` for every other method. **Provenance** (read-only mapping): `preset`, `preset_version`, `stain`, `modality`, `method`, `matrix`, `target_component`, `exclusion_component`, `polarity`, `hints` (read-only), `input_dtype`, `axes`, `colour_space` (`"srgb"`, `None` for passthrough), `skimage_version` -(`None` when scikit-image is not called), `clipped_fraction`; plus, for `custom_od`: `white`, +(`None` when scikit-image is not called), `zero_fraction`; plus, for `custom_od`: `white`, `stain_names`, `stain_matrix` (unit rows), `condition_number`, `residual_completed`, `od_definition`, `floored_fraction`, `negative_fraction`, `above_white_fraction`; for `psr_polarized`: `polarizer` (and `hints["orientation_valid"]`); for passthrough: diff --git a/docs/plan.md b/docs/plan.md index 370031d..1e34773 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -562,7 +562,11 @@ Added with the project owner on 11 November 2026, for Stage 3 (S3.2, pull reques from labelled regions (TWOMBLI). `psr_polarized` is luminance, not inverted, and requires a declared polarizer; only `"circular"` gives `orientation_valid` (linear polarizers darken fibres at some angles). Arguments a preset does not use are refused. Colour normalization - stays later, with a frozen reference only (contract §3.4). + stays later, with a frozen reference only (contract §3.4). Review fixes (13 November 2026, ChatGPT review of b152472): the P0 provenance key + `clipped_fraction` is renamed `zero_fraction` (it counts pixels whose component is exactly 0, + which includes pixels that were never clipped); masked arrays (`np.ma.MaskedArray`) and + images with no pixels are refused for every preset; `polarity = "bright"` means the + structure has higher values than its surroundings (the signal direction), not more stain. Implementation choices of S3.2 where the Spec leaves the detail open (open questions 1-7 of pull request #90), each taken on the side the Spec states or the safer one, and listed in that diff --git a/src/morphobricks_normalization/scalar_inputs.py b/src/morphobricks_normalization/scalar_inputs.py index 004a62d..15a8cb4 100644 --- a/src/morphobricks_normalization/scalar_inputs.py +++ b/src/morphobricks_normalization/scalar_inputs.py @@ -17,6 +17,17 @@ SAMPLE_COUNT = 3 +def _refuse_masked_or_empty(pixels: np.ndarray, spatial: tuple[int, int]) -> None: + """Refuse a masked array (its mask would be dropped) and an image with no pixels.""" + if isinstance(pixels, np.ma.MaskedArray): + raise TypeError( + "A masked array is refused: its mask would be dropped silently. Pass the data and " + "an explicit mask to the consumer instead." + ) + if 0 in spatial: + raise ValueError(f"The image has no pixels (Y, X = {spatial}); it is refused.") + + def rgb_yxc(pixels: np.ndarray, axes: str) -> np.ndarray: """Return a new C-contiguous ``(Y, X, 3)`` array after checking the RGB input contract. @@ -37,6 +48,7 @@ def rgb_yxc(pixels: np.ndarray, axes: str) -> np.ndarray: f"Expected exactly 3 samples (R, G, B) on the colour axis of {axes}; " f"got {pixels.shape[colour_axis]}.{extra}" ) + _refuse_masked_or_empty(pixels, tuple(np.delete(np.array(pixels.shape), colour_axis))) if pixels.dtype != np.uint8: raise TypeError( f"Only 8-bit sRGB (uint8) is supported; got {pixels.dtype}. " @@ -61,6 +73,7 @@ def plane_float32(pixels: np.ndarray, axes: str) -> tuple[np.ndarray, dict[str, raise ValueError(f"axes must be exactly 'YX' for a single plane; got {axes!r}.") if pixels.ndim != 2: raise ValueError(f"Expected a 2-D array with axes YX; got {pixels.ndim}-D.") + _refuse_masked_or_empty(pixels, pixels.shape) if pixels.dtype.kind not in "iuf": raise TypeError( f"Expected a real numeric dtype (integer or float); got {pixels.dtype}. " diff --git a/src/morphobricks_normalization/stain.py b/src/morphobricks_normalization/stain.py index 66ca8e3..48ead53 100644 --- a/src/morphobricks_normalization/stain.py +++ b/src/morphobricks_normalization/stain.py @@ -15,7 +15,7 @@ - ``fixed_deconvolution`` (``he_eosin``, ``he_hematoxylin``, ``trichrome_blue``): ``skimage.color.separate_stains`` with a fixed scikit-image matrix. It floors every pixel at 1e-6, treats 255 as white, divides by ``-log(1e-6)`` and **sets negative component values to - 0**; ``clipped_fraction`` is the fraction of pixels whose target component is exactly 0 + 0**; ``zero_fraction`` is the fraction of pixels whose target component is exactly 0 after the call (a pixel clipped from a negative value and one that was 0 cannot be told apart). - ``luminance`` (``luminance_inverted``: ``1 - rgb2gray``; ``psr_polarized``: ``rgb2gray``, not @@ -97,7 +97,9 @@ class StainPreset: exclusion_component Component returned as the optional exclusion map (haematoxylin), or ``None``. polarity - ``"bright"``: a higher value means more of the measured structure. + ``"bright"``: the structure of interest has higher values than its surroundings. It + is the signal direction for ridge detection, not a stain amount (``lab_a`` is redness, + ``luminance_inverted`` darkness). fixed_range Range for fixed intensity scaling; ``None`` until calibrated from a reference set. modality @@ -361,7 +363,7 @@ def _provenance( "axes": axes, "colour_space": None if recipe.method == "passthrough" else COLOUR_SPACE, "skimage_version": None, - "clipped_fraction": None, + "zero_fraction": None, } record.update(extra) return MappingProxyType(record) @@ -403,7 +405,7 @@ def to_scalar( ``image`` float32 ``(Y, X)``; ``exclusion`` or ``None``; read-only ``provenance`` (contract §3.4): preset id, version, stain, modality, method, matrix, components, polarity, hints, input dtype, axes, colour space, scikit-image version (``None`` when - not used), ``clipped_fraction``, and the method's own entries (``custom_od``: white, + not used), ``zero_fraction``, and the method's own entries (``custom_od``: white, stain names, matrix, condition number, floored, negative and above-white fractions; ``psr_polarized``: polarizer; ``passthrough``: ``float32_exact``). @@ -432,7 +434,7 @@ def to_scalar( else: image, exclusion, version = _convert(inputs.rgb_yxc(pixels, axes), recipe) clipped = float(np.mean(image == 0)) if recipe.method == "fixed_deconvolution" else None - extra = {"skimage_version": version, "clipped_fraction": clipped} + extra = {"skimage_version": version, "zero_fraction": clipped} if "polarizer" in recipe.needs: extra["polarizer"] = polarizer provenance = _provenance(recipe, pixels, axes, extra) diff --git a/tests/test_stain.py b/tests/test_stain.py index 8857935..b2accac 100644 --- a/tests/test_stain.py +++ b/tests/test_stain.py @@ -114,7 +114,7 @@ def test_yxc_and_cyx_give_bitwise_identical_outputs(preset): assert a.image.tobytes() == b.image.tobytes() if a.exclusion is not None: assert a.exclusion.tobytes() == b.exclusion.tobytes() - assert a.provenance["clipped_fraction"] == b.provenance["clipped_fraction"] + assert a.provenance["zero_fraction"] == b.provenance["zero_fraction"] @pytest.mark.parametrize("axes", ["yxc", "XYC", "YXCZ", "", "CXY", None, ("Y", "X", "C")]) @@ -210,7 +210,7 @@ def test_lab_a_is_signed_and_not_clipped(): pytest.importorskip("skimage") result = rgb_to_scalar(patches(EOSIN, (60, 160, 60)), axes="YXC", preset="lab_a") assert result.image.min() < 0 < result.image.max() - assert result.provenance["clipped_fraction"] is None + assert result.provenance["zero_fraction"] is None def test_he_eosin_exclusion_is_the_hematoxylin_image_and_others_have_none(): @@ -224,7 +224,7 @@ def test_he_eosin_exclusion_is_the_hematoxylin_image_and_others_have_none(): assert rgb_to_scalar(image, axes="YXC", preset=preset).exclusion is None -def test_provenance_contents_and_clipped_fraction(): +def test_provenance_contents_and_zero_fraction(): skimage = pytest.importorskip("skimage") image = patches(EOSIN, HEMATOXYLIN) result = rgb_to_scalar(image, axes="YXC", preset="he_eosin") @@ -248,10 +248,10 @@ def test_provenance_contents_and_clipped_fraction(): "axes": "YXC", "colour_space": "srgb", "skimage_version": skimage.__version__, - "clipped_fraction": float(np.mean(result.image == 0)), + "zero_fraction": float(np.mean(result.image == 0)), } assert dict(result.provenance) == expected - assert 0.0 <= result.provenance["clipped_fraction"] <= 1.0 + assert 0.0 <= result.provenance["zero_fraction"] <= 1.0 lab = rgb_to_scalar(image, axes="YXC", preset="lab_a").provenance assert lab["stain"] is None and lab["matrix"] is None and lab["method"] == "cielab" assert ( @@ -260,12 +260,12 @@ def test_provenance_contents_and_clipped_fraction(): ) -def test_the_clipped_fraction_counts_zeros_of_the_target_component(): +def test_the_zero_fraction_counts_zeros_of_the_target_component(): pytest.importorskip("skimage") # Left half pure white: every component is 0 there; right half the haematoxylin colour. image = patches((255, 255, 255), HEMATOXYLIN) result = rgb_to_scalar(image, axes="YXC", preset="he_hematoxylin") - assert result.provenance["clipped_fraction"] == pytest.approx(0.5) + assert result.provenance["zero_fraction"] == pytest.approx(0.5) def test_provenance_is_immutable(): diff --git a/tests/test_stain_catalogue.py b/tests/test_stain_catalogue.py index c84f6fa..d36a558 100644 --- a/tests/test_stain_catalogue.py +++ b/tests/test_stain_catalogue.py @@ -315,7 +315,7 @@ def test_custom_od_keeps_negative_concentrations_and_reports_the_fractions(): assert (result.image[0, :2] < 0).all() and provenance["negative_fraction"] >= 0.2 assert provenance["above_white_fraction"] == pytest.approx(2 / 10) assert provenance["floored_fraction"] == pytest.approx(2 / 10) - assert provenance["clipped_fraction"] is None + assert provenance["zero_fraction"] is None assert np.isfinite(result.image).all() and np.isfinite(result.exclusion).all() @@ -358,7 +358,7 @@ def test_trichrome_blue_is_methyl_blue_of_bpx(): expected = color.separate_stains(rgb, color.bpx_from_rgb)[..., 0].astype(np.float32) assert result.image.tobytes() == expected.tobytes() assert result.exclusion is None - assert result.provenance["clipped_fraction"] == float(np.mean(expected == 0)) + assert result.provenance["zero_fraction"] == float(np.mean(expected == 0)) assert result.provenance["matrix"] == "bpx_from_rgb" blue = to_scalar( np.full((2, 2, 3), (40, 90, 200), np.uint8), axes="YXC", preset="trichrome_blue" @@ -376,7 +376,7 @@ def test_trichrome_green_is_minus_lab_a_signed(): expected = (-color.rgb2lab(rgb)[..., 1]).astype(np.float32) assert result.image.tobytes() == expected.tobytes() assert result.image.min() < 0 < result.image.max() - assert result.provenance["clipped_fraction"] is None + assert result.provenance["zero_fraction"] is None # --- shared contract ---------------------------------------------------------------------- @@ -485,3 +485,47 @@ def test_numpy_only_presets_work_without_scikit_image(): ) assert result.returncode == 0, result.stderr assert result.stdout.split() == ["float32", "float32"] + + +# Review fixes (ChatGPT review of b152472): masked arrays, empty images, fixed reference values. +@pytest.mark.parametrize("axes", ["YXC", "CYX"]) +def test_a_masked_rgb_array_is_refused_not_unmasked(axes): + data = np.zeros((4, 5, 3) if axes == "YXC" else (3, 4, 5), dtype=np.uint8) + masked = np.ma.MaskedArray(data, mask=np.zeros(data.shape, dtype=bool)) + with pytest.raises(TypeError, match="masked array"): + to_scalar(masked, axes=axes, preset="luminance_inverted") + + +def test_a_masked_plane_is_refused(): + masked = np.ma.MaskedArray(np.ones((4, 5), dtype=np.float32), mask=False) + with pytest.raises(TypeError, match="masked array"): + to_scalar(masked, axes="YX", preset="fluorescence") + + +@pytest.mark.parametrize( + ("shape", "axes"), + [((0, 5, 3), "YXC"), ((4, 0, 3), "YXC"), ((3, 0, 5), "CYX"), ((3, 4, 0), "CYX")], +) +def test_an_rgb_image_without_pixels_is_refused(shape, axes): + with pytest.raises(ValueError, match="no pixels"): + to_scalar(np.zeros(shape, dtype=np.uint8), axes=axes, preset="he_eosin") + + +@pytest.mark.parametrize("shape", [(0, 5), (4, 0)]) +def test_a_plane_without_pixels_is_refused(shape): + with pytest.raises(ValueError, match="no pixels"): + to_scalar(np.zeros(shape, dtype=np.float32), axes="YX", preset="shg") + + +def test_fixed_reference_values_of_the_colour_projections(): + pytest.importorskip("skimage") + rgb = np.array([[[255, 255, 255], [0, 0, 0], [128, 128, 128]]], dtype=np.uint8) + lum = to_scalar(rgb, axes="YXC", preset="luminance_inverted").image + assert lum[0, 0] == pytest.approx(0.0, abs=1e-6) # white is not dark + assert lum[0, 1] == pytest.approx(1.0, abs=1e-6) # black is fully dark + lab = to_scalar(rgb, axes="YXC", preset="lab_a").image + # neutral greys have no redness; scikit-image gives |a*| ~ 0.002 from its D65 white point + assert np.abs(lab).max() < 1e-2 + eosin = to_scalar(rgb, axes="YXC", preset="he_eosin") + assert eosin.image[0, 0] == 0.0 # white: no stain + assert eosin.provenance["zero_fraction"] >= 1 / 3 From 32c48d751541f932e51a2ff18fb5067f19c6f26b Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 12:22:02 +0000 Subject: [PATCH 07/12] Add default stain vectors and an optional white for the colour projections he_custom and trichrome_custom get published default stain vectors (Ruifrok-Johnston H&E rows of rgb_from_hed; Methyl Blue and Ponceau Fuchsin rows of rgb_from_bpx), used when stain_vectors is not given and recorded as stain_vectors_source = preset_default with a citation; psr_brightfield keeps requiring user vectors. trichrome_custom's exclusion stain is now ponceau_fuchsin (the published matrix has no nuclear stain). white stays required for custom_od. luminance_inverted, lab_a, psr_polarized and trichrome_green accept an optional white: the RGB is rescaled to clip(rint(I * 255 / white), 0, 255) as uint8 before scikit-image, with the clipped fraction recorded. Without white the outputs are bitwise unchanged; white = 255 is the identity. The preset registry moves to stain_presets.py (file length). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- pyproject.toml | 1 + .../optical_density.py | 15 + src/morphobricks_normalization/stain.py | 271 +++++----------- .../stain_presets.py | 254 +++++++++++++++ tests/stain_simulation.py | 2 + tests/test_stain_catalogue.py | 22 +- tests/test_stain_defaults.py | 288 ++++++++++++++++++ 7 files changed, 649 insertions(+), 204 deletions(-) create mode 100644 src/morphobricks_normalization/stain_presets.py create mode 100644 tests/test_stain_defaults.py diff --git a/pyproject.toml b/pyproject.toml index 10dff82..d7d19cc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -124,5 +124,6 @@ baseline = "tools/code_rules_baseline.json" "morphobricks_normalization.tail_moments" = "pure" "morphobricks_normalization.sums" = "pure" "morphobricks_normalization.stain" = "pure" +"morphobricks_normalization.stain_presets" = "pure" "morphobricks_normalization.optical_density" = "pure" "morphobricks_normalization.scalar_inputs" = "pure" diff --git a/src/morphobricks_normalization/optical_density.py b/src/morphobricks_normalization/optical_density.py index 4ff5b18..78ae447 100644 --- a/src/morphobricks_normalization/optical_density.py +++ b/src/morphobricks_normalization/optical_density.py @@ -225,6 +225,21 @@ def above_white_fraction(rgb: np.ndarray, white: np.ndarray) -> float: return float(np.mean((rgb.astype(np.float64) > white).any(axis=-1))) +def rescale_to_white(rgb: np.ndarray, white: np.ndarray) -> tuple[np.ndarray, float]: + """Return ``(rescaled, clipped_fraction)``: ``rgb`` rescaled per channel so ``white`` is 255. + + ``rescaled = clip(rint(I * 255 / white), 0, 255)`` as ``uint8``, computed in float64 in + that order; ``rint`` rounds half to even. With ``white = (255, 255, 255)`` it is the + identity (``I * 255 / 255`` is exact in float64). ``clipped_fraction`` is the fraction of + pixels with at least one channel above 255 after rounding (brighter than the white by more + than half a grey level), set to 255. + """ + scaled = np.rint(rgb.astype(np.float64) * WHITE_MAX / white) + clipped = scaled > WHITE_MAX + rescaled = np.clip(scaled, 0.0, WHITE_MAX).astype(np.uint8) + return rescaled, float(np.mean(clipped.any(axis=-1))) + + def concentrations(od: np.ndarray, matrix: StainMatrix) -> np.ndarray: """Return ``OD @ inv(M)``: one float64 concentration per stain, signed, never clipped.""" return od @ np.linalg.inv(matrix.array()) diff --git a/src/morphobricks_normalization/stain.py b/src/morphobricks_normalization/stain.py index 48ead53..7ee6218 100644 --- a/src/morphobricks_normalization/stain.py +++ b/src/morphobricks_normalization/stain.py @@ -20,17 +20,23 @@ apart). - ``luminance`` (``luminance_inverted``: ``1 - rgb2gray``; ``psr_polarized``: ``rgb2gray``, not inverted, since birefringent fibres are bright on a dark field) and ``cielab`` (``lab_a``: - a*; ``trichrome_green``: minus a*), signed or bounded, not clipped. + a*; ``trichrome_green``: minus a*), signed or bounded, not clipped. They accept an optional + declared ``white`` (owner decision of 13 November 2026, N69): the RGB is then rescaled per + channel to ``clip(rint(I * 255 / white), 0, 255)`` as uint8 before the scikit-image call, and + the fraction of clipped pixels is reported. Without ``white`` nothing changes. - ``custom_od`` (``psr_brightfield``, ``he_custom``, ``trichrome_custom``): optical density - against a declared ``white`` and user-supplied ``stain_vectors``, in NumPy + against a declared ``white`` and ``stain_vectors``, in NumPy (:mod:`morphobricks_normalization.optical_density`), negative concentrations kept. + ``he_custom`` and ``trichrome_custom`` have published default vectors + (:data:`HE_DEFAULT_VECTORS`, :data:`BPX_DEFAULT_VECTORS`), used when ``stain_vectors`` is not + given; ``psr_brightfield`` has none. ``white`` is required for all three. RGB input contract. ``axes`` is exactly ``"YXC"`` or ``"CYX"`` and is never guessed; the array is 3-D with exactly three samples on the colour axis (R, G, B in that order; an alpha channel is refused, never dropped); the dtype is ``uint8`` (sRGB; 16-bit and ICC handling are later work). Nothing is reordered, rescaled or repaired silently, and the input array is never modified. -Arguments a preset does not use (``white`` for ``lab_a``, ``polarizer`` for ``he_eosin``) are +Arguments a preset does not use (``white`` for ``he_eosin``, ``polarizer`` for ``lab_a``) are refused, never ignored. ``psr_polarized`` needs ``polarizer`` (``"circular"`` or ``"linear"``): a linear polarizer darkens fibres at some angles to its axis, so the provenance hints record ``orientation_valid = False`` for it and consumers must refuse orientation features then. @@ -46,7 +52,7 @@ from __future__ import annotations from collections.abc import Mapping -from dataclasses import dataclass, field +from dataclasses import dataclass from types import MappingProxyType, ModuleType from typing import Any @@ -54,83 +60,29 @@ from morphobricks_normalization import optical_density as od from morphobricks_normalization import scalar_inputs as inputs +from morphobricks_normalization.stain_presets import BPX_CITATION as BPX_CITATION +from morphobricks_normalization.stain_presets import BPX_DEFAULT_VECTORS as BPX_DEFAULT_VECTORS +from morphobricks_normalization.stain_presets import HE_CITATION as HE_CITATION +from morphobricks_normalization.stain_presets import HE_DEFAULT_VECTORS as HE_DEFAULT_VECTORS +from morphobricks_normalization.stain_presets import MASK_HINTS as MASK_HINTS +from morphobricks_normalization.stain_presets import MODALITIES as MODALITIES +from morphobricks_normalization.stain_presets import PRESETS as PRESETS +from morphobricks_normalization.stain_presets import PRESETS_BY_MODALITY as PRESETS_BY_MODALITY +from morphobricks_normalization.stain_presets import SCALING_HINTS as SCALING_HINTS +from morphobricks_normalization.stain_presets import StainPreset as StainPreset AXES = inputs.AXES PLANE_AXES = inputs.PLANE_AXES SAMPLE_COUNT = inputs.SAMPLE_COUNT COLOUR_SPACE = "srgb" INSTALL_HINT = "pip install morphobricks-normalization[histology]" -MODALITIES = ("fluorescence", "shg", "brightfield", "polarized") -SCALING_HINTS = ("per_image", "fixed") -MASK_HINTS = ( - "exclude_nuclei_rbc_empty", - "exclude_nuclei", - "exclude_counterstain", - "background", - "none", - "optional", -) POLARIZERS = ("circular", "linear") OPTIONAL_ARGUMENTS = ("white", "stain_vectors", "polarizer") -EMPTY: Mapping[str, Any] = MappingProxyType({}) - - -@dataclass(frozen=True) -class StainPreset: - """A frozen, versioned recipe for one scalar map from one imaging unit. - - Attributes - ---------- - id - Preset identifier, the key in :data:`PRESETS`. - version - Version of the recipe; any change of numbers is a new version. - stain - Stain as proposed for the MorphoBricks Spec (``"he"``, ``"picrosirius_red"``, - ``"trichrome_blue"``, ``"trichrome_green"``), or ``None``. - method - ``passthrough``, ``fixed_deconvolution``, ``luminance``, ``cielab`` or ``custom_od``. - matrix - Name of the scikit-image stain matrix (``hed_from_rgb``, ``bpx_from_rgb``), or ``None``. - target_component - Index of the component (deconvolution) or CIELAB channel returned, or ``None``. - exclusion_component - Component returned as the optional exclusion map (haematoxylin), or ``None``. - polarity - ``"bright"``: the structure of interest has higher values than its surroundings. It - is the signal direction for ridge detection, not a stain amount (``lab_a`` is redness, - ``luminance_inverted`` darkness). - fixed_range - Range for fixed intensity scaling; ``None`` until calibrated from a reference set. - modality - One of :data:`MODALITIES`. - hints - Read-only consumer hints: ``polarity``, ``scaling`` (:data:`SCALING_HINTS`), ``mask`` - (:data:`MASK_HINTS`) and ``measures`` (short text). Recommendations, not applied here. - needs - The optional arguments of :func:`to_scalar` this preset requires; the others are - refused. - inverted - ``luminance``: ``1 - gray`` instead of ``gray``; ``cielab``: minus the channel. - target_stain, exclusion_stain - ``custom_od``: names of the target and exclusion vectors in ``stain_vectors``. - """ - - id: str - version: int = 1 - stain: str | None = None - method: str = "fixed_deconvolution" - matrix: str | None = None - target_component: int | None = None - exclusion_component: int | None = None - polarity: str = "bright" - fixed_range: tuple[float, float] | None = None - modality: str = "brightfield" - hints: Mapping[str, str] = field(default_factory=lambda: EMPTY, hash=False) - needs: frozenset[str] = frozenset() - inverted: bool = False - target_stain: str | None = None - exclusion_stain: str | None = None +PRESET_DEFAULT = "preset_default" +USER = "user" +"""Values of the provenance entry ``stain_vectors_source``.""" +WHITE_RESCALE = "clip(rint(I * 255 / white), 0, 255) as uint8, rint half to even" +"""The rescaling applied by a ``luminance`` or ``cielab`` preset given a ``white``.""" @dataclass(frozen=True) @@ -153,122 +105,6 @@ class ScalarResult: provenance: Mapping[str, Any] -def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: - """Return read-only consumer hints; every preset is ``bright``.""" - return MappingProxyType( - {"polarity": "bright", "scaling": scaling, "mask": mask, "measures": measures} - ) - - -_HE_MASK = "exclude_nuclei_rbc_empty" -_CUSTOM = frozenset({"white", "stain_vectors"}) - -PRESETS: Mapping[str, StainPreset] = MappingProxyType( - { - preset.id: preset - for preset in ( - StainPreset( - "he_eosin", - stain="he", - matrix="hed_from_rgb", - target_component=1, - exclusion_component=0, - hints=_hints("fixed", _HE_MASK, "eosin fibrillar texture, not collagen-specific"), - ), - StainPreset( - "he_hematoxylin", - stain="he", - matrix="hed_from_rgb", - target_component=0, - hints=_hints("fixed", "none", "haematoxylin, mainly nuclei"), - ), - StainPreset( - "luminance_inverted", - method="luminance", - inverted=True, - hints=_hints("fixed", _HE_MASK, "inverted brightness (baseline), nuclei-dominated"), - ), - StainPreset( - "lab_a", - method="cielab", - target_component=1, - hints=_hints("fixed", _HE_MASK, "CIELAB a* redness (comparator), low contrast"), - ), - StainPreset( - "fluorescence", - method="passthrough", - modality="fluorescence", - hints=_hints("per_image", "optional", "the labelled structure"), - ), - StainPreset( - "shg", - method="passthrough", - modality="shg", - hints=_hints("fixed", "none", "fibrillar collagen (mainly type I)"), - ), - StainPreset( - "psr_brightfield", - stain="picrosirius_red", - method="custom_od", - needs=_CUSTOM, - target_stain="psr", - exclusion_stain="counterstain", - hints=_hints("fixed", "exclude_counterstain", "collagen"), - ), - StainPreset( - "psr_polarized", - stain="picrosirius_red", - method="luminance", - modality="polarized", - needs=frozenset({"polarizer"}), - hints=_hints("fixed", "background", "birefringent packed collagen"), - ), - StainPreset( - "trichrome_blue", - stain="trichrome_blue", - matrix="bpx_from_rgb", - target_component=0, - hints=_hints("fixed", "exclude_nuclei", "collagen (methyl blue)"), - ), - StainPreset( - "trichrome_green", - stain="trichrome_green", - method="cielab", - target_component=1, - inverted=True, - hints=_hints("fixed", "exclude_nuclei", "collagen (green, minus CIELAB a*)"), - ), - StainPreset( - "he_custom", - stain="he", - method="custom_od", - needs=_CUSTOM, - target_stain="eosin", - exclusion_stain="hematoxylin", - hints=_hints("fixed", _HE_MASK, "eosin fibrillar texture, not collagen-specific"), - ), - StainPreset( - "trichrome_custom", - method="custom_od", - needs=_CUSTOM, - target_stain="collagen", - exclusion_stain="nuclei", - hints=_hints("fixed", "exclude_nuclei", "collagen"), - ), - ) - } -) -"""The preset catalogue (issue #94, owner decision of 13 November 2026), read-only.""" - -PRESETS_BY_MODALITY: Mapping[str, tuple[str, ...]] = MappingProxyType( - { - modality: tuple(key for key, preset in PRESETS.items() if preset.modality == modality) - for modality in MODALITIES - } -) -"""Preset ids by modality, read-only, in registry order.""" - - def _skimage_color() -> tuple[ModuleType, ModuleType]: """Import ``skimage`` and ``skimage.color`` lazily; name the extra when missing.""" try: @@ -279,6 +115,24 @@ def _skimage_color() -> tuple[ModuleType, ModuleType]: return skimage, color +def _white_rescaled(rgb: np.ndarray, white: object) -> tuple[np.ndarray, dict[str, Any]]: + """Rescale ``rgb`` to a declared ``white`` if given; return it and the provenance entries.""" + if white is None: + entries = {"white": None, "white_applied": False, "white_rescale": None} + entries |= {"white_clipped_fraction": None, "above_white_fraction": None} + return rgb, entries + values = od.check_white(white) + rescaled, clipped = od.rescale_to_white(rgb, values) + entries = { + "white": tuple(values.tolist()), + "white_applied": True, + "white_rescale": WHITE_RESCALE, + "white_clipped_fraction": clipped, + "above_white_fraction": od.above_white_fraction(rgb, values), + } + return rescaled, entries + + def _convert(rgb: np.ndarray, preset: StainPreset) -> tuple[np.ndarray, np.ndarray | None, str]: """Return ``(image, exclusion, skimage version)`` as float32, from a ``(Y, X, 3)`` uint8.""" skimage, color = _skimage_color() @@ -305,8 +159,16 @@ def _convert(rgb: np.ndarray, preset: StainPreset) -> tuple[np.ndarray, np.ndarr def _custom_od( rgb: np.ndarray, preset: StainPreset, white: object, stain_vectors: object ) -> tuple[np.ndarray, np.ndarray, dict[str, Any]]: - """Unmix ``rgb`` with the declared white and vectors; return image, exclusion, provenance.""" + """Unmix ``rgb`` with the declared white and vectors; return image, exclusion, provenance. + + Without ``stain_vectors`` the preset's published defaults are used (checked as user + vectors are); a preset without defaults never gets here without them. + """ white_values = od.check_white(white) + source, citation = USER, None + if stain_vectors is None: + stain_vectors = preset.default_stain_vectors + source, citation = PRESET_DEFAULT, preset.default_stain_vectors_citation matrix = od.stain_matrix(stain_vectors, preset.target_stain, preset.exclusion_stain) density, floored = od.optical_density(rgb, white_values) amounts = od.concentrations(density, matrix) @@ -314,6 +176,8 @@ def _custom_od( exclusion = np.ascontiguousarray(amounts[..., 1], dtype=np.float32) extra = { "white": tuple(white_values.tolist()), + "stain_vectors_source": source, + "stain_vectors_citation": citation, "stain_names": matrix.names, "stain_matrix": matrix.rows, "condition_number": matrix.condition, @@ -331,7 +195,8 @@ def _check_arguments(recipe: StainPreset, given: Mapping[str, Any]) -> None: for name in OPTIONAL_ARGUMENTS: if name in recipe.needs and given[name] is None: raise ValueError(f"Preset {recipe.id!r} needs {name}; it has no default.") - if name not in recipe.needs and given[name] is not None: + used = recipe.needs | recipe.accepts + if name not in used and given[name] is not None: raise ValueError( f"Preset {recipe.id!r} does not use {name}; it is refused rather than ignored." ) @@ -390,12 +255,15 @@ def to_scalar( preset A key of :data:`PRESETS`. white - ``custom_od`` only: background R, G, B of the 8-bit image, three finite values in - ``(0, 255]`` (see :func:`white_from_region`). + Background R, G, B of the 8-bit image, three finite values in ``(0, 255]`` (see + :func:`white_from_region`). Required for ``custom_od``; optional for ``luminance`` and + ``cielab`` (the RGB is then rescaled to it, :data:`WHITE_RESCALE`); refused otherwise. stain_vectors ``custom_od`` only: mapping stain name -> three R, G, B optical-density components; the preset's target and exclusion names, and optionally ``"residual"`` (otherwise the cross - product of the first two). See :func:`stain_vectors_from_regions`. + product of the first two). See :func:`stain_vectors_from_regions`. Optional for + ``he_custom`` and ``trichrome_custom`` (published defaults), required for + ``psr_brightfield``. polarizer ``psr_polarized`` only: ``"circular"`` or ``"linear"``. @@ -406,8 +274,11 @@ def to_scalar( (contract §3.4): preset id, version, stain, modality, method, matrix, components, polarity, hints, input dtype, axes, colour space, scikit-image version (``None`` when not used), ``zero_fraction``, and the method's own entries (``custom_od``: white, - stain names, matrix, condition number, floored, negative and above-white fractions; - ``psr_polarized``: polarizer; ``passthrough``: ``float32_exact``). + ``stain_vectors_source`` and citation, stain names, matrix, condition number, floored, + negative and above-white fractions; ``luminance`` and ``cielab``: ``white``, + ``white_applied``, ``white_rescale``, ``white_clipped_fraction``, + ``above_white_fraction``; ``psr_polarized``: polarizer; ``passthrough``: + ``float32_exact``). Raises ------ @@ -432,9 +303,13 @@ def to_scalar( inputs.rgb_yxc(pixels, axes), recipe, white, stain_vectors ) else: - image, exclusion, version = _convert(inputs.rgb_yxc(pixels, axes), recipe) + rgb = inputs.rgb_yxc(pixels, axes) + white_entries: dict[str, Any] = {} + if "white" in recipe.accepts: + rgb, white_entries = _white_rescaled(rgb, white) + image, exclusion, version = _convert(rgb, recipe) clipped = float(np.mean(image == 0)) if recipe.method == "fixed_deconvolution" else None - extra = {"skimage_version": version, "zero_fraction": clipped} + extra = {"skimage_version": version, "zero_fraction": clipped, **white_entries} if "polarizer" in recipe.needs: extra["polarizer"] = polarizer provenance = _provenance(recipe, pixels, axes, extra) diff --git a/src/morphobricks_normalization/stain_presets.py b/src/morphobricks_normalization/stain_presets.py new file mode 100644 index 0000000..321f92e --- /dev/null +++ b/src/morphobricks_normalization/stain_presets.py @@ -0,0 +1,254 @@ +"""The preset catalogue of :mod:`morphobricks_normalization.stain` (contract §3.4, N68, N69). + +Frozen, versioned :class:`StainPreset` records in the read-only registry :data:`PRESETS`, their +consumer hints and the published default stain vectors of the ``custom_od`` presets +(:data:`HE_DEFAULT_VECTORS`, :data:`BPX_DEFAULT_VECTORS`). Data only: the conversions are in +``stain.py``. Importing this module needs nothing but the standard library. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from dataclasses import dataclass, field +from types import MappingProxyType +from typing import Any + +MODALITIES = ("fluorescence", "shg", "brightfield", "polarized") +SCALING_HINTS = ("per_image", "fixed") +MASK_HINTS = ( + "exclude_nuclei_rbc_empty", + "exclude_nuclei", + "exclude_counterstain", + "background", + "none", + "optional", +) +EMPTY: Mapping[str, Any] = MappingProxyType({}) + +# Ruifrok AC, Johnston DA. Quantification of histochemical staining by color deconvolution. +# Anal Quant Cytol Histol. 2001;23(4):291-299. Rows 0 (haematoxylin) and 1 (eosin) of +# scikit-image's ``skimage.color.rgb_from_hed`` (scikit-image 0.26), copied so that the preset +# needs no scikit-image; tests/test_stain_defaults.py checks them against scikit-image. +HE_CITATION = ( + "Ruifrok & Johnston 2001, Anal Quant Cytol Histol 23:291-299; " + "skimage.color.rgb_from_hed rows 0 (hematoxylin) and 1 (eosin)" +) +HE_DEFAULT_VECTORS: Mapping[str, tuple[float, float, float]] = MappingProxyType( + {"hematoxylin": (0.65, 0.70, 0.29), "eosin": (0.07, 0.99, 0.11)} +) +# Methyl Blue and Ponceau Fuchsin (Masson's trichrome), rows 0 and 1 of scikit-image's +# ``skimage.color.rgb_from_bpx`` (scikit-image 0.26), the values of the "Masson Trichrome" +# vectors of G. Landini's Colour Deconvolution plugin for ImageJ/Fiji, with the method of +# Ruifrok & Johnston 2001. Copied here; checked against scikit-image by the tests. +BPX_CITATION = ( + "Ruifrok & Johnston 2001 method; Masson Trichrome vectors of Landini's Colour " + "Deconvolution (Fiji); skimage.color.rgb_from_bpx rows 0 (methyl blue) and 1 " + "(ponceau fuchsin)" +) +BPX_DEFAULT_VECTORS: Mapping[str, tuple[float, float, float]] = MappingProxyType( + { + "collagen": (0.7995107, 0.5913521, 0.10528667), + "ponceau_fuchsin": (0.09997159, 0.73738605, 0.6680326), + } +) + + +@dataclass(frozen=True) +class StainPreset: + """A frozen, versioned recipe for one scalar map from one imaging unit. + + Attributes + ---------- + id + Preset identifier, the key in :data:`PRESETS`. + version + Version of the recipe; any change of numbers is a new version. + stain + Stain as proposed for the MorphoBricks Spec (``"he"``, ``"picrosirius_red"``, + ``"trichrome_blue"``, ``"trichrome_green"``), or ``None``. + method + ``passthrough``, ``fixed_deconvolution``, ``luminance``, ``cielab`` or ``custom_od``. + matrix + Name of the scikit-image stain matrix (``hed_from_rgb``, ``bpx_from_rgb``), or ``None``. + target_component + Index of the component (deconvolution) or CIELAB channel returned, or ``None``. + exclusion_component + Component returned as the optional exclusion map (haematoxylin), or ``None``. + polarity + ``"bright"``: the structure of interest has higher values than its surroundings. It + is the signal direction for ridge detection, not a stain amount (``lab_a`` is redness, + ``luminance_inverted`` darkness). + fixed_range + Range for fixed intensity scaling; ``None`` until calibrated from a reference set. + modality + One of :data:`MODALITIES`. + hints + Read-only consumer hints: ``polarity``, ``scaling`` (:data:`SCALING_HINTS`), ``mask`` + (:data:`MASK_HINTS`) and ``measures`` (short text). Recommendations, not applied here. + needs + The optional arguments of :func:`to_scalar` this preset requires. + accepts + The optional arguments this preset uses when given but does not require (``white`` + for ``luminance`` and ``cielab``, ``stain_vectors`` when there are defaults). Every + argument in neither ``needs`` nor ``accepts`` is refused. + inverted + ``luminance``: ``1 - gray`` instead of ``gray``; ``cielab``: minus the channel. + target_stain, exclusion_stain + ``custom_od``: names of the target and exclusion vectors in ``stain_vectors``. + default_stain_vectors + ``custom_od``: read-only published stain vectors (target and exclusion names), used + when ``stain_vectors`` is not given, or ``None`` (then ``stain_vectors`` is needed). + default_stain_vectors_citation + Where the default vectors come from, recorded in the provenance. + """ + + id: str + version: int = 1 + stain: str | None = None + method: str = "fixed_deconvolution" + matrix: str | None = None + target_component: int | None = None + exclusion_component: int | None = None + polarity: str = "bright" + fixed_range: tuple[float, float] | None = None + modality: str = "brightfield" + hints: Mapping[str, str] = field(default_factory=lambda: EMPTY, hash=False) + needs: frozenset[str] = frozenset() + accepts: frozenset[str] = frozenset() + inverted: bool = False + target_stain: str | None = None + exclusion_stain: str | None = None + default_stain_vectors: Mapping[str, tuple[float, float, float]] | None = field( + default=None, hash=False + ) + default_stain_vectors_citation: str | None = None + + +def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: + """Return read-only consumer hints; every preset is ``bright``.""" + return MappingProxyType( + {"polarity": "bright", "scaling": scaling, "mask": mask, "measures": measures} + ) + + +_HE_MASK = "exclude_nuclei_rbc_empty" +_CUSTOM = frozenset({"white", "stain_vectors"}) +_WHITE = frozenset({"white"}) + +PRESETS: Mapping[str, StainPreset] = MappingProxyType( + { + preset.id: preset + for preset in ( + StainPreset( + "he_eosin", + stain="he", + matrix="hed_from_rgb", + target_component=1, + exclusion_component=0, + hints=_hints("fixed", _HE_MASK, "eosin fibrillar texture, not collagen-specific"), + ), + StainPreset( + "he_hematoxylin", + stain="he", + matrix="hed_from_rgb", + target_component=0, + hints=_hints("fixed", "none", "haematoxylin, mainly nuclei"), + ), + StainPreset( + "luminance_inverted", + method="luminance", + inverted=True, + accepts=_WHITE, + hints=_hints("fixed", _HE_MASK, "inverted brightness (baseline), nuclei-dominated"), + ), + StainPreset( + "lab_a", + method="cielab", + target_component=1, + accepts=_WHITE, + hints=_hints("fixed", _HE_MASK, "CIELAB a* redness (comparator), low contrast"), + ), + StainPreset( + "fluorescence", + method="passthrough", + modality="fluorescence", + hints=_hints("per_image", "optional", "the labelled structure"), + ), + StainPreset( + "shg", + method="passthrough", + modality="shg", + hints=_hints("fixed", "none", "fibrillar collagen (mainly type I)"), + ), + StainPreset( + "psr_brightfield", + stain="picrosirius_red", + method="custom_od", + needs=_CUSTOM, + target_stain="psr", + exclusion_stain="counterstain", + hints=_hints("fixed", "exclude_counterstain", "collagen"), + ), + StainPreset( + "psr_polarized", + stain="picrosirius_red", + method="luminance", + modality="polarized", + needs=frozenset({"polarizer"}), + accepts=_WHITE, + hints=_hints("fixed", "background", "birefringent packed collagen"), + ), + StainPreset( + "trichrome_blue", + stain="trichrome_blue", + matrix="bpx_from_rgb", + target_component=0, + hints=_hints("fixed", "exclude_nuclei", "collagen (methyl blue)"), + ), + StainPreset( + "trichrome_green", + stain="trichrome_green", + method="cielab", + target_component=1, + inverted=True, + accepts=_WHITE, + hints=_hints("fixed", "exclude_nuclei", "collagen (green, minus CIELAB a*)"), + ), + StainPreset( + "he_custom", + stain="he", + method="custom_od", + needs=_WHITE, + accepts=frozenset({"stain_vectors"}), + target_stain="eosin", + exclusion_stain="hematoxylin", + default_stain_vectors=HE_DEFAULT_VECTORS, + default_stain_vectors_citation=HE_CITATION, + hints=_hints("fixed", _HE_MASK, "eosin fibrillar texture, not collagen-specific"), + ), + StainPreset( + "trichrome_custom", + method="custom_od", + needs=_WHITE, + accepts=frozenset({"stain_vectors"}), + target_stain="collagen", + # Ponceau Fuchsin (cytoplasm, muscle), not nuclei: the published Masson + # matrix has no nuclear stain (iron haematoxylin is near black). The mask hint + # stays exclude_nuclei: a nuclear mask is built in preprocessing (N69). + exclusion_stain="ponceau_fuchsin", + default_stain_vectors=BPX_DEFAULT_VECTORS, + default_stain_vectors_citation=BPX_CITATION, + hints=_hints("fixed", "exclude_nuclei", "collagen"), + ), + ) + } +) +"""The preset catalogue (issue #94, owner decision of 13 November 2026), read-only.""" + +PRESETS_BY_MODALITY: Mapping[str, tuple[str, ...]] = MappingProxyType( + { + modality: tuple(key for key, preset in PRESETS.items() if preset.modality == modality) + for modality in MODALITIES + } +) +"""Preset ids by modality, read-only, in registry order.""" diff --git a/tests/stain_simulation.py b/tests/stain_simulation.py index 5ab7d35..309105c 100644 --- a/tests/stain_simulation.py +++ b/tests/stain_simulation.py @@ -16,6 +16,8 @@ DAB = (0.27, 0.57, 0.78) # Ruifrok-Johnston PSR = (0.10, 0.80, 0.55) # picrosirius red: absorbs green and blue, passes red COUNTERSTAIN = HEMATOXYLIN +METHYL_BLUE = (0.80, 0.59, 0.11) # Masson trichrome, about rgb_from_bpx row 0 +PONCEAU_FUCHSIN = (0.10, 0.74, 0.67) # about rgb_from_bpx row 1 def unit(vector): diff --git a/tests/test_stain_catalogue.py b/tests/test_stain_catalogue.py index d36a558..8e9516a 100644 --- a/tests/test_stain_catalogue.py +++ b/tests/test_stain_catalogue.py @@ -23,6 +23,8 @@ DAB, EOSIN, HEMATOXYLIN, + METHYL_BLUE, + PONCEAU_FUCHSIN, PSR, WHITE, beer_lambert, @@ -35,7 +37,7 @@ CUSTOM = { "psr_brightfield": ("psr", "counterstain", PSR, COUNTERSTAIN), "he_custom": ("eosin", "hematoxylin", EOSIN, HEMATOXYLIN), - "trichrome_custom": ("collagen", "nuclei", (0.80, 0.59, 0.11), HEMATOXYLIN), + "trichrome_custom": ("collagen", "ponceau_fuchsin", METHYL_BLUE, PONCEAU_FUCHSIN), } SKIMAGE_RGB = ["he_eosin", "he_hematoxylin", "luminance_inverted", "lab_a"] SKIMAGE_RGB += ["psr_polarized", "trichrome_blue", "trichrome_green"] @@ -191,12 +193,17 @@ def test_passthrough_refusals(plane, axes, error, message): # --- arguments ---------------------------------------------------------------------------- UNUSED = [ - ("lab_a", "white"), + ("lab_a", "polarizer"), + ("lab_a", "stain_vectors"), + ("luminance_inverted", "polarizer"), + ("he_eosin", "white"), + ("he_hematoxylin", "white"), ("he_eosin", "polarizer"), ("he_eosin", "stain_vectors"), ("fluorescence", "white"), + ("shg", "white"), ("shg", "polarizer"), - ("psr_polarized", "white"), + ("trichrome_blue", "white"), ("psr_brightfield", "polarizer"), ("trichrome_blue", "stain_vectors"), ] @@ -212,9 +219,12 @@ def test_unused_arguments_are_refused_not_ignored(preset, name): to_scalar(pixels, axes=axes, preset=preset, **given) -@pytest.mark.parametrize("preset", list(CUSTOM)) -@pytest.mark.parametrize("missing", ["white", "stain_vectors"]) -def test_custom_presets_need_white_and_vectors(preset, missing): +@pytest.mark.parametrize( + ("preset", "missing"), + [(preset, "white") for preset in CUSTOM] + [("psr_brightfield", "stain_vectors")], +) +def test_custom_presets_need_white_and_psr_needs_vectors(preset, missing): + # he_custom and trichrome_custom have default vectors (N69); white is never assumed. given = arguments(preset) del given[missing] with pytest.raises(ValueError, match=f"needs {missing}"): diff --git a/tests/test_stain_defaults.py b/tests/test_stain_defaults.py new file mode 100644 index 0000000..ae3818d --- /dev/null +++ b/tests/test_stain_defaults.py @@ -0,0 +1,288 @@ +"""Default stain vectors and the optional white of the colour projections (plan N69).""" + +import subprocess +import sys +import textwrap + +import numpy as np +import pytest + +from morphobricks_normalization import optical_density as od +from morphobricks_normalization.stain import ( + BPX_CITATION, + BPX_DEFAULT_VECTORS, + HE_CITATION, + HE_DEFAULT_VECTORS, + PRESET_DEFAULT, + PRESETS, + USER, + WHITE_RESCALE, + to_scalar, +) +from stain_simulation import ( + EOSIN, + HEMATOXYLIN, + WHITE, + beer_lambert, + concentration_bound, + concentrations, + to_uint8, + unit_matrix, +) + +WHITE_PRESETS = ["luminance_inverted", "lab_a", "psr_polarized", "trichrome_green"] +P0_WHITE = ["luminance_inverted", "lab_a"] +DEFAULTS = { + "he_custom": (HE_DEFAULT_VECTORS, HE_CITATION, "eosin", "hematoxylin"), + "trichrome_custom": (BPX_DEFAULT_VECTORS, BPX_CITATION, "collagen", "ponceau_fuchsin"), +} + + +def random_rgb(seed=0, shape=(7, 9)): + return np.random.default_rng(seed).integers(0, 256, (*shape, 3), dtype=np.uint8) + + +def extra(preset): + return {"polarizer": "circular"} if preset == "psr_polarized" else {} + + +# --- default stain vectors ---------------------------------------------------------------- + + +def test_he_defaults_are_the_hematoxylin_and_eosin_rows_of_rgb_from_hed(): + color = pytest.importorskip("skimage.color") + assert list(HE_DEFAULT_VECTORS) == ["hematoxylin", "eosin"] + assert np.array_equal(HE_DEFAULT_VECTORS["hematoxylin"], color.rgb_from_hed[0]) + assert np.array_equal(HE_DEFAULT_VECTORS["eosin"], color.rgb_from_hed[1]) + + +def test_trichrome_defaults_are_the_methyl_blue_and_ponceau_rows_of_rgb_from_bpx(): + color = pytest.importorskip("skimage.color") + assert list(BPX_DEFAULT_VECTORS) == ["collagen", "ponceau_fuchsin"] + assert np.array_equal(BPX_DEFAULT_VECTORS["collagen"], color.rgb_from_bpx[0]) + assert np.array_equal(BPX_DEFAULT_VECTORS["ponceau_fuchsin"], color.rgb_from_bpx[1]) + # scikit-image completes the third row by the cross product, as custom_od does. + assert np.array_equal( + color.rgb_from_bpx[2], np.cross(color.rgb_from_bpx[0], color.rgb_from_bpx[1]) + ) + + +@pytest.mark.parametrize("preset", list(DEFAULTS)) +def test_preset_records_hold_read_only_defaults_matching_their_stain_names(preset): + vectors, citation, target, exclusion = DEFAULTS[preset] + recipe = PRESETS[preset] + assert recipe.default_stain_vectors is vectors + assert recipe.default_stain_vectors_citation == citation + assert set(vectors) == {recipe.target_stain, recipe.exclusion_stain} == {target, exclusion} + assert recipe.needs == {"white"} and recipe.accepts == {"stain_vectors"} + with pytest.raises(TypeError): + vectors["eosin"] = (1.0, 0.0, 0.0) # type: ignore[index] + assert all(isinstance(v, tuple) for v in vectors.values()) + + +def test_psr_brightfield_has_no_default_and_needs_vectors(): + recipe = PRESETS["psr_brightfield"] + assert recipe.default_stain_vectors is None and recipe.default_stain_vectors_citation is None + assert recipe.needs == {"white", "stain_vectors"} + with pytest.raises(ValueError, match="needs stain_vectors"): + to_scalar(random_rgb(), axes="YXC", preset="psr_brightfield", white=WHITE) + + +def test_only_he_custom_and_trichrome_custom_have_defaults(): + with_defaults = [k for k, p in PRESETS.items() if p.default_stain_vectors is not None] + assert with_defaults == ["he_custom", "trichrome_custom"] + + +@pytest.mark.parametrize("preset", list(DEFAULTS)) +def test_white_is_still_required_with_default_vectors(preset): + with pytest.raises(ValueError, match="needs white"): + to_scalar(random_rgb(), axes="YXC", preset=preset) + + +@pytest.mark.parametrize("preset", list(DEFAULTS)) +def test_default_vs_user_vectors_provenance(preset): + vectors, citation, target, exclusion = DEFAULTS[preset] + rgb = random_rgb(3) + default = to_scalar(rgb, axes="YXC", preset=preset, white=WHITE) + assert default.provenance["stain_vectors_source"] == PRESET_DEFAULT + assert default.provenance["stain_vectors_citation"] == citation + np.testing.assert_allclose( + default.provenance["stain_matrix"], + unit_matrix(vectors[target], vectors[exclusion]), + rtol=0, + atol=1e-15, + ) + # The same numbers given by the user: same outputs, other source. + same = to_scalar(rgb, axes="YXC", preset=preset, white=WHITE, stain_vectors=dict(vectors)) + assert same.image.tobytes() == default.image.tobytes() + assert same.exclusion.tobytes() == default.exclusion.tobytes() + assert same.provenance["stain_vectors_source"] == USER + assert same.provenance["stain_vectors_citation"] is None + # A user override replaces the defaults. + override = {target: (0.2, 0.9, 0.4), exclusion: (0.7, 0.6, 0.3)} + user = to_scalar(rgb, axes="YXC", preset=preset, white=WHITE, stain_vectors=override) + np.testing.assert_allclose( + user.provenance["stain_matrix"], + unit_matrix(override[target], override[exclusion]), + rtol=0, + atol=1e-15, + ) + assert user.image.tobytes() != default.image.tobytes() + + +def test_he_custom_beer_lambert_round_trip_with_the_default_vectors(): + # Rows target (eosin), exclusion (haematoxylin), residual = their normalised cross product. + m = unit_matrix(EOSIN, HEMATOXYLIN) + assert dict(HE_DEFAULT_VECTORS) == {"hematoxylin": HEMATOXYLIN, "eosin": EOSIN} + amounts = concentrations(21, (24, 24), 2) + rgb = to_uint8(beer_lambert(amounts, m)) + result = to_scalar(rgb, axes="YXC", preset="he_custom", white=WHITE) + bound = concentration_bound(rgb, m) + assert (np.abs(result.image - amounts[..., 0]) <= bound[..., 0]).all() + assert (np.abs(result.exclusion - amounts[..., 1]) <= bound[..., 1]).all() + assert result.provenance["stain_vectors_source"] == PRESET_DEFAULT + assert result.provenance["residual_completed"] is True + + +def test_trichrome_custom_beer_lambert_round_trip_with_the_default_vectors(): + vectors = BPX_DEFAULT_VECTORS + m = unit_matrix(vectors["collagen"], vectors["ponceau_fuchsin"]) + amounts = concentrations(22, (24, 24), 2) + rgb = to_uint8(beer_lambert(amounts, m)) + result = to_scalar(rgb, axes="YXC", preset="trichrome_custom", white=WHITE) + bound = concentration_bound(rgb, m) + assert (np.abs(result.image - amounts[..., 0]) <= bound[..., 0]).all() + assert (np.abs(result.exclusion - amounts[..., 1]) <= bound[..., 1]).all() + assert result.provenance["stain_names"] == ("collagen", "ponceau_fuchsin", "residual") + + +def test_default_vectors_work_without_scikit_image(): + code = textwrap.dedent( + """ + import sys + sys.modules["skimage"] = None + sys.modules["skimage.color"] = None + import numpy as np + from morphobricks_normalization.stain import to_scalar + rgb = np.full((2, 2, 3), 120, np.uint8) + for preset in ("he_custom", "trichrome_custom"): + result = to_scalar(rgb, axes="YXC", preset=preset, white=(240, 235, 245)) + assert result.provenance["stain_vectors_source"] == "preset_default" + print("ok") + """ + ) + result = subprocess.run( + [sys.executable, "-c", code], capture_output=True, text=True, check=False + ) + assert result.returncode == 0, result.stderr + assert result.stdout.split() == ["ok"] + + +# --- white for the colour projections ----------------------------------------------------- + + +def test_rescale_to_white_rounds_half_to_even_and_clips(): + rgb = np.array([[[1, 3, 170], [200, 201, 255]]], dtype=np.uint8) + white = np.array([170.0, 170.0, 200.0]) + rescaled, clipped = od.rescale_to_white(rgb, white) + # 1 * 255 / 170 = 1.5 -> 2; 3 * 255 / 170 = 4.5 -> 4 (half to even); 170 * 255 / 200 + # = 216.75 -> 217; 200 * 255 / 170 = 300 -> 255; 201 * 255 / 170 = 301.5 -> 255; + # 255 * 255 / 200 = 325.125 -> 255. + assert rescaled.dtype == np.uint8 + assert rescaled.tolist() == [[[2, 4, 217], [255, 255, 255]]] + assert clipped == 0.5 + + +def test_rescale_to_white_counts_only_values_rounded_above_255(): + # 255 * 255 / 254.6 = 255.40 rounds to 255: above the white, but not clipped. + rgb = np.array([[[255, 255, 255], [254, 254, 254]]], dtype=np.uint8) + white = np.array([254.6, 254.6, 254.6]) + rescaled, clipped = od.rescale_to_white(rgb, white) + assert rescaled.tolist() == [[[255, 255, 255], [254, 254, 254]]] + assert clipped == 0.0 + assert od.above_white_fraction(rgb, white) == 0.5 + + +def test_rescale_to_white_255_is_the_identity_for_every_uint8_value(): + values = np.arange(256, dtype=np.uint8) + rgb = np.stack([values, values[::-1], values], axis=-1)[None] + rescaled, clipped = od.rescale_to_white(rgb, np.full(3, 255.0)) + assert np.array_equal(rescaled, rgb) and clipped == 0.0 + + +@pytest.mark.parametrize("preset", WHITE_PRESETS) +def test_white_rescaling_equals_the_manual_formula(preset): + color = pytest.importorskip("skimage.color") + rgb = random_rgb(4, shape=(11, 13)) + white = (170.0, 200.0, 230.5) + manual = np.clip(np.rint(rgb.astype(np.float64) * 255 / np.array(white)), 0, 255) + manual = manual.astype(np.uint8) + expected = { + "luminance_inverted": lambda x: 1.0 - color.rgb2gray(x), + "lab_a": lambda x: color.rgb2lab(x)[..., 1], + "psr_polarized": color.rgb2gray, + "trichrome_green": lambda x: -color.rgb2lab(x)[..., 1], + }[preset](manual) + result = to_scalar(rgb, axes="YXC", preset=preset, white=white, **extra(preset)) + assert result.image.tobytes() == expected.astype(np.float32).tobytes() + provenance = result.provenance + assert provenance["white"] == white and provenance["white_applied"] is True + assert provenance["white_rescale"] == WHITE_RESCALE + clipped = (np.rint(rgb.astype(np.float64) * 255 / np.array(white)) > 255).any(axis=-1) + assert provenance["white_clipped_fraction"] == float(np.mean(clipped)) > 0 + above = (rgb.astype(np.float64) > np.array(white)).any(axis=-1) + assert provenance["above_white_fraction"] == float(np.mean(above)) + + +@pytest.mark.parametrize("preset", WHITE_PRESETS) +def test_white_255_gives_bitwise_the_outputs_without_white(preset): + pytest.importorskip("skimage") + rgb = random_rgb(5, shape=(13, 17)) + without = to_scalar(rgb, axes="YXC", preset=preset, **extra(preset)) + with_white = to_scalar(rgb, axes="YXC", preset=preset, white=(255, 255, 255), **extra(preset)) + assert with_white.image.tobytes() == without.image.tobytes() + assert with_white.provenance["white_applied"] is True + assert with_white.provenance["white_clipped_fraction"] == 0.0 + assert with_white.provenance["above_white_fraction"] == 0.0 + + +@pytest.mark.parametrize("preset", WHITE_PRESETS) +def test_without_white_the_provenance_says_so(preset): + pytest.importorskip("skimage") + provenance = to_scalar(random_rgb(), axes="YXC", preset=preset, **extra(preset)).provenance + assert provenance["white"] is None and provenance["white_applied"] is False + assert provenance["white_rescale"] is None + assert provenance["white_clipped_fraction"] is None + assert provenance["above_white_fraction"] is None + + +@pytest.mark.parametrize("preset", P0_WHITE) +def test_a_background_at_the_declared_white_becomes_neutral(preset): + pytest.importorskip("skimage") + # A grey-blue background of (170, 180, 190) reads as dark and bluish against 255, and as + # white (no darkness, no redness) against its declared white. + rgb = np.full((3, 4, 3), (170, 180, 190), np.uint8) + without = to_scalar(rgb, axes="YXC", preset=preset).image + declared = to_scalar(rgb, axes="YXC", preset=preset, white=(170, 180, 190)).image + assert np.abs(declared).max() < 1e-2 + assert np.abs(without).min() > 0.2 + + +@pytest.mark.parametrize("preset", WHITE_PRESETS) +@pytest.mark.parametrize("white", [(240, 235), (0, 235, 245), (240, 235, 300), (np.nan, 1, 1)]) +def test_an_invalid_white_is_refused_for_the_projections(preset, white): + with pytest.raises(ValueError, match="white"): + to_scalar(random_rgb(), axes="YXC", preset=preset, white=white, **extra(preset)) + + +@pytest.mark.parametrize("preset", P0_WHITE) +def test_yxc_and_cyx_with_white_are_bitwise_identical(preset): + pytest.importorskip("skimage") + yxc = random_rgb(6) + cyx = np.moveaxis(yxc, 2, 0).copy() + before = yxc.copy() + a = to_scalar(yxc, axes="YXC", preset=preset, white=(200, 210, 220)) + b = to_scalar(cyx, axes="CYX", preset=preset, white=(200, 210, 220)) + assert a.image.tobytes() == b.image.tobytes() + assert np.array_equal(yxc, before) + assert {k for k in a.provenance if a.provenance[k] != b.provenance[k]} == {"axes"} From ee72fe37790461b2acfa32aa69f773d30b5bf9ae Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 12:24:01 +0000 Subject: [PATCH 08/12] =?UTF-8?q?Document=20default=20stain=20vectors=20an?= =?UTF-8?q?d=20the=20declared=20white=20(contract=20=C2=A73.4,=20plan=20N6?= =?UTF-8?q?9)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- .github/workflows/ci.yml | 2 +- CHANGELOG.md | 10 ++++++ docs/contract.md | 69 ++++++++++++++++++++++++++++++++-------- docs/plan.md | 20 ++++++++++++ 4 files changed, 86 insertions(+), 15 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f3ef0cf..723dc0a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -69,4 +69,4 @@ jobs: with: python-version: "3.12" - run: python -m pip install -e ".[dev,histology]" - - run: python -m pytest tests/test_stain.py tests/test_stain_catalogue.py tests/test_optical_density.py -rs + - run: python -m pytest tests/test_stain.py tests/test_stain_catalogue.py tests/test_stain_defaults.py tests/test_optical_density.py -rs diff --git a/CHANGELOG.md b/CHANGELOG.md index a223636..4b19e67 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,16 @@ The decisions are in [docs/plan.md](docs/plan.md) §12. `target_stain` and `exclusion_stain`; the provenance gains `modality` and `hints`. New helpers `white_from_region`, `stain_vectors_from_regions`; new pure modules `optical_density` and `scalar_inputs`; `PRESETS_BY_MODALITY`. +- **Default stain vectors and declared white** (issue #94, plan N69, contract §3.4): + `he_custom` (Ruifrok & Johnston H&E) and `trichrome_custom` (Methyl Blue / Ponceau Fuchsin of + `rgb_from_bpx`) use published default vectors when `stain_vectors` is not given; + `psr_brightfield` still needs them. New `StainPreset` fields `accepts`, + `default_stain_vectors`, `default_stain_vectors_citation`; provenance `stain_vectors_source`, + `stain_vectors_citation`. `trichrome_custom`'s exclusion stain is renamed `nuclei` -> + `ponceau_fuchsin`. `luminance_inverted`, `lab_a`, `psr_polarized`, `trichrome_green` accept an + optional `white` (rescale to `clip(rint(I * 255 / white), 0, 255)`, `white_clipped_fraction` + reported); without it outputs are bitwise unchanged. The registry moves to `stain_presets.py` + (re-exported by `stain.py`). ## v0.5.0 (2026-11-11) diff --git a/docs/contract.md b/docs/contract.md index f561c1b..db6a13b 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -488,7 +488,7 @@ y_c = (x_c − offset) / scale for c = R, G, B --- -### 3.4 Image-to-scalar presets (`stain.py`, issue #94, plan N67, N68) +### 3.4 Image-to-scalar presets (`stain.py`, issue #94, plan N67, N68, N69) A method family of its own, separate from §3.3 (which is unchanged): `to_scalar(pixels, *, axes, preset, white=None, stain_vectors=None, polarizer=None)` returns `ScalarResult(image, @@ -496,8 +496,8 @@ exclusion, provenance)`, one float32 `(Y, X)` map of one declared imaging unit, that need one channel (FiberFeatures first). normalization owns the conversion; consumers call it and never reimplement it. `rgb_to_scalar(pixels, *, axes, preset)` (N67) is `to_scalar` without the optional arguments. No file reading, tiling, segmentation or WSI handling is done -here. Modules: `stain.py` (API, registry), `optical_density.py` (OD maths, NumPy), -`scalar_inputs.py` (input checks). +here. Modules: `stain.py` (API), `stain_presets.py` (registry, hints, default vectors; re-exported +by `stain.py`), `optical_density.py` (OD maths, NumPy), `scalar_inputs.py` (input checks). **Catalogue** (every preset version 1, `polarity = "bright"`, `fixed_range = None`): @@ -505,16 +505,16 @@ here. Modules: `stain.py` (API, registry), `optical_density.py` (OD maths, NumPy |---|---|---|---|---|---|---|---|---| | `he_eosin` | brightfield | `fixed_deconvolution` | `separate_stains(rgb, hed_from_rgb)[..., 1]` | component 0 (haematoxylin) | — | fixed | `exclude_nuclei_rbc_empty` | eosin fibrillar texture, not collagen-specific | | `he_hematoxylin` | brightfield | `fixed_deconvolution` | the same call, component 0 | `None` | — | fixed | `none` | haematoxylin, mainly nuclei | -| `luminance_inverted` | brightfield | `luminance` | `1 - rgb2gray(rgb)` (baseline) | `None` | — | fixed | `exclude_nuclei_rbc_empty` | inverted brightness, nuclei-dominated | -| `lab_a` | brightfield | `cielab` | `rgb2lab(rgb)[..., 1]` (a*, signed; comparator) | `None` | — | fixed | `exclude_nuclei_rbc_empty` | a* redness, low contrast | +| `luminance_inverted` | brightfield | `luminance` | `1 - rgb2gray(rgb)` (baseline) | `None` | — (white optional) | fixed | `exclude_nuclei_rbc_empty` | inverted brightness, nuclei-dominated | +| `lab_a` | brightfield | `cielab` | `rgb2lab(rgb)[..., 1]` (a*, signed; comparator) | `None` | — (white optional) | fixed | `exclude_nuclei_rbc_empty` | a* redness, low contrast | | `fluorescence` | fluorescence | `passthrough` | the `"YX"` plane as float32 | `None` | — | per_image | `optional` | the labelled structure | | `shg` | shg | `passthrough` | the `"YX"` plane as float32 | `None` | — | fixed | `none` | fibrillar collagen (mainly type I) | | `psr_brightfield` | brightfield | `custom_od` | concentration of `"psr"` | concentration of `"counterstain"` | white, stain_vectors | fixed | `exclude_counterstain` | collagen | -| `psr_polarized` | polarized | `luminance` | `rgb2gray(rgb)`, **not inverted** | `None` | polarizer | fixed | `background` | birefringent packed collagen | +| `psr_polarized` | polarized | `luminance` | `rgb2gray(rgb)`, **not inverted** | `None` | polarizer (white optional) | fixed | `background` | birefringent packed collagen | | `trichrome_blue` | brightfield | `fixed_deconvolution` | `separate_stains(rgb, bpx_from_rgb)[..., 0]` (Methyl Blue) | `None` | — | fixed | `exclude_nuclei` | collagen (methyl blue) | -| `trichrome_green` | brightfield | `cielab` | `-rgb2lab(rgb)[..., 1]` (minus a*, signed) | `None` | — | fixed | `exclude_nuclei` | collagen (green) | -| `he_custom` | brightfield | `custom_od` | concentration of `"eosin"` | concentration of `"hematoxylin"` | white, stain_vectors | fixed | `exclude_nuclei_rbc_empty` | eosin fibrillar texture, not collagen-specific | -| `trichrome_custom` | brightfield | `custom_od` | concentration of `"collagen"` | concentration of `"nuclei"` | white, stain_vectors | fixed | `exclude_nuclei` | collagen | +| `trichrome_green` | brightfield | `cielab` | `-rgb2lab(rgb)[..., 1]` (minus a*, signed) | `None` | — (white optional) | fixed | `exclude_nuclei` | collagen (green) | +| `he_custom` | brightfield | `custom_od` | concentration of `"eosin"` | concentration of `"hematoxylin"` | white (stain_vectors optional: defaults) | fixed | `exclude_nuclei_rbc_empty` | eosin fibrillar texture, not collagen-specific | +| `trichrome_custom` | brightfield | `custom_od` | concentration of `"collagen"` | concentration of `"ponceau_fuchsin"` | white (stain_vectors optional: defaults) | fixed | `exclude_nuclei` | collagen | `stain` is `"he"`, `"picrosirius_red"`, `"trichrome_blue"`, `"trichrome_green"` (the values proposed in MorphoBricks-Spec #36) or `None` (colour projections, passthrough, @@ -529,9 +529,10 @@ amount), `exclude_counterstain`, `background`, `none`, `optional`) and `measures` (short text). They are recommendations for consumers (FiberFeatures #214), **never applied here**. -**Arguments.** A preset's `needs` lists the optional arguments it requires: each one missing is -refused, and each argument it does not use (`white` for `lab_a`, `polarizer` for `he_eosin`) is -refused, never ignored. `polarizer` is `"circular"` or `"linear"`; the provenance hints of +**Arguments.** A preset's `needs` lists the optional arguments it requires (each one missing is +refused) and its `accepts` those it uses when given; every other argument (`white` for +`he_eosin` or the passthrough presets, `polarizer` for `lab_a`, `stain_vectors` for +`trichrome_blue`) is refused, never ignored. `polarizer` is `"circular"` or `"linear"`; the provenance hints of `psr_polarized` record `orientation_valid = (polarizer == "circular")`: a linear polarizer darkens fibres at some angles to its axis, so consumers must refuse orientation features then. Hue is not used (later, as a separate feature). @@ -587,13 +588,53 @@ told apart); it is `None` for every other method. background region; an empty or wrongly shaped region, or a median of 0 or 255 (dark or saturated, so the true white is unknown), is refused. +**Default stain vectors** (N69, owner decision of 13 November 2026). Where published values +exist, a `custom_od` preset carries them as read-only data (`default_stain_vectors`, +`default_stain_vectors_citation`); they are used when `stain_vectors` is not given, and a +given `stain_vectors` (for example from `stain_vectors_from_regions`) overrides them. They are +checked like user vectors (unit length, cross-product residual, angle and condition limits). +The numbers are copied into the code (no scikit-image needed) and tested equal to scikit-image's +rows when it is installed. `white` stays required: it is never assumed to be 255. + +| `preset` | target | exclusion | R, G, B OD (before normalisation) | source | +|---|---|---|---|---| +| `he_custom` | `eosin` | `hematoxylin` | eosin (0.07, 0.99, 0.11); haematoxylin (0.65, 0.70, 0.29) | Ruifrok & Johnston, *Anal Quant Cytol Histol* 23:291-299 (2001); `skimage.color.rgb_from_hed` rows 1 and 0 | +| `trichrome_custom` | `collagen` (Methyl Blue) | `ponceau_fuchsin` | methyl blue (0.7995107, 0.5913521, 0.10528667); ponceau fuchsin (0.09997159, 0.73738605, 0.6680326) | Masson Trichrome vectors of Landini's Colour Deconvolution (Fiji), Ruifrok & Johnston method; `skimage.color.rgb_from_bpx` rows 0 and 1 | +| `psr_brightfield` | `psr` | `counterstain` | none: no reliable published picrosirius red vectors; `stain_vectors` is required | — | + +The residual is the normalised cross product of the two (scikit-image's `rgb_from_bpx` is +completed the same way; `rgb_from_hed` has DAB as its third row, so `he_custom` with defaults is +not `he_eosin`: other residual, OD in log10 against the declared white, no clipping). +`trichrome_custom`'s exclusion is Ponceau Fuchsin (cytoplasm, muscle, the red counterstain), not +nuclei: the published Masson matrix has no nuclear stain (iron haematoxylin is near black and +absorbs in all channels). The mask hint stays `exclude_nuclei`: a nuclear mask is built in +preprocessing, not taken from this exclusion map. User vectors for `trichrome_custom` name +`collagen` and `ponceau_fuchsin`. Green trichrome stays user-measured (no default). + +**Declared white for the colour projections** (N69). `luminance_inverted`, `lab_a`, +`psr_polarized` and `trichrome_green` accept an optional `white` (three finite values in +`(0, 255]`, checked as for `custom_od`). When given, the RGB is rescaled per channel before the +scikit-image call to `clip(rint(I * 255 / white), 0, 255)` as `uint8`, computed in float64 in +that order; `rint` rounds half to even, and a value above 255 after rounding (brighter than the +white by more than half a grey level) is set to 255. `white_clipped_fraction` is the fraction of +pixels with at least one such clipped channel; `above_white_fraction` (as for `custom_od`) the +fraction with a channel strictly above the white, which includes pixels within half a grey level +of it that round to 255 without clipping. `white = (255, 255, 255)` is the identity, so the +outputs are bitwise those without `white`; without `white` nothing changes. The +`fixed_deconvolution` presets (`he_eosin`, `he_hematoxylin`, `trichrome_blue`) do not accept +`white` (scikit-image's `separate_stains` fixes white at 255): use `he_custom` or +`trichrome_custom` with a declared white instead. + **Provenance** (read-only mapping): `preset`, `preset_version`, `stain`, `modality`, `method`, `matrix`, `target_component`, `exclusion_component`, `polarity`, `hints` (read-only), `input_dtype`, `axes`, `colour_space` (`"srgb"`, `None` for passthrough), `skimage_version` (`None` when scikit-image is not called), `zero_fraction`; plus, for `custom_od`: `white`, -`stain_names`, `stain_matrix` (unit rows), `condition_number`, `residual_completed`, +`stain_vectors_source` (`"preset_default"` or `"user"`), `stain_vectors_citation` (the +default's citation, `None` for user vectors), `stain_names`, `stain_matrix` (unit rows), `condition_number`, `residual_completed`, `od_definition`, `floored_fraction`, `negative_fraction`, `above_white_fraction`; for -`psr_polarized`: `polarizer` (and `hints["orientation_valid"]`); for passthrough: +`luminance` and `cielab`: `white` (`None` when not given), `white_applied`, `white_rescale` (the +formula, or `None`), `white_clipped_fraction` and `above_white_fraction` (`None` without white); +for `psr_polarized`: `polarizer` (and `hints["orientation_valid"]`); for passthrough: `float32_exact`. **Dependency.** scikit-image is the optional extra `histology` diff --git a/docs/plan.md b/docs/plan.md index 1e34773..e717346 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -567,6 +567,26 @@ Added with the project owner on 11 November 2026, for Stage 3 (S3.2, pull reques which includes pixels that were never clipped); masked arrays (`np.ma.MaskedArray`) and images with no pixels are refused for every preset; `polarity = "bright"` means the structure has higher values than its surroundings (the signal direction), not more stain. +63. **N69** (13 November 2026, issue #94, owner decisions "default stain vectors per preset") + Custom-vector presets carry published **default stain vectors** where they exist, as frozen + read-only data in the preset record (`default_stain_vectors`, with a citation): `he_custom` + the Ruifrok & Johnston haematoxylin and eosin rows of `rgb_from_hed`, `trichrome_custom` the + Methyl Blue (`collagen`) and Ponceau Fuchsin (`ponceau_fuchsin`) rows of `rgb_from_bpx` (the + values of Fiji's Masson Trichrome). The numbers are copied into the code and tested equal to + scikit-image's. `psr_brightfield` and green trichrome stay user-measured (no reliable + published vectors). A given `stain_vectors` overrides the default; the provenance records + `stain_vectors_source` (`"preset_default"` with the citation, or `"user"`). User vectors are + measured in preprocessing from marked regions with a declared white and registered with the + experiment description (MorphoBricks-Spec #36). `white` stays required for every + `custom_od` preset. `trichrome_custom`'s exclusion stain becomes `ponceau_fuchsin` (it was + `nuclei`): the published matrix has no nuclear stain, so the exclusion map is the red + counterstain and the `exclude_nuclei` mask hint is met by a mask built in preprocessing. + `luminance_inverted`, `lab_a`, `psr_polarized` and `trichrome_green` accept an optional + declared **white**: the RGB is rescaled to `clip(rint(I * 255 / white), 0, 255)` as uint8 + before scikit-image, `white_clipped_fraction` is reported, and without `white` (or with + 255) outputs are bitwise unchanged. The `fixed_deconvolution` presets do not take `white` + (`custom_od` is the path for deconvolution against a declared white). Masks are built in + preprocessing, never in an analysis brick (contract §3.4). Implementation choices of S3.2 where the Spec leaves the detail open (open questions 1-7 of pull request #90), each taken on the side the Spec states or the safer one, and listed in that From 734c808e2681c5c0669e2719c68b1fd0035143b3 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 14:49:55 +0000 Subject: [PATCH 09/12] Name each stain preset's image CellProfiler style (#94) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit StainPreset.image_name (recorded in the provenance as image_name) names the image a preset makes, as CellProfiler's UnmixColors names its outputs: he_eosin -> Eosin, he_hematoxylin -> Hematoxylin, trichrome_blue -> MethylBlue, psr_brightfield -> SiriusRed, and so on; unique tokens, which consumers use as the derived channel's name in feature columns (owner decision: CellProfiler-style names). Contract §3.4. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- CHANGELOG.md | 3 +++ docs/contract.md | 10 ++++++++++ src/morphobricks_normalization/stain.py | 1 + .../stain_presets.py | 19 +++++++++++++++++++ tests/test_stain.py | 1 + tests/test_stain_catalogue.py | 14 ++++++++++++++ 6 files changed, 48 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b19e67..229a366 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,9 @@ The decisions are in [docs/plan.md](docs/plan.md) §12. `target_stain` and `exclusion_stain`; the provenance gains `modality` and `hints`. New helpers `white_from_region`, `stain_vectors_from_regions`; new pure modules `optical_density` and `scalar_inputs`; `PRESETS_BY_MODALITY`. +- **CellProfiler-style image names** (issue #94, contract §3.4): each stain preset names its image + (`StainPreset.image_name`, e.g. `he_eosin` -> `Eosin`, `trichrome_blue` -> `MethylBlue`), recorded in the + provenance; consumers use it as the derived channel's name in feature columns. - **Default stain vectors and declared white** (issue #94, plan N69, contract §3.4): `he_custom` (Ruifrok & Johnston H&E) and `trichrome_custom` (Methyl Blue / Ponceau Fuchsin of `rgb_from_bpx`) use published default vectors when `stain_vectors` is not given; diff --git a/docs/contract.md b/docs/contract.md index db6a13b..d68527b 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -499,6 +499,16 @@ without the optional arguments. No file reading, tiling, segmentation or WSI han here. Modules: `stain.py` (API), `stain_presets.py` (registry, hints, default vectors; re-exported by `stain.py`), `optical_density.py` (OD maths, NumPy), `scalar_inputs.py` (input checks). +**Image names.** Every preset names the image it makes, CellProfiler style (`StainPreset.image_name`, +recorded in the provenance as `image_name`), as CellProfiler's `UnmixColors` names its outputs: +`he_eosin` → `Eosin`, `he_hematoxylin` → `Hematoxylin`, `luminance_inverted` → `InvertedLuminance`, +`lab_a` → `LabA`, `fluorescence` → `Fluorescence`, `shg` → `SHG`, `psr_brightfield` → `SiriusRed`, +`psr_polarized` → `SiriusRedPolarized`, `trichrome_blue` → `MethylBlue`, `trichrome_green` → +`TrichromeGreen`, `he_custom` → `EosinCustom`, `trichrome_custom` → `TrichromeCustom`. Names are +unique tokens; consumers use them as the derived image's channel name in feature columns +(`Eosin__Cells__ridge_support_mean`). A name says which image was measured, not what it means +biologically (`Eosin` is eosin-derived texture, not collagen). + **Catalogue** (every preset version 1, `polarity = "bright"`, `fixed_range = None`): | `preset` | modality | method | `image` | `exclusion` | needs | scaling | mask | measures | diff --git a/src/morphobricks_normalization/stain.py b/src/morphobricks_normalization/stain.py index 7ee6218..859a1d5 100644 --- a/src/morphobricks_normalization/stain.py +++ b/src/morphobricks_normalization/stain.py @@ -216,6 +216,7 @@ def _provenance( record = { "preset": recipe.id, "preset_version": recipe.version, + "image_name": recipe.image_name, "stain": recipe.stain, "modality": recipe.modality, "method": recipe.method, diff --git a/src/morphobricks_normalization/stain_presets.py b/src/morphobricks_normalization/stain_presets.py index 321f92e..8f00aca 100644 --- a/src/morphobricks_normalization/stain_presets.py +++ b/src/morphobricks_normalization/stain_presets.py @@ -100,6 +100,12 @@ class StainPreset: when ``stain_vectors`` is not given, or ``None`` (then ``stain_vectors`` is needed). default_stain_vectors_citation Where the default vectors come from, recorded in the provenance. + image_name + The CellProfiler-style name of the image the preset makes (as CellProfiler's + ``UnmixColors`` names its outputs ``Hematoxylin``, ``Eosin``, ``Methyl blue``): a token of + letters and digits, unique among the presets, which consumers use as the derived + image's channel name in feature columns (``Eosin__Cells__ridge_support_mean``). It names + the image, not what it measures: ``Eosin`` is eosin-derived texture, not collagen. """ id: str @@ -122,6 +128,7 @@ class StainPreset: default=None, hash=False ) default_stain_vectors_citation: str | None = None + image_name: str | None = None def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: @@ -141,6 +148,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: for preset in ( StainPreset( "he_eosin", + image_name="Eosin", stain="he", matrix="hed_from_rgb", target_component=1, @@ -149,6 +157,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "he_hematoxylin", + image_name="Hematoxylin", stain="he", matrix="hed_from_rgb", target_component=0, @@ -156,6 +165,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "luminance_inverted", + image_name="InvertedLuminance", method="luminance", inverted=True, accepts=_WHITE, @@ -163,6 +173,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "lab_a", + image_name="LabA", method="cielab", target_component=1, accepts=_WHITE, @@ -170,18 +181,21 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "fluorescence", + image_name="Fluorescence", method="passthrough", modality="fluorescence", hints=_hints("per_image", "optional", "the labelled structure"), ), StainPreset( "shg", + image_name="SHG", method="passthrough", modality="shg", hints=_hints("fixed", "none", "fibrillar collagen (mainly type I)"), ), StainPreset( "psr_brightfield", + image_name="SiriusRed", stain="picrosirius_red", method="custom_od", needs=_CUSTOM, @@ -191,6 +205,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "psr_polarized", + image_name="SiriusRedPolarized", stain="picrosirius_red", method="luminance", modality="polarized", @@ -200,6 +215,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "trichrome_blue", + image_name="MethylBlue", stain="trichrome_blue", matrix="bpx_from_rgb", target_component=0, @@ -207,6 +223,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "trichrome_green", + image_name="TrichromeGreen", stain="trichrome_green", method="cielab", target_component=1, @@ -216,6 +233,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "he_custom", + image_name="EosinCustom", stain="he", method="custom_od", needs=_WHITE, @@ -228,6 +246,7 @@ def _hints(scaling: str, mask: str, measures: str) -> Mapping[str, str]: ), StainPreset( "trichrome_custom", + image_name="TrichromeCustom", method="custom_od", needs=_WHITE, accepts=frozenset({"stain_vectors"}), diff --git a/tests/test_stain.py b/tests/test_stain.py index b2accac..02d6faf 100644 --- a/tests/test_stain.py +++ b/tests/test_stain.py @@ -231,6 +231,7 @@ def test_provenance_contents_and_zero_fraction(): expected = { "preset": "he_eosin", "preset_version": 1, + "image_name": "Eosin", "stain": "he", "modality": "brightfield", "method": "fixed_deconvolution", diff --git a/tests/test_stain_catalogue.py b/tests/test_stain_catalogue.py index 8e9516a..0e49cac 100644 --- a/tests/test_stain_catalogue.py +++ b/tests/test_stain_catalogue.py @@ -539,3 +539,17 @@ def test_fixed_reference_values_of_the_colour_projections(): eosin = to_scalar(rgb, axes="YXC", preset="he_eosin") assert eosin.image[0, 0] == 0.0 # white: no stain assert eosin.provenance["zero_fraction"] >= 1 / 3 + + +def test_every_preset_names_its_image_cellprofiler_style(): + """Each preset names its image as CellProfiler's UnmixColors does: unique tokens.""" + import re + + names = {preset_id: preset.image_name for preset_id, preset in PRESETS.items()} + assert all( + isinstance(name, str) and re.fullmatch(r"[A-Z][A-Za-z0-9]*", name) + for name in names.values() + ), names + assert len(set(names.values())) == len(names) + assert names["he_eosin"] == "Eosin" and names["he_hematoxylin"] == "Hematoxylin" + assert names["trichrome_blue"] == "MethylBlue" From 8cc45da642fa27a9c1665dc3adb1df9f4fdc5919 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 15:03:14 +0000 Subject: [PATCH 10/12] Refuse masked auxiliary inputs and floored reference regions (review of 3a7af3e) - check_white and stain vectors refuse a masked array, or a list holding a masked value, before np.asarray could drop the mask (an all-false mask too). - Region helpers refuse a masked region array. - stain_vectors_from_regions refuses a reference region holding a pixel with a channel at 0 (its OD would be taken at the floor and bend the vector); zeros outside the regions are allowed. - MAX_CONDITION documented as a heuristic norm-wise threshold, not a per-map error guarantee; psr_polarized's green-weighted greyscale bias documented; white_from_region's median-255 refusal says a declared white of 255 is accepted. - Regression tests through to_scalar and the public region helpers. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- CHANGELOG.md | 6 ++ docs/contract.md | 4 ++ .../optical_density.py | 58 +++++++++++++++---- tests/test_stain_defaults.py | 57 ++++++++++++++++++ 4 files changed, 113 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 229a366..fc61349 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,12 @@ The decisions are in [docs/plan.md](docs/plan.md) §12. `target_stain` and `exclusion_stain`; the provenance gains `modality` and `hints`. New helpers `white_from_region`, `stain_vectors_from_regions`; new pure modules `optical_density` and `scalar_inputs`; `PRESETS_BY_MODALITY`. +- **Review fixes (ChatGPT review of `3a7af3e`)**: a masked `white`, a masked stain-vector component (or a masked + value inside a list) and a masked region are refused (`TypeError`), never unmasked by `np.asarray`; + `stain_vectors_from_regions` refuses a reference region holding a pixel with a channel at 0 (its optical density + would be taken at the floor), while zeros outside the regions are allowed; the condition-number limit is + documented as a heuristic norm-wise threshold; the `psr_polarized` green-weighted greyscale bias is documented; + the median-255 refusal of `white_from_region` says a declared white of 255 is still accepted. - **CellProfiler-style image names** (issue #94, contract §3.4): each stain preset names its image (`StainPreset.image_name`, e.g. `he_eosin` -> `Eosin`, `trichrome_blue` -> `MethylBlue`), recorded in the provenance; consumers use it as the derived channel's name in feature columns. diff --git a/docs/contract.md b/docs/contract.md index d68527b..a0c2aff 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -545,6 +545,10 @@ refused) and its `accepts` those it uses when given; every other argument (`whit `trichrome_blue`) is refused, never ignored. `polarizer` is `"circular"` or `"linear"`; the provenance hints of `psr_polarized` record `orientation_valid = (polarizer == "circular")`: a linear polarizer darkens fibres at some angles to its axis, so consumers must refuse orientation features then. +`psr_polarized` is scikit-image's `rgb2gray`, a weighted sum `0.2125 R + 0.7154 G + 0.0721 B` +dominated by green: birefringent fibres of different colours (red-orange thick, green-yellow thin +fibres) that look equally bright can get different values, which can bias detection towards some +fibre colours even with a circular polarizer. Hue is not used (later, as a separate feature). **RGB inputs (refused, never repaired).** `axes` is exactly `"YXC"` or `"CYX"`, with no default; diff --git a/src/morphobricks_normalization/optical_density.py b/src/morphobricks_normalization/optical_density.py index 78ae447..174adce 100644 --- a/src/morphobricks_normalization/optical_density.py +++ b/src/morphobricks_normalization/optical_density.py @@ -44,13 +44,17 @@ vectors and :data:`MAX_CONDITION` the one for three nearly coplanar vectors. """ MAX_CONDITION = 50.0 -"""Largest accepted 2-norm condition number of the unit stain matrix. - -The relative error of the concentrations is at most ``cond(M)`` times the relative error of -the optical densities. One 8-bit grey level at a mid-range value (``I`` about 80, ``OD`` about -0.5) is an optical-density error of about 1 %, so 50 bounds the worst-case amplification at -about 50 %: beyond it the unmixed maps are dominated by quantisation noise. Ruifrok-Johnston -H&E-DAB has 3.8 and scikit-image's Methyl Blue/Ponceau Fuchsin matrix 2.0. +"""Largest accepted 2-norm condition number of the unit stain matrix: a heuristic threshold. + +The norm-wise relative error of a pixel's concentration vector is at most ``cond(M)`` times the +norm-wise relative error of its optical-density vector. One 8-bit grey level at a mid-range +value (``I`` about 80, ``OD`` about 0.5) is an optical-density error of about 1 %, so at 50 the +worst-case norm-wise amplification is of order 50 %. This is a provisional refusal of nearly +singular matrices, not a guarantee on any one stain map or pixel (the error of a single +component can be smaller or larger than the norm-wise bound suggests); concentration errors are +to be benchmarked on realistic noise and backgrounds. The measured condition number is recorded +in the provenance. Ruifrok-Johnston H&E-DAB has 3.8 and scikit-image's Methyl Blue/Ponceau +Fuchsin matrix 2.0. """ MIN_REGION_OD = 0.01 """Smallest optical-density norm of a region's mean for it to define a stain vector. @@ -86,6 +90,16 @@ def array(self) -> np.ndarray: return np.array(self.rows, dtype=np.float64) +def _refuse_masked(value: object, what: str) -> None: + """Refuse a masked array, or a sequence holding a masked value: ``np.asarray`` drops masks.""" + items = value if isinstance(value, (list, tuple)) else () + if isinstance(value, np.ma.MaskedArray) or any(isinstance(v, np.ma.MaskedArray) for v in items): + raise TypeError( + f"{what} is a masked array (or holds a masked value): its mask would be dropped " + "silently. Pass plain numbers." + ) + + def check_white(white: object) -> np.ndarray: """Return a declared white as a float64 ``(3,)`` array after checking it. @@ -104,6 +118,7 @@ def check_white(white: object) -> np.ndarray: ValueError Not three finite numbers in ``(0, 255]``. """ + _refuse_masked(white, "white") values = np.asarray(white) if values.dtype.kind not in "iuf" or values.shape != (CHANNELS,): raise ValueError( @@ -121,6 +136,7 @@ def check_white(white: object) -> np.ndarray: def _unit_vector(name: str, values: object) -> np.ndarray: """Return the stain vector ``name`` normalised to unit length, or refuse it.""" + _refuse_masked(values, f"Stain vector {name!r}") vector = np.asarray(values) if vector.dtype.kind not in "iuf" or vector.shape != (CHANNELS,): raise ValueError( @@ -246,7 +262,12 @@ def concentrations(od: np.ndarray, matrix: StainMatrix) -> np.ndarray: def _check_region(region: object, shape: tuple[int, int], name: str) -> np.ndarray: - """Return a boolean ``(Y, X)`` region, or refuse a wrong type, shape or empty region.""" + """Return a boolean ``(Y, X)`` region; refuse a masked, mistyped, misshaped or empty one.""" + if isinstance(region, np.ma.MaskedArray): + raise TypeError( + f"Region {name!r} is a masked array: which pixels belong to it would be ambiguous. " + "Pass a plain boolean array." + ) if not isinstance(region, np.ndarray) or region.dtype != np.bool_: raise TypeError(f"Region {name!r} must be a boolean NumPy array (Y, X).") if region.shape != shape: @@ -272,8 +293,9 @@ def region_white(rgb: np.ndarray, region: object) -> tuple[float, float, float]: if (median <= 0).any() or (median >= WHITE_MAX).any(): raise ValueError( f"The background median {tuple(median.tolist())} has a channel at 0 or 255 " - "(dark or saturated); the white cannot be measured from it. Lower the exposure " - "or choose another background region." + "(dark or saturated): an unsaturated white cannot be inferred from this region. " + "Lower the exposure or choose another background region; a white declared " + "explicitly (for example 255 for a scanner that maps background to 255) is accepted." ) return tuple(float(v) for v in median) @@ -283,11 +305,16 @@ def region_vectors( ) -> Mapping[str, tuple[float, float, float]]: """Return the unit mean optical-density vector of each labelled region. + A selected pixel with a channel at 0 would enter at the floor (:data:`FLOOR`, intensity 1) + and bend the region's direction, so such a region is refused; zeros outside the regions do + not matter. + Raises ------ ValueError - No regions, an empty or wrongly shaped region, or a region whose mean optical density - is below :data:`MIN_REGION_OD` (background, not stain). + No regions, an empty or wrongly shaped region, a region holding a pixel with a channel at + 0, or a region whose mean optical density is below :data:`MIN_REGION_OD` (background, + not stain). TypeError ``regions`` is not a mapping, or a region is not a boolean array. """ @@ -297,6 +324,13 @@ def region_vectors( vectors = {} for name, region in regions.items(): mask = _check_region(region, rgb.shape[:2], name) + floored = int((rgb[mask] < FLOOR).any(axis=-1).sum()) + if floored: + raise ValueError( + f"Region {name!r} holds {floored} pixel(s) with a channel at 0: their optical " + "density would be taken at the floor (intensity 1) and bend the stain vector. " + "Choose a reference region without saturated-dark pixels." + ) mean = od[mask].mean(axis=0) norm = float(np.linalg.norm(mean)) if norm < MIN_REGION_OD: diff --git a/tests/test_stain_defaults.py b/tests/test_stain_defaults.py index ae3818d..bc1605a 100644 --- a/tests/test_stain_defaults.py +++ b/tests/test_stain_defaults.py @@ -286,3 +286,60 @@ def test_yxc_and_cyx_with_white_are_bitwise_identical(preset): assert a.image.tobytes() == b.image.tobytes() assert np.array_equal(yxc, before) assert {k for k in a.provenance if a.provenance[k] != b.provenance[k]} == {"axes"} + + +# --- Review of 3a7af3e: masked auxiliary inputs and floored reference regions --------- + + +def _plain_rgb() -> np.ndarray: + rng = np.random.default_rng(3) + return rng.integers(60, 230, (8, 8, 3), dtype=np.uint8) + + +@pytest.mark.parametrize( + "white", + [ + np.ma.array([240.0, 235.0, 245.0], mask=[True, False, False]), + np.ma.array([240.0, 235.0, 245.0], mask=False), # an all-false mask is refused too + [240.0, np.ma.masked, 245.0], + ], +) +def test_a_masked_white_is_refused_through_to_scalar(white): + with pytest.raises(TypeError, match="white is a masked array"): + to_scalar(_plain_rgb(), axes="YXC", preset="he_custom", white=white) + + +def test_a_masked_stain_vector_component_is_refused_through_to_scalar(): + vectors = { + "eosin": np.ma.array([0.07, 0.99, 0.11], mask=[False, True, False]), + "hematoxylin": (0.65, 0.70, 0.29), + } + with pytest.raises(TypeError, match="Stain vector 'eosin' is a masked array"): + to_scalar(_plain_rgb(), axes="YXC", preset="he_custom", white=WHITE, stain_vectors=vectors) + + +def test_masked_regions_are_refused_by_both_region_helpers(): + from morphobricks_normalization.stain import stain_vectors_from_regions, white_from_region + + region = np.zeros((8, 8), bool) + region[:2, :2] = True + masked = np.ma.array(region, mask=np.zeros((8, 8), bool)) + masked.mask[0, 0] = True # one selected pixel marked invalid, three valid ones + with pytest.raises(TypeError, match="masked array"): + white_from_region(_plain_rgb(), axes="YXC", region=masked) + with pytest.raises(TypeError, match="masked array"): + stain_vectors_from_regions(_plain_rgb(), axes="YXC", white=WHITE, regions={"eosin": masked}) + + +def test_a_reference_region_with_a_floored_pixel_is_refused_and_zeros_outside_do_not_matter(): + from morphobricks_normalization.stain import stain_vectors_from_regions + + rgb = _plain_rgb() + region = np.zeros((8, 8), bool) + region[2:5, 2:5] = True + rgb[7, 7, 1] = 0 # outside the region: allowed + vectors = stain_vectors_from_regions(rgb, axes="YXC", white=WHITE, regions={"eosin": region}) + assert set(vectors) == {"eosin"} + rgb[3, 3, 0] = 0 # inside the region: its OD would be taken at the floor + with pytest.raises(ValueError, match="holds 1 pixel"): + stain_vectors_from_regions(rgb, axes="YXC", white=WHITE, regions={"eosin": region}) From 2506e51ac2582964e187114315d493761684d146 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 15:28:55 +0000 Subject: [PATCH 11/12] State he_custom as a two-stain model and test it without scikit-image (review of 8cc45da) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Contract §3.4: he_custom is deliberately not the three-stain HED conversion (cross-product residual instead of DAB, log10 OD against the declared white, no clipping); the wrong "completion of rgb_from_hed" sentence is corrected; the declared-white path's uint8 re-quantisation (at most half a grey level) is documented. - Tests: literal published vectors independent of scikit-image; he_custom on a reference pixel computed by hand; its difference from scikit-image's HED eosin. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- CHANGELOG.md | 6 +++++ docs/contract.md | 18 ++++++++++--- tests/test_stain_defaults.py | 50 ++++++++++++++++++++++++++++++++++++ 3 files changed, 70 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fc61349..7ccc81c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,12 @@ The decisions are in [docs/plan.md](docs/plan.md) §12. `target_stain` and `exclusion_stain`; the provenance gains `modality` and `hints`. New helpers `white_from_region`, `stain_vectors_from_regions`; new pure modules `optical_density` and `scalar_inputs`; `PRESETS_BY_MODALITY`. +- **`he_custom` stated as a two-stain model** (review of `8cc45da`): contract §3.4 says plainly that it is not the + three-stain HED conversion (cross-product residual instead of DAB, log10 OD against the declared white, no + clipping); the wrong "completion of `rgb_from_hed`" sentence is corrected; the uint8 re-quantisation of the + declared-white path is documented (at most half a grey level per channel). New tests: the literal published + vectors independent of scikit-image, `he_custom` on a reference pixel computed by hand, and its difference from + scikit-image's HED eosin. - **Review fixes (ChatGPT review of `3a7af3e`)**: a masked `white`, a masked stain-vector component (or a masked value inside a list) and a masked region are refused (`TypeError`), never unmasked by `np.asarray`; `stain_vectors_from_regions` refuses a reference region holding a pixel with a channel at 0 (its optical density diff --git a/docs/contract.md b/docs/contract.md index a0c2aff..b6e5e36 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -585,7 +585,8 @@ told apart); it is `None` for every other method. the preset's target and exclusion names, plus optionally `"residual"`. Missing or other names are refused. Each vector is normalised to unit length (non-finite values, a wrong length or a norm below 1e-6 are refused). Without `"residual"` the third row is the normalised cross - product of the first two (scikit-image's completion of `rgb_from_hed`). + product of the first two (as scikit-image completes `rgb_from_bpx`; `rgb_from_hed` has an explicit + DAB third row instead). - Rows of `M`: target, exclusion, residual. Two vectors closer than 5 degrees (either direction) are refused as collinear; for two vectors completed by the cross product `cond(M) = cot(angle / 2)`, so 5 degrees gives 22.9. A 2-norm condition number above 50 is @@ -617,8 +618,14 @@ rows when it is installed. `white` stays required: it is never assumed to be 255 | `psr_brightfield` | `psr` | `counterstain` | none: no reliable published picrosirius red vectors; `stain_vectors` is required | — | The residual is the normalised cross product of the two (scikit-image's `rgb_from_bpx` is -completed the same way; `rgb_from_hed` has DAB as its third row, so `he_custom` with defaults is -not `he_eosin`: other residual, OD in log10 against the declared white, no clipping). +completed the same way). **`he_custom` is a two-stain model, deliberately not the three-stain HED +conversion**: it shares Ruifrok & Johnston's haematoxylin and eosin vectors, but its third row is +the cross product, not DAB (0.27, 0.57, 0.78); its optical density is log10 against the declared +white, not scikit-image's natural log divided by `log(1e-6)` with white 255; and negative +concentrations are kept, not clipped. Its eosin map therefore differs from `he_eosin`'s in scale +and in value (tested on a reference pixel: 0.294 against scikit-image's 0.046). Use `he_eosin` for +the published HED conversion, `he_custom` for a declared white or measured vectors. The literal +published values are tested independently of scikit-image, as well as against it. `trichrome_custom`'s exclusion is Ponceau Fuchsin (cytoplasm, muscle, the red counterstain), not nuclei: the published Masson matrix has no nuclear stain (iron haematoxylin is near black and absorbs in all channels). The mask hint stays `exclude_nuclei`: a nuclear mask is built in @@ -630,7 +637,10 @@ preprocessing, not taken from this exclusion map. User vectors for `trichrome_cu `(0, 255]`, checked as for `custom_od`). When given, the RGB is rescaled per channel before the scikit-image call to `clip(rint(I * 255 / white), 0, 255)` as `uint8`, computed in float64 in that order; `rint` rounds half to even, and a value above 255 after rounding (brighter than the -white by more than half a grey level) is set to 255. `white_clipped_fraction` is the fraction of +white by more than half a grey level) is set to 255. Re-quantising to `uint8` keeps scikit-image's +colour conversions on the 8-bit input they are defined for and adds at most half a grey level of +rounding per channel (plus the reported clipping); scikit-image would also accept float RGB, which +would avoid this rounding, but the method is fixed as stated here. `white_clipped_fraction` is the fraction of pixels with at least one such clipped channel; `above_white_fraction` (as for `custom_od`) the fraction with a channel strictly above the white, which includes pixels within half a grey level of it that round to 255 without clipping. `white = (255, 255, 255)` is the identity, so the diff --git a/tests/test_stain_defaults.py b/tests/test_stain_defaults.py index bc1605a..51b9e87 100644 --- a/tests/test_stain_defaults.py +++ b/tests/test_stain_defaults.py @@ -343,3 +343,53 @@ def test_a_reference_region_with_a_floored_pixel_is_refused_and_zeros_outside_do rgb[3, 3, 0] = 0 # inside the region: its OD would be taken at the floor with pytest.raises(ValueError, match="holds 1 pixel"): stain_vectors_from_regions(rgb, axes="YXC", white=WHITE, regions={"eosin": region}) + + +# --- Review of 8cc45da: he_custom is a two-stain model, tested without scikit-image ------------ + +# Literal reference values (Ruifrok & Johnston 2001, p. 293, normalised OD rows; Landini's +# Masson Trichrome), independent of scikit-image, so an upstream change cannot move them. +LITERAL_HE = {"hematoxylin": (0.65, 0.70, 0.29), "eosin": (0.07, 0.99, 0.11)} +LITERAL_BPX = { + "collagen": (0.7995107, 0.5913521, 0.10528667), + "ponceau_fuchsin": (0.09997159, 0.73738605, 0.6680326), +} +# he_custom's unit matrix rows: eosin, haematoxylin, their normalised cross product (computed once +# by hand from the literal rows above; not scikit-image's DAB row 0.27, 0.57, 0.78). +HE_CUSTOM_ROWS = np.array( + [ + [0.07010172, 0.99143863, 0.11015985], + [0.65110783, 0.70119304, 0.29049426], + [0.33211686, 0.08093471, -0.93975952], + ] +) + + +def test_default_vectors_equal_the_literal_published_values(): + assert dict(HE_DEFAULT_VECTORS) == LITERAL_HE + assert dict(BPX_DEFAULT_VECTORS) == LITERAL_BPX + + +def test_he_custom_is_the_two_stain_model_on_a_reference_pixel(): + """One pixel by hand: log10 OD against the declared white, then OD @ inv(rows above).""" + pixel = np.array([[[180, 90, 170]]], np.uint8) + white = (250.0, 250.0, 250.0) + od = -np.log10(np.array([180.0, 90.0, 170.0]) / 250.0) + expected = od @ np.linalg.inv(HE_CUSTOM_ROWS) + result = to_scalar(pixel, axes="YXC", preset="he_custom", white=white) + np.testing.assert_allclose(result.image.ravel(), expected[0], rtol=1e-6) + np.testing.assert_allclose(result.exclusion.ravel(), expected[1], rtol=1e-6) + np.testing.assert_allclose(result.provenance["stain_matrix"], HE_CUSTOM_ROWS, atol=1e-8) + + +def test_he_custom_is_not_the_three_stain_hed_conversion(): + """Same H and E vectors, different model: scikit-image's HED has a DAB row, natural-log OD + scaled by log(1e-6), white 255 and clipping; its eosin differs from he_custom's.""" + color = pytest.importorskip("skimage.color") + pixel = np.array([[[180, 90, 170]]], np.uint8) + ours = to_scalar(pixel, axes="YXC", preset="he_custom", white=(255.0, 255.0, 255.0)).image + hed = color.separate_stains(pixel, color.hed_from_rgb)[..., 1] + assert not np.allclose(ours, hed, rtol=0.1) + assert not np.allclose( + color.rgb_from_hed[2], np.cross(color.rgb_from_hed[0], color.rgb_from_hed[1]) + ) From a5ee1648c44ea41a88c6d9740679a2604a79bf8f Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 13 Nov 2026 19:17:24 +0000 Subject: [PATCH 12/12] Release v0.6.0: version, release notes and known limitations pyproject 0.6.0; CHANGELOG 'Unreleased' becomes v0.6.0 (2026-11-13), an additive release (stain presets, optional extra histology; nothing that existed changes) with the known limitations the release review asked to state; README release line, install command and a stain-presets bullet; plan status; AGENTS.md decision range N1-N69 (release review FILE 1). No code change. Checked: full check 2026 passed; a clean NumPy-only install reports 0.6.0, runs he_custom and refuses he_eosin naming the histology extra; a clean [histology] install runs he_eosin (scikit-image 0.26.0). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VbytvF3Gsa85ix8Ne2yMDW --- AGENTS.md | 4 ++-- CHANGELOG.md | 24 +++++++++++++++++++++++- README.md | 8 +++++--- docs/plan.md | 2 +- pyproject.toml | 2 +- 5 files changed, 32 insertions(+), 8 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 3a7c9fc..19594b4 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -6,7 +6,7 @@ Read it, and the documents below, before changing anything. | File | What it fixes | |---|---| -| `docs/plan.md` | Goals, methods, scopes, milestones, the owner's decisions N1–N64 (§12) and what is still to be measured (§13) | +| `docs/plan.md` | Goals, methods, scopes, milestones, the owner's decisions N1–N69 (§12) and what is still to be measured (§13) | | `docs/contract.md` | Rules other bricks rely on: methods and formulas, interface, scopes and repeat groups, image types, the statistics file, intensity checks, parameter columns, validation | | `CHANGELOG.md` | Release notes and migration notes | | `docs/archive/` | Historical record: the two done follow-ups for the calling bricks | @@ -127,7 +127,7 @@ range check. ## 3. How we work 1. **Read first:** the plan, the contract and the modules you touch. -2. **Ask before deciding.** Decisions N1–N64 are final (those marked withdrawn in plan §12 no longer apply); do not re-open them. A choice the +2. **Ask before deciding.** Decisions N1–N69 are final (those marked withdrawn in plan §12 no longer apply); do not re-open them. A choice the plan and the contract leave open, and any change of product behaviour or file format, is the owner's: ask with 2–4 options, the recommended one first. 3. **Outputs stay bit-identical** unless a change says otherwise. `tests/test_fit_goldens.py` keeps the fitted diff --git a/CHANGELOG.md b/CHANGELOG.md index 7ccc81c..bfd390c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,29 @@ Release notes and migration notes of `morphobricks-normalization`. Numbers are p issues of [CellMigrationLab/morphobricks-normalization](https://github.com/CellMigrationLab/morphobricks-normalization). The decisions are in [docs/plan.md](docs/plan.md) §12. -## Unreleased +## v0.6.0 (2026-11-13) + +Additive release (issue #94, PR #95; plan N67-N69, contract §3.4): a new `stain` module turns one +declared imaging unit into one scalar map with a named preset; scikit-image is the new optional extra +`histology`. Nothing that existed changes: `rgb.py`, `NormalizationParameters`, +`plane_histograms.csv`, the CSBDeep parity path and every public name and output of v0.5.0 are +untouched, and the package still imports with NumPy only. Reviewed independently at `3a7af3e`, +`8cc45da` and `2506e51` (all findings fixed). + +**Known limitations** (state them with results): +- RGB inputs are declared 8-bit sRGB only, three samples and no alpha; 16-bit, ICC profiles and + scanner colour management are not handled. +- `he_custom` is a two-stain model (cross-product residual, declared white), not scikit-image's + three-stain HED deconvolution; use `he_eosin` for the published conversion. +- The fixed HED/BPX deconvolution keeps scikit-image's floor and zero-clipping; `zero_fraction` + cannot tell naturally zero components from clipped negatives. +- With a declared white, the luminance and CIELAB projections re-quantise to uint8 (clipping + reported); `custom_od` floors zero intensity to 1 and keeps negative concentrations. +- The trichrome vectors are one Masson profile, not every protocol; a linear polarizer makes + orientation unreliable; eosin-derived maps are not collagen-specific. +- Beer–Lambert on scanner sRGB values is an approximation; preset scaling hints are recorded, never + applied. + - **RGB-to-scalar stain presets** (issue #94, plan N67, contract §3.4): new `stain` module with `rgb_to_scalar(pixels, *, axes, preset)` and the preset registry `PRESETS` (`he_eosin`, diff --git a/README.md b/README.md index c0bc02c..da602b0 100644 --- a/README.md +++ b/README.md @@ -29,12 +29,14 @@ out of scope for now. ## Status -Pre-alpha; the latest release is v0.5.0. The library offers: +Pre-alpha; the latest release is v0.6.0. The library offers: - supported pixel types, bin counts and per-plane histograms (`dtypes`, `histogram`), with a robust range and a range check; - plane records, `plane_histograms.csv` reading and writing and pooling of group statistics (`records`, `csv_io`, `pooling`): schema version 5 since v0.5.0, with channels named by canonical channel ID (ADR 0003; version 4 tables are refused, see the CHANGELOG migration notes); - settings, scopes, repeat groups and `fit`, with `Parameters` (`apply`, `to_columns`, `from_columns`), for fluorescence and RGB images; -- per-image summaries and intensity checks, returned as `Issue` warnings. +- per-image summaries and intensity checks, returned as `Issue` warnings; +- since v0.6.0, image-to-scalar stain presets (`stain`: H&E, picrosirius red, trichrome, colour + projections and passthrough; optional extra `histology` for scikit-image; contract §3.4). Release notes and migration notes are in [CHANGELOG.md](CHANGELOG.md). @@ -56,7 +58,7 @@ for key, rows in groups.items(): ## Installation ```bash -pip install "git+https://github.com/CellMigrationLab/morphobricks-normalization.git@v0.5.0" +pip install "git+https://github.com/CellMigrationLab/morphobricks-normalization.git@v0.6.0" ``` For development, from a clone: diff --git a/docs/plan.md b/docs/plan.md index e717346..254e993 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -2,7 +2,7 @@ > See also: [Preprocessing follow-up](archive/followup_preprocessing.md) (done) · [CellPatcher follow-up](archive/followup_cellpatcher.md) (done) -**Status:** implemented through v0.5.0 (Milestones 0 to 5 done; later changes in [CHANGELOG.md](../CHANGELOG.md)); design brief simplified on 26 October 2026 (decisions N29–N40, §12), with +**Status:** implemented through v0.6.0 (Milestones 0 to 5 done; later changes in [CHANGELOG.md](../CHANGELOG.md)); design brief simplified on 26 October 2026 (decisions N29–N40, §12), with channels moved into columns on 27 October 2026 (N41–N45). The earlier design with six modalities, two purposes, colour and stain methods and reference targets (N4–N28) is withdrawn where §12 says so. A robust histogram range was added on 28 October diff --git a/pyproject.toml b/pyproject.toml index d7d19cc..134b8da 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "morphobricks-normalization" -version = "0.5.0" +version = "0.6.0" description = "MorphoBricks normalization brick: per-plane statistics, their file formats, and normalization methods for microscopy images." readme = "README.md" license = "MIT"