"""Minimal resting-state preprocessing for OpenNeuro ds002748 (one subject per call). Steps: drop the first 4 volumes; rigid-body realignment of every volume to the mean volume (DIPY, mutual information); framewise displacement (Power et al. 2012, 50 mm radius); 12-parameter affine registration of the mean volume to the MNI152 template bundled with nilearn (mutual information); resampling of the realigned series to MNI space at 3 mm; nuisance regression of 6 motion parameters, their temporal derivatives and the first 5 principal components of the white-matter signal, with linear detrending and a 0.01-0.1 Hz band-pass; averaging within each of the 100 Schaefer-2018 (7-network) parcels. Usage: python code/preprocess.py Writes _parcels.npy (time x 100) and _qc.json. """ import json import sys import nibabel as nib import numpy as np from dipy.align.imaffine import AffineRegistration, MutualInformationMetric, AffineMap from dipy.align.transforms import RigidTransform3D, TranslationTransform3D, AffineTransform3D from nilearn import datasets, image, masking, signal TR = 2.5 DROP = 4 FD_RADIUS = 50.0 def register(static, static_aff, moving, moving_aff, kind, starting=None, level_iters=(100, 30, 10)): k = len(level_iters) reg = AffineRegistration(metric=MutualInformationMetric(32, None), level_iters=list(level_iters), sigmas=[3.0, 1.0, 0.0][-k:], factors=[4, 2, 1][-k:], verbosity=0) params0 = None if kind == "rigid": return reg.optimize(static, moving, RigidTransform3D(), params0, static_aff, moving_aff, starting_affine=starting) tr = reg.optimize(static, moving, TranslationTransform3D(), params0, static_aff, moving_aff, starting_affine="mass") rg = reg.optimize(static, moving, RigidTransform3D(), params0, static_aff, moving_aff, starting_affine=tr.affine) return reg.optimize(static, moving, AffineTransform3D(), params0, static_aff, moving_aff, starting_affine=rg.affine) def rigid_params(mat): """Translations (mm) and rotations (rad, x-y-z Euler) from a 4x4 rigid matrix.""" R = mat[:3, :3] t = mat[:3, 3] ry = np.arcsin(-np.clip(R[2, 0], -1, 1)) rx = np.arctan2(R[2, 1], R[2, 2]) rz = np.arctan2(R[1, 0], R[0, 0]) return np.r_[t, rx, ry, rz] def main(): bold_path, atlas_path, out = sys.argv[1], sys.argv[2], sys.argv[3] bold = nib.load(bold_path) data = np.asarray(bold.dataobj, dtype=np.float32)[..., DROP:] aff = bold.affine n = data.shape[-1] ref = data.mean(-1) # realignment: each volume to the mean volume (one pass) mats = [register(ref, aff, data[..., i], aff, "rigid", starting=np.eye(4), level_iters=(25, 10)).affine for i in range(n)] realigned = np.stack([AffineMap(mats[i], ref.shape, aff, ref.shape, aff).transform(data[..., i]) for i in range(n)], -1) ref = realigned.mean(-1) params = np.array([rigid_params(m) for m in mats]) d = np.abs(np.diff(params, axis=0)) fd = np.r_[0.0, d[:, :3].sum(1) + FD_RADIUS * d[:, 3:].sum(1)] # mean volume to MNI (12-parameter affine) tmpl = datasets.load_mni152_template(resolution=2) tdata = tmpl.get_fdata() ref_img = nib.Nifti1Image(ref, aff) brain = masking.compute_epi_mask(ref_img) ref_brain = ref * brain.get_fdata() to_mni = register(tdata, tmpl.affine, ref_brain, aff, "affine") grid = image.resample_img(tmpl, target_affine=np.diag([3.0, 3.0, 3.0]), interpolation="nearest", force_resample=True, copy_header=True) gshape, gaff = grid.shape[:3], grid.affine m = AffineMap(to_mni.affine, gshape, gaff, ref.shape, aff) mni = np.stack([m.transform(realigned[..., i]) for i in range(n)], -1) mni_img = nib.Nifti1Image(mni.astype(np.float32), gaff) # coverage: voxels with signal in every volume covered = (mni > 0).all(-1) # nuisance: motion, derivatives, white-matter PCA wm = image.resample_to_img(datasets.load_mni152_wm_mask(resolution=2, threshold=0.5), mni_img, interpolation="nearest", force_resample=True, copy_header=True).get_fdata() > 0 wm &= covered wm_ts = mni[wm].T wm_ts = (wm_ts - wm_ts.mean(0)) / (wm_ts.std(0) + 1e-6) u, s, vt = np.linalg.svd(wm_ts - wm_ts.mean(0), full_matrices=False) conf = np.column_stack([params, np.vstack([np.zeros(6), np.diff(params, axis=0)]), u[:, :5]]) # parcels atlas = image.resample_to_img(nib.load(atlas_path), mni_img, interpolation="nearest", force_resample=True, copy_header=True).get_fdata().astype(int) ts = np.zeros((n, 100)) cover = np.zeros(100) for p in range(1, 101): vox = atlas == p cover[p - 1] = (vox & covered).sum() / max(vox.sum(), 1) sel = vox & covered ts[:, p - 1] = mni[sel].mean(0) if sel.sum() else np.nan good = np.isfinite(ts).all(0) clean = np.full_like(ts, np.nan) clean[:, good] = signal.clean(ts[:, good], detrend=True, standardize="zscore_sample", confounds=conf, low_pass=0.1, high_pass=0.01, t_r=TR) np.save(out + "_parcels.npy", clean) # registration check: Dice overlap of the warped EPI brain mask with the template brain mask warped_mask = AffineMap(to_mni.affine, domain_grid_shape=gshape, domain_grid2world=gaff, codomain_grid_shape=ref.shape, codomain_grid2world=aff).transform( brain.get_fdata().astype(np.float32), interpolation="nearest") > 0.5 tmask = image.resample_to_img(datasets.load_mni152_brain_mask(resolution=2), mni_img, interpolation="nearest", force_resample=True, copy_header=True).get_fdata() > 0 reg_r = float(2 * (warped_mask & tmask).sum() / (warped_mask.sum() + tmask.sum())) nib.save(nib.Nifti1Image(mni.mean(-1).astype(np.float32), gaff), out + "_meanmni.nii.gz") qc = {"volumes": int(n), "mean_fd_mm": float(fd.mean()), "max_fd_mm": float(fd.max()), "share_fd_over_0_5": float((fd > 0.5).mean()), "max_translation_mm": float(np.abs(params[:, :3]).max()), "registration_dice": reg_r, "parcel_coverage": [round(float(c), 3) for c in cover]} with open(out + "_qc.json", "w") as fh: json.dump(qc, fh) print(json.dumps({k: v for k, v in qc.items() if k != "parcel_coverage"})) if __name__ == "__main__": main()