diff --git a/scripts/qa_label_interpolation.py b/scripts/qa_label_interpolation.py index e856ea1..a23d1cb 100644 --- a/scripts/qa_label_interpolation.py +++ b/scripts/qa_label_interpolation.py @@ -9,20 +9,25 @@ For one M6 case (default: 2BBN5XA4 / 20250411): ants_rigid : Rigid[0.25] (Mattes, 500x500x200, 4x2x1) ants_affine : Affine[0.25] (Mattes, 2100x1200x1200x10, 6x4x2x1) ants_syn : SyN (Affine[0.25] + SyN[0.2,3,0] stages) - - the binary TV is then warped to the T1c grid by nearest-neighbor (nn) and - by SignedMaurerDistanceMap (mm, inside<0) -> {linear, B-spline order-3} - interp -> threshold (smd <= thr, default 0.5) + - the binary TV is then warped ONTO THE 0.5mm RESAMPLED T1c GRID by + nearest-neighbor (nn) and by SignedMaurerDistanceMap (mm, inside<0) + -> {linear, B-spline order-3} interp -> threshold (smd <= thr), where the + candidate thresholds are derived per case from the CT grid spacing + (normal-distribution theory, see smd_threshold()): the mean spacing + halved, thr = (sx + sy + sz) / 6; the registration transform is + grid-independent, so it applies unchanged - one PNG per registration: 3 rows (axial / coronal / sagittal) x 4 columns (leftmost = bare MR, then nn + each SMD method at its threshold); all - columns show the SAME slice - (the nn warp's largest-TV cross-section per view), ZOOMED to the TV +columns show the SAME slice (the nn warp's largest-TV cross-section per + view), ZOOMED to the TV, rendered on a 0.5mm display grid (MR linear, + masks NN; metrics stay on the native grid) bounding box + margin (the NN box, shared by all columns); only the 1-voxel TV contour is drawn (yellow, no fill), plus the volume + % of the original CT-grid TV volume in the column header. Outputs (default results/label_interp_qa/_/): _T.h5 registration transform (cached; re-runs skip it) - ___tv.nii.gz warped TV on the T1c grid (warp_key includes thr) + ___tv05.nii.gz warped TV on the 0.5mm T1c grid (warp_key includes thr) __label_qa_.png one 3x5 figure per registration summary.json per-combination metrics (nvox, mm3, dice, shifts) + the SMD thresholding volume-preservation check @@ -57,15 +62,36 @@ from src import registration_pipelines as RP M6 = RP._load_m6() REGS = ["ants_rigid", "ants_affine", "ants_syn"] -SMD_THRESHOLDS = [0.5] # SMD mask = (resampled smd <= thr); 0.0 / 1.0 were -# tested earlier (36% under / 40% over volume) and dropped from the figure METHOD_LABELS = {"nn": "NN", "smd_linear": "SMD + linear", "smd_bspline": "SMD + B-spline(3)"} VIEWS = ["axial", "coronal", "sagittal"] SMD_DEFAULT = 1e5 # out-of-grid fill for the SMD volume: far OUTSIDE (inside < 0) -# figure / warp columns: NN, then each SMD method at each threshold -COLS = ([("nn", None)] - + [(m, t) for m in ("smd_linear", "smd_bspline") for t in SMD_THRESHOLDS]) -ZOOM_MARGIN = 12 # voxels of context around the TV bbox in the displayed slices +ZOOM_MARGIN = 24 # voxels of context around the TV bbox (0.5mm grid = 12 mm) + + +def smd_threshold(spacing): + """Final SMD threshold from the normal-distribution theory. + + The SMD zero plane is inset from the true boundary by half a voxel + projected onto the boundary normal: s(n)/2 with + s(n) = sx|nx| + sy|ny| + sz|nz|, which ranges over the bounds + thr_lo = min(sx, sy, sz) / 2 (normal along the smallest spacing) + thr_hi = max(sx, sy, sz) / 2 (normal along the largest spacing) + for the axis-aligned normals a voxelized surface actually has. The + final single scalar threshold (best mean-abs-volume error on the 9-case + batch) is the mean spacing halved: + thr = (sx + sy + sz) / 6 (e.g. (0.586+0.586+1.25)/6 = 0.404) + Returns (thr, thr_lo, thr_hi).""" + s = [abs(float(x)) for x in spacing] + lo, hi = 0.5 * min(s), 0.5 * max(s) + return round(sum(s) / 6.0, 3), lo, hi + + +def cols(thr): + """Figure / warp columns: NN, then each SMD method at the threshold.""" + return [("nn", None)] + [(m, thr) for m in ("smd_linear", "smd_bspline")] +DISP_SPACING = 0.5 # display-grid spacing (mm): MR/masks are resampled to this +# for the PNG so the zoomed slices render smoothly; all metrics stay on the +# native T1c grid def warp_key(method, thr): @@ -166,6 +192,37 @@ def best_slices(mask, shape): "sagittal": (int(np.argmax(ax)), int(ax.max()))} +def resample_grid(img, spacing=0.5, nn=False): + """Resample ``img`` onto an isotropic ``spacing`` grid (same origin and + direction). Linear for scalar images, nearest for label/mask images.""" + size = [max(1, int(round(s * sp / spacing))) + for s, sp in zip(img.GetSize(), img.GetSpacing())] + ref = sitk.Image(size, img.GetPixelID()) + ref.SetSpacing((spacing,) * 3) + ref.SetOrigin(img.GetOrigin()) + ref.SetDirection(img.GetDirection()) + rs = sitk.ResampleImageFilter() + rs.SetReferenceImage(ref) + rs.SetTransform(sitk.Transform(3, sitk.sitkIdentity)) + rs.SetInterpolator(sitk.sitkNearestNeighbor if nn else sitk.sitkLinear) + rs.SetDefaultPixelValue(0) + rs.SetOutputPixelType(img.GetPixelID()) + return rs.Execute(img) + + +def largest_cc(mask): + """The largest connected component of the mask (for slice selection and + the zoom box when the TV has multiple separated targets). Returns + (cc_mask, n_components); an empty mask gives (empty, 0).""" + lab, n = ndimage.label(mask) + if n == 0: + return mask.astype(np.uint8), 0 + if n == 1: + return (lab > 0).astype(np.uint8), 1 + sizes = ndimage.sum(mask, lab, index=range(1, n + 1)) + return (lab == (int(np.argmax(sizes)) + 1)).astype(np.uint8), n + + def tv_crop_box(mask, shape, margin): """(lo, hi) 3-D index box around the mask (per axis of (z, y, x)) with a margin; centered box when the mask is empty.""" @@ -204,16 +261,16 @@ def centroid_mm(mask, spacing): return c -def check_volume(tv_bin, t1c, source_vol_mm3): +def check_volume(tv_bin, t1c, source_vol_mm3, colspec): """Does thresholding the SMD preserve the TV volume? Resampled with the IDENTITY transform (CT grid -> MR grid, no registration), which isolates the interpolation + thresholding step from registration error. One row - per figure column (NN + each SMD method x threshold).""" + per figure column (NN + each SMD method at its threshold).""" T = sitk.Transform(3, sitk.sitkIdentity) sp = [float(s) for s in t1c.GetSpacing()] vv = sp[0] * sp[1] * sp[2] out = {} - for method, thr in COLS: + for method, thr in colspec: mask = warp_tv(tv_bin, t1c, T, method, thr or 0.0) out[warp_key(method, thr)] = round(float(mask.sum() * vv), 1) print(f"volume preservation (identity resample CT grid -> MR grid): " @@ -232,7 +289,7 @@ def make_png(mr, combos, out_png, title, spacing, window, source_vol): bare-MR column, "% of it" under each warp column.""" lo, hi = window ncol, nrow = len(combos) + 1, len(VIEWS) - cw, ch, gap, m = 3.6, 3.4, 0.15, 0.62 + cw, ch, gap, m = 3.0, 3.4, 0.15, 0.62 top = 1.5 fig_w = ncol * cw + (ncol - 1) * gap + 2 * m fig_h = nrow * ch + (nrow - 1) * gap + m + top @@ -293,24 +350,42 @@ def make_png(mr, combos, out_png, title, spacing, window, source_vol): plt.close(fig) -def main(): - ap = argparse.ArgumentParser(description=__doc__, - formatter_class=argparse.RawDescriptionHelpFormatter) - ap.add_argument("--subject", default="2BBN5XA4") - ap.add_argument("--date", default="20250411") - ap.add_argument("--regs", nargs="*", default=REGS, - help="subset of: " + " ".join(REGS)) - ap.add_argument("--tv", default=None, - help="explicit TV file (default: first TV-*.nii.gz in RT)") - ap.add_argument("--out", default=None, - help="output dir (default results/label_interp_qa/_)") - ap.add_argument("--no-cache", action="store_true", - help="re-register and re-warp even if cached outputs exist") - args = ap.parse_args() +def ct_complete(ct): + """Production guard (copied from scripts/bench_registration.py): skip + corrupt CTs (upper brain missing).""" + vz = M6._ct_vertex_slice(ct) + if vz is None: + return False + sp = ct.GetSpacing()[2] + per = (sitk.GetArrayFromImage(ct) > -150).sum((1, 2)) + iz = vz - int(round(100 / sp)) + return bool(0 <= iz < len(per) and per[iz] > 5000) - row = find_case(args.subject, args.date) - case_dir = args.out or os.path.join(ROOT, "results", "label_interp_qa", - f"{args.subject}_{args.date}") + +def select_rows(args): + """Cases to run: an explicit --subject/--date, or the first / a seeded + random N labeled M6 cases.""" + if args.n_cases: + labeled, _unlabeled = M6.discover() + pool = list(labeled) + print(f"[cases] m6 labeled candidates (CT+T1c+TV): {len(pool)}") + if args.sample == "random": + rng = np.random.default_rng(args.seed) + idx = rng.choice(len(pool), min(args.n_cases, len(pool)), replace=False) + sel = [pool[i] for i in sorted(idx.tolist())] + else: + sel = pool[:args.n_cases] + print(f"[cases] selected {len(sel)} (sample={args.sample} seed={args.seed})") + for r in sel: + print(f" - {r['key']}") + return sel + return [find_case(args.subject, args.date)] + + +def run_case(row, args, out_base): + ap_case_dir = args.out + case_dir = ap_case_dir or os.path.join(out_base, + f"{row['subject']}_{row['date']}") os.makedirs(case_dir, exist_ok=True) ct_path = row["ct"] t1c_path = row["img"] @@ -322,45 +397,62 @@ def main(): t1c = read_nii_img(t1c_path) tv = read_nii_img(tv_path) tv_bin = make_tv_bin(tv) - mr = sitk.GetArrayFromImage(t1c).astype("float32") - spacing = [float(s) for s in t1c.GetSpacing()] # (x, y, z) - print(f"t1c {t1c.GetSize()} @ {spacing}mm, tv nvox={int((sitk.GetArrayFromImage(tv_bin) > 0).sum())} " + # the TV interpolation targets the 0.5mm resampled MR grid (same origin / + # direction; the registration transform T is grid-independent, so it + # applies unchanged) + t1c05 = resample_grid(t1c, DISP_SPACING) + mr05 = sitk.GetArrayFromImage(t1c05).astype("float32") + spacing = [DISP_SPACING] * 3 # warp-grid spacing (x, y, z) + print(f"t1c {t1c.GetSize()} @ {[float(s) for s in t1c.GetSpacing()]}mm -> warp grid " + f"{t1c05.GetSize()} @ {DISP_SPACING}mm, tv nvox={int((sitk.GetArrayFromImage(tv_bin) > 0).sum())} " f"({int((sitk.GetArrayFromImage(tv_bin) > 0).sum()) * ct.GetSpacing()[0] * ct.GetSpacing()[1] * ct.GetSpacing()[2]:.0f} mm3 on the CT grid)", flush=True) # MR display window from the head - pos = mr > 0 - thr = np.percentile(mr[pos], 2) if pos.any() else 0.0 - head = mr > max(thr, 1e-6) - lo, hi = (np.percentile(mr[head], (2, 98)) if head.any() else (0.0, 1.0)) + pos = mr05 > 0 + thr = np.percentile(mr05[pos], 2) if pos.any() else 0.0 + head = mr05 > max(thr, 1e-6) + lo, hi = (np.percentile(mr05[head], (2, 98)) if head.any() else (0.0, 1.0)) - # volume-preservation check for the SMD thresholding (no registration) + # volume-preservation check for the SMD thresholding (no registration); + # the SMD threshold is derived per case from the CT grid spacing sp_ct = [float(s) for s in tv_bin.GetSpacing()] source_vol = float((sitk.GetArrayFromImage(tv_bin) > 0).sum()) * sp_ct[0] * sp_ct[1] * sp_ct[2] - vol_check = check_volume(tv_bin, t1c, source_vol) + THR, THR_LO, THR_HI = smd_threshold(sp_ct) + colspec = cols(THR) + print(f"smd threshold from CT spacing {['%.3f' % s for s in sp_ct]}: " + f"bounds lo={THR_LO:.3f} hi={THR_HI:.3f} -> used (sx+sy+sz)/6 = {THR:.3f}", + flush=True) + vol_check = check_volume(tv_bin, t1c05, source_vol, colspec) combos = [] reg_meta = {} for reg in args.regs: - T, meta = get_transform(reg, ct, t1c, case_dir, use_cache=not args.no_cache) + try: + T, meta = get_transform(reg, ct, t1c, case_dir, use_cache=not args.no_cache) + except Exception as e: + reg_meta[reg] = {"error": repr(e)} + print(f"[{reg}] FAILED: {e!r}", flush=True) + continue reg_meta[reg] = meta reg_combos = [] - # warp every column (NN + each SMD method x threshold) - for method, thr in COLS: + # warp every column (NN + each SMD method at the case threshold) -> 0.5mm grid + for method, thr in colspec: key = warp_key(method, thr) - p = os.path.join(case_dir, f"{reg}__{key}_tv.nii.gz") + p = os.path.join(case_dir, f"{reg}__{key}_tv05.nii.gz") if not args.no_cache and os.path.exists(p): print(f"[{reg}/{key}] loading cached warp {p}", flush=True) mask = (sitk.GetArrayFromImage(read_nii_img(p)) > 0).astype(np.uint8) dt = 0.0 else: t0 = time.time() - mask = warp_tv(tv_bin, t1c, T, method, thr or 0.0) + mask = warp_tv(tv_bin, t1c05, T, method, thr or 0.0) dt = time.time() - t0 - mimg = sitk.GetImageFromArray(mask) - mimg.CopyInformation(t1c) - sitk.WriteImage(mimg, p, True) nvox = int(mask.sum()) + mimg = sitk.GetImageFromArray(mask) + mimg.CopyInformation(t1c05) + if not os.path.exists(p) or args.no_cache: + sitk.WriteImage(mimg, p, True) reg_combos.append({"reg": reg, "method": method, "thr": thr, "key": key, "ilabel": col_label(method, thr), "mask": mask, "nvox": nvox, @@ -368,14 +460,21 @@ def main(): "seconds_warp": round(dt, 1)}) print(f"[{reg}/{key}] nvox={nvox} ({reg_combos[-1]['vol_mm3']} mm3) ({dt:.1f}s)", flush=True) - # all columns show the SAME slices (the NN warp's largest-area slice - # per view) and the SAME zoom box (NN mask bbox + margin) + # all columns show the SAME slices and the SAME zoom box, selected from the + # LARGEST connected component of the NN mask (when the TV has several + # separated targets, only the biggest one drives the view), on the + # 0.5mm warp grid nn = next(c for c in reg_combos if c["method"] == "nn") - sl = best_slices(nn["mask"], mr.shape) - lo_b, hi_b = tv_crop_box(nn["mask"], mr.shape, ZOOM_MARGIN) + nn_cc, n_ccs = largest_cc(nn["mask"]) + if n_ccs > 1: + print(f"[{reg}] TV has {n_ccs} components; slices/zoom follow the " + f"largest ({int(nn_cc.sum())} of {int(nn['mask'].sum())} voxels)", + flush=True) + sl = best_slices(nn_cc, mr05.shape) + lo_b, hi_b = tv_crop_box(nn_cc, mr05.shape, ZOOM_MARGIN) for c in reg_combos: c["slices"] = [sl[k][0] for k in VIEWS] - c["views"] = {kind: {"mr": view_crop(mr, kind, sl[kind][0], lo_b, hi_b), + c["views"] = {kind: {"mr": view_crop(mr05, kind, sl[kind][0], lo_b, hi_b), "mask": view_crop(c["mask"], kind, sl[kind][0], lo_b, hi_b), "idx": sl[kind][0], "area": sl[kind][1]} for kind in VIEWS} c["pct_of_source"] = 100.0 * c["vol_mm3"] / source_vol @@ -397,21 +496,31 @@ def main(): pngs = {} for reg in args.regs: - title = (f"{args.subject} {args.date} -- {reg}: TV warping QA " + reg_combos = [c for c in combos if c["reg"] == reg] + if not reg_combos: + continue + title = (f"{row['subject']} {row['date']} -- {reg}: TV warping QA " f"(same slices = NN largest-area per view)\n" f"fixed=T1c ({os.path.basename(t1c_path)[:40]}...) moving=CT | " f"yellow = TV contour | MR window {lo:.0f}..{hi:.0f}") - png = os.path.join(case_dir, f"{args.subject}_{args.date}_label_qa_{reg}.png") - make_png(mr, [c for c in combos if c["reg"] == reg], png, title, spacing, - (float(lo), float(hi)), source_vol) + png = os.path.join(case_dir, + f"{row['subject']}_{row['date']}_label_qa_{reg}.png") + make_png(mr05, reg_combos, png, title, + [DISP_SPACING] * 3, (float(lo), float(hi)), source_vol) pngs[reg] = png print(f"png: {png}", flush=True) summary = { - "case": row["key"], "subject": args.subject, "date": args.date, + "case": row["key"], "subject": row["subject"], "date": row["date"], "inputs": {"ct": ct_path, "t1c": t1c_path, "tv": tv_path}, - "t1c_geometry": {"size": list(t1c.GetSize()), "spacing": spacing, + "t1c_geometry": {"size": list(t1c.GetSize()), + "spacing": [float(s) for s in t1c.GetSpacing()], "origin": [float(o) for o in t1c.GetOrigin()]}, + "warp_grid": {"spacing": DISP_SPACING, "size": list(t1c05.GetSize()), + "note": "the TV interpolation (all methods) runs on this grid"}, + "smd_threshold": {"ct_spacing": [round(s, 3) for s in sp_ct], + "lo": round(THR_LO, 3), "hi": round(THR_HI, 3), + "formula": "(sx+sy+sz)/6", "used": round(THR, 3)}, "regs": {k: {kk: vv for kk, vv in v.items() if kk != "stages"} for k, v in reg_meta.items()}, "volume_check": {"source_tv_mm3": round(source_vol, 1), @@ -429,6 +538,98 @@ def main(): with open(sf, "w") as f: json.dump(summary, f, indent=1) print(f"summary: {sf}", flush=True) + return summary + + +def print_batch_table(rows, regs): + """Compact cross-case table of the warp volumes (% of the CT TV).""" + short = {"nn": "nn", "smd_linear": "lin", "smd_bspline": "bspl"} + + def key_label(key): + if key == "nn": + return "nn" + m, _, t = key.partition("_thr") + return short[m] + f"<={t}" + + def reg_keys(s, r): + return [k.split("/", 1)[1] for k in s.get("combos", {}) + if k.startswith(r + "/")] + + width = 28 + 7 + 49 * len(regs) + print("\n" + "=" * width) + first = next((s for s in rows if not s.get('error')), None) + print(f"{'case':28s} {'CT TV':>6s} | " + + " | ".join(f"{r:48s}" for r in regs) + + (" (% of CT TV; per reg: " + + " ".join(key_label(k) for k in reg_keys(first, regs[0])) + ")" + if first else "")) + print("-" * width) + for s in rows: + if s.get("error"): + print(f"{s['case']:28s} {'':>6s} | ERROR {s['error'][:70]}") + continue + cells = [] + for r in regs: + cells.append(" ".join( + f"{s['combos'].get(f'{r}/{k}', {}).get('pct_of_source_tv', '-'):>4}" + for k in reg_keys(s, r))) + print(f"{s['case']:28s} {s['volume_check']['source_tv_mm3']:6.0f} | " + + " | ".join(f"{c:48s}" for c in cells)) + print("=" * width) + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--subject", default="2BBN5XA4", + help="case (ignored when --n-cases is set)") + ap.add_argument("--date", default="20250411", + help="case (ignored when --n-cases is set)") + ap.add_argument("--n-cases", type=int, default=0, + help="QA the first (or --sample random) N labeled M6 cases " + "instead of --subject/--date") + ap.add_argument("--sample", choices=["first", "random"], default="first", + help="batch selection: first N by discover order, or a " + "seeded random sample") + ap.add_argument("--seed", type=int, default=0, help="seed for --sample random") + ap.add_argument("--regs", nargs="*", default=REGS, + help="subset of: " + " ".join(REGS)) + ap.add_argument("--tv", default=None, + help="explicit TV file (single-case mode only)") + ap.add_argument("--out", default=None, + help="single output dir (single-case mode only)") + ap.add_argument("--no-cache", action="store_true", + help="re-register and re-warp even if cached outputs exist") + args = ap.parse_args() + if args.n_cases and args.tv: + raise SystemExit("--tv is only valid for a single case") + + out_base = os.path.join(ROOT, "results", "label_interp_qa") + rows = select_rows(args) + results = [] + for i, row in enumerate(rows): + print(f"\n[{i + 1}/{len(rows)}] ===== case {row['key']} =====", flush=True) + try: + ct = read_nii_img(row["ct"]) + if not ct_complete(ct): + print(f"skip {row['key']}: CT head incomplete (data issue)", flush=True) + results.append({"case": row["key"], "error": "ct_incomplete"}) + continue + results.append(run_case(row, args, out_base)) + except Exception as e: + print(f"case {row['key']}: ERROR {e!r}", flush=True) + results.append({"case": row["key"], "error": repr(e)}) + + if args.n_cases: + ok = [s for s in results if not s.get("error")] + print_batch_table(results, args.regs) + bs = os.path.join(out_base, "batch_summary.json") + with open(bs, "w") as f: + json.dump({"n_cases": len(rows), + "selection": f"{args.sample} seed={args.seed}", + "regs": args.regs, + "cases": results}, f, indent=1) + print(f"batch summary: {bs} ({len(ok)}/{len(rows)} cases ok)") print("done", flush=True)