diff --git a/brainprep/interfaces/__init__.py b/brainprep/interfaces/__init__.py index b76e9fcd..f689b867 100644 --- a/brainprep/interfaces/__init__.py +++ b/brainprep/interfaces/__init__.py @@ -96,6 +96,7 @@ copyfiles, htmlmin, movedir, + regridmask, ungzfile, write_uuid_mapping, ) @@ -147,6 +148,7 @@ "reconall_summary", "reconall_surfreg", "reconall_tissues", + "regridmask", "reorient", "scale", "subject_level_qa", diff --git a/brainprep/interfaces/utils.py b/brainprep/interfaces/utils.py index a782fe52..9a677de2 100644 --- a/brainprep/interfaces/utils.py +++ b/brainprep/interfaces/utils.py @@ -15,6 +15,7 @@ import shutil import socket +import nibabel import pandas as pd from ..decorators import ( @@ -84,6 +85,106 @@ def copyfiles( shutil.copy(src_path, dest_path) +@step( + hooks=[ + CoerceparamsHook(), + OutputdirHook(), + LogRuntimeHook( + bunched=False + ), + PythonWrapperHook(), + SignatureHook(), + ] +) +def regridmask( + image_file: File, + mask_file: File, + output_dir: Directory, + entities: dict, + dryrun: bool = False, + ) -> tuple[File]: + """ + Rewrite a mask file onto a reference image's exact grid and give it a + self-consistent qform/sform, using ``nibabel``. + + A mask produced by a different tool than the one that wrote the + reference image (e.g. FreeSurfer's `mri_synthstrip` vs. FSL's + `fslreorient2std`) can end up on a grid that is numerically close but + not bit-for-bit identical to it. Most FSL tools tolerate this, but + ANTs' `N4BiasFieldCorrection` does not: it raises + ``itk::ExceptionObject: Inputs do not occupy the same physical space`` + when handed a mask (``-x``) that isn't exactly co-registered with the + image it's correcting (``-i``), down to a ~1e-7 tolerance on the + origin. + + A NIfTI file also carries both a qform and an sform, meant to agree + but not always written that way — e.g. `fslreorient2std` can leave + them a float32 ULP or two apart — and N4's ITK reader has been + observed to read one of the two forms for its main image (``-i``) + but the other for its mask (``-x``), so two files that are each + internally self-consistent can still trip that check if regridded + through an external tool that only aligns one form (e.g. FSL's + `flirt`, which only touches sform, or `fslcpgeom`, which round-trips + the copied values through the same float32 fields and can leave + residual noise right at the tolerance's edge). + + Recomputing either form from an affine matrix doesn't reliably close + that gap either — qform is quaternion-based, and any fresh matrix -> + quaternion decomposition (nibabel's own `set_qform()` included) can + land on a different last-bit rounding than whatever decomposition + originally produced the reference image's own qform (e.g. FSL's). + The only way to get a genuinely identical header is to not recompute + anything at all: the mask's raw qform/sform NIfTI fields (quaternion + components, offsets, sform rows, codes) are transcribed byte-for-byte + from the reference image's own header. + + Parameters + ---------- + image_file : File + Path to the reference image file whose grid will be matched. + mask_file : File + Path to the mask file to regrid. + output_dir : Directory + Directory where the regridded mask will be saved. + entities : dict + A dictionary of parsed BIDS entities including modality. + dryrun : bool + If True, skip actual computation and file writing. + Default False. + + Returns + ------- + outputs : tuple[File] + - regridded_mask_file : File - the mask file, rewritten onto + image_file's exact grid with a self-consistent qform/sform. + """ + basename = "sub-{sub}_ses-{ses}_run-{run}_mod-{mod}_regridmask".format( + **entities) + regridded_mask_file = output_dir / f"{basename}.nii.gz" + if dryrun: + return (regridded_mask_file, ) + + image_img = nibabel.load(image_file) + mask_img = nibabel.load(mask_file) + regridded_mask_img = nibabel.Nifti1Image( + mask_img.get_fdata().astype(mask_img.get_data_dtype()), + image_img.affine, + header=mask_img.header, + ) + # The constructor above only syncs the mask's sform from the passed + # affine, leaving its qform untouched — so both forms are transcribed + # here directly from the reference image's own header fields (see + # the docstring for why nothing here is recomputed). + for field in ("qform_code", "sform_code", "pixdim", + "quatern_b", "quatern_c", "quatern_d", + "qoffset_x", "qoffset_y", "qoffset_z", + "srow_x", "srow_y", "srow_z"): + regridded_mask_img.header[field] = image_img.header[field] + nibabel.save(regridded_mask_img, regridded_mask_file) + + return (regridded_mask_file, ) + + @step( hooks=[ CoerceparamsHook(), diff --git a/brainprep/workflow/quasiraw.py b/brainprep/workflow/quasiraw.py index a42b1866..547d0e8d 100644 --- a/brainprep/workflow/quasiraw.py +++ b/brainprep/workflow/quasiraw.py @@ -177,6 +177,20 @@ def brainprep_quasiraw( workspace_dir / "02-brainmask", entities, ) + # mri_synthstrip (FreeSurfer) and fslreorient2std (FSL) each write their + # own NIfTI header independently, which can leave the mask on a grid + # that's numerically equivalent but not bit-for-bit identical to the + # reoriented image. ANTs' N4BiasFieldCorrection (unlike most FSL tools) + # enforces exact geometry equality between -i and -x and raises + # "Inputs do not occupy the same physical space" otherwise, so the + # mask is resampled onto the image's exact grid before it's used + # further. + mask_file = interfaces.regridmask( + reoriented_anatomical_file, + mask_file, + workspace_dir / "02-regridmask", + entities, + ) bc_anatomical_file, _ = interfaces.biasfield( reoriented_anatomical_file, mask_file,