diff --git a/docs/data_dictionary.md b/docs/data_dictionary.md index 684e3d01..080764f1 100644 --- a/docs/data_dictionary.md +++ b/docs/data_dictionary.md @@ -56,7 +56,8 @@ Produced by `rbc functional`. These are functional MRI (BOLD) processing results | `*_space-MNI152NLin6Asym_desc-preproc_bold.nii.gz` | `bold` | BOLD timeseries resampled to MNI152NLin6Asym template space in a single interpolation step (before denoising) | ANTs resampling | 4D NIfTI | | `*_space-MNI152NLin6Asym_desc-bold_mask.nii.gz` | `mask` | Brain mask warped to template space at the BOLD resolution | ANTs resampling | 3D NIfTI, binary mask | | `*_space-MNI152NLin6Asym_reg-{regressor}_desc-preproc_bold.nii.gz` | `bold` | Denoised BOLD timeseries in template space after nuisance regression and bandpass filtering. `{regressor}` is `36parameter` or `aCompCor` | Nuisance regression | 4D NIfTI | -| `*_desc-{regressor}_regressors.1D` | `regressors` | Nuisance regressor matrix used for denoising. `{regressor}` is `36parameter` or `aCompCor` | Computed from motion parameters and tissue masks | Text, multi-column 1D file | +| `*_desc-{regressor}_regressors.1D` | `regressors` | Raw (unfiltered) nuisance regressor matrix. `{regressor}` is `36parameter` or `aCompCor`. Carried forward for longitudinal regression reuse | Computed from motion parameters and tissue masks | Text, multi-column 1D file | +| `*_desc-{regressor}Filtered_regressors.1D` | `regressors` | Bandpass-filtered nuisance regressor matrix matching what `3dTproject -bandpass` applied. For provenance only | FFT-based bandpass filter | Text, multi-column 1D file | --- @@ -133,3 +134,9 @@ Produced by the `rbc longitudinal` subcommand group (`template`, `anatomical`, ` | `*_space-longitudinal_desc-T1w_mask.nii.gz` | `mask` | Brain mask in longitudinal template space | ANTs registration to longitudinal template | 3D NIfTI, binary mask | | `*_from-T1w_to-longitudinal_mode-image_xfm.nii.gz` | `xfm` | Warp field mapping subject anatomy to the longitudinal template | ANTs registration | 3D NIfTI, displacement field | | `*_from-longitudinal_to-T1w_mode-image_xfm.nii.gz` | `xfm` | Inverse warp field mapping longitudinal template back to subject anatomy | ANTs registration | 3D NIfTI, displacement field | +| `*_space-longitudinal_sbref.nii.gz` | `sbref` | Motion reference volume warped to longitudinal template space | ANTs warping (composed BOLD-to-longitudinal) | 3D NIfTI | +| `*_space-longitudinal_desc-preproc_bold.nii.gz` | `bold` | Preprocessed BOLD warped to longitudinal template space | ANTs warping (composed BOLD-to-longitudinal) | 4D NIfTI | +| `*_space-longitudinal_desc-brain_mask.nii.gz` | `mask` | BOLD brain mask in longitudinal template space | ANTs warping (nearest-neighbor) | 3D NIfTI, binary mask | +| `*_from-bold_to-longitudinal_...xfm.nii.gz` | `xfm` | Composite BOLD-to-longitudinal-template warp field | ANTs compose transforms | 3D NIfTI, displacement field | +| `*_space-longitudinal_desc-regressed_reg-_bold.nii.gz` | `bold` | Nuisance-regressed BOLD (no bandpass) in longitudinal space, per regressor strategy | AFNI 3dTproject | 4D NIfTI | +| `*_space-longitudinal_desc-preproc_reg-_bold.nii.gz` | `bold` | Nuisance-regressed + bandpass-filtered BOLD in longitudinal space, per regressor strategy | AFNI 3dTproject -bandpass | 4D NIfTI | diff --git a/src/rbc/bids/functional.py b/src/rbc/bids/functional.py index 1991fd18..17ece5b2 100644 --- a/src/rbc/bids/functional.py +++ b/src/rbc/bids/functional.py @@ -158,6 +158,12 @@ def export_functional( desc=bids_safe_label(reg), extension=".1D", ) + func.save( + outputs.bpf_regressor_file[reg], + suffix="regressors", + desc=f"{bids_safe_label(reg)}Filtered", + extension=".1D", + ) mni = func.derive(space=TemplateSpace.MNI152NLIN6ASYM) for reg in regressors: diff --git a/src/rbc/bids/longitudinal/functional.py b/src/rbc/bids/longitudinal/functional.py index b02c4247..0ac9d457 100644 --- a/src/rbc/bids/longitudinal/functional.py +++ b/src/rbc/bids/longitudinal/functional.py @@ -7,6 +7,7 @@ from rbc.bids import Suffix, bids_safe_label if TYPE_CHECKING: + from collections.abc import Sequence from pathlib import Path import polars as pl @@ -22,7 +23,8 @@ def resolve_longitudinal_func( tpl_df: pl.DataFrame, *, ses: str, -) -> dict[str, Path | None]: + regressors: Sequence[str] = ("36-parameter",), +) -> dict[str, Path | dict[str, Path]]: """Resolve inputs for longitudinal functional processing. Args: @@ -31,10 +33,23 @@ def resolve_longitudinal_func( func_df: DataFrame of functional derivatives. tpl_df: DataFrame of longitudinal template files. ses: Session label (used for template xfm lookup). + regressors: Regressor strategy names to resolve raw regressor + files for. Returns: - Dict with keys matching ``longitudinal_process`` parameters. + Dict with keys matching ``longitudinal_process`` parameters, + including ``regressor_files`` keyed by strategy name. """ + regressor_files: dict[str, Path] = {} + for reg in regressors: + regressor_files[reg] = func_q.expect( + func_df, + suffix="regressors", + desc=bids_safe_label(reg), + extension=".1D", + without=["space"], + ) + return { "template": tpl_q.expect(tpl_df, suffix="T1w"), "anat_to_template_xfm": tpl_q.expect( @@ -54,18 +69,25 @@ def resolve_longitudinal_func( "bold": func_q.expect( func_df, suffix=Suffix.BOLD, desc="preproc", without=["space"] ), - "bold_mask": func_q.find( + "bold_mask": func_q.expect( func_df, suffix=Suffix.MASK, desc="brain", without=["space"] ), + "regressor_files": regressor_files, } -def export_longitudinal_func(fex: Bids, outputs: FunctionalLongOutputs) -> None: +def export_longitudinal_func( + fex: Bids, + outputs: FunctionalLongOutputs, + *, + regressors: Sequence[str], +) -> None: """Export longitudinal functional outputs. Args: fex: Bids builder with ``space="longitudinal"`` and identity entities. outputs: Results from the longitudinal functional workflow. + regressors: Regressor strategy names that were applied. """ fex.save(outputs.sbref, suffix=Suffix.SBREF) fex.save(outputs.bold, suffix=Suffix.BOLD, desc="preproc") @@ -75,5 +97,18 @@ def export_longitudinal_func(fex: Bids, outputs: FunctionalLongOutputs) -> None: desc="composite", extra={"from": "bold", "to": "longitudinal", "mode": "image"}, ) - if outputs.bold_mask: - fex.save(outputs.bold_mask, suffix=Suffix.MASK, desc="brain") + fex.save(outputs.bold_mask, suffix=Suffix.MASK, desc="brain") + + for reg in regressors: + fex.save( + outputs.regressed_bold[reg], + suffix=Suffix.BOLD, + desc="regressed", + extra={"reg": bids_safe_label(reg)}, + ) + fex.save( + outputs.cleaned_bold[reg], + suffix=Suffix.BOLD, + desc="preproc", + extra={"reg": bids_safe_label(reg)}, + ) diff --git a/src/rbc/cli/longitudinal/functional.py b/src/rbc/cli/longitudinal/functional.py index 5eea0cce..c1214557 100644 --- a/src/rbc/cli/longitudinal/functional.py +++ b/src/rbc/cli/longitudinal/functional.py @@ -3,7 +3,7 @@ from __future__ import annotations from dataclasses import dataclass -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Literal from rbc.cli.base import _validate_task from rbc.cli.longitudinal._base import LongitudinalBaseArgs, add_fs_license_argument @@ -20,6 +20,7 @@ class FunctionalLongArgs(LongitudinalBaseArgs): """Arguments for ``rbc longitudinal functional``.""" task: str | None + regressor: Sequence[Literal["36-parameter", "aCompCor"]] @classmethod def validate_namespace(cls, ns: argparse.Namespace) -> FunctionalLongArgs: @@ -28,6 +29,7 @@ def validate_namespace(cls, ns: argparse.Namespace) -> FunctionalLongArgs: return cls( **LongitudinalBaseArgs.validate_namespace(ns).__dict__, task=ns.task, + regressor=ns.regressor, ) @@ -41,6 +43,7 @@ def main(args: FunctionalLongArgs) -> int: session_label=args.session_label, task=args.task, ), + regressors=args.regressor, runner_config=RunnerConfig( runner=args.runner, verbose=bool(args.verbose), @@ -61,7 +64,7 @@ def register_command( parents=parents, description=( "Warp preprocessed BOLD derivatives into each subject's " - "longitudinal template space." + "longitudinal template space and re-run nuisance regression." ), help="Longitudinal functional stage", usage=( @@ -75,6 +78,13 @@ def register_command( default=None, help="Task label to filter BOLD runs (without 'task-' prefix).", ) + parser.add_argument( + "--regressor", + nargs="+", + choices=["36-parameter", "aCompCor"], + default=["36-parameter"], + help="Space-delimited nuisance regression method(s) to apply.", + ) parser.set_defaults( func=lambda args: main(FunctionalLongArgs.validate_namespace(args)) diff --git a/src/rbc/orchestration/__init__.py b/src/rbc/orchestration/__init__.py index 48b00121..a3a1904e 100644 --- a/src/rbc/orchestration/__init__.py +++ b/src/rbc/orchestration/__init__.py @@ -63,7 +63,7 @@ def apply(self, df: pl.DataFrame, *base_exprs: pl.Expr) -> pl.DataFrame: if len(self.session_label) > 0: exprs.append(pl.col("ses").is_in(self.session_label)) if self.task is not None: - exprs.append(pl.col("task") == self.task) + exprs.append(pl.col("task").is_null() | (pl.col("task") == self.task)) if not exprs: return df return df.filter(pl.all_horizontal(exprs)) diff --git a/src/rbc/orchestration/longitudinal/functional.py b/src/rbc/orchestration/longitudinal/functional.py index 724634f9..d5145be1 100644 --- a/src/rbc/orchestration/longitudinal/functional.py +++ b/src/rbc/orchestration/longitudinal/functional.py @@ -20,6 +20,7 @@ if TYPE_CHECKING: from collections.abc import Sequence from pathlib import Path + from typing import Literal import polars as pl @@ -34,6 +35,8 @@ def process_func( pipe_ctx: RunContext, func_df: pl.DataFrame, tpl_df: pl.DataFrame, + *, + regressors: Sequence[Literal["36-parameter", "aCompCor"]] = ("36-parameter",), ) -> None: """Handle functional longitudinal processing for one BOLD run. @@ -41,6 +44,7 @@ def process_func( pipe_ctx: RunContext bound to this subject/session. func_df: Functional derivative DataFrame for this run. tpl_df: Longitudinal template DataFrame. + regressors: Regressor strategies to apply in longitudinal space. """ row = func_df.filter(suffix=Suffix.BOLD).row(0, named=True) ents = extract_entities(row, ["task", "run"]) @@ -54,10 +58,11 @@ def process_func( func_df, tpl_df, ses=pipe_ctx.ses, # type: ignore[arg-type] + regressors=regressors, ) func_outputs = functional_longitudinal(**resolved) # type: ignore[arg-type] fex = func_q.derive(space="longitudinal") - export_longitudinal_func(fex, func_outputs) + export_longitudinal_func(fex, func_outputs, regressors=regressors) def run( @@ -65,6 +70,7 @@ def run( output_dir: Path, *, filters: Filters, + regressors: Sequence[Literal["36-parameter", "aCompCor"]] = ("36-parameter",), runner_config: RunnerConfig | None = None, ) -> None: """Run longitudinal functional processing for all matching subjects/sessions. @@ -75,6 +81,7 @@ def run( and longitudinal templates). output_dir: Output directory for derivatives. filters: Participant/session/task filters applied before grouping. + regressors: Regressor strategies to apply in longitudinal space. runner_config: Execution backend configuration. """ config = runner_config or RunnerConfig() @@ -91,7 +98,12 @@ def run( input_dirs, output_dir, filters=filters, verbose=verbose ): for func_df, _ in iter_session_files(session, groupby=FUNC_GROUP_ENTITIES): - process_func(pipe_ctx=pipe_ctx, func_df=func_df, tpl_df=tpl_df) + process_func( + pipe_ctx=pipe_ctx, + func_df=func_df, + tpl_df=tpl_df, + regressors=regressors, + ) pipe_ctx.ensure_dataset_description() _logger.info("RBC longitudinal functional workflow complete") diff --git a/src/rbc/workflows/functional.py b/src/rbc/workflows/functional.py index 21b779c8..11f44655 100644 --- a/src/rbc/workflows/functional.py +++ b/src/rbc/workflows/functional.py @@ -76,7 +76,12 @@ class FunctionalOutputs(NamedTuple): template_bold: BOLD resampled to template space. regressed_bold: Nuisance-regressed & non-bandpassed BOLD. cleaned_bold: Nuisance-regressed & bandpass-filtered BOLD. - regressor_file: Bandpass-filtered nuisance regressor ``.1D`` file. + regressor_file: Raw (unfiltered) nuisance regressor ``.1D`` file, + as computed from native-space BOLD. Carried forward so + longitudinal regression can reuse it without recomputation. + bpf_regressor_file: Bandpass-filtered nuisance regressor ``.1D`` + file, matching what ``3dTproject -bandpass`` actually applied. + For BIDS export only. template_brain_mask: Brain mask warped to template space. """ @@ -100,6 +105,7 @@ class FunctionalOutputs(NamedTuple): regressed_bold: dict[str, Path] cleaned_bold: dict[str, Path] regressor_file: dict[str, Path] + bpf_regressor_file: dict[str, Path] template_brain_mask: Path @@ -329,6 +335,7 @@ def single_session_preprocess( regression: dict[str, ApplyRegressionOutputs] = {} cleaned: dict[str, ApplyRegressionOutputs] = {} + raw_regressors: dict[str, Path] = {} filtered_regressors: dict[str, Path] = {} for regressor in regressor_set: # 15. Nuisance regression without bandpass (pre-bandpass residuals @@ -350,8 +357,11 @@ def single_session_preprocess( regressor_file=regressors[regressor].regressor_file, ) - # 17. Export bandpass-filtered regressors (matches what 3dTproject - # actually applied; raw regressors still in compute_regressors output) + # 17a. Carry raw (unfiltered) regressors forward for longitudinal reuse + raw_regressors[regressor] = regressors[regressor].regressor_file + + # 17b. Export bandpass-filtered regressors (matches what 3dTproject + # actually applied) filtered_regressors[regressor] = bandpass_regressor_file( regressors[regressor].regressor_file, tr=metadata.tr, @@ -379,6 +389,7 @@ def single_session_preprocess( template_bold=template_bold, regressed_bold={r: regression[r].regressed_bold for r in regressor_set}, cleaned_bold={r: cleaned[r].regressed_bold for r in regressor_set}, - regressor_file=filtered_regressors, + regressor_file=raw_regressors, + bpf_regressor_file=filtered_regressors, template_brain_mask=tmpl_brain, ) diff --git a/src/rbc/workflows/longitudinal/functional.py b/src/rbc/workflows/longitudinal/functional.py index da8dc0ce..fdf7b669 100644 --- a/src/rbc/workflows/longitudinal/functional.py +++ b/src/rbc/workflows/longitudinal/functional.py @@ -1,14 +1,17 @@ """Longitudinal functional processing workflow. Transforms preprocessed functional outputs to a pre-computed longitudinal -template space and returns all output paths as a +template space, then re-runs nuisance regression on the warped BOLD using +raw regressors from the cross-sectional run. Returns all output paths as a :class:`FunctionalLongOutputs` named tuple. """ from __future__ import annotations +import logging from typing import TYPE_CHECKING, NamedTuple +from rbc.core.functional import apply_regression, apply_regression_bandpass from rbc.core.longitudinal.transform import ( compose_transform, func_transform, @@ -18,6 +21,8 @@ if TYPE_CHECKING: from pathlib import Path +_logger = logging.getLogger("rbc") + class FunctionalLongOutputs(NamedTuple): """Outputs from the longitudinal functional preprocessing pipeline. @@ -26,14 +31,19 @@ class FunctionalLongOutputs(NamedTuple): bold_to_long_xfm: BOLD-to-longitudinal-template composite warp. sbref: Motion reference volume warped to longitudinal template space. bold: Preprocessed BOLD warped to longitudinal template space. - bold_mask: Brain mask warped to longitudinal template space, - or *None* if no mask was provided. + bold_mask: Brain mask warped to longitudinal template space. + regressed_bold: Per-regressor nuisance-regressed BOLD (no bandpass) + in longitudinal template space, keyed by strategy name. + cleaned_bold: Per-regressor nuisance-regressed + bandpass-filtered + BOLD in longitudinal template space, keyed by strategy name. """ bold_to_long_xfm: Path sbref: Path bold: Path - bold_mask: Path | None = None + bold_mask: Path + regressed_bold: dict[str, Path] + cleaned_bold: dict[str, Path] def longitudinal_process( @@ -43,13 +53,19 @@ def longitudinal_process( bold_to_anat_itk: Path, sbref: Path, bold: Path, - bold_mask: Path | None, + bold_mask: Path, + regressor_files: dict[str, Path], ) -> FunctionalLongOutputs: """Transform preprocessed functional outputs to longitudinal template space. - Assumes a longitudinal template has been generated, the subject-to-template - composite warp is available, and anatomical data has already been processed - to longitudinal template space. + After warping the BOLD timeseries, re-runs nuisance regression using the + raw (unfiltered) regressors produced by the cross-sectional pipeline. + No regressor recomputation is performed: the same regressor matrix is + applied in the new target space. + + Regression is applied for every strategy present in *regressor_files*. + The caller controls which strategies to run by passing only the desired + keys (the resolve layer filters to the requested ``--regressor`` set). Args: template: Longitudinal template image. @@ -57,11 +73,14 @@ def longitudinal_process( bold_to_anat_itk: BOLD-to-T1w affine in ITK format. sbref: Motion reference (single-band reference) volume. bold: Preprocessed bold image. - bold_mask: Bold brain mask, if available. + bold_mask: Bold brain mask. + regressor_files: Raw (unfiltered) regressor ``.1D`` files from + the cross-sectional run, keyed by strategy name. Regression + is applied for every key in this dict. Returns: - :class:`FunctionalLongOutputs` with all non-null inputs transformed to template - space. + :class:`FunctionalLongOutputs` with all inputs transformed to + longitudinal template space and per-regressor regression outputs. """ bold_to_tpl_xfm = compose_transform( ref=template, @@ -69,15 +88,36 @@ def longitudinal_process( anat_to_tpl_xfm=anat_to_template_xfm, ) + long_sbref = func_transform( + in_file=sbref, template=template, xfm=bold_to_tpl_xfm, strategy="single" + ) + long_bold = func_transform( + in_file=bold, template=template, xfm=bold_to_tpl_xfm, strategy="chunked" + ) + long_mask = mask_transform(mask=bold_mask, template=template, xfm=bold_to_tpl_xfm) + + regressed: dict[str, Path] = {} + cleaned: dict[str, Path] = {} + for reg, reg_file in regressor_files.items(): + _logger.info("Longitudinal %s nuisance regression (no bandpass)", reg) + regressed[reg] = apply_regression( + bold_file=long_bold, + brain_mask_file=long_mask, + regressor_file=reg_file, + ).regressed_bold + + _logger.info("Longitudinal %s nuisance regression + bandpass filtering", reg) + cleaned[reg] = apply_regression_bandpass( + bold_file=long_bold, + brain_mask_file=long_mask, + regressor_file=reg_file, + ).regressed_bold + return FunctionalLongOutputs( - sbref=func_transform( # 3D volume - in_file=sbref, template=template, xfm=bold_to_tpl_xfm, strategy="single" - ), - bold=func_transform( - in_file=bold, template=template, xfm=bold_to_tpl_xfm, strategy="chunked" - ), - bold_mask=mask_transform(mask=bold_mask, template=template, xfm=bold_to_tpl_xfm) - if bold_mask - else None, bold_to_long_xfm=bold_to_tpl_xfm, + sbref=long_sbref, + bold=long_bold, + bold_mask=long_mask, + regressed_bold=regressed, + cleaned_bold=cleaned, ) diff --git a/tests/integration/longitudinal/conftest.py b/tests/integration/longitudinal/conftest.py index daaf1bc4..89b53759 100644 --- a/tests/integration/longitudinal/conftest.py +++ b/tests/integration/longitudinal/conftest.py @@ -25,6 +25,7 @@ ) _SUB = "01" +_TASK = "fingerfootlips" def _rbc_exe() -> str: @@ -143,3 +144,73 @@ def longitudinal_template_output( ], ) return ds000114_anat_derivatives + + +@pytest.fixture(scope="session") +def ds000114_func_derivatives( + ds000114_dataset: Path, + ds000114_anat_derivatives: Path, + _runner: str, +) -> Path: + """Run ``rbc functional`` on ds000114 sub-01 ses-test. + + Produces cross-sectional functional derivatives (including raw + regressor ``.1D`` files) that the longitudinal functional stage + consumes. Writes into the same derivatives tree as the anatomical + stage so all outputs are visible to downstream fixtures. + + Only ses-test is processed (one session is sufficient to exercise + the longitudinal functional chain). + """ + _run_rbc( + [ + "functional", + str(ds000114_dataset), + str(ds000114_anat_derivatives), + "-o", + str(ds000114_anat_derivatives), + "--runner", + _runner, + "--participant-label", + _SUB, + "--session-label", + "test", + "--task", + _TASK, + ], + ) + return ds000114_anat_derivatives + + +@pytest.fixture(scope="session") +def longitudinal_func_output( + ds000114_func_derivatives: Path, + longitudinal_template_output: Path, # noqa: ARG001 — fixture dep for ordering + _runner: str, +) -> Path: + """Run ``rbc longitudinal functional`` on ds000114 sub-01 ses-test. + + Produces longitudinal functional derivatives (warped BOLD, + per-regressor regressed/cleaned BOLD) by consuming the + cross-sectional functional outputs and the longitudinal template. + Both ``ds000114_func_derivatives`` and ``longitudinal_template_output`` + write into the same derivatives directory, so all files are visible. + """ + _run_rbc( + [ + "longitudinal", + "functional", + str(ds000114_func_derivatives), + "-o", + str(ds000114_func_derivatives), + "--runner", + _runner, + "--participant-label", + _SUB, + "--session-label", + "test", + "--task", + _TASK, + ], + ) + return ds000114_func_derivatives diff --git a/tests/integration/longitudinal/test_regression_reuse.py b/tests/integration/longitudinal/test_regression_reuse.py new file mode 100644 index 00000000..2879a757 --- /dev/null +++ b/tests/integration/longitudinal/test_regression_reuse.py @@ -0,0 +1,106 @@ +"""Integration test for longitudinal functional regression reuse. + +Exercises the full ``rbc longitudinal functional`` CLI on ds000114, +which chains: composed BOLD-to-longitudinal warp, BOLD resampling, +and re-application of cross-sectional regressors in longitudinal space. + +Tier-2 integration test for Stage 5 of the longitudinal refactor +(tracker: #301). +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import nibabel as nib +import numpy as np +import pytest + +if TYPE_CHECKING: + from pathlib import Path + + +@pytest.mark.slow +def test_longitudinal_functional_produces_expected_tree( + longitudinal_func_output: Path, +) -> None: + """``rbc longitudinal functional`` writes per-regressor BOLD derivatives. + + Checks that the ``ses-test`` func directory under ``space-longitudinal`` + contains the expected files: warped BOLD, sbref, mask, composite xfm, + and per-regressor regressed + cleaned BOLD. + """ + func_dir = longitudinal_func_output / "sub-01" / "ses-test" / "func" + stem = "sub-01_ses-test_task-fingerfootlips" + + expected_fragments = [ + f"{stem}_space-longitudinal_sbref.nii.gz", + f"{stem}_space-longitudinal_desc-preproc_bold.nii.gz", + f"{stem}_space-longitudinal_desc-brain_mask.nii.gz", + f"{stem}_space-longitudinal_reg-36parameter_desc-regressed_bold.nii.gz", + f"{stem}_space-longitudinal_reg-36parameter_desc-preproc_bold.nii.gz", + ] + tree = sorted( + str(p.relative_to(longitudinal_func_output)) + for p in longitudinal_func_output.rglob("*") + if p.is_file() + ) + for name in expected_fragments: + assert (func_dir / name).is_file(), ( + f"Missing: {name}\n--- file tree ---\n" + "\n".join(tree) + ) + + +@pytest.mark.slow +def test_regressed_bold_non_degenerate( + longitudinal_func_output: Path, +) -> None: + """Regressed BOLD in longitudinal space has non-zero variance.""" + path = ( + longitudinal_func_output + / "sub-01" + / "ses-test" + / "func" + / "sub-01_ses-test_task-fingerfootlips" + "_space-longitudinal_reg-36parameter_desc-regressed_bold.nii.gz" + ) + img = nib.nifti1.load(path) + data = img.get_fdata() + assert data.var() > 0, "Regressed BOLD has zero variance" + + +@pytest.mark.slow +def test_cleaned_bold_non_degenerate( + longitudinal_func_output: Path, +) -> None: + """Cleaned (bandpassed) BOLD in longitudinal space has non-zero variance.""" + path = ( + longitudinal_func_output + / "sub-01" + / "ses-test" + / "func" + / "sub-01_ses-test_task-fingerfootlips" + "_space-longitudinal_reg-36parameter_desc-preproc_bold.nii.gz" + ) + img = nib.nifti1.load(path) + data = img.get_fdata() + assert data.var() > 0, "Cleaned BOLD has zero variance" + + +@pytest.mark.slow +def test_bold_mask_is_binary( + longitudinal_func_output: Path, +) -> None: + """Warped bold mask in longitudinal space is binary.""" + path = ( + longitudinal_func_output + / "sub-01" + / "ses-test" + / "func" + / "sub-01_ses-test_task-fingerfootlips" + "_space-longitudinal_desc-brain_mask.nii.gz" + ) + img = nib.nifti1.load(path) + data = img.get_fdata() + unique = np.unique(data) + assert set(unique).issubset({0, 1}), f"Mask has non-binary values: {unique}" diff --git a/tests/unit/bids/test_exports.py b/tests/unit/bids/test_exports.py index 1c96faab..33ec4852 100644 --- a/tests/unit/bids/test_exports.py +++ b/tests/unit/bids/test_exports.py @@ -74,6 +74,7 @@ def _make_func_outputs(w: Path, regressors: list[str]) -> FunctionalOutputs: regressed_bold={r: _dummy(w, f"regressed_{r}.nii.gz") for r in regressors}, cleaned_bold={r: _dummy(w, f"cleaned_{r}.nii.gz") for r in regressors}, regressor_file={r: _dummy(w, f"regressors_{r}.1D") for r in regressors}, + bpf_regressor_file={r: _dummy(w, f"regressors_bpf_{r}.1D") for r in regressors}, template_brain_mask=_dummy(w, "template_mask.nii.gz"), ) @@ -195,26 +196,27 @@ def test_file_count_single_regressor( ) -> None: """Correct file count with one regressor. - 8 native-space fixed + 1 regressor file + 2 per-regressor MNI - + 2 fixed MNI = 13. + 8 native-space fixed + 2 regressor files (raw + filtered) + + 2 per-regressor MNI + 2 fixed MNI = 14. """ outputs = _make_func_outputs(workdir, ["36-parameter"]) export_functional(func_bids, outputs, regressors=["36-parameter"]) saved = list(pipe_ctx.output_dir.rglob("*.*")) - assert len(saved) == 13 + assert len(saved) == 14 def test_file_count_two_regressors( self, func_bids: Bids, workdir: Path, pipe_ctx: RunContext ) -> None: """Correct file count with two regressors. - 8 fixed + 2 regressor files + 4 per-regressor MNI + 2 fixed MNI = 16. + 8 fixed + 4 regressor files (2x raw + 2x filtered) + + 4 per-regressor MNI + 2 fixed MNI = 18. """ regs = ["36-parameter", "aCompCor"] outputs = _make_func_outputs(workdir, regs) export_functional(func_bids, outputs, regressors=regs) saved = list(pipe_ctx.output_dir.rglob("*.*")) - assert len(saved) == 16 + assert len(saved) == 18 # --------------------------------------------------------------------------- diff --git a/tests/unit/bids/test_longitudinal_functional.py b/tests/unit/bids/test_longitudinal_functional.py new file mode 100644 index 00000000..98a17d41 --- /dev/null +++ b/tests/unit/bids/test_longitudinal_functional.py @@ -0,0 +1,417 @@ +"""Unit tests for ``rbc.bids.longitudinal.functional``.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import polars as pl +import pytest + +from rbc.bids.longitudinal.functional import ( + export_longitudinal_func, + resolve_longitudinal_func, +) +from rbc.context import RunContext +from rbc.workflows.longitudinal.functional import FunctionalLongOutputs + +if TYPE_CHECKING: + from pathlib import Path + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _func_row( + *, + sub: str, + ses: str, + suffix: str, + desc: str | None = None, + ext: str = ".nii.gz", + space: str | None = None, + extra: list[dict[str, str]] | None = None, +) -> dict[str, object]: + """Build a single BIDS-like row for functional derivatives.""" + desc_part = f"_desc-{desc}" if desc else "" + space_part = f"_space-{space}" if space else "" + path = ( + f"sub-{sub}/ses-{ses}/func/" + f"sub-{sub}_ses-{ses}{space_part}{desc_part}_{suffix}{ext}" + ) + return { + "datatype": "func", + "suffix": suffix, + "ext": ext, + "sub": sub, + "ses": ses, + "space": space, + "desc": desc, + "root": "/data", + "path": path, + "extra_entities": extra or [], + } + + +def _anat_row( + *, + sub: str, + ses: str, + suffix: str, + desc: str | None = None, + ext: str = ".nii.gz", + extra: list[dict[str, str]] | None = None, +) -> dict[str, object]: + """Build a single BIDS-like row for anatomical/template derivatives.""" + desc_part = f"_desc-{desc}" if desc else "" + path = f"sub-{sub}/ses-{ses}/anat/sub-{sub}_ses-{ses}{desc_part}_{suffix}{ext}" + return { + "datatype": "anat", + "suffix": suffix, + "ext": ext, + "sub": sub, + "ses": ses, + "space": None, + "desc": desc, + "root": "/data", + "path": path, + "extra_entities": extra or [], + } + + +def _df(*rows: dict[str, object]) -> pl.DataFrame: + return pl.DataFrame(list(rows)) + + +def _make_long_outputs(workdir: Path) -> FunctionalLongOutputs: + """Build a populated FunctionalLongOutputs pointing at dummy files.""" + + def _dummy(name: str) -> Path: + p = workdir / name + p.write_bytes(b"\x00") + return p + + return FunctionalLongOutputs( + bold_to_long_xfm=_dummy("bold_to_long_xfm.nii.gz"), + sbref=_dummy("sbref.nii.gz"), + bold=_dummy("bold.nii.gz"), + bold_mask=_dummy("bold_mask.nii.gz"), + regressed_bold={ + "36-parameter": _dummy("regressed_36p.nii.gz"), + "aCompCor": _dummy("regressed_acompcor.nii.gz"), + }, + cleaned_bold={ + "36-parameter": _dummy("cleaned_36p.nii.gz"), + "aCompCor": _dummy("cleaned_acompcor.nii.gz"), + }, + ) + + +# --------------------------------------------------------------------------- +# resolve_longitudinal_func +# --------------------------------------------------------------------------- + + +class TestResolveLongitudinalFunc: + """Tests for :func:`resolve_longitudinal_func`.""" + + def test_resolves_single_regressor(self, tmp_path: Path) -> None: + """Single regressor resolves raw regressor file from derivatives.""" + func_df = _df( + _func_row(sub="01", ses="baseline", suffix="sbref"), + _func_row(sub="01", ses="baseline", suffix="bold", desc="preproc"), + _func_row(sub="01", ses="baseline", suffix="mask", desc="brain"), + _func_row( + sub="01", + ses="baseline", + suffix="xfm", + desc="linearITK", + ext=".txt", + extra=[ + {"key": "from", "value": "bold"}, + {"key": "to", "value": "T1w"}, + {"key": "mode", "value": "image"}, + ], + ), + _func_row( + sub="01", + ses="baseline", + suffix="regressors", + desc="36parameter", + ext=".1D", + ), + ) + tpl_df = _df( + _anat_row(sub="01", ses="longitudinal", suffix="T1w"), + _anat_row( + sub="01", + ses="longitudinal", + suffix="xfm", + ext=".txt", + extra=[ + {"key": "from", "value": "baseline"}, + {"key": "to", "value": "longitudinal"}, + ], + ), + ) + + ctx = RunContext(sub="01", ses="baseline", output_dir=tmp_path) + func_q = ctx.bids(datatype="func") + tpl_q = ctx.bids(datatype="anat").derive(ses="longitudinal") + + resolved = resolve_longitudinal_func( + func_q, + tpl_q, + func_df, + tpl_df, + ses="baseline", + regressors=["36-parameter"], + ) + + assert set(resolved) == { + "template", + "anat_to_template_xfm", + "bold_to_anat_itk", + "sbref", + "bold", + "bold_mask", + "regressor_files", + } + assert isinstance(resolved["regressor_files"], dict) + assert "36-parameter" in resolved["regressor_files"] + assert str(resolved["regressor_files"]["36-parameter"]).endswith( + "regressors.1D" + ) + + def test_resolves_multiple_regressors(self, tmp_path: Path) -> None: + """Multiple regressors each get their own raw regressor file resolved.""" + func_df = _df( + _func_row(sub="01", ses="baseline", suffix="sbref"), + _func_row(sub="01", ses="baseline", suffix="bold", desc="preproc"), + _func_row(sub="01", ses="baseline", suffix="mask", desc="brain"), + _func_row( + sub="01", + ses="baseline", + suffix="xfm", + desc="linearITK", + ext=".txt", + extra=[ + {"key": "from", "value": "bold"}, + {"key": "to", "value": "T1w"}, + {"key": "mode", "value": "image"}, + ], + ), + _func_row( + sub="01", + ses="baseline", + suffix="regressors", + desc="36parameter", + ext=".1D", + ), + _func_row( + sub="01", + ses="baseline", + suffix="regressors", + desc="aCompCor", + ext=".1D", + ), + ) + tpl_df = _df( + _anat_row(sub="01", ses="longitudinal", suffix="T1w"), + _anat_row( + sub="01", + ses="longitudinal", + suffix="xfm", + ext=".txt", + extra=[ + {"key": "from", "value": "baseline"}, + {"key": "to", "value": "longitudinal"}, + ], + ), + ) + + ctx = RunContext(sub="01", ses="baseline", output_dir=tmp_path) + func_q = ctx.bids(datatype="func") + tpl_q = ctx.bids(datatype="anat").derive(ses="longitudinal") + + resolved = resolve_longitudinal_func( + func_q, + tpl_q, + func_df, + tpl_df, + ses="baseline", + regressors=["36-parameter", "aCompCor"], + ) + + reg_files = resolved["regressor_files"] + assert isinstance(reg_files, dict) + assert set(reg_files) == {"36-parameter", "aCompCor"} + + def test_missing_regressor_raises(self, tmp_path: Path) -> None: + """Requesting a regressor not present in derivatives raises.""" + func_df = _df( + _func_row(sub="01", ses="baseline", suffix="sbref"), + _func_row(sub="01", ses="baseline", suffix="bold", desc="preproc"), + _func_row(sub="01", ses="baseline", suffix="mask", desc="brain"), + _func_row( + sub="01", + ses="baseline", + suffix="xfm", + desc="linearITK", + ext=".txt", + extra=[ + {"key": "from", "value": "bold"}, + {"key": "to", "value": "T1w"}, + {"key": "mode", "value": "image"}, + ], + ), + _func_row( + sub="01", + ses="baseline", + suffix="regressors", + desc="36parameter", + ext=".1D", + ), + ) + tpl_df = _df( + _anat_row(sub="01", ses="longitudinal", suffix="T1w"), + _anat_row( + sub="01", + ses="longitudinal", + suffix="xfm", + ext=".txt", + extra=[ + {"key": "from", "value": "baseline"}, + {"key": "to", "value": "longitudinal"}, + ], + ), + ) + + ctx = RunContext(sub="01", ses="baseline", output_dir=tmp_path) + func_q = ctx.bids(datatype="func") + tpl_q = ctx.bids(datatype="anat").derive(ses="longitudinal") + + with pytest.raises(FileNotFoundError): + resolve_longitudinal_func( + func_q, + tpl_q, + func_df, + tpl_df, + ses="baseline", + regressors=["aCompCor"], + ) + + def test_bold_mask_mandatory(self, tmp_path: Path) -> None: + """bold_mask is now resolved with expect(), so missing raises.""" + func_df = _df( + _func_row(sub="01", ses="baseline", suffix="sbref"), + _func_row(sub="01", ses="baseline", suffix="bold", desc="preproc"), + _func_row( + sub="01", + ses="baseline", + suffix="xfm", + desc="linearITK", + ext=".txt", + extra=[ + {"key": "from", "value": "bold"}, + {"key": "to", "value": "T1w"}, + {"key": "mode", "value": "image"}, + ], + ), + _func_row( + sub="01", + ses="baseline", + suffix="regressors", + desc="36parameter", + ext=".1D", + ), + ) + tpl_df = _df( + _anat_row(sub="01", ses="longitudinal", suffix="T1w"), + _anat_row( + sub="01", + ses="longitudinal", + suffix="xfm", + ext=".txt", + extra=[ + {"key": "from", "value": "baseline"}, + {"key": "to", "value": "longitudinal"}, + ], + ), + ) + + ctx = RunContext(sub="01", ses="baseline", output_dir=tmp_path) + func_q = ctx.bids(datatype="func") + tpl_q = ctx.bids(datatype="anat").derive(ses="longitudinal") + + with pytest.raises(FileNotFoundError): + resolve_longitudinal_func( + func_q, + tpl_q, + func_df, + tpl_df, + ses="baseline", + regressors=["36-parameter"], + ) + + +# --------------------------------------------------------------------------- +# export_longitudinal_func +# --------------------------------------------------------------------------- + + +class TestExportLongitudinalFunc: + """Tests for :func:`export_longitudinal_func`.""" + + def test_writes_expected_files(self, tmp_path: Path) -> None: + """Exports BOLD, sbref, mask, xfm, plus per-regressor regressed/cleaned.""" + workdir = tmp_path / "work" + workdir.mkdir() + out_dir = tmp_path / "out" + + ctx = RunContext(sub="01", ses="baseline", output_dir=out_dir) + fex = ctx.bids(datatype="func").derive(space="longitudinal") + + outputs = _make_long_outputs(workdir) + export_longitudinal_func(fex, outputs, regressors=["36-parameter", "aCompCor"]) + + saved = sorted(p.name for p in out_dir.rglob("*.*")) + # 4 fixed (sbref, bold, xfm, mask) + # + 2 regressors x 2 (regressed + cleaned) = 4 per-regressor + # = 8 total + assert len(saved) == 8 + + def test_regressor_entity_in_filenames(self, tmp_path: Path) -> None: + """Per-regressor outputs carry the reg- entity.""" + workdir = tmp_path / "work" + workdir.mkdir() + out_dir = tmp_path / "out" + + ctx = RunContext(sub="01", ses="baseline", output_dir=out_dir) + fex = ctx.bids(datatype="func").derive(space="longitudinal") + + outputs = _make_long_outputs(workdir) + export_longitudinal_func(fex, outputs, regressors=["36-parameter", "aCompCor"]) + + names = [p.name for p in out_dir.rglob("*.*")] + reg_files = [n for n in names if "reg-" in n] + assert len(reg_files) == 4 + + assert any("reg-36parameter" in n for n in reg_files) + assert any("reg-aCompCor" in n for n in reg_files) + + def test_single_regressor_file_count(self, tmp_path: Path) -> None: + """Single regressor produces 6 files: 4 fixed + 2 per-regressor.""" + workdir = tmp_path / "work" + workdir.mkdir() + out_dir = tmp_path / "out" + + ctx = RunContext(sub="01", ses="baseline", output_dir=out_dir) + fex = ctx.bids(datatype="func").derive(space="longitudinal") + + outputs = _make_long_outputs(workdir) + export_longitudinal_func(fex, outputs, regressors=["36-parameter"]) + + saved = list(out_dir.rglob("*.*")) + assert len(saved) == 6 diff --git a/tests/unit/cli/test_longitudinal.py b/tests/unit/cli/test_longitudinal.py index 4608cd01..66f97918 100644 --- a/tests/unit/cli/test_longitudinal.py +++ b/tests/unit/cli/test_longitudinal.py @@ -80,15 +80,24 @@ class TestFunctionalLongArgs: def test_valid_task(self, base_ns: argparse.Namespace) -> None: """Alphanumeric task labels pass validation.""" base_ns.task = "rest" + base_ns.regressor = ["36-parameter"] args = FunctionalLongArgs.validate_namespace(base_ns) assert args.task == "rest" def test_invalid_task_rejected(self, base_ns: argparse.Namespace) -> None: """Task labels with special characters are rejected.""" base_ns.task = "rest/invalid" + base_ns.regressor = ["36-parameter"] with pytest.raises(ValueError, match="Task"): FunctionalLongArgs.validate_namespace(base_ns) + def test_regressor_preserved(self, base_ns: argparse.Namespace) -> None: + """Regressor choices round-trip through validation.""" + base_ns.task = None + base_ns.regressor = ["36-parameter", "aCompCor"] + args = FunctionalLongArgs.validate_namespace(base_ns) + assert list(args.regressor) == ["36-parameter", "aCompCor"] + class TestMetricsLongArgs: """Tests for the metrics longitudinal subcommand validator.""" diff --git a/tests/unit/orchestration/test_filters.py b/tests/unit/orchestration/test_filters.py index dd94e708..ca711f89 100644 --- a/tests/unit/orchestration/test_filters.py +++ b/tests/unit/orchestration/test_filters.py @@ -70,10 +70,13 @@ def test_session_filter(self, bids_df: pl.DataFrame) -> None: result = Filters(session_label=["baseline"]).apply(bids_df) assert set(result["ses"].unique().to_list()) == {"baseline"} - def test_task_filter(self, bids_df: pl.DataFrame) -> None: - """Task filter keeps only matching tasks (nulls excluded).""" + def test_task_filter_preserves_anat(self, bids_df: pl.DataFrame) -> None: + """Task filter keeps matching tasks AND rows with null task (anat).""" result = Filters(task="rest").apply(bids_df) - assert all(t == "rest" for t in result["task"].to_list()) + tasks = result["task"].to_list() + assert all(t in ("rest", None) for t in tasks) + # Anat rows (task=null) must survive + assert result.filter(pl.col("datatype") == "anat").height > 0 def test_combined_filters(self, bids_df: pl.DataFrame) -> None: """Participant + session + task filters compose correctly.""" @@ -82,7 +85,8 @@ def test_combined_filters(self, bids_df: pl.DataFrame) -> None: session_label=["baseline"], task="rest", ).apply(bids_df) - assert len(result) == 3 # raw bold + MNI preproc + T1w preproc + # anat T1w (task=null) + raw bold + MNI preproc + T1w preproc = 4 + assert len(result) == 4 assert all(s == "01" for s in result["sub"].to_list()) assert all(s == "baseline" for s in result["ses"].to_list()) diff --git a/tests/unit/orchestration/test_longitudinal.py b/tests/unit/orchestration/test_longitudinal.py index 43b321f8..4195992f 100644 --- a/tests/unit/orchestration/test_longitudinal.py +++ b/tests/unit/orchestration/test_longitudinal.py @@ -67,6 +67,8 @@ def _mock_func_outputs(*, with_bold_mask: bool = True) -> Mock: m.bold = fake / "bold.nii.gz" m.bold_to_long_xfm = fake / "bold_to_long_xfm.nii.gz" m.bold_mask = (fake / "bold_mask.nii.gz") if with_bold_mask else None + m.regressed_bold = {"36-parameter": fake / "regressed_36.nii.gz"} + m.cleaned_bold = {"36-parameter": fake / "cleaned_36.nii.gz"} return m @@ -382,24 +384,23 @@ def test_missing_required_file_raises( ): process_func(pipe_ctx=pipe_ctx, func_df=func_df, tpl_df=tpl_df) - def test_optional_bold_mask_file_not_found( + def test_missing_bold_mask_raises( self, func_df: pl.DataFrame, tpl_df: pl.DataFrame, tmp_path: Path ) -> None: - """Optional bold_mask not found is caught; 3 exports emitted.""" + """bold_mask is now mandatory; missing file raises FileNotFoundError.""" pipe_ctx = RunContext(sub="01", ses="baseline", output_dir=tmp_path) with ( patch( "rbc.orchestration.longitudinal.functional.functional_longitudinal", - return_value=_mock_func_outputs(with_bold_mask=False), + return_value=_mock_func_outputs(), ), patch( "rbc.bids.query.find_file", side_effect=_none_for(suffix="mask", desc="brain"), ), - patch("rbc.bids.builder.shutil.copy2") as mock_copy, + pytest.raises(FileNotFoundError), ): process_func(pipe_ctx=pipe_ctx, func_df=func_df, tpl_df=tpl_df) - assert mock_copy.call_count == 3 class TestLongitudinalAnatomicalRun: