diff --git a/BENCHMARKING.md b/BENCHMARKING.md index 4889866..1e75559 100644 --- a/BENCHMARKING.md +++ b/BENCHMARKING.md @@ -22,6 +22,37 @@ onto the T1c grid with SMD+linear. After the repair pass: volumes, 3893 subjects; unlabeled pool 13280 - tumor volume stats: median 3116 mm3, p95 24.8k, max 3.2M, zero_frac 0.12% +## Dataset exclusions (2026-10-04) + +Row-level curation lives in `src/exclusions.py` and is applied at the **first +stage** by `scripts/make_tv_slices.py` (right after the rows are fetched from +the sources, before any split): it classifies every labeled volume and writes +its TV-slice QA figure to `data/tv_slices/inclusion/` (kept) or +`data/tv_slices/exclusion/` (excluded, reason in the title). `05_build_splits` +re-applies the same rules **after** split assignment (kept rows retain their +split; the subject-level shuffle boundaries are unchanged): + +- total TV > 1e6 mm3: **1** row (ntuh_6019528, label outlines the full skull) +- single target (1 connected component) with TV < 100 mm3: **109** rows +- non-brain images: **33** rows (all ntuh — spine/neck/torso CyberKnife + cases whose series names are indistinguishable from brain T1c; identified + by visual review of the TV-slice figures, now in `data/tv_slices/exclusion/`) + +Audit trail: `data/manifests/excluded.jsonl` (143 rows, reason + split); +pre-exclusion manifests in `data/manifests/backup_20261004/`. New dataset: +**4057** labeled rows → train 3447 / val 205 / test 405 (3771 subjects); +unlabeled pool unchanged (13280); volume stats now median 3212 mm3, p95 +23.4k, max 163k (m6 meningioma), zero_frac still 0.12% (5 rows, e.g. +ntuh_3436472 ×3). + +After the initial 39-row non-brain list, expert inspection overturned 6 +(skull-base / posterior-fossa brain targets that look like cervical spine at +the TV-max slice: ntuh_5902230, ntuh_5231609, ntuh_5083686, ntuh_5157194, +ntuh_5303758, ntuh_6373987). The remaining 33 were re-audited at the +mid-cranial slice (widest head-mask axial) — all still non-brain. Review +protocol now recorded in `src/exclusions.py`: classify at full resolution +on the mid-cranial slice, not just the TV-max slice. + ## Smoke run (2026-10-04) Small/fast end-to-end check that all three pipelines run on the rebuilt data @@ -52,6 +83,36 @@ ntuh). Artifacts: `results/benchmark_pipelines_smoke.json`, ModuleNotFoundError) — installed `pip install monai==1.6.1 --no-deps` (torch 2.14+cu126 untouched); `import monai` verified. +## Full-dataset 10-epoch run (2026-10-04) — interim + +`--train-epochs 10`, default budgets (2856 train rows = 80% of split_train, +42 eval rows = 10% of split_test, TTA on), A/C/GPU0·GPU2 finished, B still +preprocessing at time of writing: + +``` +python scripts/benchmark_pipelines.py --pipelines A --gpu 0 --train-epochs 10 \ + --out results/benchmark_full_A.json # done +python scripts/benchmark_pipelines.py --pipelines C --gpu 2 --train-epochs 10 \ + --out results/benchmark_full_C.json # done +python scripts/benchmark_pipelines.py --pipelines B --gpu 1 --train-epochs 10 \ + --out results/benchmark_full_B.json # in progress (nnU-Net prep) +``` + +| pipeline | train s/epoch | wall | tr VRAM | s/vol (TTA) | in VRAM | params | weights | dice @10ep | +|---|---|---|---|---|---|---|---|---| +| A Unet3D(16, depth 4) | 258.8 (952 steps) | 46 min | 2.0 GB | 3.51 ± 0.93 | 0.5 GB | 1.46 M | 17.6 MB | **0.0008** | +| C MONAI UNet(16→128) | 284.3 (952 steps) | 48 min | 0.4 GB | 0.38 ± 0.06 | 0.5 GB | 1.19 M | 14.3 MB | **0.2119** | +| B nnU-Net 3d_fullres | (pending) | (pending) | — | — | — | 30.8 M | 246 MB | (pending) | + +Dice distribution over 41 timed test vols: A p25/p75 = 0.0000/0.0000 (98% ≈ 0); +C p25/p75 = 0.0000/0.4048 (29% ≈ 0). Training loss (epoch 1 → 10): A +0.335 → 0.209, C 0.403 → 0.188. 10 epochs is far below convergence for these +sparse targets (median tumor 3.1 k mm³); the A/C numbers here are +intermediate, and C is already clearly ahead at this budget. Checkpoints at +`runs_bench/{A,C}/best.pt` (extend with `--resume --train-epochs `). +Note: B's first `plan_and_preprocess` attempt died and the built-in 60 s retry +re-ran it successfully. + ## Resuming training (`--resume`) All three pipelines can continue from their latest benchmark checkpoint: @@ -83,8 +144,9 @@ python scripts/benchmark_pipelines.py --pipelines "A B C" --resume \ ## Readiness notes for a full run - Budget knobs: `--train-epochs` (default 100), `--max-rows` (default 80% of - split_train ≈ 2856 for A/C and for B's `--nnu-cases`), `--max-eval-rows` - (default 10% of split_test ≈ 42), `--no-tta`. + split_train ≈ 2753 for A/C and for B's `--nnu-cases`, post-exclusion; + ≈ 2856 pre-exclusion), `--max-eval-rows` (default 10% of split_test ≈ 40 + post-exclusion), `--no-tta`. - `--skip-train` is only usable for **A** (`runs/round0/best.pt`, trained 2026-09-25 on the *previous* split — a few patients moved between train/test in the rebuild, so its test dice carries a small leakage diff --git a/README.md b/README.md index 279f5b6..cc91ff6 100644 --- a/README.md +++ b/README.md @@ -40,6 +40,8 @@ data/ proc/ .nii.gz (1mm, cropped, normalized) proc/_label.nii.gz procmeta/.json geometry (origin/direction/crop_vox) + normalization + tv_slices/inclusion/ TV-slice QA figures, kept rows (_.png) + tv_slices/exclusion/ TV-slice QA figures, excluded rows (reason in title) vols.json labeled tumor volume stats (mm3 percentiles) pseudo/roundK/ in-house pseudo-labels (rows.jsonl, masks, summary.json) pseudo_nnu/roundK/ nnU-Net pseudo-labels (rows.jsonl, masks, summary.json) @@ -69,12 +71,13 @@ DDP, per-sample weighted loss. | 02 | `02_build_m6_dataset.py` | Build M6-2025 manifests (GTV CT→T1c registration/warp) | | 03 | `03_scan_lee_t1c.py` | Scan lee for brain T1c series → raw + selected manifests | | 04 | `04_reconstruct_lee.py` | Reconstruct lee T1c niftis from JPG slices + txt metadata | -| 05 | `05_build_splits.py` | Patient-level train/val/test splits + unlabeled pool + volume stats | +| 05 | `05_build_splits.py` | Patient-level train/val/test splits + unlabeled pool + volume stats (re-applies `src/exclusions` post-split) | | 06 | `06_pseudo_label.py` | Pseudo-label the pool with a round model (sliding-window + TTA) | | 07 | `07_train.py` | DDP training entrypoint (labeled + weighted pseudo rows) | | 08 | `08_eval.py` | Held-out test evaluation (Dice from probability map) | | 09 | `09_run_iterative.py` | Orchestrates rounds 0…K and produces the summary table/plot | | — | `preprocess.py` | Crop/normalize source volumes → `data/proc/` | +| — | `make_tv_slices.py` | First-stage curation (right after source fetch): classifies labeled rows via `src/exclusions.py` + TV-slice QA figures → `data/tv_slices/{inclusion,exclusion}/` | | — | `scan_procs.py` | Flag corrupt/non-finite processed volumes | | — | `test_dataloader.py` | Smoke-test the training dataloader | diff --git a/scripts/05_build_splits.py b/scripts/05_build_splits.py index 4f4dce5..b0f5e13 100644 --- a/scripts/05_build_splits.py +++ b/scripts/05_build_splits.py @@ -11,6 +11,11 @@ normalized, used by Pipelines A and C) and native fields (img/label, the raw source NIfTIs — used by Pipeline B, which feeds nnU-Net its own plan_and_preprocess). Row membership still requires the processed volume to exist, so the patient-level splits stay identical across pipelines. + +Before splitting, every labeled row is screened with the row-level exclusion +rules (src/exclusions.py: total TV > 1e6 mm3, single target < 100 mm3, and +the non-brain ntuh rows in NON_BRAIN_KEYS). Excluded rows are written to +data/manifests/excluded.jsonl with their reason and never enter a split. """ import os import sys @@ -19,20 +24,16 @@ import json import argparse import random import numpy as np -from src.common import ROOT, DATA, d, load_jsonl, save_jsonl, read_nii_arr +from src.common import ROOT, DATA, d, load_jsonl, save_jsonl, proc_row +from src.exclusions import label_stats, exclusion_reason -def proc_row(r, prefix): - key = r["key"] - p = os.path.join(d("data/proc"), key + ".nii.gz") - if not os.path.exists(p): - return None - row = {"key": key, "subject": f"{prefix}_{r['subject']}", "pimg": p, - "img": r["img"], "label": r.get("label"), "date": r.get("date")} - lab = os.path.join(d("data/proc"), key + "_label.nii.gz") - if r.get("label") and os.path.exists(lab): - row["plabel"] = lab - return row +def _label_stats_row(r): + """(key, label stats) for the ProcessPoolExecutor in main().""" + try: + return r["key"], label_stats(r["plabel"]) + except Exception: + return r["key"], {"total": 0, "n_cc": 0} def main(): @@ -53,6 +54,31 @@ def main(): labeled_rows.append(p) print(f"labeled processed rows: {len(labeled_rows)}") + # row-level curation (src/exclusions.py): label stats + exclusion reason + # per row. Exclusions are applied AFTER the split assignment, so kept rows + # retain the split of the 2026-10-04 rebuild and excluded.jsonl is the + # audit trail (excluded subjects still count in the shuffle boundaries). + from collections import Counter + from concurrent.futures import ProcessPoolExecutor + + tv_vols, reason_of = {}, {} + with ProcessPoolExecutor(max_workers=32) as ex: + for k, stats in ex.map(_label_stats_row, labeled_rows, chunksize=16): + tv_vols[k] = stats["total"] + reason_of[k] = exclusion_reason(k, stats) + n_excl = sum(1 for v in reason_of.values() if v) + + def _cat(reason): + return (reason.split(" image", 1)[0] + " image" + if "non-brain" in reason else + reason.split(" TV ", 1)[0] + " TV < 100 mm3" + if "single target" in reason else + "total TV > 1e6 mm3") + rc = Counter(_cat(v) for v in reason_of.values() if v) + print(f"excluded rows: {n_excl}") + for reason, n in sorted(rc.items()): + print(f" excluded {n:4d} {reason}") + subjects = sorted({r["subject"] for r in labeled_rows}) random.shuffle(subjects) n_test = max(1, int(len(subjects) * args.frac_test)) @@ -62,18 +88,24 @@ def main(): train_subj = set(subjects[n_test + n_val:]) print(f"subjects: {len(subjects)} train={len(train_subj)} val={len(val_subj)} test={len(test_subj)}") - splits = {"train": [], "val": [], "test": []} + splits, excluded = {"train": [], "val": [], "test": []}, [] for r in labeled_rows: r["w"] = 1.0 if r["subject"] in test_subj: - splits["test"].append(r) + sk = "test" elif r["subject"] in val_subj: - splits["val"].append(r) + sk = "val" else: - splits["train"].append(r) + sk = "train" + if reason_of[r["key"]]: + excluded.append({**r, "split": sk, "exclude_reason": reason_of[r["key"]]}) + else: + splits[sk].append(r) for k, v in splits.items(): save_jsonl(v, os.path.join(DATA, "manifests", f"split_{k}.jsonl")) print(f" split_{k}: {len(v)} volumes / {len(set(r['subject'] for r in v))} subjects") + save_jsonl(excluded, os.path.join(DATA, "manifests", "excluded.jsonl")) + kept_rows = [r for v in splits.values() for r in v] labeled_subjects = train_subj | val_subj | test_subj @@ -114,15 +146,12 @@ def main(): print(f"unlabeled pool: {len(pool)} volumes ({sum(1 for r in pool if r['source']=='m6')} m6, " f"{sum(1 for r in pool if r['source']=='lee')} lee) from {len(set(r['subject'] for r in pool))} subjects") - # tumor volume stats (1mm voxels == mm3) + # tumor volume stats of the kept dataset rows (1mm voxels == mm3; totals from +# the curation pass) vols = [] - for r in labeled_rows: - try: - l = read_nii_arr(r["plabel"]) - n = int((l > 0.5).sum()) - vols.append({"key": r["key"], "subject": r["subject"], "vol_mm3": n}) - except Exception: - pass + for r in kept_rows: + vols.append({"key": r["key"], "subject": r["subject"], + "vol_mm3": tv_vols[r["key"]]}) v = np.array([x["vol_mm3"] for x in vols]) stats = {"n": int(v.size), "p2": float(np.percentile(v, 2)), "p5": float(np.percentile(v, 5)), diff --git a/scripts/make_tv_slices.py b/scripts/make_tv_slices.py new file mode 100644 index 0000000..6499cd9 --- /dev/null +++ b/scripts/make_tv_slices.py @@ -0,0 +1,230 @@ +"""First-stage curation + TV-slice QA figures for labeled rows. + +Runs on the labeled source rows (data/manifests/{ntuh,m6_labeled}.jsonl, 1mm +processed) — the first place rows exist after fetching from the sources, +before any split (05). For each labeled volume: + + 1. classify with the exclusion rules in src/exclusions.py (total TV > 1e6 + mm3, single target < 100 mm3, curated non-brain list), + 2. render one PNG with the three orthogonal slices (axial | coronal | + sagittal) that carry the largest TV (target-volume) area, with the TV + contour overlaid; when the label has several connected targets + (26-connectivity), only the largest one is shown, + 3. write the figure to data/tv_slices/inclusion/ (kept) or + data/tv_slices/exclusion/ (excluded, reason in the title). + +Slices are picked per body axis (S/P/L) from the procmeta direction, so the +panel order follows radiological conventions (same mapping as +scripts/make_screenshots.py); volumes that are not body-aligned fall back to +the image-axis slices with a note in the title. Filenames are +_.png, so each folder lists small TV -> large TV. + +The split files (05) are the dataset of record; this stage only curates and +documents. 05 applies the same src/exclusions rules after split assignment, +so the two always agree on kept/excluded. + +Usage: + python scripts/make_tv_slices.py # all source rows, 32 workers + python scripts/make_tv_slices.py --from-splits # rows/tag from split files + python scripts/make_tv_slices.py --keys k1 k2 # specific rows + python scripts/make_tv_slices.py --limit 50 # first N rows + python scripts/make_tv_slices.py --force # overwrite existing PNGs +""" +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import glob +import json +import argparse +from collections import Counter + +import numpy as np +import SimpleITK as sitk +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from scipy import ndimage + +from src.common import DATA, d, load_jsonl, proc_row, labeled_source_rows +from src.exclusions import label_stats, exclusion_reason +from scripts.make_screenshots import vol_to_body + + +def out_name(vol_total, key): + """Flat-folder filename: zero-padded TOTAL TV volume (mm3, all targets) + first so a plain listing sorts small -> large, then the row key.""" + return f"{int(vol_total):07d}_{key}.png" + + +def largest_cc(mask, structure=None): + """Largest connected component of a boolean mask (26-conn by default).""" + if structure is None: + structure = np.ones((3, 3, 3), dtype=int) + lab, n = ndimage.label(mask, structure=structure) + if n == 0: + return mask, 0 + sizes = ndimage.sum(mask, lab, range(1, n + 1)) + return lab == (int(np.argmax(sizes)) + 1), int(n) + + +def max_area_slices(cc): + """Per-axis index of the slice containing the largest CC area.""" + return tuple(int(np.argmax(cc.sum(axis=tuple(j for j in range(cc.ndim) if j != i)))) + for i in range(cc.ndim)) + + +def render(key, tag, mr, cc, n_cc, vol_total, direction, reason, out, dpi=110): + """1x3 figure: axial | coronal | sagittal, max-CC-area slice each, with the + largest-CC contour in red. Writes to `out`. Returns the path.""" + body, ok, _ = vol_to_body(mr, direction) + if ok: + c, _, _ = vol_to_body(cc.astype(np.float32), direction) + si, pi, li = max_area_slices(c) + panels = [ + (body[si], c[si], "axial", "S", si), + (body[:, pi, :][::-1], c[:, pi, :][::-1], "coronal", "P", pi), + (body[:, :, li][::-1], c[:, :, li][::-1], "sagittal", "L", li), + ] + else: # not body-aligned: image-axis fallback + zi, yi, xi = max_area_slices(cc) + panels = [ + (mr[zi], cc[zi], "axial*", "z", zi), + (mr[:, yi, :][::-1], cc[:, yi, :][::-1], "coronal*", "y", yi), + (mr[:, :, xi][::-1], cc[:, :, xi][::-1], "sagittal*", "x", xi), + ] + + lo, hi = np.percentile(mr[mr > 0], [1, 99.5]) + panel_specs = [] + for im, c2, title, axn, idx in panels: + nr, nc = im.shape + a = 1.0 # processed volumes are 1mm isotropic + panel_specs.append((im, c2, f"{title} {axn}={idx}", a, nr, nc)) + + max_in, gap_in, m_in, top_in, bot_in = 4.4, 0.25, 0.3, 1.7, 0.2 + w_in, h_in = [], [] + for _, _, _, a, nrow, ncol in panel_specs: + pw, ph = ncol, nrow * a + m = max(pw, ph) + w_in.append(max_in * pw / m) + h_in.append(max_in * ph / m) + plot_h = max_in + fig_w = sum(w_in) + 2 * gap_in + 2 * m_in + fig_h = plot_h + top_in + bot_in + fig = plt.figure(figsize=(fig_w, fig_h)) + left = m_in / fig_w + for (im, c2, title, a, _, _), w, h in zip(panel_specs, w_in, h_in): + axi = fig.add_axes([left, (bot_in + (plot_h - h) / 2) / fig_h, w / fig_w, h / fig_h]) + axi.imshow(im, cmap="gray", vmin=lo, vmax=hi, origin="upper", aspect=a) + if c2.any(): + axi.contour(c2, levels=[0.5], colors="r", linewidths=1.2) + axi.set_title(title, color="w", fontsize=11) + axi.set_xticks([]) + axi.set_yticks([]) + for spine in axi.spines.values(): + spine.set_edgecolor("0.35") + left += (w + gap_in) / fig_w + tkey = key.split("_", 1)[1] if "_" in key else key + align = "body-aligned" if ok else "NOT body-aligned (image axes)" + status = f"EXCLUDED: {reason}" if reason else "kept" + fig.suptitle( + f"{tkey} [{tag}] TV total={int(vol_total)} mm3, " + f"largest={int(cc.sum())} mm3 ({n_cc} targets) {status} {align}", + color="w", fontsize=12) + fig.patch.set_facecolor("k") + fig.savefig(out, facecolor="k", dpi=dpi) + plt.close(fig) + return out + + +def worker(args): + key, tag, pimg, plabel, out_base, force, dpi = args + try: + meta = json.load(open(os.path.join(d("data/procmeta"), key + ".json"))) + stats = label_stats(plabel) + except Exception: + return key, "error: unreadable label/procmeta" + reason = exclusion_reason(key, stats) + cat = reason.split(" (", 1)[0] if reason else None + # reconcile a stale figure for this key left in the other folder by an + # earlier run under a different exclusion status (e.g. a key newly added to + # NON_BRAIN_KEYS would otherwise keep its old inclusion/ figure) + folder = "exclusion" if reason else "inclusion" + other = "inclusion" if reason else "exclusion" + for stale in glob.glob(os.path.join(out_base, other, f"*_{key}.png")): + os.remove(stale) + if stats["total"] == 0: + return key, "empty-label", cat + status = "excluded" if reason else "included" + out = os.path.join(out_base, folder, out_name(stats["total"], key)) + if os.path.exists(out) and not force: + return key, f"skip-{status}", cat + try: + mr = np.asarray(sitk.GetArrayFromImage(sitk.ReadImage(pimg)), dtype=np.float32) + lab = sitk.GetArrayFromImage(sitk.ReadImage(plabel)) + cc, n_cc = largest_cc(lab > 0.5) + try: + render(key, tag, mr, cc, n_cc, int(stats["total"]), + np.array(meta["direction"]), reason, out, dpi) + finally: + plt.close("all") + return key, status, cat + except Exception as e: # noqa + return key, f"error: {e!r}", cat + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--from-splits", action="store_true", + help="rows/tags from split_{train,val,test}.jsonl instead of " + "the source manifests (ntuh + m6_labeled)") + ap.add_argument("--keys", nargs="*", default=[]) + ap.add_argument("--limit", type=int, default=0) + ap.add_argument("--workers", type=int, default=32) + ap.add_argument("--force", action="store_true") + ap.add_argument("--out", default=d("data/tv_slices")) + ap.add_argument("--dpi", type=int, default=110) + args = ap.parse_args() + os.makedirs(os.path.join(args.out, "inclusion"), exist_ok=True) + os.makedirs(os.path.join(args.out, "exclusion"), exist_ok=True) + + if args.from_splits: + rows = [] + for s in ("train", "val", "test"): + for r in load_jsonl(os.path.join(DATA, "manifests", f"split_{s}.jsonl")): + rows.append((r["key"], s, r["pimg"], r.get("plabel"))) + else: + rows = [(r["key"], + "ntuh" if r["key"].startswith("ntuh_") else "m6", + r["pimg"], r.get("plabel")) for r in labeled_source_rows()] + if args.keys: + keep = set(args.keys) + rows = [r for r in rows if r[0] in keep] + if args.limit: + rows = rows[: args.limit] + print(f"rows: {len(rows)}") + + from concurrent.futures import ProcessPoolExecutor, as_completed + tasks = [(k, tag, pimg, plabel, args.out, args.force, args.dpi) + for k, tag, pimg, plabel in rows] + stats = {} + reasons = Counter() + done = 0 + with ProcessPoolExecutor(max_workers=args.workers) as ex: + futs = {ex.submit(worker, t): t for t in tasks} + for fu in as_completed(futs): + k, status, cat = fu.result() + stats[status] = stats.get(status, 0) + 1 + if status == "excluded" and cat: + reasons[cat] += 1 + done += 1 + if done % 100 == 0 or done == len(tasks): + print(f" {done}/{len(tasks)} {stats}", flush=True) + print(f"done: {stats}") + if reasons: + print("excluded by rule:") + for r, n in sorted(reasons.items()): + print(f" {n:4d} {r}") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/common.py b/src/common.py index e128407..5bacac6 100644 --- a/src/common.py +++ b/src/common.py @@ -223,6 +223,45 @@ def write_arr(arr, path, itk_img=None): sitk.WriteImage(itk_img, path, True) +def proc_row(r, prefix): + """Processed-volume view of a source manifest row, or None when the + processed volume does not exist. + + Source rows come from data/manifests/{ntuh,m6_labeled,m6_unlabeled}.jsonl + ({key, subject, img, label?, date?} — bare subjects, no source prefix). + Returns {key, subject (prefixed), pimg, img, label, date[, plabel]} where + pimg/plabel point at data/proc/{,_label}.nii.gz (1mm processed). + """ + key = r["key"] + proc_root = d("data/proc") + p = os.path.join(proc_root, key + ".nii.gz") + if not os.path.exists(p): + return None + row = {"key": key, "subject": f"{prefix}_{r['subject']}", "pimg": p, + "img": r["img"], "label": r.get("label"), "date": r.get("date")} + lab = os.path.join(proc_root, key + "_label.nii.gz") + if r.get("label") and os.path.exists(lab): + row["plabel"] = lab + return row + + +def labeled_source_rows(): + """All labeled rows fetched from the sources (ntuh + m6_labeled), as + processed rows — the first dataset stage, before any split (05).""" + rows = [] + for name, prefix in (("ntuh", "ntuh"), ("m6_labeled", "m6")): + p = os.path.join(DATA, "manifests", name + ".jsonl") + if not os.path.exists(p): + continue + for r in load_jsonl(p): + if not r.get("label"): + continue + pr = proc_row(r, prefix) + if pr: + rows.append(pr) + return rows + + def dice(a, b, eps=1e-7): a = a > 0 b = b > 0 diff --git a/src/exclusions.py b/src/exclusions.py new file mode 100644 index 0000000..b466a37 --- /dev/null +++ b/src/exclusions.py @@ -0,0 +1,117 @@ +"""Row-level exclusion rules for the labeled dataset. + +Used at two places, both via exclusion_reason() (single source of truth): + - scripts/make_tv_slices.py — first-stage curation right after the rows are + fetched from the sources (ntuh + m6_labeled): classifies each labeled + volume and writes its QA figure to data/tv_slices/inclusion/ or + data/tv_slices/exclusion/; + - scripts/05_build_splits.py — applied after the patient-level split + assignment so kept rows retain their split; excluded rows are recorded + in data/manifests/excluded.jsonl. + +Rules (2026-10-04 curation): + 1. total TV volume > MAX_TOTAL_VOL mm3 (implausible targets, e.g. a + full-skull outline instead of a lesion) + 2. a single target (1 connected component, 26-conn) with total volume + < MIN_SINGLE_TARGET_VOL mm3 (sub-target fragments) + 3. non-brain images (spine/neck CyberKnife cases in + ntuh; their series names are indistinguishable from brain T1c, so the + rows were identified by visual review of the TV-slice figures, + scripts/make_tv_slices.py — 43 rows, see NON_BRAIN_KEYS) + +label_stats() / exclusion_reason() are the single source of truth; 05 calls +them for every labeled row and records the excluded rows (with reason) in +data/manifests/excluded.jsonl. +""" +import numpy as np +from scipy import ndimage + +MAX_TOTAL_VOL = 1_000_000 # mm3 (1mm isotropic voxels == mm3) +MIN_SINGLE_TARGET_VOL = 100 # mm3 + +# Non-brain rows identified by visual review of the TV-slice figures on +# 2026-10-04 (spine / neck / torso T1c with vertebral-slab or neck targets; +# series filenames carry no anatomical marker). All ntuh. +# +# Review protocol: classify at FULL resolution (data/tv_slices/exclusion/), +# using the mid-cranial slice (widest head-mask axial), not just the +# TV-max slice — skull-base / posterior-fossa brain targets look like +# cervical spine at the TV-max level (6 such false positives were found and +# removed on 2026-10-04 after expert inspection: ntuh_5902230, ntuh_5231609, +# ntuh_5083686, ntuh_5157194, ntuh_5303758, ntuh_6373987). TV-size gating +# also misses small vertebral targets (e.g. ntuh_5022161, 16k mm3); the +# mid-cranial slice is the only reliable discriminator, so the full ntuh +# population was audited slice-by-slice on 2026-10-05. +NON_BRAIN_KEYS = frozenset({ + # cervical / thoracic spine slabs + "ntuh_5237879_case2009.12.10.13.25.32_20091210120310", + "ntuh_5584533_case2011.03.18.13.24.12_20110318115157", + "ntuh_4500992_case2009.08.13.13.51.21_20090813120210", + "ntuh_2922239_case2010.08.23.15.43.52_20100823123647", + "ntuh_0915393_case2012.12.14.17.58.27_20121214121137", + "ntuh_5085147_case2013.05.31.12.51.10_20130531121922", + "ntuh_2486132_case2013.04.24.14.32.34_20130424122008", + "ntuh_4493086_case2009.08.06.13.37.10_20090806114937", + "ntuh_5185487_case2009.07.27.14.18.34_20090727114447", + "ntuh_2177312_case2012.09.10.15.19.55_20120910111101", + "ntuh_5302085_case2014.01.06.14.04.54_20140106123022", + "ntuh_6069590_case2014.08.28.14.06.40_20140828105732", + "ntuh_4408551_case2012.06.28.16.21.19_20120316123153", + "ntuh_4408551_case2012.06.28.16.21.19_20120628114905", + "ntuh_4408551_case2012.03.16.15.05.30_20120316123153", + "ntuh_2444262_case2009.07.07.13.36.39_20090707113235", + "ntuh_5406976_case2009.02.12.14.05.41_20090212112402", + "ntuh_8009538_case2009.11.04.12.46.23_20091104123842", + "ntuh_3563834_case2010.04.20.13.15.45_20100420113442", + "ntuh_5712878_case2011.12.05.14.11.00_20111205120244", + "ntuh_4272908_case2016.06.24.13.33.59_20160624113223", + "ntuh_4494910_case2008.06.03.13.40.34_20080603114804", + "ntuh_2167997_case2016.05.23.13.44.36_20160523114640", + "ntuh_4572226_case2017.11.06.14.11.35_20171106115826", + "ntuh_3420930_case2012.07.30.15.01.30_20120730111910", + "ntuh_4219462_case2010.05.31.13.50.31_20100601120540", + "ntuh_4572226_case2017.11.06.14.12.33_20171106115826", + "ntuh_5022161_case2011.03.02.13.24.52_20110302113819", + "ntuh_5650997_case2011.06.29.13.21.10_20110629113315", + "ntuh_6264497_case2016.09.21.14.48.29_20160921170338", + "ntuh_6182961_case2016.08.08.13.07.21_20160808123711", + "ntuh_2325341_case2010.12.10.13.47.32_20101210122902", + "ntuh_5552572_case2010.04.13.13.27.12_20100413120757", + "ntuh_0711619_case2011.05.10.13.47.39_20110509161659", + "ntuh_5745927_case2013.09.17.13.27.40_20130917115939", + "ntuh_4615006_case2010.07.29.13.43.30_20100722121347", + "ntuh_5332969_case2008.07.21.13.40.07_20080721140145", + "ntuh_5808597_case2012.07.02.15.22.10_20120702115701", + "ntuh_5886224_case2013.12.26.14.59.16_20131226115857", + # byte-identical duplicate images of reviewed non-brain rows above (same + # T1c file registered under a second case-folder name; matched by MD5 on + # 2026-10-05), so they escape the row-level list and the TV-size gates + "ntuh_2177312_case2012.09.10.16.05.26_20120910111101", + "ntuh_2922239_case2010.08.23.14.46.37_20100823123647", + "ntuh_4408551_case2012.06.28.16.14.15_20120628114905", + "ntuh_5085147_case2013.05.31.12.56.43_20130531121922", +}) + +_CC_STRUCTURE = np.ones((3, 3, 3), dtype=int) + + +def label_stats(label_path): + """{total: voxel count, n_cc: connected-component count (26-conn)} for a + processed label (1mm isotropic, so voxels == mm3).""" + import SimpleITK as sitk + m = sitk.GetArrayFromImage(sitk.ReadImage(label_path)) > 0.5 + total = int(m.sum()) + n_cc = int(ndimage.label(m, structure=_CC_STRUCTURE)[1]) + return {"total": total, "n_cc": n_cc} + + +def exclusion_reason(key, stats): + """Return the exclusion reason for a labeled row, or None to keep it.""" + if stats["total"] > MAX_TOTAL_VOL: + return f"total TV {stats['total']} mm3 > {MAX_TOTAL_VOL} mm3" + if stats["n_cc"] == 1 and stats["total"] < MIN_SINGLE_TARGET_VOL: + return (f"single target, TV {stats['total']} mm3 " + f"< {MIN_SINGLE_TARGET_VOL} mm3") + if key in NON_BRAIN_KEYS: + return "non-brain image (spine/neck, reviewed 2026-10-04)" + return None \ No newline at end of file