From 12143e5bb8b8d9264ae1ea313529b56e720afc03 Mon Sep 17 00:00:00 2001 From: Furen Xiao Date: Sun, 4 Oct 2026 14:59:39 +0800 Subject: [PATCH] fix(label_qa): resample TV warps to the 0.5 mm T1c grid Warp the binary TV mask onto the resampled 0.5 mm T1c grid before visual QA, derive SMD thresholds from the native CT spacing using normal-distribution theory, and widen the zoom context to 24 voxels. Keep metrics on the native grid while rendering 0.5 mm slices, and rename warped outputs to *_tv05.nii.gz. --- scripts/qa_label_interpolation.py | 323 ++++++++++++++++++++++++------ 1 file changed, 262 insertions(+), 61 deletions(-) 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)