Resample BOLD with nitransforms (closes #112)

pyproject.toml

@@ -9,6 +9,7 @@ requires-python = ">=3.12"
 dependencies = [
   "bids2table>=2.1.2",
   "nibabel>=5.3.3",
+  "nitransforms>=24.0.0",
   "niwrap>=0.9.3",
   "numpy>=2.4.2",
   "polars>=1.38.1",

src/rbc/core/common.py

@@ -2,7 +2,7 @@
 
 Currently provides:
 - Deobliquing and RPI reorientation (initial preprocessing for T1w and BOLD).
-- 4D NIfTI splitting and merging utilities.
+- 4D NIfTI merging utility.
 """
 
 from __future__ import annotations
@@ -19,9 +19,8 @@
 
 from rbc.core.fileops import file_tmp_copy
 from rbc.core.nifti import strip_afni_volatile_metadata
-from rbc.core.niwrap import generate_exec_folder
 
-__all__ = ["deoblique_and_reorient", "merge_3d_to_4d", "split_4d"]
+__all__ = ["deoblique_and_reorient", "merge_3d_to_4d"]
 
 
 def deoblique_and_reorient(
@@ -55,29 +54,6 @@ def deoblique_and_reorient(
         )
 
 
-def split_4d(img_4d: Path) -> list[Path]:
-    """Split a 4D NIfTI timeseries into individual 3D volumes.
-
-    Volumes are written as uncompressed NIfTI (.nii) to avoid gzip
-    overhead on float intermediates that are read back immediately.
-
-    Args:
-        img_4d: Path to a 4D NIfTI image.
-
-    Returns:
-        Sorted list of paths to the individual 3D volume files.
-    """
-    img = nib.nifti1.load(img_4d)
-    volumes = nib.four_to_three(img)
-    out_dir = generate_exec_folder(suffix="split4d")
-    paths: list[Path] = []
-    for idx, vol in enumerate(volumes):
-        out_path = out_dir / f"vol_{idx:04d}.nii"
-        nib.save(vol, out_path)
-        paths.append(out_path)
-    return paths
-
-
 def merge_3d_to_4d(volumes: Sequence[Path], output: Path) -> Path:
     """Merge a sequence of 3D NIfTI volumes into a single 4D timeseries.
 

src/rbc/core/functional/distortion.py

@@ -23,7 +23,7 @@
 from niwrap import fsl
 from scipy.ndimage import binary_erosion
 
-from rbc.core.functional.resampling import merge_3d_to_4d
+from rbc.core.common import merge_3d_to_4d
 from rbc.core.niwrap import generate_exec_folder
 
 if TYPE_CHECKING:
@@ -187,7 +187,9 @@ def _shiftmap_to_itk_warp(
 
     FUGUE shift maps contain per-voxel shifts (in voxels) along the PE
     axis.  This converts to a 3-component mm displacement field in LPS
-    (ANTs/ITK convention).
+    (ANTs/ITK convention), written in the canonical ANTs 5D shape
+    ``(X, Y, Z, 1, 3)`` so it loads cleanly via
+    :class:`nitransforms.io.itk.ITKDisplacementsField`.
 
     Args:
         shiftmap: FUGUE shift map NIfTI (voxel units, single component).
@@ -202,17 +204,15 @@ def _shiftmap_to_itk_warp(
     voxel_sizes = np.array(shift_img.header.get_zooms()[:3])
 
     axis = _pe_axis_index(pe_direction)
+    warp = np.zeros((*shift_data.shape[:3], 1, 3), dtype=np.float32)
+    warp[..., 0, axis] = shift_data * voxel_sizes[axis]
 
-    # Build 3-component displacement field (mm, RAS)
-    warp = np.zeros((*shift_data.shape[:3], 3), dtype=np.float32)
-    warp[..., axis] = shift_data * voxel_sizes[axis]
-
-    # RAS → LPS for ANTs
-    warp[..., 0] *= -1
-    warp[..., 1] *= -1
+    # RAS -> LPS for ANTs
+    warp[..., 0, 0] *= -1
+    warp[..., 0, 1] *= -1
 
     header = shift_img.header.copy()
-    header.set_intent(1007)  # NIFTI_INTENT_VECTOR
+    header.set_intent("vector")
     header.set_data_dtype(np.float32)
 
     nib.save(nib.Nifti1Image(warp, shift_img.affine, header), output)
@@ -249,15 +249,15 @@ def _fieldmap_hz_to_itk_warp(
     axis = _pe_axis_index(pe_direction)
     sign = _pe_sign(pe_direction)
 
-    warp = np.zeros((*fmap_data.shape[:3], 3), dtype=np.float32)
-    warp[..., axis] = fmap_data * readout_time * voxel_sizes[axis] * sign
+    warp = np.zeros((*fmap_data.shape[:3], 1, 3), dtype=np.float32)
+    warp[..., 0, axis] = fmap_data * readout_time * voxel_sizes[axis] * sign
 
     # RAS -> LPS for ANTs
-    warp[..., 0] *= -1
-    warp[..., 1] *= -1
+    warp[..., 0, 0] *= -1
+    warp[..., 0, 1] *= -1
 
     header = fmap_img.header.copy()
-    header.set_intent(1007)
+    header.set_intent("vector")
     header.set_data_dtype(np.float32)
 
     nib.save(nib.Nifti1Image(warp, fmap_img.affine, header), output)