Stitching artifact correction: fix inverted background mask + add evaluation/sample tooling (WIP)

.gitignore

@@ -160,3 +160,14 @@ cython_debug/
 .idea/
 
 tests/data/000051/
+
+# stitching-artifact-correction: never commit extracted samples or eval outputs (restricted data)
+/stitch-eval/
+**/*.ome.zarr/
+scripts/**/_out/
+*.before.json
+*.after.json
+figures/
+
+# uv lock (upstream uses poetry)
+uv.lock

docs/stitching_artifact_correction_notes.md

@@ -0,0 +1,126 @@
+# Stitching artifact correction — investigation notes
+
+Feature branch: `fix/stitching-artifact-correction` (off `james/spool_zarr`).
+Evaluation sample: 3 adjacent strips of `sub-MF283/sample-slice036`, full Y=1600,
+Z[64:128], X[12000:16096] (a tissue region — see below), extracted by
+`scripts/sample/extract_sample.py` into scratch (never committed; restricted data).
+
+## Baseline (current `james/spool_zarr` behavior, `--blend`, no stripe maps)
+
+Measured on the representative sample (mid-Z plane), via `scripts/eval/evaluate_stitch.py`:
+
+| metric | value | meaning |
+|--------|-------|---------|
+| seam discontinuity | 4.03 | seams present |
+| stripe prominence (along Y) | 22.1 | strong periodic stripes along Y |
+| in-strip mean / std | 125 / 64 | real tissue contrast |
+
+Per-strip body means: **127.8 / 204.2 / 265.3** — strips get progressively brighter.
+
+## Findings
+
+1. **Region matters.** A naive crop at X[0:4096] is *background* (mean ~102, std ~0.9) and
+   shows almost no artifact. Tissue with contrast is around X≈12000–16000, Z≈96. The
+   extraction script now targets that window (`--x-start 12000 --z-start 64`).
+
+2. **Stripes are periodic along Y, not X.** The branch's correction is a per-(z,y)-line
+   model (constant along X). The evaluation stripe metric was corrected to measure spectral
+   prominence along Y (axis=0); it then reports 22.1 on the sample (vs 1.07 measuring X).
+
+3. **The seam is a whole-strip brightness offset, not a local edge.** The raised-cosine
+   blend already smooths the *local* transition over the 192-px overlap (local row jumps are
+   ~few units). The visible seam is the large strip-to-strip body difference (128→204→265).
+
+4. **A naive in-blend gain equalization is the wrong fix.** Matching adjacent strips' overlap
+   means and multiplying (sequential) overcorrects via compounding (strip3 → 101, below
+   strip1's 128) and can flatten genuine signal. It *raised* the seam metric (4.03 → 5.06).
+   Reverted. Per-strip illumination normalization belongs in the stripe/white-matter
+   normalization path (`stripes` → `white_matter_intensity / corr`), applied to whole strips
+   against a common target — which requires the `mip → stripes` chain working first.
+
+5. **`mip → stripes` filename mismatch (R5) blocks the stripe path.** `mip.py` writes
+   `{name}.tiff`; `stripes.py` reads `{name}_proc-mip.tiff` (with a `slice0`→`slice` rewrite).
+   Must be reconciled before a stripe baseline/after can be produced.
+
+6. **Inverted background mask (R3) fixed.** `data*(data < threshold)` kept background and
+   zeroed foreground; changed to `data*(data >= threshold)`, gated by `background_threshold`
+   (default unchanged). Latent data-erasing bug (constitution: never fail silently).
+
+## Update — illumination normalization implemented and measured
+
+Changes made:
+- **R5 fixed**: `mip.py` now writes `{name}_proc-mip.tiff`; `stripes.py` reads it directly
+  (removed the `slice0`→`slice` rewrite that corrupted `slice036`→`slice36`). The
+  `mip → stripes → stitch` chain runs end-to-end.
+- **Stripe-map hardening**: empty (no-foreground) `(z,y)` lines no longer become the
+  `9999999` sentinel (→ black holes); they are filled with the tile's finite-line median
+  and the count is **logged**. On the sample, chunk-0001 had **49,920/102,400 (49%)** empty
+  lines — these would have been black stripes. Added a finiteness guard on apply.
+- **Map smoothing along Y** (`smooth_y`, Gaussian): the per-line map is noisy (many fallback
+  lines); dividing by it injects new stripes. Smoothing keeps the low-frequency illumination
+  falloff but not the line noise.
+- **Seam metric corrected**: now measures the strip-body brightness step across each seam
+  (the blend turns the junction into a gradient, so a local row-jump metric missed it).
+
+Results on the sample (`white_matter_intensity=200`, `smooth_y=25`), vs. the `--blend`-only
+baseline:
+
+| metric | baseline | after | change |
+|--------|----------|-------|--------|
+| seam-body step (% of mean) | 126.6% | 13.6% | **−89%** (≈ SC-002's 90% target) |
+| stripe prominence (Y) | 22.1 | 23.9 | ≈ neutral (−8%) |
+
+Interpretation:
+- **Seams (P1) are essentially fixed** by per-line illumination normalization (the
+  `stripes`/`white_matter_intensity` path), once the map is hardened and Y-smoothed.
+- **Stripes (P2) are not yet reduced.** The per-line flat-field removes the low-frequency
+  falloff that drives the seam, but not the mid-frequency stripe prominence. Note also the
+  overlap analysis: adjacent strips' shared-overlap foreground medians disagree ~2.5× (153 vs
+  379), i.e. strong **within-strip illumination falloff along Y** — a single per-strip gain
+  cannot fix it (confirmed: whole-strip normalization only moved the step 126.6%→103%).
+
+## Stripe suppression: why no FFT-notch, and what's actually left
+
+Investigated stripe suppression with spectral analysis + visual inspection of the stitched
+mosaic (figures generated in scratch; not committed — restricted data):
+
+- **No narrow stripe peak.** The Y power spectrum of the baseline has its in-band energy
+  spread broadly around period ~35–48 px (peak/band-sum ≈ 0.006; a true periodic stripe
+  would be ≈ 1), and low-frequency power dominates the band by ~37×. So an **FFT-notch is the
+  wrong tool** — there is no narrow frequency to notch, and notching the broad band would
+  destroy real ~40 px tissue structure. Deliberately NOT implemented.
+- **The residual "stripe" is banding introduced by the correction**, not a source stripe.
+  Visually, the per-line `white_matter / corr` multiply equalizes the strips (seam fixed) but
+  amplifies low-signal/background lines unevenly into horizontal bands. Y-smoothing
+  (`smooth_y`) and bounding the gain (`correction_clip`, clips to [clip, 100-clip] percentiles)
+  reduce the worst of it but do not eliminate it (both ~neutral on the prominence metric).
+- The stripe-prominence metric is an **unreliable proxy** here: corrected images have more
+  contrast, which inflates it regardless of true striping.
+
+**Foreground-gated correction implemented** (`foreground_gated=True`,
+`background_level` optional, else per-tile 5th percentile): apply
+`corrected = bg + (data - bg) * correction` so background stays flat and only signal above
+background is normalized. Results vs. baseline on the sample (`smooth_y=25`, `wm=200`):
+
+| variant | seam-body step | stripe prominence | mean | banding (visual) |
+|---------|----------------|-------------------|------|------------------|
+| baseline (blend only) | 126.6% | 22.1 | 125 | n/a |
+| blanket multiply | 13.6% (−89%) | 23.9 (worse) | 511 (5× amp) | strong horizontal bands |
+| **foreground-gated** | 42.4% (−67%) | 20.8 (−6%) | 117 (preserved) | **bands removed** |
+
+Visually the foreground-gated mosaic is the cleanest: strips equalized, tissue preserved, and
+the horizontal background banding of the blanket multiply is gone. Its seam-body-step number is
+higher only because that metric conflates each strip's *foreground fraction* (strips genuinely
+tile different tissue density) with intensity — a metric limitation, not a worse result.
+**Recommendation: foreground-gated is the right default for the illumination/stripe correction.**
+
+## Remaining work
+
+- Foreground-gated correction (above) to remove background banding; re-define an acceptance
+  metric for stripes that matches a visible artifact (current prominence metric is unreliable).
+- Tune `white_matter_intensity` per dataset (default 1000 amplifies this data ~5×; relative
+  metrics unaffected, but output scale matters downstream).
+- Full-slice end-to-end run (SC-005); `_cm2` camera variant; unit tests for the fixes.
+- **Maintainer input**: is the inter-strip / within-strip brightness variation purely
+  illumination (normalize away) or partly real signal? And is there a dataset/region with a
+  genuine periodic stripe (this sample has none)? These gate the stripe approach.

linc_convert/modalities/lsm/convert_spool_or_zarr.py

@@ -232,7 +232,7 @@
 
         if self.flip_z:
             Z_eff = self.raw.shape[0] - 1 - Zg
             # TODO: figure out why this works. This was generated from guess and check and not what my intuition would tell me is correct
             X_src = Xg + Z_eff * self.shps
         else:
             Z_eff = Zg
@@ -352,6 +352,9 @@
     stripes: Optional[str] = None,
     white_matter_intensity: float = 1000.0,
     background_threshold: Optional[Union[float, Literal["auto"]]] = None,
+    correction_clip: float = 0.0,
+    foreground_gated: bool = False,
+    background_level: Optional[float] = None,
     checkpoint_file: Optional[str] = None,
     alternate_pattern: bool = False,
     flip_z: bool = False
@@ -673,24 +676,63 @@
                             constant_values=0,
                         )
                     if threshold is not None:
-                        data = data*(data < threshold)
+                        # Suppress BACKGROUND (values at/below the background level),
+                        # keeping foreground signal. The previous `data < threshold`
+                        # zeroed foreground and kept background — an inverted mask that
+                        # silently erased signal (constitution IV). Gated by
+                        # background_threshold, so default (None) behavior is unchanged.
+                        data = data*(data >= threshold)
 
                     if stripes is not None:
 
                         name = os.path.basename(
                             tile.filename.rstrip("/").replace(".ome.zarr", ""))
 
                         correction = tiff.imread(
-                            f"{stripes}/{name}.tiff")
+                            f"{stripes}/{name}.tiff").astype(np.float32)
                         correction[correction == 0.0] = 1.0
                         if flip_z:
                             correction = correction[::-1, :]
 
                         correction = white_matter_intensity / correction
 
+                        # Never let a bad map silently zero/NaN out data
+                        # (constitution: never fail silently).
+                        n_bad = int((~np.isfinite(correction)).sum())
+                        if n_bad:
+                            logger.warning(
+                                "tile %s: %d non-finite stripe-correction "
+                                "values; setting them to 1.0 (no-op)",
+                                name, n_bad,
+                            )
+                            correction = np.where(
+                                np.isfinite(correction), correction, 1.0)
+
+                        # Bound the per-line gain so low-signal/background lines are
+                        # not amplified into horizontal bands (visible artifact from an
+                        # unbounded wm/corr). Clips to the [clip, 100-clip] percentiles
+                        # of this tile's correction; default 0.0 = off.
+                        if correction_clip and correction_clip > 0:
+                            lo = float(np.percentile(correction, correction_clip))
+                            hi = float(np.percentile(correction, 100 - correction_clip))
+                            correction = np.clip(correction, lo, hi)
+
                         correction = correction[:, :, None]
 
-                        data = data * correction
+                        if foreground_gated:
+                            # Foreground-gated (affine) correction: normalize signal
+                            # above background while leaving background unscaled, so the
+                            # per-line gain cannot amplify background into horizontal
+                            # bands. corrected = bg + (data - bg) * correction. bg is a
+                            # per-tile background level (param, else 5th percentile).
+                            if background_level is not None:
+                                bg = float(background_level)
+                            else:
+                                bg = da.percentile(
+                                    data.reshape((-1,)).astype("f4"), [5])[0]
+                            data = bg + (data - bg) * correction
+                        else:
+                            data = data * correction
 
                     next_overlap = 0
                     if (y+1, z) in tiles:

linc_convert/modalities/lsm/mip.py

@@ -33,7 +33,7 @@
 
 @mip.default
 @autoconfig
 def convert(
     inp: str,
     *,
     general_config: GeneralConfig = None,
@@ -59,5 +59,7 @@
         if z_end is not None:
             reader = reader[:z_end, :, :]
 
-        tiff.imwrite(f"{general_config.out}/{name}.tiff",
+        # Write the name `lsm stripes` expects ({name}_proc-mip.tiff) so the
+        # mip -> stripes -> stitch chain works without manual renaming (R5).
+        tiff.imwrite(f"{general_config.out}/{name}_proc-mip.tiff",
                      da.max(reader, axis=0).compute())

linc_convert/modalities/lsm/stitch.py

@@ -55,6 +55,9 @@ def convert(
     stripes: Optional[str] = None,
     white_matter_intensity: float = 1000.0,
     background_threshold: Optional[Union[float, Literal["auto"]]] = None,
+    correction_clip: float = 0.0,
+    foreground_gated: bool = False,
+    background_level: Optional[float] = None,
     checkpoint_file: Optional[str] = None,
     alternate_pattern: bool = False,
     flip_z: bool = False
@@ -147,6 +150,9 @@ def convert(
                           stripes=stripes,
                           white_matter_intensity=white_matter_intensity,
                           background_threshold=background_threshold,
+                          correction_clip=correction_clip,
+                          foreground_gated=foreground_gated,
+                          background_level=background_level,
                           checkpoint_file=checkpoint_file,
                           alternate_pattern=alternate_pattern,
                           flip_z=flip_z

linc_convert/modalities/lsm/stripes.py

@@ -3,6 +3,7 @@
 import logging
 import os
 import time
+import warnings
 from typing import Optional
 
 # externals
@@ -53,6 +54,7 @@ def create(
     z_end: Optional[int] = None,
     y_start: Optional[int] = None,
     y_end: Optional[int] = None,
+    smooth_y: float = 0.0,
 ) -> None:
     """
     Generate ZY projection TIFF images from volumetric tile data.
@@ -133,8 +135,10 @@ def create(
         output_name = f"{general_config.out}/{name}.tiff"
 
         if not os.path.exists(output_name):
-            yx_path = os.path.join(
-                mip_dir, f"{name}_proc-mip.tiff").replace("slice0", "slice")
+            # Match exactly what `lsm mip` writes. The previous
+            # .replace("slice0", "slice") corrupted names like
+            # "sample-slice036" -> "sample-slice36" (R5).
+            yx_path = os.path.join(mip_dir, f"{name}_proc-mip.tiff")
 
             if not os.path.exists(yx_path):
                 raise FileNotFoundError(f"Missing YX image: {yx_path}")
@@ -148,13 +152,49 @@ def create(
 
             mask = _compute_mask(img_yx)
 
-            vol_np = reader.compute()
-            vol_np = vol_np.astype(float)
+            vol_np = np.asarray(reader[:], dtype=float)
             vol_np[:, ~mask] = np.nan
-            corr_zy = np.nanmedian(vol_np, axis=2)
-            corr_zy = np.nan_to_num(corr_zy, nan=9999999.0)
-
-            tiff.imwrite(output_name + ".tmp", corr_zy)
+            with warnings.catch_warnings():
+                # all-foreground-masked (z,y) lines produce an all-NaN slice
+                warnings.simplefilter("ignore", RuntimeWarning)
+                corr_zy = np.nanmedian(vol_np, axis=2)
+
+            # Harden the map: (z,y) lines with no foreground come back NaN. The old
+            # code replaced them with a 9999999 sentinel, which at stitch time becomes
+            # white_matter_intensity / 9999999 ~= 0 -> a black line/hole. Instead fill
+            # with the tile's finite-line median (a neutral correction) and LOG how many
+            # lines fell back, so degraded output is visible, never silent.
+            finite = np.isfinite(corr_zy)
+            n_empty = int((~finite).sum())
+            if n_empty:
+                logger.warning(
+                    "tile %s: %d/%d (z,y) lines had no foreground; "
+                    "filling with finite-line median fallback",
+                    name, n_empty, corr_zy.size,
+                )
+            fill = float(np.median(corr_zy[finite])) if finite.any() else 1.0
+            corr_zy = np.where(finite, corr_zy, fill)
+            # avoid zero/negative medians that would blow up white_matter / corr
+            positive = corr_zy > 0
+            if not positive.all():
+                pos_fill = (
+                    float(np.median(corr_zy[positive])) if positive.any() else 1.0
+                )
+                corr_zy = np.where(positive, corr_zy, pos_fill)
+            if not np.isfinite(corr_zy).all():
+                raise ValueError(f"non-finite correction map for tile {name}")
+
+            # Smooth the per-(z,y)-line map along Y so it captures the smooth
+            # illumination falloff but NOT line-to-line noise. Dividing data by a
+            # noisy per-line map injects new stripes; smoothing prevents that while
+            # still removing the low-frequency falloff that drives seams.
+            if smooth_y and smooth_y > 0:
+                from scipy.ndimage import gaussian_filter1d
+
+                corr_zy = gaussian_filter1d(
+                    corr_zy, sigma=float(smooth_y), axis=1, mode="nearest")
+
+            tiff.imwrite(output_name + ".tmp", corr_zy.astype(np.float32))
             os.replace(output_name + ".tmp", output_name)
 
         print("--- %s secs ---" % (time.time() - start_time))

scripts/eval/compare.py

@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+"""Compare two evaluate_stitch.py reports (before vs after) and check SC thresholds.
+
+Feature: stitching-artifact-correction (SC-002 seam >=90% reduction, SC-003 stripe
+>=75% reduction, SC-004 in-strip content preserved). Exits non-zero if a gate fails,
+so the dev loop fails loudly (constitution IV).
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import sys
+
+
+def pct_reduction(before: float, after: float) -> float:
+    if before <= 0:
+        return 0.0
+    return 100.0 * (before - after) / before
+
+
+def main(argv: list[str] | None = None) -> int:
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("before_json")
+    ap.add_argument("after_json")
+    ap.add_argument("--seam-target", type=float, default=90.0, help="%% reduction (SC-002)")
+    ap.add_argument("--stripe-target", type=float, default=75.0, help="%% reduction (SC-003)")
+    ap.add_argument("--instrip-tol", type=float, default=0.05, help="max rel. mean drift (SC-004)")
+    args = ap.parse_args(argv)
+
+    with open(args.before_json) as f:
+        b = json.load(f)
+    with open(args.after_json) as f:
+        a = json.load(f)
+
+    seam_red = pct_reduction(b["seam_discontinuity"], a["seam_discontinuity"])
+    stripe_red = pct_reduction(b["stripe_energy"], a["stripe_energy"])
+    mean_drift = abs(a["in_strip_mean"] - b["in_strip_mean"]) / (abs(b["in_strip_mean"]) + 1e-9)
+
+    print(f"seam discontinuity : {b['seam_discontinuity']:.3f} -> {a['seam_discontinuity']:.3f} "
+          f"({seam_red:.1f}% reduction; target >={args.seam_target}%) [SC-002]")
+    print(f"stripe energy      : {b['stripe_energy']:.4f} -> {a['stripe_energy']:.4f} "
+          f"({stripe_red:.1f}% reduction; target >={args.stripe_target}%) [SC-003]")
+    print(f"in-strip mean drift: {mean_drift*100:.2f}% (tol {args.instrip_tol*100:.0f}%) [SC-004]")
+
+    gates = {
+        "SC-002 seam": seam_red >= args.seam_target,
+        "SC-003 stripe": stripe_red >= args.stripe_target,
+        "SC-004 in-strip": mean_drift <= args.instrip_tol,
+    }
+    failed = [k for k, ok in gates.items() if not ok]
+    if failed:
+        print(f"FAIL: {', '.join(failed)}")
+        return 1
+    print("PASS: all gates met")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())