From 4db61229f1d5e8afae729ba25d5fcaf8db2a882e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Hidalgo?= <47150974+IvanHCenalmor@users.noreply.github.com> Date: Tue, 27 Oct 2026 20:22:26 +0200 Subject: [PATCH] Check the images' intensities from the table, with plots in the notebook (P5) After every written run, and in step 4 for any table, the normalization brick's image_summaries and intensity_checks give the four intensity warnings from the table alone (intensity module). The notebook shows them in steps 3 and 4 with a collapsed Check settings panel and per-channel box plots and cumulative distributions per biological repeat; nothing is written. The command line prints them, with --outlier-limit, --saturation-limit and --min-images; the manifest records the limits as the optional key intensity_checks (version 1). Plan question 30. Co-Authored-By: Claude Opus 5.5 --- AGENTS.md | 21 +- README.md | 6 +- docs/contract.md | 19 +- docs/plan.md | 29 ++ docs/user_guide.md | 55 ++- notebooks/MorphoBricks_Preprocessing.ipynb | 34 +- src/morphobricks_preprocessing/__init__.py | 16 + src/morphobricks_preprocessing/cli.py | 44 +- src/morphobricks_preprocessing/intensity.py | 277 ++++++++++++ src/morphobricks_preprocessing/manifest.py | 16 +- src/morphobricks_preprocessing/notebook_ui.py | 422 +++++++++++++++++- src/morphobricks_preprocessing/run.py | 19 +- tests/test_cli.py | 37 ++ tests/test_intensity.py | 177 ++++++++ tests/test_manifest.py | 14 + tests/test_notebook_ui.py | 81 ++++ tests/test_run.py | 32 ++ 17 files changed, 1280 insertions(+), 19 deletions(-) create mode 100644 src/morphobricks_preprocessing/intensity.py create mode 100644 tests/test_intensity.py diff --git a/AGENTS.md b/AGENTS.md index d137f545..584b614d 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -76,6 +76,9 @@ src/morphobricks_preprocessing/ image_type.py image type (P4): FLUORESCENCE/RGB, review_image_type (suggestion from TiffMetadata.rgb_samples and the confirmed roles, whether RGB is possible via the brick's check_rgb), image_type_issues (rgb_not_possible) + intensity.py intensity checks (P5): check_intensities/check_table (the brick's image_summaries + and intensity_checks, IntensityReport), repeat_boxes and repeat_distributions for + the notebook's plots, intensity_text for the CLI validation.py ValidationIssue (severity, code, message, paths, action), error()/warning() builders axes.py AxisMapping (role per dimension), prefill, check_mapping, review_axes (one mapping applied to every file, missing axes as singletons only when confirmed) @@ -386,6 +389,18 @@ functions the tests use. - **Manifest:** optional top-level `image_type` (null without one) and `channels_from_rgb_samples`; changing the type recomputes nothing. +### Intensity checks (contract §14, follow-up P5) + +- **From the table only:** `run_experiment` checks its rows in memory (`RunResult.intensity`, + apart from `RunResult.issues`, never blocking); `check_table` reads any table. The four + codes are the normalization brick's; comparisons stay within a condition. +- **Limits:** `IntensityCheckSettings` (3.0, 0.01, 5), recorded in the manifest + (`intensity_checks`, optional); changing them recomputes nothing. +- **Notebook:** `IntensityPanel` in step 3 (after a successful run, `run_checks`) and step 4 + (`table_checks`): status, tiles, warnings table with actions, collapsed **Check + settings**, channel buttons, box plots (`Issue.channel` marks images) and cumulative + distributions per biological repeat. Nothing written. CLI prints them, no plots. + ### Validation messages (contract §8) Every expected problem is a `ValidationIssue` with: @@ -426,6 +441,7 @@ Technical details go to `logging` (shown in the notebook's "Technical details" p | Dataset consistency (Q21, P2) | Step 1 "Dataset consistency" panel: balance table (images and technical-repeat folders per condition x biological repeat), `unbalanced_images` (< 0.5 or > 2 x median), `repeat_missing_conditions`, `different_pixel_sizes`, `missing_pixel_size` (OME-Zarr scale included), `different_plane_counts`; `different_sizes` narrowed to width/height; consistency issues shown only in that panel, not in Messages. Warnings only; outputs unchanged | | Experiment design (Q22, P3) | Step 1 after the labels: pairing dropdowns, one yellow "Yes, they are comparable" question per repeat group (single whole-experiment question without biological-repeat folders), control checkboxes; `experiment_design.csv` written after every answer; reruns load it and fold the questions behind a Change button; step 3 gated. CLI `run` reads it or `--design`, else writes a template and exits 1. `run_experiment(design=...)` optional; manifest `experiment_design_sha256` (optional key, version 1). Two folders of one condition in a group block that group | | Image type (Q29, P4) | Step 2 card after the axis roles: Fluorescence or RGB, suggestion preselected, confirmed with **Yes, this is the right image type**; RGB suggested for TIFF photometric RGB with three samples on the C axis; RGB needs three integer channels (`check_rgb`), else disabled with the reason and `rgb_not_possible` in a run. Manifest keys `image_type`, `channels_from_rgb_samples` (optional, version 1); type changes recompute nothing. CLI `--image-type`, default the printed suggestion | +| Intensity checks (Q30, P5) | Warnings from the table alone with the brick's `image_summaries`/`intensity_checks`, shown automatically after step 3 and in step 4 for any table; limits under a collapsed **Check settings** panel, recorded in the manifest (`intensity_checks`); per-channel box plots with flagged images and cumulative distributions per biological repeat, notebook only; the CLI prints the warnings (`--outlier-limit`, `--saturation-limit`, `--min-images`) | | Position identity (Q20, D02) | `position_id` = image name without extension, in every CSV row and manifest file entry; FOV key = labels + `position_id`; duplicates in one leaf folder are the error `duplicate_position_id`. CSV columns and schema 2 unchanged; manifest version 1 with an optional key | | Notebook style | Follow the lab's NucleiSky and LabChronicle notebooks (see §7) | | Widget widths | **99%**, never 100% | @@ -714,9 +730,8 @@ jupyter lab notebooks/MorphoBricks_Preprocessing.ipynb - **Released as v0.1.0** (27 October 2026), with everything above; the normalization brick is pinned to its v0.2.0 (question 28). The next tag, v0.2.0, follows follow-ups P4 and P5. -- **Since v0.1.0:** follow-up P4 of the normalization plan, the image type - (question 29). -- **Next:** follow-up P5 of the normalization plan (intensity checks and plots). +- **Since v0.1.0:** follow-ups P4 (the image type, question 29) and P5 (intensity checks + and plots, question 30) of the normalization plan. The owner tags v0.2.0 next. - **Not yet verified:** - the exact look of the notebook in a real JupyterLab; - the optional `jl-hidecode` JupyterLab extension, which keeps code cells collapsed. It diff --git a/README.md b/README.md index 37ad4649..3d64fac0 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,9 @@ Pre-alpha. Implemented and tested: confirmed by the user and recorded in the run manifest (`image_type`) - Running an experiment end to end and updating earlier output incrementally, with a run manifest (`run`, `manifest`) +- Intensity checks of the written table: images far from the others of their biological + repeat, shifted biological repeats, saturated images and images without spread, with + plots per channel in the notebook (`intensity`) - A step-by-step notebook ([notebooks/MorphoBricks_Preprocessing.ipynb](notebooks/MorphoBricks_Preprocessing.ipynb)) and a command line (`morphobricks-preprocessing check` / `run`) @@ -96,7 +99,8 @@ conditions, and that they are comparable) and write `experiment_design.csv`; in pass it with `design=` (see the [user guide](docs/user_guide.md) §9). Later bricks read it with `read_experiment_design`. They also ask for the image type, `fluorescence` or `rgb` (suggested by `review_image_type`), and record it in the manifest; in Python, pass it with -`image_type=`. +`image_type=`. After writing, the run checks the images' intensities from the table +(`result.intensity`); `check_table(path)` checks any table. ## Example diff --git a/docs/contract.md b/docs/contract.md index b85d83f1..1a6ad8fd 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -108,7 +108,7 @@ method is an offset and a scale, `(pixels − offset) / scale`. | Severity | Situations | Effect | |---|---|---| | Error | Experiment design not confirmed (`design_not_confirmed`) or a repeat group holding two folders of one condition (`repeat_group_conflict`), when a design is used (§12); image type RGB for images that are not three channels of an integer type (`rgb_not_possible`, §13); mixed dtypes; images with different numbers of channels (`different_channel_counts`, §2); unreadable TIFF (header or pixel data) or OME-Zarr folder in a leaf folder that is not excluded; OME-Zarr images, or converting TIFF images to OME-Zarr, without the `zarr` extra; OME-Zarr images whose declared axes do not fit the experiment's roles; every image excluded; more than one image series in a file; unsupported sample format or axis; missing/duplicate X or Y, duplicate roles, unassigned non-singleton axis, inconsistent mapping; plane that cannot be read within the memory limit, not even in blocks (for example a plane stored compressed in one piece); folder named `single` in another case, or beside other folders (§1); sibling folders differing only in case; images in one leaf folder with the same position identifier (§1); folder below the technical-repeat level; output folder inside the experiment; no usable planes | Processing is blocked | -| Warning | Files that are not images; empty folders; images outside leaf folders; images excluded by the user, or excluded paths that match no image; channel names differing between files; non-finite pixels; different image widths or heights; different numbers of depth slices, time points or other planes; physical pixel sizes that differ or are known for only some images; unbalanced numbers of images; biological repeats missing from some conditions (§11); planes without finite pixels while others remain; planes over the memory limit that are read in blocks | Processing continues; the item is reported | +| Warning | Intensity checks of the written table (`image_intensity_outlier`, `repeat_intensity_shift`, `saturated_image`, `no_intensity_spread`, §14), reported apart from the run's issues; files that are not images; empty folders; images outside leaf folders; images excluded by the user, or excluded paths that match no image; channel names differing between files; non-finite pixels; different image widths or heights; different numbers of depth slices, time points or other planes; physical pixel sizes that differ or are known for only some images; unbalanced numbers of images; biological repeats missing from some conditions (§11); planes without finite pixels while others remain; planes over the memory limit that are read in blocks | Processing continues; the item is reported | Every message states what happened in plain language, lists the affected paths, and suggests an action. Technical details go to the diagnostic log. @@ -152,6 +152,7 @@ file's schema or statistical interpretation. | `experiment_design_sha256` | SHA-256 of the `experiment_design.csv` written with this run (§12), so later bricks can check that the design they read is the one used; null when the run had no design, or absent in manifests written before it. Optional key: the manifest stays version 1. | | `image_type` | `fluorescence` or `rgb`, the image type the user confirmed (§13); null when the run was given none, or absent in manifests written before it. Optional key: the manifest stays version 1. Changing it recomputes nothing. | | `channels_from_rgb_samples` | `true` when the C role is the RGB samples axis of images stored as RGB colour (§13), `false` otherwise; absent in manifests written before it. Optional key. | +| `intensity_checks` | The limits of the intensity checks run on the table (§14): `IntensityCheckSettings.to_dict()` of the normalization brick (`library_version`, `outlier_limit`, `saturation_limit`, `min_images`); absent in manifests written before it. Optional key. Changing the limits recomputes nothing. | | `versions` | Versions of `morphobricks-preprocessing`, `morphobricks-normalization` (which computes the histograms and writes the CSV; absent in manifests written before it), `numpy`, `tifffile` and Python. Information only: a different version never causes a rebuild. | Channel names and significant bits come from OME metadata and are for information only. @@ -257,3 +258,19 @@ plan N34); this brick suggests it, asks for it after the axis roles, and records | Confirmation | Never applied without the user's confirmation. The notebook asks in step 2 once the axis roles are confirmed, with the suggestion preselected and a **Yes, this is the right image type** button; `rgb` cannot be chosen, with the reason shown, when it is not possible. Changing the axis roles or the folder asks again. Step 3 runs only once the type is confirmed. The command line takes `--image-type`, or uses the suggestion and prints it. | | Run | `run_experiment(..., image_type=...)` records the type and `channels_from_rgb_samples` in the manifest (§9). `rgb` for images that cannot be RGB is the error `rgb_not_possible`; any other value raises `ValueError`. Without a type the manifest records null. | | Outputs | `plane_histograms.csv` and `experiment_design.csv` do not depend on the type: a plane of an RGB image is one row with channels C0, C1 and C2, as for any image. Changing the type on a later run reads no image again. | + +## 14. Intensity checks + +After the table is written, the images' intensities are checked from its rows alone, with +the normalization brick's `image_summaries` and `intensity_checks` ([normalization +contract][normalization contract] §7, plan N38): no image is read and nothing is written +apart from the limits in the manifest. + +| Rule | Detail | +|---|---| +| Summaries | One per image, all Z and T planes pooled per channel: median, 1st and 99.8th percentiles, saturated fraction (normalization contract §7.1). | +| Warnings | `image_intensity_outlier` (an image's median or 99.8th percentile in a channel more than `outlier_limit` robust deviations from the other images of its condition and biological repeat), `repeat_intensity_shift` (a biological repeat's median of image medians far from the other repeats' images of its condition), `saturated_image` (more than `saturation_limit` of the pixels at the type's maximum), `no_intensity_spread` (1st and 99.8th percentiles equal). One warning per image or repeat, channel and compared value, naming the channel `C`, the image or `/`, the value and the group value or limit. Comparisons stay within a condition, so treatment effects are never reported. | +| Limits | `IntensityCheckSettings`: `outlier_limit` 3.0, `saturation_limit` 0.01, `min_images` 5 (provisional in the normalization brick). `run_experiment(..., intensity_settings=...)` uses them and records them in the manifest (`intensity_checks`, §9); they change no row. | +| Run | `RunResult.intensity` holds the summaries and warnings (`IntensityReport`) of a written table, None otherwise. The warnings never block and are kept apart from `RunResult.issues`. | +| Notebook | Step 3 shows the checks of the table just written after every successful run, step 4 those of any chosen table: a status line with the limits, counts per kind, a table of the warnings (check, channel, image or repeat, what was found) and what to do about each kind, a collapsed **Check settings** panel (changing a limit updates the warnings at once; in step 3 the next run records it), and per channel (chosen with buttons) box plots of the image medians and high percentiles per condition and biological repeat, with images that have a warning in that channel (`Issue.channel`, N52) as red diamonds and shifted repeats as amber boxes, and the cumulative intensity distribution of each biological repeat, one panel per condition (shifted repeats dashed; at most 12 conditions). The plots are drawn in the notebook only; no file is written. | +| Command line | `run` prints the warnings after the run, grouped by kind with what to do once per kind, and draws no plots; `--outlier-limit`, `--saturation-limit` (a fraction) and `--min-images` set the limits. | diff --git a/docs/plan.md b/docs/plan.md index 0046ad45..62486590 100644 --- a/docs/plan.md +++ b/docs/plan.md @@ -346,6 +346,12 @@ API, CLI). The notebook and CLI only call the API. needs three channels of an integer type. The manifest records the type; the command line takes `--image-type` or prints the suggestion. Outputs are otherwise unchanged. *Implemented, see question 29.* +- **Intensity checks (follow-up P5 of the normalization plan):** after every successful + run, and in step 4 for any table, the normalization brick's `image_summaries` and + `intensity_checks` give the four intensity warnings from the table alone (`intensity`); + the notebook shows them with adjustable limits and per-channel box plots and + distributions, the command line prints them, and the manifest records the limits. + *Implemented, see question 30.* ## 11. Tests and acceptance criteria @@ -719,3 +725,26 @@ Raised during review; to be decided before the milestones that need them: `--image-type fluorescence|rgb`, or uses the suggestion and prints it with the reason; `check` prints the suggestion. `plane_histograms.csv` and `experiment_design.csv` do not change.* +30. How are the intensity checks shown (normalization plan N38, N50–N53, its contract §7, + follow-up P5)? + *Decided (27 October 2026): `run_experiment` checks the rows it wrote, in memory, with + the normalization brick's `image_summaries` and `intensity_checks` + (`intensity.check_intensities`, `RunResult.intensity`), and `check_table` checks any + table; no image is read. The four warnings are kept apart from the run's issues and + never block. Agreed with the owner: the warnings and the plots are shown both in step 3, + automatically after every successful run, and in step 4 for the chosen table; the + manifest records the limits used (optional key `intensity_checks`, + `IntensityCheckSettings.to_dict()`; version 1), changing them afterwards only updates + the warnings on screen and is recorded by the next run; the distributions are + cumulative curves, like step 4's existing plot. Each panel has a status line, tiles per + kind, a table (check, channel `C`, image or repeat, what was found) with what to do + per kind, a collapsed **Check settings** panel (`outlier_limit` 3.0, `saturation_limit` + 1 %, `min_images` 5, and **Restore the defaults**), channel buttons, box plots of the + image medians and high percentiles per condition and biological repeat (images with a + warning in that channel, from `Issue.channel`, as red diamonds; shifted repeats as amber + boxes), and the cumulative distribution of each biological repeat, one panel per + condition (shifted repeats dashed; at most 12 conditions). Plots are drawn in the + notebook only and nothing is written. The command line prints the warnings grouped by + kind with the action once per kind, takes `--outlier-limit`, `--saturation-limit` and + `--min-images`, and draws no plots. Comparisons stay within a condition, as the brick + does. `plane_histograms.csv` and `experiment_design.csv` do not change.* diff --git a/docs/user_guide.md b/docs/user_guide.md index f3315569..496be4c4 100644 --- a/docs/user_guide.md +++ b/docs/user_guide.md @@ -250,6 +250,28 @@ results or comparing another grouping. then summarise planes written, images read, kept, new, changed, removed and excluded. Nothing is written unless the whole run succeeds. **Technical details (log)** holds the log for troubleshooting. + - **Intensity checks** then compares the images, using only the table just written (no + image is read again, and nothing else is written). Each image is compared with the + other images of its condition and biological repeat, and each biological repeat with + the other repeats of its condition, so differences between conditions (treatment + effects) are never reported: + - Tiles count the images checked and each kind of warning; a table lists every + warning with the check, the channel (**C0**, **C1**, …), the image or biological + repeat, and what was found (the values compared and the limit), followed by what to + do about each kind. The four kinds are images far from the others of their + biological repeat, biological repeats shifted from the others of their condition, + saturated images, and images without intensity spread (see §8). + - **Check settings** (collapsed) holds the limits: how many robust deviations count as + far (default 3), the share of saturated pixels (default 1 %), and the fewest images a + comparison needs (default 5). Changing them updates the warnings and plots at once; + the run records the limits in **run_manifest.json**, so the next run records the + new ones. **Restore the defaults** sets them back. + - Plots, one channel at a time (choose it above the plots): box plots of the image + medians and high percentiles for each condition and biological repeat, one point per + image, with images that have a warning in that channel as red diamonds and shifted + repeats as amber boxes; and one panel per condition with the intensity distribution + of each biological repeat (comparable repeats lie close together; a shifted one is + dashed). The plots are only shown here, never saved. 4. **Calculate the statistics (optional).** The table from step 3 is all that later bricks need; this step only reads it to compare groups and writes no files. - **Histogram table** is filled in with the table from step 3. To look at another @@ -264,6 +286,9 @@ results or comparing another grouping. before; its **planes** column counts each channel of a plane as one plane. - The plot compares the groups' intensity distributions: curves further right are brighter. + - **Intensity checks** shows the same warnings, settings and plots as step 3 for the + chosen table, so tables of other experiments can be checked too. It needs only the + table, not the images. In JupyterLab, the optional `jl-hidecode` extension keeps the code cells collapsed and adds a run button to each; without it, click the "…" of a cell to see its code. @@ -316,9 +341,13 @@ the design is confirmed. Other options: | `--chunk-size PX` | With `--convert-tiff`: tile size of the copies (16 to 8192, default 1024) | | `--compression lz4\|zstd\|none`, `--compression-level N` | With `--convert-tiff`: lossless compression (default lz4, level 5) | | `--pyramid-levels auto\|none\|N` | With `--convert-tiff`: smaller levels for viewing (default auto) | +| `--outlier-limit N`, `--saturation-limit FRACTION`, `--min-images N` | Limits of the intensity checks printed after the run (default 3, 0.01 and 5) | | `--log FILE` | Write technical details to a file | -Both commands exit with status 1 when an error blocks processing. +After writing, `run` prints the intensity checks of the table, grouped by kind, each with +the channel, the image or biological repeat and the values compared, and what to do once +per kind; it draws no plots (the notebook does). Both commands exit with status 1 when an +error blocks processing; intensity warnings never do. `python -m morphobricks_preprocessing` works the same way. ## 6. Output files @@ -345,8 +374,8 @@ The default output folder is next to the experiment folder: [contract.md](contract.md) §12. - **`run_manifest.json`:** what the run read and with which settings: labels, axis roles, the image type (and whether the channels are RGB colour samples), bin count, memory - limit, excluded images, package versions, the fingerprint of the experiment design, and - each image's position, size and modification time. Channel names and significant bits are recorded when the metadata + limit, excluded images, the limits of the intensity checks, package versions, the + fingerprint of the experiment design, and each image's position, size and modification time. Channel names and significant bits are recorded when the metadata provides them, and whether the image's planes were too large to read whole and are read in blocks. Format in [contract.md](contract.md) §9. - **`ome_zarr/`** (only when TIFF images are converted): one `.ome.zarr` folder per TIFF @@ -432,6 +461,10 @@ continues. | `design_not_confirmed` | Error | Some repeat groups are not confirmed as comparable. Confirm them in step 1, or set `confirmed` to `true` in `experiment_design.csv`. | | `rgb_not_possible` | Error | The image type RGB was chosen, but the images are not three channels of a whole-number pixel type. Choose fluorescence, or check the axis roles: the RGB samples must have the role C. | | `repeat_group_conflict` | Error | A repeat group holds two folders of the same condition. Move one into another group. | +| `image_intensity_outlier` | Warning (intensity check) | An image's median or high percentile in a channel is far from the other images of its condition and biological repeat. Look at it: it may be out of focus, contain debris or bubbles, or have been acquired differently. Exclude it if it is faulty. | +| `repeat_intensity_shift` | Warning (intensity check) | A biological repeat is brighter or darker than the other repeats of its condition. Check whether it was prepared or acquired differently; normalising per biological repeat in later bricks removes such differences. | +| `saturated_image` | Warning (intensity check) | More than the limit (default 1 %) of an image's pixels in a channel are at the largest value of its pixel type, so their true intensity is lost. Check the exposure, gain or laser power. | +| `no_intensity_spread` | Warning (intensity check) | An image's channel has the same 1st and 99.8th percentile (empty or blank), so it cannot be normalised on its own. Check the channel, or exclude the image. | | `channel_names_differ` | Warning | Channel names differ between images at the same position. Check the channel order. | | `non_finite_pixels` | Warning | NaN or infinite pixels are left out of the histograms and counted. | | `no_finite_pixels` | Warning | Some planes have a channel without finite pixels; those channels are left out of group statistics. | @@ -469,6 +502,22 @@ pairing = design.repeat_groups() # {(condition, biological_repeat): repeat_grou Later bricks read the file with `read_experiment_design(path)`. +The intensity checks of a written table are in `result.intensity` (summaries and the +normalization brick's warnings, each with its channel); `check_table(path)` checks any +table: + +```python +from morphobricks_normalization import IntensityCheckSettings +from morphobricks_preprocessing import check_table, repeat_boxes + +limits = IntensityCheckSettings(outlier_limit=4) +result = run_experiment("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/data/exp1", mapping, intensity_settings=limits) +for issue in result.intensity.issues: + print(issue.code, issue.channel, issue.items, issue.message) +records, report = check_table(result.csv_path, limits) +boxes = repeat_boxes(report, channel=0) # what the notebook's box plots show +``` + The image type is optional too: `run_experiment(..., image_type="rgb")` records it in the manifest, and without it the manifest records none. `review_image_type` gives the suggestion and whether RGB is possible, from the headers and the axis roles: diff --git a/notebooks/MorphoBricks_Preprocessing.ipynb b/notebooks/MorphoBricks_Preprocessing.ipynb index 2369e3c2..3c49c601 100644 --- a/notebooks/MorphoBricks_Preprocessing.ipynb +++ b/notebooks/MorphoBricks_Preprocessing.ipynb @@ -25,6 +25,9 @@ " colour images.\n", "- **Computes the histograms** and writes them, with a record of what was done, to a separate\n", " output folder.\n", + "- **Checks the intensities** of every image from the histograms: images far from the others of\n", + " their biological repeat, shifted biological repeats, saturated images and images without\n", + " spread, with plots per channel.\n", "- **Shows pooled statistics** per condition, repeat or image.\n", "\n", "### Quick workflow\n", @@ -36,8 +39,10 @@ "3. **Step 2:** check the role of each axis, press **Show preview**, press\n", " **Yes, the preview looks right** and then **Confirm axis roles**. Check the suggested image\n", " type and press **Yes, this is the right image type**.\n", - "4. **Step 3:** press **Compute histograms**. This writes the table the next bricks need.\n", - "5. **Step 4 (optional):** press **Calculate statistics** to compare your groups.\n", + "4. **Step 3:** press **Compute histograms**. This writes the table the next bricks need. Read\n", + " the **Intensity checks** that follow.\n", + "5. **Step 4 (optional):** press **Calculate statistics** to compare your groups and check the\n", + " intensities of any table.\n", "\n", "### Before you start\n", "\n", @@ -349,7 +354,23 @@ "> the images a second time, so the space they occupy is duplicated; the card shows how much space is\n", "> needed and lets you choose the chunk size, compression and pyramid levels. It needs the zarr package.\n", ">\n", - "> 🛡️ Nothing is written unless the whole run succeeds." + "> 🛡️ Nothing is written unless the whole run succeeds.\n", + "\n", + "After the run, **Intensity checks** compares the images using the table just written (no image is\n", + "read again, nothing more is written):\n", + "\n", + "- warnings about images far from the other images of their condition and biological repeat,\n", + " biological repeats shifted from the other repeats of their condition, saturated images and\n", + " images without intensity spread, each with the channel (**C0**, **C1**, …), the image or\n", + " repeat, and the values compared, and what to do;\n", + "- box plots of the image medians and high percentiles per condition and biological repeat\n", + " (red diamonds: images with a warning; amber boxes: shifted repeats), and the intensity\n", + " distribution of each biological repeat within each condition. Choose the channel above the\n", + " plots.\n", + "\n", + "> ℹ️ Images are compared only within their condition, so a treatment effect is never reported.\n", + "> The limits (3 robust deviations, 1 % saturated pixels, at least 5 images) are under\n", + "> **Check settings**; changing them updates the warnings at once, and the next run records them." ] }, { @@ -391,7 +412,8 @@ "2. Choose how to group the planes (**Condition**, **Biological repeat**, **Technical repeat** or\n", " **Image**) and the percentiles to report.\n", "3. Press **Calculate statistics**. The table lists pooled statistics per group, and the plot\n", - " compares their intensity distributions." + " compares their intensity distributions. **Intensity checks** shows the same warnings and\n", + " plots as step 3 for the chosen table, with its own **Check settings**." ] }, { @@ -463,6 +485,10 @@ "**\"RGB images need exactly three colour channels of a whole-number pixel type\"**: choose\n", "**Fluorescence** in step 2, or check the axis roles: the RGB samples must have the role C.\n", "\n", + "**\"Image …, channel C0: its median intensity … is … robust deviations from …\"**: an\n", + "intensity check, never an error. Look at the image: it may be out of focus, contain debris or\n", + "have been acquired with other settings. If it is faulty, exclude it and run again.\n", + "\n", "**The preview looks like stripes**: X or Y is on the wrong dimension; change the roles in step 2.\n", "\n", "More help: the [user guide](../docs/user_guide.md)." diff --git a/src/morphobricks_preprocessing/__init__.py b/src/morphobricks_preprocessing/__init__.py index f536eb4f..9e673220 100644 --- a/src/morphobricks_preprocessing/__init__.py +++ b/src/morphobricks_preprocessing/__init__.py @@ -66,6 +66,15 @@ image_type_issues, review_image_type, ) +from .intensity import ( + INTENSITY_CODES, + IntensityReport, + RepeatBox, + check_intensities, + check_table, + repeat_boxes, + repeat_distributions, +) from .manifest import ( MANIFEST_FILENAME, FileEntry, @@ -151,6 +160,7 @@ def __dir__() -> list[str]: "DESIGN_SCHEMA_VERSION", "FLUORESCENCE", "IMAGE_TYPES", + "INTENSITY_CODES", "MANIFEST_FILENAME", "OUTPUT_FILENAME", "RGB", @@ -174,6 +184,7 @@ def __dir__() -> list[str]: "GroupSummary", "ImageMetadata", "ImageTypeReview", + "IntensityReport", "LeafFolder", "ManifestFormatError", "MetadataReport", @@ -192,6 +203,7 @@ def __dir__() -> list[str]: "PooledHistogram", "PooledStatistics", "Preflight", + "RepeatBox", "RunManifest", "RunResult", "RunSettings", @@ -207,10 +219,12 @@ def __dir__() -> list[str]: "biological_repeats", "check_dataset", "check_image_type", + "check_intensities", "check_mapping", "check_metadata", "check_output_folder", "check_planes", + "check_table", "compute_block_histogram", "compute_plane_histogram", "declared_mapping", @@ -256,6 +270,8 @@ def __dir__() -> list[str]: "read_sample_plane", "reconcile_design", "reference_axes", + "repeat_boxes", + "repeat_distributions", "review_axes", "review_image_type", "robust_z_score", diff --git a/src/morphobricks_preprocessing/cli.py b/src/morphobricks_preprocessing/cli.py index f3ff9a66..58c0e02f 100644 --- a/src/morphobricks_preprocessing/cli.py +++ b/src/morphobricks_preprocessing/cli.py @@ -7,7 +7,9 @@ experiment design confirmed in ``experiment_design.csv``: without it, ``run`` writes a template to the output folder for the user to fill in, and stops. The image type (fluorescence or RGB) is given with ``--image-type``, or the -suggestion from the image headers is used and printed. +suggestion from the image headers is used and printed. After writing, ``run`` +prints the intensity checks of the table (warnings only; plots are drawn in the +notebook only). Both exit with status 1 when an error blocks processing. """ @@ -21,6 +23,8 @@ from datetime import date from pathlib import Path +from morphobricks_normalization.checks import IntensityCheckSettings + from .axes import AxisMapping, reference_axes, review_axes from .design import ( DESIGN_FILENAME, @@ -34,6 +38,7 @@ ) from .discovery import check_output_folder, default_output_folder from .image_type import IMAGE_TYPES, RGB, ImageTypeReview, review_image_type +from .intensity import intensity_text from .planes import DEFAULT_MEMORY_LIMIT, format_bytes from .report import ( balance_summary, @@ -196,6 +201,30 @@ def build_parser() -> argparse.ArgumentParser: help="smaller levels for viewing: auto (default, halve until a plane fits 512 px), " f"none, or a number from 1 to {MAX_LEVELS - 1}", ) + checks = run.add_argument_group( + "intensity checks", "limits of the warnings printed after the run; the defaults suit most" + ) + checks.add_argument( + "--outlier-limit", + type=float, + default=3.0, + metavar="N", + help="robust deviations beyond which an image or biological repeat is reported (default 3)", + ) + checks.add_argument( + "--saturation-limit", + type=float, + default=0.01, + metavar="FRACTION", + help="fraction of saturated pixels above which an image is reported (default 0.01, 1 %%)", + ) + checks.add_argument( + "--min-images", + type=int, + default=5, + metavar="N", + help="fewest images a comparison needs (default 5)", + ) run.add_argument("--log", metavar="FILE", help="write technical details to this file") return parser @@ -296,6 +325,15 @@ def _run(args: argparse.Namespace) -> int: except ValueError as exc: print(exc) return 1 + try: + intensity_settings = IntensityCheckSettings( + outlier_limit=args.outlier_limit, + saturation_limit=args.saturation_limit, + min_images=args.min_images, + ) + except ValueError as exc: + print(exc) + return 1 design = None # Folder problems are reported by the run itself; the design is only asked about # once the folders can be processed. @@ -339,11 +377,15 @@ def progress(done: int, total: int, path: str) -> None: conversion=conversion, design=design, image_type=image_type, + intensity_settings=intensity_settings, ) if result.issues: print(format_issues(result.issues)) print() print(run_summary(result)) + if result.intensity is not None: + print() + print(intensity_text(result.intensity)) return 0 if result.written else 1 diff --git a/src/morphobricks_preprocessing/intensity.py b/src/morphobricks_preprocessing/intensity.py new file mode 100644 index 00000000..5ecc503b --- /dev/null +++ b/src/morphobricks_preprocessing/intensity.py @@ -0,0 +1,277 @@ +"""Intensity checks on the histogram table (``docs/contract.md`` §14, follow-up P5). + +After the histograms are computed, the normalization brick's +:func:`~morphobricks_normalization.summaries.image_summaries` pools every image's +Z and T planes per channel, and +:func:`~morphobricks_normalization.checks.intensity_checks` compares the images +(normalization contract §7, N38): images far from the others of their condition +and biological repeat, biological repeats shifted from the other repeats of their +condition, saturated images and images without intensity spread. Comparisons +stay within a condition, so a treatment effect is never reported. Everything is +computed from the table's rows; no image is read and nothing is written. + +This module also prepares what the notebook plots, channel by channel: the +image medians and high percentiles per condition and biological repeat, with the +flagged images, and the pooled intensity distribution of each biological repeat. +The drawing itself is in :mod:`.notebook_ui`; the command line prints the +warnings only. +""" + +from __future__ import annotations + +import os +from collections.abc import Iterable +from dataclasses import dataclass, field + +import numpy as np +from morphobricks_normalization.checks import IntensityCheckSettings, intensity_checks +from morphobricks_normalization.csv_io import read_plane_histograms +from morphobricks_normalization.issues import ( + IMAGE_INTENSITY_OUTLIER, + NO_INTENSITY_SPREAD, + REPEAT_INTENSITY_SHIFT, + SATURATED_IMAGE, + Issue, +) +from morphobricks_normalization.pooling import NoValidPixelsError, PooledHistogram +from morphobricks_normalization.records import PlaneRecord +from morphobricks_normalization.summaries import ImageSummary, image_summaries + +from .validation import ValidationIssue, plural, warning + +INTENSITY_CODES = ( + IMAGE_INTENSITY_OUTLIER, + REPEAT_INTENSITY_SHIFT, + SATURATED_IMAGE, + NO_INTENSITY_SPREAD, +) +# Findings about one image, which the plots mark on that image. +IMAGE_CODES = (IMAGE_INTENSITY_OUTLIER, SATURATED_IMAGE, NO_INTENSITY_SPREAD) + + +@dataclass(frozen=True) +class IntensityReport: + """The intensity checks of one histogram table. + + Attributes + ---------- + summaries + One per image, sorted by condition, biological repeat, technical repeat + and path. + issues + The normalization brick's warnings, in its order: image outliers, repeat + shifts, saturated images, images without spread. Each has the channel it + is about (``Issue.channel``). + settings + The limits used. + """ + + summaries: tuple[ImageSummary, ...] = () + issues: tuple[Issue, ...] = () + settings: IntensityCheckSettings = field(default_factory=IntensityCheckSettings) + + @property + def channel_count(self) -> int: + """Number of channels of the images; 0 without images.""" + return self.summaries[0].channel_count if self.summaries else 0 + + @property + def high_percentile(self) -> float: + """The high percentile of the summaries, in percent (99.8 by default).""" + return self.summaries[0].percentile[1] if self.summaries else 99.8 + + def flagged_images(self, channel: int) -> dict[str, tuple[str, ...]]: + """Return the codes of the warnings about each image in ``channel``, by path.""" + found: dict[str, list[str]] = {} + for issue in self.issues: + if issue.channel == channel and issue.code in IMAGE_CODES: + for path in issue.items: + codes = found.setdefault(path, []) + if issue.code not in codes: + codes.append(issue.code) + return {path: tuple(codes) for path, codes in found.items()} + + def shifted_repeats(self, channel: int) -> set[tuple[str, str]]: + """Return the ``(condition, biological_repeat)`` reported as shifted in ``channel``.""" + shifted = set() + for issue in self.issues: + if issue.channel == channel and issue.code == REPEAT_INTENSITY_SHIFT: + for item in issue.items: + condition, _, repeat = item.partition("/") + shifted.add((condition, repeat)) + return shifted + + def validation_issues(self) -> list[ValidationIssue]: + """Return the warnings as this brick's validation issues, for printing and showing.""" + return [warning(i.code, i.message, i.items, i.action) for i in self.issues] + + +@dataclass(frozen=True) +class RepeatBox: + """The image values of one condition and biological repeat in one channel, for a box plot. + + Attributes + ---------- + condition, biological_repeat + The folder labels. + paths + The images with a value in the channel. + medians, highs + Their median and high percentile, in the order of ``paths``. + flagged + For each image, the codes of the warnings about it in the channel. + shifted + True when the biological repeat is reported as shifted in the channel. + """ + + condition: str + biological_repeat: str + paths: tuple[str, ...] + medians: tuple[float, ...] + highs: tuple[float, ...] + flagged: tuple[tuple[str, ...], ...] + shifted: bool = False + + @property + def label(self) -> str: + """``condition / biological repeat``, for the plot's axis.""" + return f"{self.condition} / {self.biological_repeat}" + + +def check_intensities( + records: Iterable[PlaneRecord], settings: IntensityCheckSettings | None = None +) -> IntensityReport: + """Summarise every image of the table and run the four intensity checks. + + Parameters + ---------- + records + Rows of ``plane_histograms.csv``. + settings + The limits; None uses the defaults (``outlier_limit`` 3.0, + ``saturation_limit`` 0.01, ``min_images`` 5). + + Returns + ------- + IntensityReport + The summaries and warnings. Images are never read. + """ + settings = IntensityCheckSettings() if settings is None else settings + summaries = tuple(image_summaries(records)) + return IntensityReport(summaries, tuple(intensity_checks(summaries, settings)), settings) + + +def check_table( + path: str | os.PathLike, settings: IntensityCheckSettings | None = None +) -> tuple[list[PlaneRecord], IntensityReport]: + """Read a ``plane_histograms.csv`` and check its intensities. + + Raises + ------ + morphobricks_normalization.csv_io.CsvFormatError + If the file is not a histogram table this version can read. + """ + records = read_plane_histograms(path) + return records, check_intensities(records, settings) + + +def repeat_boxes(report: IntensityReport, channel: int) -> list[RepeatBox]: + """Return the image values of every condition and biological repeat in ``channel``. + + Images without finite pixels in the channel are left out; a repeat with no + such image is left out too. Sorted by condition, then biological repeat. + """ + flagged = report.flagged_images(channel) + shifted = report.shifted_repeats(channel) + groups: dict[tuple[str, str], list[ImageSummary]] = {} + for summary in report.summaries: + if summary.median[channel] is None: + continue + groups.setdefault((summary.condition, summary.biological_repeat), []).append(summary) + return [ + RepeatBox( + condition, + repeat, + tuple(s.source_relative_path for s in members), + tuple(s.median[channel] for s in members), + tuple(s.p_high[channel] for s in members), + tuple(flagged.get(s.source_relative_path, ()) for s in members), + (condition, repeat) in shifted, + ) + for (condition, repeat), members in sorted(groups.items()) + ] + + +def repeat_distributions( + records: Iterable[PlaneRecord], channel: int +) -> dict[str, dict[str, PooledHistogram]]: + """Pool ``channel`` of every biological repeat, condition by condition. + + Returns + ------- + dict + ``{condition: {biological_repeat: PooledHistogram}}``, sorted; repeats + without finite pixels in the channel are left out. + """ + groups: dict[tuple[str, str], list[PlaneRecord]] = {} + for record in records: + groups.setdefault((record.condition, record.biological_repeat), []).append(record) + curves: dict[str, dict[str, PooledHistogram]] = {} + for (condition, repeat), members in sorted(groups.items()): + try: + pooled = PooledHistogram(members, channel=channel) + except NoValidPixelsError: + continue + curves.setdefault(condition, {})[repeat] = pooled + return curves + + +def cumulative_share(pooled: PooledHistogram) -> np.ndarray: + """Share of the pooled pixels at or below each bin's value.""" + return np.cumsum(pooled.counts) / pooled.pixel_count + + +CHECK_TITLES = { + IMAGE_INTENSITY_OUTLIER: "Images far from the others of their biological repeat", + REPEAT_INTENSITY_SHIFT: "Biological repeats shifted from the others of their condition", + SATURATED_IMAGE: "Saturated images", + NO_INTENSITY_SPREAD: "Images without intensity spread", +} + + +def describe_settings(settings: IntensityCheckSettings) -> str: + """Describe the limits in plain words, for the command line and the notebook.""" + return ( + f"more than {settings.outlier_limit:g} robust deviations from the others, more than " + f"{100 * settings.saturation_limit:g} % saturated pixels, comparisons with at least " + f"{settings.min_images} images" + ) + + +def intensity_text(report: IntensityReport) -> str: + """Describe the checks as plain text: the summary, the limits, then each kind of warning. + + Each kind lists its findings (image or repeat, channel ``C`` and the + compared values, as the normalization brick words them) and says once what + to do. + """ + lines = [ + f"Intensity checks: {intensity_summary(report)}", + f" Limits: {describe_settings(report.settings)}.", + ] + for code in INTENSITY_CODES: + found = [issue for issue in report.issues if issue.code == code] + if not found: + continue + lines += ["", f"{CHECK_TITLES[code]} ({code}, {len(found)}):"] + lines += [f" - {issue.message}" for issue in found] + lines.append(f" What to do: {found[0].action}") + return "\n".join(lines) + + +def intensity_summary(report: IntensityReport) -> str: + """Describe the checks in one line, e.g. ``"12 images checked: 2 intensity warnings."``.""" + images = plural(len(report.summaries), "image") + if not report.issues: + return f"{images} checked: no intensity warnings." + return f"{images} checked: {plural(len(report.issues), 'intensity warning')}." diff --git a/src/morphobricks_preprocessing/manifest.py b/src/morphobricks_preprocessing/manifest.py index f90ee47a..7684a286 100644 --- a/src/morphobricks_preprocessing/manifest.py +++ b/src/morphobricks_preprocessing/manifest.py @@ -10,7 +10,8 @@ The manifest also records, for information only, each file's channel names and significant bits when its metadata provides them, the SHA-256 of the experiment design used, so later bricks can check that it was not changed, and -the image type the user confirmed, which later bricks normalize by. +the image type the user confirmed, which later bricks normalize by, and the +limits of the intensity checks run on the table. """ from __future__ import annotations @@ -27,6 +28,7 @@ import morphobricks_normalization import numpy as np import tifffile +from morphobricks_normalization.checks import IntensityCheckSettings from morphobricks_normalization.csv_io import SCHEMA_VERSION from .axes import AxisMapping @@ -173,6 +175,9 @@ class RunManifest: channels_from_rgb_samples True when the channel axis is the RGB samples axis of images stored as RGB colour; None in manifests written before it was recorded. + intensity_checks + The limits of the intensity checks run on the table (:mod:`.intensity`); + None in manifests written before they were recorded. schema_version Version of the CSV format. manifest_version @@ -190,6 +195,7 @@ class RunManifest: experiment_design_sha256: str | None = None image_type: str | None = None channels_from_rgb_samples: bool | None = None + intensity_checks: IntensityCheckSettings | None = None schema_version: str = SCHEMA_VERSION manifest_version: str = MANIFEST_VERSION @@ -236,6 +242,9 @@ def to_dict(self) -> dict: "experiment_design_sha256": self.experiment_design_sha256, "image_type": self.image_type, "channels_from_rgb_samples": self.channels_from_rgb_samples, + "intensity_checks": ( + None if self.intensity_checks is None else self.intensity_checks.to_dict() + ), "versions": dict(self.versions), } @@ -283,6 +292,11 @@ def from_dict(cls, data: dict) -> RunManifest: if data.get("channels_from_rgb_samples") is None else bool(data["channels_from_rgb_samples"]) ), + intensity_checks=( + None + if data.get("intensity_checks") is None + else IntensityCheckSettings.from_dict(data["intensity_checks"]) + ), schema_version=str(data["schema_version"]), ) except ManifestFormatError: diff --git a/src/morphobricks_preprocessing/notebook_ui.py b/src/morphobricks_preprocessing/notebook_ui.py index 27db62bc..0d477543 100644 --- a/src/morphobricks_preprocessing/notebook_ui.py +++ b/src/morphobricks_preprocessing/notebook_ui.py @@ -10,8 +10,10 @@ 2. **Axes** - review the role of every image axis, look at planes with the preview, and confirm; then confirm the image type (fluorescence or RGB). 3. **Run** - check that everything is ready, see what will happen, and compute - or update the histogram table. -4. **Results** - pooled statistics and intensity distributions per group. + or update the histogram table; then the intensity checks of the table, with + their plots. +4. **Results** - pooled statistics and intensity distributions per group, and the + intensity checks of any histogram table. Each step is a page of cards (see :mod:`.notebook_style`): a blue box says what the step is for, every control has a help line, the main button is at the @@ -38,12 +40,15 @@ try: import ipywidgets as widgets from matplotlib.figure import Figure + from matplotlib.lines import Line2D + from matplotlib.patches import Patch except ImportError as exc: # pragma: no cover - depends on the environment raise ImportError( "The notebook interface needs ipywidgets and matplotlib. Install them with: " 'pip install "morphobricks-preprocessing[notebook]"' ) from exc +from morphobricks_normalization.checks import IntensityCheckSettings from morphobricks_normalization.csv_io import OUTPUT_FILENAME, CsvFormatError, read_plane_histograms from morphobricks_normalization.dtypes import ( MAX_BIN_COUNT, @@ -52,7 +57,7 @@ is_integer, ) from morphobricks_normalization.pooling import NoValidPixelsError, PooledHistogram -from morphobricks_normalization.records import group_records +from morphobricks_normalization.records import PlaneRecord, group_records from . import notebook_style as ui from .axes import ( @@ -74,6 +79,18 @@ ) from .discovery import check_output_folder, default_output_folder from .image_type import FLUORESCENCE, RGB, ImageTypeReview, review_image_type +from .intensity import ( + CHECK_TITLES, + INTENSITY_CODES, + IntensityReport, + RepeatBox, + check_intensities, + cumulative_share, + describe_settings, + intensity_summary, + repeat_boxes, + repeat_distributions, +) from .manifest import MANIFEST_FILENAME, ManifestFormatError, read_manifest from .planes import ( DEFAULT_MEMORY_LIMIT, @@ -147,6 +164,13 @@ ("Choose the number of smaller copies", "fixed"), ] _MAX_PLOTTED_GROUPS = 12 +# Short names of the intensity checks, for the tiles and the table's first column. +_CHECK_LABELS = { + "image_intensity_outlier": "Image outlier", + "repeat_intensity_shift": "Shifted repeat", + "saturated_image": "Saturated", + "no_intensity_spread": "No spread", +} _PACKAGE_LOGGER = "morphobricks_preprocessing" _WIDE = {"width": "99%"} _PAIRING_COLUMNS = widgets.Layout( @@ -371,6 +395,7 @@ def run_step(self) -> widgets.Widget: self.run_status, self.run_tiles, self.run_issues, + self.run_checks.card, self.details_box, self.run_next, ], @@ -395,6 +420,7 @@ def results_step(self) -> widgets.Widget: self.stats_tiles, self.statistics_card, self.distribution_card, + self.table_checks.card, self.stats_next, ], ) @@ -1608,6 +1634,14 @@ def _build_run_step(self) -> None: self.details_box = widgets.Accordion( children=[self.details], titles=("Technical details (log)",), selected_index=None ) + self.run_checks = IntensityPanel( + "Intensity checks", + "Computed from the table just written, without reading the images again. Each " + "image is compared only with the others of its condition, so differences between " + "conditions (treatment effects) are never reported. The run records the limits " + "below in run_manifest.json; changing them here updates the warnings and " + "plots at once and is recorded by the next run.", + ) self.recommended.observe(lambda _: self._recommended_changed(), names="value") for control in (self.bins, self.memory_gib, self.force_rebuild): control.observe(lambda _: self._refresh_run_step(), names="value") @@ -2030,6 +2064,7 @@ def run(self) -> RunResult | None: return None self._busy(self.run_button, "Computing…") self.run_next.value = "" + self.run_checks.clear() ui.set_visible(self.progress, True) self.progress.value, self.progress.max, self.progress.bar_style = 0, 1, "info" self.run_status.value = ui.working("Preparing: checking everything again…") @@ -2050,6 +2085,7 @@ def run(self) -> RunResult | None: conversion=self._conversion_settings(), design=self.design, image_type=self.confirmed_image_type, + intensity_settings=self.run_checks.settings(), ) except Exception as exc: # unexpected failures still get a readable message self.log.lines.append(traceback.format_exc()) @@ -2067,6 +2103,8 @@ def run(self) -> RunResult | None: finally: self._idle(self.run_button, " Compute histograms", "play") self._show_run_result(self.result) + if self.result.written: + self._show_run_checks(self.result.csv_path) self._show_log() if self.result.written: self.csv_file.value = str(self.result.csv_path) @@ -2121,6 +2159,24 @@ def _show_run_result(self, result: RunResult) -> None: tiles.append(("⚠️", issue_counts(result.issues), "checks")) self.run_tiles.value = ui.tiles(tiles) + def _show_run_checks(self, path: Path) -> None: + """Show the intensity checks and plots of the table the run wrote.""" + try: + with self.log.attached(): + self.run_checks.show(read_plane_histograms(path)) + except Exception as exc: # the table is written; only the checks failed + self.log.lines.append(traceback.format_exc()) + self.run_checks.fail(exc) + return + report = self.run_checks.report + if report is not None and report.issues: + self.run_next.value = ui.next_steps( + "Read the intensity checks above, then go to Step 4 (optional).", + "Each warning names the image or biological repeat, the channel and the values " + "compared, and says what to do. The histogram table is ready for later " + "MorphoBricks bricks either way.", + ) + def _convert_tiffs(self) -> bool: return self.conversion.value == "convert" @@ -2217,6 +2273,12 @@ def _build_results_step(self) -> None: ], visible=False, ) + self.table_checks = IntensityPanel( + "Intensity checks", + "The same checks as after step 3, for the chosen table. Each image is compared only " + "with the others of its condition, so differences between conditions are never " + "reported. Nothing is written.", + ) self.statistics_button.on_click(lambda _: self.show_statistics()) self.csv_browse.on_click(lambda _: self._browse_csv()) self.csv_file.observe(lambda _: self._csv_changed(), names="value") @@ -2241,6 +2303,7 @@ def _csv_changed(self) -> None: self.stats_next.value = "" ui.set_visible(self.statistics_card, False) ui.set_visible(self.distribution_card, False) + self.table_checks.clear() def _statistics_option_changed(self) -> None: if self.statistics.value: @@ -2328,13 +2391,15 @@ def _calculate_statistics(self, path: Path, percentiles: Sequence[float]) -> Non png = distribution_png(groups) self.distribution.value = png or b"" ui.set_visible(self.distribution_card, bool(png)) + self.stats_status.value = ui.working("Checking the intensities of every image…") + self.table_checks.show(records) self.stats_status.value = ui.status( "ok", f"Statistics of {plural(len(summaries), LEVEL_NAMES[level])} calculated.", [f"From the table {path}"], ) self.stats_next.value = ui.next_steps( - "Review the group statistics and intensity distributions above.", + "Review the group statistics, intensity distributions and intensity checks above.", "Choose another grouping to compare repeats or images. The histogram table " "is ready for later MorphoBricks bricks; this optional step writes no files.", ) @@ -2370,9 +2435,358 @@ def _show_log(self) -> None: self.details.value = "\n".join(self.log.lines[-500:]) +class IntensityPanel: + """The intensity checks of one histogram table: warnings, limits and plots (follow-up P5). + + Shown in step 3 for the table just written and in step 4 for any table. + The checks come from :func:`.intensity.check_intensities`; the plots are + drawn per channel and nothing is written. + + Parameters + ---------- + title + Card header. + explanation + Help line under it (HTML). + """ + + def __init__(self, title: str, explanation: str): + self.records: list[PlaneRecord] = [] + self.report: IntensityReport | None = None + defaults = IntensityCheckSettings() + self.outlier_limit = widgets.BoundedFloatText( + value=defaults.outlier_limit, + min=0.5, + max=100, + step=0.5, + description="Outlier limit", + style=_label(), + ) + self.saturation_percent = widgets.BoundedFloatText( + value=100 * defaults.saturation_limit, + min=0, + max=99.9, + step=0.5, + description="Saturated pixels (%)", + style=_label(), + ) + self.min_images = widgets.BoundedIntText( + value=defaults.min_images, + min=2, + max=100_000, + description="Fewest images", + style=_label(), + ) + self.defaults_button = ui.button( + " Restore the defaults", + style="warning", + icon="undo", + tooltip="3 robust deviations, 1 % saturated pixels, at least 5 images", + ) + self.settings_box = widgets.Accordion( + children=[ + widgets.VBox( + [ + self.outlier_limit, + ui.help_text( + "How far, in robust deviations (the spread of the others, " + "insensitive to single outliers), an image's median or high " + "percentile may be from the other images of its condition and " + "biological repeat, or a biological repeat from the other repeats " + "of its condition. Default 3; larger values report fewer." + ), + self.saturation_percent, + ui.help_text( + "Share of an image's pixels at the largest value of its pixel type " + "(for example 255 in 8-bit images) above which it is reported as " + "saturated. Default 1 %." + ), + self.min_images, + ui.help_text( + "Comparisons with fewer images are skipped, since a few images say " + "little about what is usual. Default 5." + ), + self.defaults_button, + ] + ) + ], + titles=("Check settings",), + selected_index=None, + ) + self.status = widgets.HTML() + self.tiles = widgets.HTML() + self.findings = widgets.HTML() + self.channel = widgets.ToggleButtons(options=[("C0", 0)], style={"button_width": "64px"}) + self.channel_row = widgets.VBox( + [ + widgets.HTML("Channel"), + self.channel, + ui.help_text("The plots show one channel at a time."), + ] + ) + self.boxes = widgets.Image(format="png", layout=widgets.Layout(max_width="760px")) + self.curves = widgets.Image(format="png", layout=widgets.Layout(max_width="760px")) + self.plots = widgets.VBox( + [ + self.channel_row, + widgets.HTML("Image medians and high percentiles"), + ui.help_text( + "One box per condition and biological repeat, one point per image. Red " + "diamonds are images with a warning in this channel; an amber box is a " + "biological repeat shifted from the others of its condition." + ), + self.boxes, + widgets.HTML("Intensity distributions per biological repeat"), + ui.help_text( + "One panel per condition, one curve per biological repeat: the share of " + "pixels at or below each intensity. Curves of comparable repeats lie close " + "together; a dashed curve is a shifted repeat." + ), + self.curves, + ] + ) + self.card = ui.card( + title, + [ + ui.help_text(explanation), + self.status, + self.tiles, + self.findings, + self.settings_box, + self.plots, + ], + visible=False, + ) + for control in (self.outlier_limit, self.saturation_percent, self.min_images): + control.observe(lambda _: self._settings_changed(), names="value") + self.defaults_button.on_click(lambda _: self.restore_defaults()) + self.channel.observe(lambda _: self._draw(), names="value") + + def settings(self) -> IntensityCheckSettings: + """Return the limits as currently set.""" + return IntensityCheckSettings( + outlier_limit=self.outlier_limit.value, + saturation_limit=self.saturation_percent.value / 100, + min_images=self.min_images.value, + ) + + def restore_defaults(self) -> None: + """Set the default limits again.""" + defaults = IntensityCheckSettings() + self.outlier_limit.value = defaults.outlier_limit + self.saturation_percent.value = 100 * defaults.saturation_limit + self.min_images.value = defaults.min_images + + def show(self, records: Sequence[PlaneRecord]) -> None: + """Check the intensities of ``records`` and show the warnings and plots.""" + self.records = list(records) + channels = self.records[0].channel_count if self.records else 1 + with self.channel.hold_trait_notifications(): + self.channel.options = [(f"C{c}", c) for c in range(channels)] + self.channel.value = 0 + ui.set_visible(self.channel_row, channels > 1) + self._check() + ui.set_visible(self.card, True) + + def clear(self) -> None: + """Hide the panel and forget the table.""" + self.records, self.report = [], None + self.status.value = self.tiles.value = self.findings.value = "" + self.boxes.value = self.curves.value = b"" + ui.set_visible(self.card, False) + + def fail(self, exc: Exception) -> None: + """Show that the checks could not be made.""" + self.clear() + self.status.value = ui.status( + "error", + "The intensity checks could not be made.", + [f"{type(exc).__name__}: {exc}", "Open Technical details for the full report."], + ) + ui.set_visible(self.plots, False) + ui.set_visible(self.card, True) + + def _settings_changed(self) -> None: + if self.records: + self._check() + + def _check(self) -> None: + self.report = check_intensities(self.records, self.settings()) + report = self.report + lines = [ + intensity_summary(report), + f"Limits: {describe_settings(report.settings)} (change them under Check settings).", + ] + if report.issues: + self.status.value = ui.status( + "warn", + f"{plural(len(report.issues), 'intensity warning')} to read.", + lines, + ) + else: + self.status.value = ui.status("ok", "No intensity warnings.", lines) + counts = {code: 0 for code in INTENSITY_CODES} + for issue in report.issues: + counts[issue.code] += 1 + self.tiles.value = ui.tiles( + [ + ("🖼️", f"{len(report.summaries):,}", "images checked"), + ("🎯", f"{counts['image_intensity_outlier']:,}", "image outliers"), + ("📈", f"{counts['repeat_intensity_shift']:,}", "shifted repeats"), + ("☀️", f"{counts['saturated_image']:,}", "saturated"), + ("➖", f"{counts['no_intensity_spread']:,}", "without spread"), + ] + ) + self.findings.value = intensity_findings_html(report) + ui.set_visible(self.plots, bool(report.summaries)) + self._draw() + + def _draw(self) -> None: + report = self.report + if report is None or not report.summaries: + self.boxes.value = self.curves.value = b"" + return + channel = self.channel.value or 0 + boxes = repeat_boxes(report, channel) + self.boxes.value = intensity_boxes_png(boxes, report.high_percentile, channel) or b"" + curves = repeat_distributions(self.records, channel) + self.curves.value = ( + repeat_curves_png(curves, report.shifted_repeats(channel), channel) or b"" + ) + + # Rendering helpers --------------------------------------------------------------------------- +def intensity_findings_html(report: IntensityReport) -> str: + """Render the intensity warnings as a table, then what to do about each kind.""" + if not report.issues: + return "" + rows = [["Check", "Channel", "Image or biological repeat", "What was found"]] + for issue in report.issues: + rows.append( + [_CHECK_LABELS[issue.code], f"C{issue.channel}", ", ".join(issue.items), issue.message] + ) + actions = [] + for code in INTENSITY_CODES: + found = [issue for issue in report.issues if issue.code == code] + if found: + actions.append( + f"{html.escape(CHECK_TITLES[code])}: {html.escape(found[0].action)}" + ) + return ui.table(rows) + ui.status("info", "What to do", actions, raw=True) + + +def intensity_boxes_png( + boxes: Sequence[RepeatBox], high_percentile: float, channel: int +) -> bytes | None: + """Box plots of the image medians and high percentiles per condition and biological repeat. + + Images with a warning in the channel are red diamonds; a biological repeat + reported as shifted has an amber box. + """ + if not boxes: + return None + n = len(boxes) + figure = Figure(figsize=(max(7.2, 0.6 * n + 2.5), 5.6)) + pair = figure.subplots(2, 1, sharex=True) + rng = np.random.default_rng(0) # the same jitter on every redraw + jitter = [rng.uniform(-0.15, 0.15, len(box.paths)) for box in boxes] + statistics = ( + ("Image median", lambda box: box.medians), + (f"Image {high_percentile:g}th percentile", lambda box: box.highs), + ) + for axes, (name, values_of) in zip(pair, statistics, strict=True): + parts = axes.boxplot( + [values_of(box) for box in boxes], + positions=list(range(n)), + widths=0.6, + patch_artist=True, + showfliers=False, + # Dark, so the median line is not taken for the amber of a shifted repeat. + medianprops={"color": "#333333"}, + ) + for patch, box in zip(parts["boxes"], boxes, strict=True): + patch.set_facecolor("#f6c453" if box.shifted else "#dbe7f5") + patch.set_edgecolor("#555555") + for i, box in enumerate(boxes): + values = np.asarray(values_of(box), dtype=float) + flagged = np.array([bool(codes) for codes in box.flagged]) + x = i + jitter[i] + axes.scatter(x[~flagged], values[~flagged], s=12, color="#3b6ea8", zorder=3) + axes.scatter(x[flagged], values[flagged], s=40, color="#d62728", marker="D", zorder=4) + if i and box.condition != boxes[i - 1].condition: + axes.axvline(i - 0.5, color="#999999", linestyle="--", linewidth=0.8) + axes.set_ylabel(name, fontsize=9) + axes.grid(axis="y", alpha=0.3) + pair[1].set_xticks( + list(range(n)), [box.label for box in boxes], rotation=30, ha="right", fontsize=8 + ) + pair[0].set_title(f"Channel C{channel}: one point per image", fontsize=10) + pair[0].legend( + handles=[ + Line2D([], [], color="#d62728", marker="D", linestyle="", label="image with a warning"), + Patch(facecolor="#f6c453", edgecolor="#555555", label="shifted biological repeat"), + ], + fontsize=8, + frameon=False, + loc="best", + ) + figure.tight_layout() + return _png(figure) + + +def repeat_curves_png(curves: dict, shifted: set[tuple[str, str]], channel: int) -> bytes | None: + """Overlay the cumulative intensity distribution of each biological repeat, per condition. + + ``curves`` is :func:`.intensity.repeat_distributions`; shifted repeats are + dashed. At most 12 conditions are drawn. + """ + conditions = list(curves)[:_MAX_PLOTTED_GROUPS] + if not conditions: + return None + pooled = [p for condition in conditions for p in curves[condition].values()] + low = min(p.percentile(0.1) for p in pooled) + high = max(p.percentile(99.9) for p in pooled) + columns = min(3, len(conditions)) + rows = math.ceil(len(conditions) / columns) + figure = Figure(figsize=(3.7 * columns + 0.5, 2.9 * rows + 0.6)) + for index, condition in enumerate(conditions): + axes = figure.add_subplot(rows, columns, index + 1) + for repeat, histogram in curves[condition].items(): + dashed = (condition, repeat) in shifted + axes.step( + histogram.values, + cumulative_share(histogram), + where="post", + linewidth=1.4, + linestyle="--" if dashed else "-", + label=f"{repeat} (shifted)" if dashed else repeat, + ) + if high > low: + axes.set_xlim(low, high) + axes.set_ylim(0, 1) + axes.set_title(condition, fontsize=9) + axes.set_xlabel("Intensity", fontsize=8) + if index % columns == 0: + axes.set_ylabel("Share of pixels at or below", fontsize=8) + axes.tick_params(labelsize=7) + axes.grid(alpha=0.3) + axes.legend(fontsize=7, frameon=False, loc="lower right") + title = f"Channel C{channel}: one curve per biological repeat" + if len(curves) > _MAX_PLOTTED_GROUPS: + title += f" (first {_MAX_PLOTTED_GROUPS} of {len(curves)} conditions)" + figure.suptitle(title, fontsize=10) + figure.tight_layout() + return _png(figure) + + +def _png(figure: Figure) -> bytes: + buffer = io.BytesIO() + figure.savefig(buffer, format="png", dpi=100) + return buffer.getvalue() + + def issues_html(issues: Iterable[ValidationIssue]) -> str: """Render issues as status cards: errors first, each with its paths and what to do.""" issues = sorted(issues, key=lambda issue: not issue.is_error) diff --git a/src/morphobricks_preprocessing/run.py b/src/morphobricks_preprocessing/run.py index 0e36de44..c4603918 100644 --- a/src/morphobricks_preprocessing/run.py +++ b/src/morphobricks_preprocessing/run.py @@ -7,7 +7,9 @@ run, and writes ``plane_histograms.csv`` and ``run_manifest.json`` to the output folder (``docs/plan.md`` §5–§7). Given the confirmed experiment design, it also writes ``experiment_design.csv`` and records its hash in the manifest; given the -confirmed image type, it records it there too. +confirmed image type, it records it there too. After writing, it checks the +images' intensities from the table's rows (:mod:`.intensity`) and records the +limits used. Each row depends only on its own plane, so rows of unchanged files are kept as they are. Everything is recomputed when there is no usable earlier output, when @@ -28,6 +30,7 @@ from datetime import UTC, datetime from pathlib import Path +from morphobricks_normalization.checks import IntensityCheckSettings from morphobricks_normalization.csv_io import ( OUTPUT_FILENAME, SCHEMA_VERSION, @@ -49,6 +52,7 @@ ) from .discovery import DiscoveryResult, default_output_folder, discover_experiment from .image_type import check_image_type, image_type_issues, review_image_type +from .intensity import IntensityReport, check_intensities from .manifest import ( MANIFEST_FILENAME, FileEntry, @@ -164,6 +168,9 @@ class RunResult: image_type The image type recorded in the manifest, or None when the run was given none. + intensity + The intensity checks of the written table (warnings, never errors), or + None when nothing was written. """ output_folder: Path @@ -180,6 +187,7 @@ class RunResult: converted: tuple[str, ...] = () design: ExperimentDesign | None = None image_type: str | None = None + intensity: IntensityReport | None = None @property def ok(self) -> bool: @@ -312,6 +320,7 @@ def run_experiment( conversion: ConversionSettings | None = None, design: ExperimentDesign | None = None, image_type: str | None = None, + intensity_settings: IntensityCheckSettings | None = None, ) -> RunResult: """Pass 2: write or update the histogram table and the run manifest. @@ -363,6 +372,10 @@ def run_experiment( are not three channels of an integer type is the error ``rgb_not_possible``. It changes no row, so changing it recomputes nothing. Without it the manifest records none. + intensity_settings + Limits of the intensity checks run on the written table + (:func:`.intensity.check_intensities`); None uses the defaults. They are + recorded in the manifest and change no row. Returns ------- @@ -542,6 +555,8 @@ def stop(**counts) -> RunResult: return stop(read=tuple(to_read), **counts) records.sort(key=_row_order) + # From the rows in memory: no image is read again. + intensity = check_intensities(records, intensity_settings) per_file: dict[str, int] = {} for record in records: per_file[record.source_relative_path] = per_file.get(record.source_relative_path, 0) + 1 @@ -581,6 +596,7 @@ def stop(**counts) -> RunResult: experiment_design_sha256=design_sha256, image_type=image_type, channels_from_rgb_samples=type_review.channels_from_rgb_samples, + intensity_checks=intensity.settings, ) write_manifest(manifest, staged / MANIFEST_FILENAME) logger.info( @@ -600,6 +616,7 @@ def stop(**counts) -> RunResult: converted=tuple(to_convert), design=design, image_type=image_type, + intensity=intensity, **counts, ) diff --git a/tests/test_cli.py b/tests/test_cli.py index aba3de5d..474f1be1 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -132,6 +132,43 @@ def test_image_type_must_be_known(experiment): assert excinfo.value.code == 2 +def bright_tree(root): + """Six 2-channel images in ctrl/rep1, one far brighter in C0, and drug/rep1 brighter.""" + rng = np.random.default_rng(3) + second = np.clip(rng.normal(60, 10, (16, 16)), 0, 254).astype(np.uint8) + levels = {"ctrl/rep1": [100, 102, 98, 101, 99, 180], "drug/rep1": [200, 202, 198]} + for folder, values in levels.items(): + for i, level in enumerate(values): + first = np.clip(rng.normal(level, 10, (16, 16)), 0, 254).astype(np.uint8) + write_data(root / folder / f"img_{i}.tif", np.stack([first, second]), axes="CYX") + return root + + +def test_run_prints_the_intensity_checks(tmp_path, capsys): + root = bright_tree(tmp_path / "exp") + confirm_design(root) + assert main(["run", str(root), "--axes", "CYX"]) == 0 + out = capsys.readouterr().out + assert "Intensity checks: 9 images checked: " in out + assert "Limits: more than 3 robust deviations from the others, more than 1 % saturated" in out + assert "Images far from the others of their biological repeat (image_intensity_outlier" in out + assert " - Image ctrl/rep1/img_5.tif, channel C0: its " in out + assert out.count("What to do: Look at the image") == 1 + # drug is brighter than ctrl, but conditions are never compared. + assert "drug/rep1" not in out.split("Intensity checks:")[1] + args = ["run", str(root), "--axes", "CYX", "--outlier-limit", "50"] + args += ["--saturation-limit", "0.2", "--min-images", "3"] + assert main(args) == 0 + assert "no intensity warnings" in capsys.readouterr().out + limits = read_manifest(default_output_folder(root) / "run_manifest.json").intensity_checks + assert (limits.outlier_limit, limits.saturation_limit, limits.min_images) == (50, 0.2, 3) + + +def test_invalid_check_limits_are_refused(experiment, capsys): + assert main(["run", str(experiment), "--axes", "ZCYX", "--min-images", "1"]) == 1 + assert "min_images must be at least 2" in capsys.readouterr().out + + def test_run_options_reach_the_run(experiment, tmp_path, capsys): out_dir = tmp_path / "results" log = tmp_path / "run.log" diff --git a/tests/test_intensity.py b/tests/test_intensity.py new file mode 100644 index 00000000..4d8e0153 --- /dev/null +++ b/tests/test_intensity.py @@ -0,0 +1,177 @@ +import numpy as np +import pytest +from morphobricks_normalization.checks import IntensityCheckSettings, intensity_checks +from morphobricks_normalization.csv_io import write_plane_histograms +from morphobricks_normalization.histogram import compute_plane_histogram +from morphobricks_normalization.records import PlaneRecord +from morphobricks_normalization.summaries import image_summaries + +from morphobricks_preprocessing.intensity import ( + INTENSITY_CODES, + check_intensities, + check_table, + cumulative_share, + intensity_summary, + repeat_boxes, + repeat_distributions, +) + + +def record(path, condition, repeat, channels, z=0): + """One table row: ``channels`` are 2D arrays, one per channel.""" + return PlaneRecord( + source_relative_path=path, + condition=condition, + biological_repeat=repeat, + technical_repeat="single", + histograms=[compute_plane_histogram(c) for c in channels], + z=z, + ) + + +def plane(level, rng, spread=10, shape=(16, 16)): + return np.clip(rng.normal(level, spread, shape), 0, 254).astype(np.uint8) + + +def rows(levels_by_group, seed=0, saturated=(), blank=()): + """Two channels per image: channel 0 at the given level, channel 1 the same in every image.""" + rng = np.random.default_rng(seed) + second = plane(60, np.random.default_rng(99)) + records = [] + for (condition, repeat), levels in levels_by_group.items(): + for i, level in enumerate(levels): + path = f"{condition}/{repeat}/img_{i}.tif" + first = plane(level, rng) + if path in saturated: + first[:4] = 255 + if path in blank: + first[:] = 7 + for z in range(2): + records.append(record(path, condition, repeat, [first, second], z)) + return records + + +def codes(issues): + return [issue.code for issue in issues] + + +@pytest.fixture +def table(): + # ctrl/rep1 has an outlier image in C0; drug is brighter overall (a treatment + # effect, never reported); ctrl/rep3 is shifted from the other repeats of ctrl. + return rows( + { + ("ctrl", "rep1"): [100, 102, 98, 101, 99, 180], + ("ctrl", "rep2"): [100, 103, 97, 101, 99], + ("ctrl", "rep3"): [150, 152, 148, 151, 149], + ("drug", "rep1"): [200, 202, 198, 201, 199], + ("drug", "rep2"): [200, 203, 197, 201, 199], + }, + saturated=("drug/rep1/img_0.tif",), + blank=("drug/rep2/img_4.tif",), + ) + + +def test_the_checks_are_the_normalization_bricks(table): + report = check_intensities(table) + summaries = image_summaries(table) + assert report.summaries == tuple(summaries) + assert report.issues == tuple(intensity_checks(summaries)) + assert set(codes(report.issues)) == set(INTENSITY_CODES) + assert report.channel_count == 2 and report.high_percentile == 99.8 + + +def test_each_finding_names_the_channel(table): + report = check_intensities(table) + outliers = [i for i in report.issues if i.code == "image_intensity_outlier"] + assert {i.items for i in outliers} >= {("ctrl/rep1/img_5.tif",)} + assert {i.channel for i in outliers} == {0} + assert all("channel C0" in i.message for i in outliers) + shifts = [i for i in report.issues if i.code == "repeat_intensity_shift"] + assert ("ctrl/rep3",) in {i.items for i in shifts} + # Conditions are never compared with each other. + assert not any(i.items[0].startswith("drug/") for i in shifts) + + +def test_settings_change_the_findings(table): + loose = check_intensities(table, IntensityCheckSettings(outlier_limit=1000)) + assert not {"image_intensity_outlier", "repeat_intensity_shift"} & set(codes(loose.issues)) + strict = check_intensities(table, IntensityCheckSettings(saturation_limit=0.5)) + assert "saturated_image" not in codes(strict.issues) + few = check_intensities(table, IntensityCheckSettings(min_images=7)) + assert "image_intensity_outlier" not in codes(few.issues) + assert few.settings.min_images == 7 + + +def test_flagged_images_and_shifted_repeats_per_channel(table): + report = check_intensities(table) + flagged = report.flagged_images(0) + assert "image_intensity_outlier" in flagged["ctrl/rep1/img_5.tif"] + assert "saturated_image" in flagged["drug/rep1/img_0.tif"] + assert "no_intensity_spread" in flagged["drug/rep2/img_4.tif"] + assert report.flagged_images(1) == {} + assert ("ctrl", "rep3") in report.shifted_repeats(0) + assert report.shifted_repeats(1) == set() + + +def test_validation_issues_are_warnings_with_paths(table): + report = check_intensities(table) + converted = report.validation_issues() + assert len(converted) == len(report.issues) + assert all(not issue.is_error for issue in converted) + first = converted[0] + assert first.code == report.issues[0].code and first.paths == report.issues[0].items + assert first.action == report.issues[0].action + + +def test_repeat_boxes(table): + report = check_intensities(table) + boxes = repeat_boxes(report, 0) + assert [b.label for b in boxes] == [ + "ctrl / rep1", + "ctrl / rep2", + "ctrl / rep3", + "drug / rep1", + "drug / rep2", + ] + rep1 = boxes[0] + assert len(rep1.paths) == len(rep1.medians) == len(rep1.highs) == 6 + assert rep1.flagged[5] and not any(rep1.flagged[:5]) + assert [b.shifted for b in boxes] == [False, False, True, False, False] + assert not any(any(b.flagged) for b in repeat_boxes(report, 1)) + + +def test_repeat_distributions_pool_one_channel_per_repeat(table): + curves = repeat_distributions(table, 0) + assert list(curves) == ["ctrl", "drug"] and list(curves["ctrl"]) == ["rep1", "rep2", "rep3"] + pooled = curves["ctrl"]["rep2"] + # Five images, two planes each, 256 pixels per plane. + assert pooled.pixel_count == 5 * 2 * 256 + share = cumulative_share(pooled) + assert share[-1] == pytest.approx(1.0) and np.all(np.diff(share) >= 0) + + +def test_channels_without_finite_pixels_are_left_out(): + rng = np.random.default_rng(1) + empty = np.full((8, 8), np.nan, np.float32) + good = rng.normal(5, 1, (8, 8)).astype(np.float32) + table = [ + record("a/r/x.tif", "a", "r", [empty, good]), + record("a/r/y.tif", "a", "r", [good, good]), + ] + report = check_intensities(table) + assert [b.paths for b in repeat_boxes(report, 0)] == [("a/r/y.tif",)] + assert list(repeat_distributions(table, 0)["a"]) == ["r"] + + +def test_check_table_reads_the_csv(table, tmp_path): + path = write_plane_histograms(table, tmp_path / "plane_histograms.csv") + records, report = check_table(path) + assert len(records) == len(table) + assert report.issues == check_intensities(table).issues + + +def test_intensity_summary(table): + assert intensity_summary(check_intensities(table)).startswith("26 images checked: ") + quiet = check_intensities(rows({("a", "r"): [100, 100]})) + assert intensity_summary(quiet) == "2 images checked: no intensity warnings." diff --git a/tests/test_manifest.py b/tests/test_manifest.py index 7ac957c5..806b04be 100644 --- a/tests/test_manifest.py +++ b/tests/test_manifest.py @@ -156,6 +156,20 @@ def test_image_type_round_trips_and_is_none_when_absent(tmp_path): assert read.image_type is None and read.channels_from_rgb_samples is None +def test_intensity_check_limits_round_trip_and_are_none_when_absent(tmp_path): + from morphobricks_normalization.checks import IntensityCheckSettings + + limits = IntensityCheckSettings(outlier_limit=4.5, saturation_limit=0.02, min_images=3) + path = write_manifest(manifest(intensity_checks=limits), tmp_path / "m.json") + assert read_manifest(path).intensity_checks == limits + data = json.loads(path.read_bytes()) + assert data["intensity_checks"]["outlier_limit"] == 4.5 + assert "library_version" in data["intensity_checks"] + del data["intensity_checks"] # written before P5 + path.write_text(json.dumps(data), encoding="utf-8") + assert read_manifest(path).intensity_checks is None + + def test_file_entry_for_reads_size_and_time(tmp_path): path = write_tiff(tmp_path / "a.tif") found = file_entry_for(tmp_path, "a.tif", ("single",) * 3, prefill_mapping("YX")) diff --git a/tests/test_notebook_ui.py b/tests/test_notebook_ui.py index f2498121..64058bfd 100644 --- a/tests/test_notebook_ui.py +++ b/tests/test_notebook_ui.py @@ -718,6 +718,87 @@ def test_changing_the_axis_roles_asks_for_the_image_type_again(experiment): assert nb.confirmed_image_type is None and not ui.is_visible(nb.image_type_card) +# Intensity checks (follow-up P5) ---------------------------------------------------------- + + +def bright_tree(root): + """Six 2-channel images in ctrl/rep1, one far brighter in C0, and drug/rep1 brighter.""" + rng = np.random.default_rng(3) + second = np.clip(rng.normal(60, 10, (16, 16)), 0, 254).astype(np.uint8) + levels = {"ctrl/rep1": [100, 102, 98, 101, 99, 180], "drug/rep1": [200, 202, 198]} + for folder, values in levels.items(): + for i, level in enumerate(values): + first = np.clip(rng.normal(level, 10, (16, 16)), 0, 254).astype(np.uint8) + write_data(root / folder / f"img_{i}.tif", np.stack([first, second]), axes="CYX") + return root + + +def test_intensity_checks_follow_a_successful_run(tmp_path): + out = tmp_path / "out" + nb = confirmed(bright_tree(tmp_path / "exp"), output_folder=str(out)) + assert not ui.is_visible(nb.run_checks.card) + nb.run() + # The checks and plots write no file. + assert sorted(p.name for p in out.iterdir()) == [ + "experiment_design.csv", + "plane_histograms.csv", + "run_manifest.json", + ] + panel = nb.run_checks + assert ui.is_visible(panel.card) and "intensity warning" in panel.status.value + assert "9 images checked" in panel.status.value + assert "ctrl/rep1/img_5.tif" in panel.findings.value and "C0" in panel.findings.value + assert "What to do" in panel.findings.value + assert "drug/rep1" not in panel.findings.value # conditions are never compared + assert bytes(panel.boxes.value).startswith(b"\x89PNG") + assert bytes(panel.curves.value).startswith(b"\x89PNG") + assert ui.is_visible(panel.channel_row) and [o[1] for o in panel.channel.options] == [0, 1] + assert "Read the intensity checks above" in nb.run_next.value + assert panel.report.issues == nb.result.intensity.issues + # Channel C1 has no warning; its plots are drawn too. + first = bytes(panel.boxes.value) + panel.channel.value = 1 + assert bytes(panel.boxes.value) != first + # Changing a limit updates the warnings at once, without writing anything. + stamp = (out / "run_manifest.json").stat().st_mtime_ns + panel.outlier_limit.value = 100 + assert "No intensity warnings" in panel.status.value and panel.findings.value == "" + assert (out / "run_manifest.json").stat().st_mtime_ns == stamp + panel.defaults_button.click() + assert panel.outlier_limit.value == 3 and "intensity warning" in panel.status.value + # The next run records the limits shown. + panel.min_images.value = 3 + nb.run() + assert read_manifest(out / "run_manifest.json").intensity_checks.min_images == 3 + + +def test_step_4_checks_any_table_without_its_images(tmp_path): + root = bright_tree(tmp_path / "exp") + table = confirmed(root, output_folder=str(tmp_path / "out")).run().csv_path + for image in root.rglob("*.tif"): # the checks and plots need only the table + image.unlink() + nb = PreprocessingNotebook() + nb.results_step() + nb.csv_file.value = str(table) + nb.show_statistics() + panel = nb.table_checks + assert ui.is_visible(panel.card) and "ctrl/rep1/img_5.tif" in panel.findings.value + assert bytes(panel.curves.value).startswith(b"\x89PNG") + assert "intensity checks" in nb.stats_next.value + nb.csv_file.value = "" # another table: the checks of this one go + assert not ui.is_visible(panel.card) and panel.report is None + + +def test_intensity_plots_mark_flagged_images_and_skip_empty_input(): + from morphobricks_preprocessing.intensity import RepeatBox + from morphobricks_preprocessing.notebook_ui import intensity_boxes_png, repeat_curves_png + + box = RepeatBox("a", "r", ("x", "y"), (1.0, 5.0), (2.0, 9.0), ((), ("saturated_image",)), True) + assert intensity_boxes_png([box], 99.8, 0).startswith(b"\x89PNG") + assert intensity_boxes_png([], 99.8, 0) is None + assert repeat_curves_png({}, set(), 0) is None + + def test_conflicting_axes_explain_why_no_roles_are_offered(tmp_path): write_data(tmp_path / "exp/a.tif", np.zeros((2, 8, 8), np.uint8), axes="ZYX") write_data(tmp_path / "exp/b.tif", np.zeros((2, 8, 8), np.uint8), axes="CYX") diff --git a/tests/test_run.py b/tests/test_run.py index f55cfadc..d5d8fa45 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -109,6 +109,38 @@ def csv_bytes(result): # First run ---------------------------------------------------------------------------- +def test_the_run_checks_the_intensities_of_its_table(tmp_path, reads): + from morphobricks_normalization.checks import IntensityCheckSettings + + from morphobricks_preprocessing import check_table + + root = tmp_path / "exp" + rng = np.random.default_rng(5) + for i, level in enumerate([100, 101, 99, 102, 98, 190]): + image = np.clip(rng.normal(level, 8, (2, 16, 16)), 0, 254).astype(np.uint8) + write_data(root / f"ctrl/rep1/img_{i}.tif", image, axes="ZYX") + result = run_experiment(root, ZYX) + assert result.written and result.intensity is not None + _, from_table = check_table(result.csv_path) + assert result.intensity.issues == from_table.issues + assert {i.items for i in result.intensity.issues} == {("ctrl/rep1/img_5.tif",)} + # The checks are warnings about the table, not issues of the run. + assert "image_intensity_outlier" not in codes(result.issues) + manifest = read_manifest(result.manifest_path) + assert manifest.intensity_checks == IntensityCheckSettings() + # New limits change the warnings and the manifest, and read no image again. + reads.clear() + loose = IntensityCheckSettings(outlier_limit=100) + again = run_experiment(root, ZYX, intensity_settings=loose) + assert reads == [] and again.intensity.issues == () + assert read_manifest(again.manifest_path).intensity_checks == loose + + +def test_a_failed_run_has_no_intensity_checks(experiment): + result = run_experiment(experiment, prefill_mapping("YX")) + assert not result.written and result.intensity is None + + def test_first_run_writes_table_and_manifest(experiment): result = run_experiment(experiment, ZYX) assert result.written and result.ok