Functional initialization & motion correction

pyproject.toml

@@ -8,8 +8,9 @@ readme = "README.md"
 requires-python = ">=3.12"
 dependencies = [
   "bids2table>=2.1.2",
+  "nibabel>=5.3.3",
   "niwrap>=0.8.3",
-  "niwrap-helper>=0.7.0"
+  "niwrap-helper>=0.7.0",
 ]
 
 [dependency-groups]
@@ -27,7 +28,7 @@ docs = ["pdoc>=15.0.0"]
 pythonpath = ["src"]
 testpaths = ["tests"]
 markers = [
-  "unit: Fast (<1s) unit tests (e.g. utilities, BIDS parsing, file operations",
+  "unit: Fast (<1s) unit tests (e.g. utilities, BIDS parsing, file operations)",
   "integration: Medium (1-5 min) integration tests with small/downsampled data",
   "slow: Long (30+ min) tests; skip on CI",
   "full_pipeline: End-to-end workflow tests; skip on CI"

src/rbc/core/functional/__init__.py

@@ -6,3 +6,13 @@
 fMRI contains motion artifacts, timing differences, and distortions that
 must be corrected before analysis.
 """
+
+from .initialization import scale_bold, truncate_trs
+from .motion import generate_motion_reference, motion_correction
+
+__all__ = [
+    "generate_motion_reference",
+    "motion_correction",
+    "scale_bold",
+    "truncate_trs",
+]

src/rbc/core/functional/initialization.py

@@ -0,0 +1,43 @@
+"""RBC functional initialization."""
+
+from pathlib import Path
+
+from niwrap import afni
+
+
+def truncate_trs(
+    in_file: Path, output_fname: str, start_tr: int
+) -> afni.V3dcalcOutputs:
+    """Remove first N TRs from BOLD timeseries using AFNI 3dcalc.
+
+    Args:
+        in_file: Path to input BOLD timeseries to be truncated.
+        output_fname: Name of output file.
+        start_tr: Number of TRs to remove from beginning.
+
+    Returns:
+        AFNI 3dcalc output object.
+    """
+    return afni.v_3dcalc(
+        dataset_a=afni.v_3dcalc_dataset_a_file(
+            file=in_file, selectors_=f"[{start_tr}..$]"
+        ),
+        expression="a",
+        prefix=output_fname,
+    )
+
+
+def scale_bold(in_file: Path, scale_factor: float = 0.1) -> afni.V3drefitOutputs:
+    """Scale BOLD voxel dimensions using AFNI 3drefit.
+
+    Args:
+        in_file: Path to input BOLD timeseries to be scaled.
+        scale_factor: Factor to scale voxel dimensions (default: 0.1 to divide by 10).
+
+    Returns:
+        AFNI 3drefit output object.
+    """
+    return afni.v_3drefit(
+        in_file=in_file,
+        xyzscale=scale_factor,
+    )

src/rbc/core/functional/motion.py

@@ -0,0 +1,76 @@
+"""RBC motion reference & correction."""
+
+from pathlib import Path
+from typing import NamedTuple, cast
+
+import nibabel as nib
+from niwrap import afni, fsl
+
+
+def generate_motion_reference(in_file: Path, output_fname: str) -> afni.V3dcalcOutputs:
+    """Creates reference volume for motion correction by extracting middle volume.
+
+    Args:
+        in_file: Path to input BOLD timeseries.
+        output_fname: Name of output file.
+
+    Returns:
+        AFNI 3dcalc output object.
+    """
+    img = cast("nib.Nifti1Image", nib.load(in_file))
+    total_vols = img.shape[3]
+
+    mid_vol = total_vols // 2
+
+    return afni.v_3dcalc(
+        dataset_a=afni.v_3dcalc_dataset_a_file(file=in_file, selectors_=f"[{mid_vol}]"),
+        expression="a",
+        prefix=output_fname,
+    )
+
+
+class MotionCorrectedOutputs(NamedTuple):
+    """NamedTuple for motion correction outputs."""
+
+    bold: Path
+    par: Path
+    rms_rel: Path
+    rms_abs: Path
+    mat_dir: Path
+
+
+def motion_correction(
+    in_file: Path, ref_file: Path, output_prefix: str
+) -> MotionCorrectedOutputs:
+    """Estimate and correct head motion using FSL mcflirt.
+
+    Args:
+        in_file: Path to input BOLD timeseries to correct.
+        ref_file: Path to reference volume for motion correction.
+        output_prefix: Prefix for output files.
+
+    Returns:
+        NamedTuple with paths to motion corrected outputs and matrices.
+    """
+    mc_result = fsl.mcflirt(
+        in_file=in_file,
+        ref_file=ref_file,
+        save_mats=True,
+        save_plots=True,
+        save_rmsrel=True,
+        save_rmsabs=True,
+        out_file=output_prefix,
+    )
+
+    motion_mat_dir = Path(mc_result.root) / f"{output_prefix}.mat"
+
+    if not motion_mat_dir.exists():
+        raise FileNotFoundError(f"Missing .mat directory at {motion_mat_dir}")
+
+    return MotionCorrectedOutputs(
+        bold=Path(mc_result.out_file),
+        par=Path(mc_result.par_file),
+        rms_rel=Path(mc_result.rmsrel_files),
+        rms_abs=Path(mc_result.rmsabs_files),
+        mat_dir=motion_mat_dir,
+    )

src/rbc/core/utils.py

@@ -45,10 +45,12 @@ def get_base_entities(
     return {k: v for k, v in file_entities.items() if k in base_entities}
 
 
-def rename(in_file: str | Path, new_name: str | Path) -> Path:
+def rename(in_file: str | Path, new_name: str) -> Path:
     """Rename a file, keeping it in the same directory."""
     in_file = Path(in_file)
-    new_path = in_file.with_name(Path(new_name).name)
-    if new_path.exists():
-        raise FileExistsError(f"Target file already exists: {new_path}")
-    return in_file.rename(new_path)
+    target = in_file.parent / Path(new_name).name
+    if in_file == target:
+        return in_file
+    if target.exists():
+        raise FileExistsError(f"Target file {target} already exists.")
+    return Path(shutil.move(str(in_file), str(target)))

src/rbc/workflows/functional.py

@@ -0,0 +1,78 @@
+"""Functional workflows."""
+
+import shutil
+from functools import partial
+from pathlib import Path
+from typing import Any
+
+from rbc.core.bids import Datatype, Suffix, bids_name, bids_path, parse_bids_name
+from rbc.core.common import reorient
+from rbc.core.functional import (
+    generate_motion_reference,
+    motion_correction,
+    truncate_trs,
+)
+
+
+def single_session(in_bold: Path, output_dir: Path, start_tr: int = 2) -> None:
+    """Workflow for preprocessing functional data.
+
+    Args:
+        in_bold: Input BOLD timeseries to process.
+        output_dir: Parent output directory to save data to.
+        start_tr: Number of initial TRs to remove (default: 2).
+    """
+    parsed = parse_bids_name(in_bold.name)
+    entities: dict[str, Any] = parsed.entities
+
+    bn = partial(bids_name, **entities)
+    bp = partial(bids_path, **entities, datatype=Datatype.FUNC)
+
+    reoriented = reorient(
+        in_file=in_bold,
+        output_fname=bn(desc="reoriented", suffix=Suffix.BOLD, extension=".nii.gz"),
+    )
+
+    truncated = truncate_trs(
+        in_file=reoriented.out_file,
+        output_fname=bn(suffix=Suffix.BOLD, extension=".nii.gz"),
+        start_tr=start_tr,
+    )
+
+    motion_ref = generate_motion_reference(
+        in_file=truncated.output_file,
+        output_fname=bn(suffix=Suffix.SBREF, extension=".nii.gz"),
+    )
+
+    motion_corrected = motion_correction(
+        in_file=truncated.output_file,
+        ref_file=motion_ref.output_file,
+        output_prefix=bn(suffix=Suffix.MOTION, extension=""),
+    )
+
+    func_out_dir = output_dir / bp(suffix=Suffix.BOLD, extension=".nii.gz").parent
+    func_out_dir.mkdir(parents=True, exist_ok=True)
+
+    # Save transform .mat directory
+    mat_target = func_out_dir / bn(desc="motion", suffix="mat", extension="")
+    shutil.copytree(motion_corrected.mat_dir, mat_target, dirs_exist_ok=True)
+
+    # Output files
+    outputs = [
+        (reoriented.out_file, "reorient", Suffix.BOLD, ".nii.gz"),
+        (truncated.output_file, "truncated", Suffix.BOLD, ".nii.gz"),
+        (motion_ref.output_file, None, Suffix.SBREF, ".nii.gz"),
+        (
+            motion_corrected.bold.with_suffix(".nii.gz"),
+            "motion",
+            Suffix.BOLD,
+            ".nii.gz",
+        ),
+        (motion_corrected.par, "motionParams", Suffix.MOTION, ".txt"),
+        (motion_corrected.rms_rel, "relsDisplacement", Suffix.MOTION, ".rms"),
+        (motion_corrected.rms_abs, "maxDisplacement", Suffix.MOTION, ".rms"),
+    ]
+
+    for out_file, desc, suffix, ext in outputs:
+        dest_path = func_out_dir / bn(desc=desc, suffix=suffix, extension=ext)
+        shutil.move(str(out_file), str(dest_path))

tests/integration/test_functional.py

@@ -0,0 +1,116 @@
+"""Integration tests for functional workflow."""
+
+from types import SimpleNamespace
+from typing import cast
+
+import nibabel as nib
+import pytest
+from niwrap import afni
+
+from rbc.core.common import reorient
+from rbc.core.functional import (
+    generate_motion_reference,
+    motion_correction,
+    truncate_trs,
+)
+
+
+def test_truncate_trs(test_subject: SimpleNamespace) -> None:
+    """Test truncating initial TRs from BOLD timeseries."""
+    original_count = cast("nib.Nifti1Image", nib.load(test_subject.bold)).shape[3]
+
+    start_tr = 4
+    reoriented = reorient(
+        in_file=test_subject.bold, output_fname="test_reoriented.nii.gz"
+    )
+    truncated_bold = truncate_trs(
+        in_file=reoriented.out_file,
+        output_fname="test_truncated.nii.gz",
+        start_tr=start_tr,
+    )
+    # Test truncated BOLD file exists & volume count is reduced
+    assert truncated_bold.output_file.exists()
+    new_shape = cast("nib.Nifti1Image", nib.load(truncated_bold.output_file)).shape
+    assert new_shape[3] == original_count - start_tr
+
+
+def test_truncate_to_min_volume(test_subject: SimpleNamespace) -> None:
+    """Test truncating to minimum volume count of 1."""
+    original_count = cast("nib.Nifti1Image", nib.load(test_subject.bold)).shape[3]
+
+    start_tr = original_count - 1
+    truncated_bold = truncate_trs(
+        in_file=test_subject.bold,
+        output_fname="test_truncated_min.nii.gz",
+        start_tr=start_tr,
+    )
+    # Test truncated BOLD file volume count is 1
+    new_shape = cast("nib.Nifti1Image", nib.load(truncated_bold.output_file)).shape
+    nvols = new_shape[3] if len(new_shape) > 3 else 1
+    assert nvols == 1
+
+
+def test_motion_reference_volume_count(test_subject: SimpleNamespace) -> None:
+    """Test motion reference volume count is 1."""
+    reference = generate_motion_reference(
+        in_file=test_subject.bold, output_fname="test_motion_ref.nii.gz"
+    )
+    # Test motion reference file has 1 volume
+    ref_shape = cast("nib.Nifti1Image", nib.load(reference.output_file)).shape
+    nvols = ref_shape[3] if len(ref_shape) > 3 else 1
+    assert nvols == 1
+
+
+@pytest.mark.slow
+def test_motion_correction_10vols(test_subject: SimpleNamespace) -> None:
+    """Test motion correction on 10 volumes of BOLD timeseries."""
+    reoriented = reorient(
+        in_file=test_subject.bold, output_fname="test_reoriented.nii.gz"
+    )
+    truncated_10 = afni.v_3dcalc(
+        dataset_a=afni.v_3dcalc_dataset_a_file(
+            file=reoriented.out_file, selectors_="[0..9]"
+        ),
+        expression="a",
+        prefix="test_10vols.nii.gz",
+    )
+    motion_reference = generate_motion_reference(
+        in_file=truncated_10.output_file, output_fname="test_ref_10v.nii.gz"
+    )
+    motion_corrected = motion_correction(
+        in_file=truncated_10.output_file,
+        ref_file=motion_reference.output_file,
+        output_prefix="test_mc_10v",
+    )
+    assert motion_corrected.bold.with_suffix(".nii.gz").exists()
+    mc_shape = cast(
+        "nib.Nifti1Image", nib.load(motion_corrected.bold.with_suffix(".nii.gz"))
+    ).shape
+    assert mc_shape[3] == 10
+
+    assert motion_corrected.par.exists()
+    par_data = motion_corrected.par.read_text().splitlines()
+    assert len(par_data) == 10
+
+
+@pytest.mark.slow
+def test_motion_correction(test_subject: SimpleNamespace) -> None:
+    """Test motion correction on full BOLD timeseries."""
+    reoriented = reorient(
+        in_file=test_subject.bold, output_fname="test_reoriented.nii.gz"
+    )
+    truncated = truncate_trs(
+        in_file=reoriented.out_file, output_fname="test_truncated.nii.gz", start_tr=2
+    )
+    motion_reference = generate_motion_reference(
+        in_file=truncated.output_file, output_fname="test_full_ref.nii.gz"
+    )
+    motion_corrected = motion_correction(
+        in_file=truncated.output_file,
+        ref_file=motion_reference.output_file,
+        output_prefix="test_full_mc",
+    )
+    # Test motion corrected BOLD files exists
+    assert motion_corrected.bold.with_suffix(".nii.gz").exists()
+    assert motion_corrected.par.exists()
+    assert motion_corrected.rms_rel.exists()