Michał Szczepanik
e744431151
The focus is on visual grepping, ie. printing a few lines for each subject to get an overview of the situation. We require JSON sidecars to be present (to have been gotten).
77 lines
2.9 KiB
Python
77 lines
2.9 KiB
Python
# /// script
|
|
# dependencies = ["numpy"]
|
|
# ///
|
|
#
|
|
# This script prints whether shims are identical within a subset of
|
|
# acquisitions, and how far off they are from the first. See arguments
|
|
# (or --help) and code comments for details.
|
|
|
|
import argparse
|
|
import json
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
|
|
def load_shim(p: str | Path):
|
|
p = Path(p) if isinstance(p, str) else p
|
|
with p.open() as fp:
|
|
d = json.load(fp)
|
|
return np.array(d["ShimSetting"])
|
|
|
|
|
|
parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
|
|
parser.add_argument("ds_path", type=Path, help="Path to BIDS dataset")
|
|
parser.add_argument("--func", action="store_true", help="Include func in comparison")
|
|
parser.add_argument("--dwi", action="store_true", help="Include dwi in comparison")
|
|
parser.add_argument("--session", default="001", help="Session label")
|
|
parser.add_argument("--subject", help="Focus on this subject (label) and print more")
|
|
args = parser.parse_args()
|
|
|
|
max_rel_deltas = []
|
|
|
|
for subject_path in sorted(list(args.ds_path.glob("sub-*"))):
|
|
|
|
if args.subject is not None and subject_path.name.split("-")[-1] != args.subject:
|
|
continue
|
|
|
|
# load shims
|
|
fmap_sidecars = list((subject_path / f"ses-{args.session}" / "fmap").glob("*.json"))
|
|
func_sidecars = list((subject_path / f"ses-{args.session}" / "func").glob("*.json"))
|
|
dwi_sidecars = list((subject_path / f"ses-{args.session}" / "dwi").glob("*.json"))
|
|
if not all(p.is_file() for p in fmap_sidecars + func_sidecars + dwi_sidecars):
|
|
print(subject_path.name, "----")
|
|
continue
|
|
if len(fmap_sidecars) == 0:
|
|
print(subject_path.name, "----")
|
|
continue
|
|
fmap_shims = np.array([load_shim(p) for p in sorted(fmap_sidecars)]) # file x shim
|
|
func_shims = np.array([load_shim(p) for p in sorted(func_sidecars)])
|
|
dwi_shims = np.array([load_shim(p) for p in sorted(dwi_sidecars)])
|
|
|
|
# restrict to chosen modalities
|
|
to_stack = [fmap_shims]
|
|
if args.func:
|
|
to_stack.append(func_shims)
|
|
if args.dwi:
|
|
to_stack.append(dwi_shims)
|
|
shims = np.vstack(to_stack)
|
|
|
|
all_equal = np.all(shims == shims[0])
|
|
rel_delta = np.abs((shims - shims[0]) / shims[0])
|
|
# rel_vec_dist = np.linalg.norm(shims - shims[0], axis=1) / np.linalg.norm(shims[0])
|
|
|
|
# if showing single subject, print all shims, and biggest relative difference per dimension
|
|
if args.subject is not None:
|
|
print(shims)
|
|
print(np.round(np.max(rel_delta, axis=0), 3))
|
|
exit()
|
|
|
|
# if showing all subjects, print whether all shims are equal, followed by biggest relative difference if non-zero
|
|
if all_equal:
|
|
print(subject_path.name, all_equal)
|
|
max_rel_deltas.append(0)
|
|
else:
|
|
max_rel_deltas.append(np.max(rel_delta))
|
|
print(subject_path.name, all_equal, np.round(np.max(rel_delta), 3))
|
|
|
|
print("Median of max relative deltas", np.round(np.median(np.array(max_rel_deltas)), 3))
|