From 004ac4545f272d511caceca20dac2cdc77d4e0d1 Mon Sep 17 00:00:00 2001 From: Benoit Dufumier Date: Wed, 23 Sep 2026 11:59:23 +0200 Subject: [PATCH] quasiraw: fix N4BiasFieldCorrection "not same physical space" crash MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit brainprep_quasiraw() was passing brainmask()'s output straight to ANTs' N4BiasFieldCorrection as the -x mask. mri_synthstrip (FreeSurfer) and fslreorient2std (FSL) each write their own NIfTI header independently, so the mask can end up on a grid that's numerically equivalent to the reoriented image but not bit-for-bit identical — down to the qform and sform of a single file sometimes disagreeing with each other by a float32 ULP or two. N4's ITK reader enforces exact geometry equality between -i and -x and aborts otherwise with: "itk::ExceptionObject: Inputs do not occupy the same physical space". Observed on real data (OpenMind), this affected a meaningful share of T1w/T2w/FLAIR scans, heavily concentrated in a handful of source datasets whose particular header conventions land on the wrong side of that tolerance. Adds interfaces.regridmask(), inserted between brainmask() and biasfield(). It rewrites the mask onto the reference image's exact grid and transcribes the raw qform/sform NIfTI fields byte-for-byte from the reference image's own header — deliberately not recomputing either form from an affine matrix, since any fresh matrix -> quaternion decomposition (including nibabel's own set_qform()) can round differently than whatever decomposition originally produced the reference's qform. Verified on 30 real OpenMind scans that previously failed this way (across 3 source datasets, both FLAIR and T1w) plus 5 that already succeeded — 35/35 pass after this change. Co-Authored-By: Claude Sonnet 5 --- brainprep/interfaces/__init__.py | 2 + brainprep/interfaces/utils.py | 101 +++++++++++++++++++++++++++++++ brainprep/workflow/quasiraw.py | 14 +++++ 3 files changed, 117 insertions(+) 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,