xfr/longitudinal
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

fix(label_qa): resample TV warps to the 0.5 mm T1c grid

Browse source 
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.
This commit is contained in:
Furen Xiao 2026-10-04 14:59:39 +08:00
parent 89bb1e6331
commit 12143e5bb8
1 changed files with 262 additions and 61 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

321
scripts/qa_label_interpolation.py
View file

@ -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/<subject>_<date>/):
<reg>_T.h5 registration transform (cached; re-runs skip it)
<reg>__<warp_key>_tv.nii.gz warped TV on the T1c grid (warp_key includes thr)
<reg>__<warp_key>_tv05.nii.gz warped TV on the 0.5mm T1c grid (warp_key includes thr)
<subject>_<date>_label_qa_<reg>.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/<subject>_<date>)")
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:
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)