diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 790011e..723dc0a 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 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 + 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 tests/test_stain_catalogue.py tests/test_stain_defaults.py tests/test_optical_density.py -rs diff --git a/AGENTS.md b/AGENTS.md index da806db..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 | @@ -56,6 +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 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_) @@ -113,13 +120,14 @@ 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 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 22047e3..bfd390c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,71 @@ 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. +## 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`, + `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`. +- **`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 + 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. +- **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) Breaking release. Stage 3 S3.2 (MorphoBricks-Spec roadmap; ADR 0003, plan N61–N66, PR #90): 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/contract.md b/docs/contract.md index b27d4d2..b6e5e36 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -488,6 +488,186 @@ y_c = (x_c − offset) / scale for c = R, G, B --- +### 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, +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), `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 | +|---|---|---|---|---|---|---|---|---| +| `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` | — (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 (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` | — (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, +`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"`: 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**. + +**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. +`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; +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. 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. + +**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 `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. + +**`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 (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 + 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. + +**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). **`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 +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. 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 +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_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 +`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` +(`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 ### 4.1 Scopes diff --git a/docs/plan.md b/docs/plan.md index 1618a6c..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 @@ -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,57 @@ 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). +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). 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. +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 diff --git a/pyproject.toml b/pyproject.toml index 0b1e568..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" @@ -24,6 +24,7 @@ dependencies = [ ] [project.optional-dependencies] +histology = ["scikit-image>=0.22"] dev = [ "pytest>=7.4", "ruff>=0.6", @@ -122,3 +123,7 @@ baseline = "tools/code_rules_baseline.json" "morphobricks_normalization.bin_reachability" = "pure" "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 new file mode 100644 index 0000000..174adce --- /dev/null +++ b/src/morphobricks_normalization/optical_density.py @@ -0,0 +1,342 @@ +"""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: 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. + +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 _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. + + 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]``. + """ + _refuse_masked(white, "white") + 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.""" + _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( + 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 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()) + + +def _check_region(region: object, shape: tuple[int, int], name: str) -> np.ndarray: + """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: + 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): 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) + + +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. + + 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, 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. + """ + 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) + 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: + 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..15a8cb4 --- /dev/null +++ b/src/morphobricks_normalization/scalar_inputs.py @@ -0,0 +1,107 @@ +"""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 _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. + + 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}" + ) + _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}. " + "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.") + _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}. " + "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 new file mode 100644 index 0000000..859a1d5 --- /dev/null +++ b/src/morphobricks_normalization/stain.py @@ -0,0 +1,418 @@ +"""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**; ``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 + inverted, since birefringent fibres are bright on a dark field) and ``cielab`` (``lab_a``: + 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 ``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 ``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. + +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 types import MappingProxyType, ModuleType +from typing import Any + +import numpy as np + +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]" +POLARIZERS = ("circular", "linear") +OPTIONAL_ARGUMENTS = ("white", "stain_vectors", "polarizer") +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) +class ScalarResult: + """The scalar map of :func:`to_scalar`. + + Attributes + ---------- + image + float32 ``(Y, X)``, finite. + exclusion + 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. + """ + + image: np.ndarray + exclusion: np.ndarray | None + provenance: Mapping[str, Any] + + +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 _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() + 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": + gray = color.rgb2gray(rgb) + image = 1.0 - gray if preset.inverted else gray + elif preset.method == "cielab": + 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) + excluded = None if exclusion is None else np.ascontiguousarray(exclusion, dtype=np.float32) + 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. + + 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) + image = np.ascontiguousarray(amounts[..., 0], dtype=np.float32) + 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, + "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.") + 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." + ) + 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, + "image_name": recipe.image_name, + "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, + "zero_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 + 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`. Optional for + ``he_custom`` and ``trichrome_custom`` (published defaults), required for + ``psr_brightfield``. + 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), ``zero_fraction``, and the method's own entries (``custom_od``: white, + ``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 + ------ + 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: + 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, **white_entries} + 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 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 + ---------- + 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 + As :func:`to_scalar`. + + Raises + ------ + ValueError + 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]``). + """ + 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/src/morphobricks_normalization/stain_presets.py b/src/morphobricks_normalization/stain_presets.py new file mode 100644 index 0000000..8f00aca --- /dev/null +++ b/src/morphobricks_normalization/stain_presets.py @@ -0,0 +1,273 @@ +"""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. + 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 + 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 + image_name: 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", + image_name="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", + image_name="Hematoxylin", + stain="he", + matrix="hed_from_rgb", + target_component=0, + hints=_hints("fixed", "none", "haematoxylin, mainly nuclei"), + ), + StainPreset( + "luminance_inverted", + image_name="InvertedLuminance", + method="luminance", + inverted=True, + accepts=_WHITE, + hints=_hints("fixed", _HE_MASK, "inverted brightness (baseline), nuclei-dominated"), + ), + StainPreset( + "lab_a", + image_name="LabA", + method="cielab", + target_component=1, + accepts=_WHITE, + hints=_hints("fixed", _HE_MASK, "CIELAB a* redness (comparator), low contrast"), + ), + 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, + target_stain="psr", + exclusion_stain="counterstain", + hints=_hints("fixed", "exclude_counterstain", "collagen"), + ), + StainPreset( + "psr_polarized", + image_name="SiriusRedPolarized", + stain="picrosirius_red", + method="luminance", + modality="polarized", + needs=frozenset({"polarizer"}), + accepts=_WHITE, + hints=_hints("fixed", "background", "birefringent packed collagen"), + ), + StainPreset( + "trichrome_blue", + image_name="MethylBlue", + stain="trichrome_blue", + matrix="bpx_from_rgb", + target_component=0, + hints=_hints("fixed", "exclude_nuclei", "collagen (methyl blue)"), + ), + StainPreset( + "trichrome_green", + image_name="TrichromeGreen", + 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", + image_name="EosinCustom", + 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", + image_name="TrichromeCustom", + 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 new file mode 100644 index 0000000..309105c --- /dev/null +++ b/tests/stain_simulation.py @@ -0,0 +1,57 @@ +"""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 +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): + 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 new file mode 100644 index 0000000..02d6faf --- /dev/null +++ b/tests/test_stain.py @@ -0,0 +1,288 @@ +import subprocess +import sys +import textwrap + +import numpy as np +import pytest + +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) + + +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_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.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] + with pytest.raises(AttributeError): + PRESETS["he_eosin"].version = 2 # type: ignore[misc] + + +@pytest.mark.parametrize("preset", P0) +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["zero_fraction"] == b.provenance["zero_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", P0) +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["zero_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_zero_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, + "image_name": "Eosin", + "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", + "skimage_version": skimage.__version__, + "zero_fraction": float(np.mean(result.image == 0)), + } + assert dict(result.provenance) == expected + 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 ( + rgb_to_scalar(image, axes="YXC", preset="luminance_inverted").provenance["method"] + == "luminance" + ) + + +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["zero_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] + with pytest.raises(TypeError): + provenance["hints"]["scaling"] = "per_image" # type: ignore[index] + + +@pytest.mark.parametrize("preset", P0) +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) diff --git a/tests/test_stain_catalogue.py b/tests/test_stain_catalogue.py new file mode 100644 index 0000000..0e49cac --- /dev/null +++ b/tests/test_stain_catalogue.py @@ -0,0 +1,555 @@ +"""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, + METHYL_BLUE, + PONCEAU_FUCHSIN, + 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", "ponceau_fuchsin", METHYL_BLUE, PONCEAU_FUCHSIN), +} +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", "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"), + ("trichrome_blue", "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", "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}"): + 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["zero_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["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" + ).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["zero_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"] + + +# 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 + + +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" diff --git a/tests/test_stain_defaults.py b/tests/test_stain_defaults.py new file mode 100644 index 0000000..51b9e87 --- /dev/null +++ b/tests/test_stain_defaults.py @@ -0,0 +1,395 @@ +"""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"} + + +# --- 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}) + + +# --- 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]) + )