diff --git a/.maint/CONTRIBUTORS.md b/.maint/CONTRIBUTORS.md
index 483f44c56..829ff6ef6 100644
--- a/.maint/CONTRIBUTORS.md
+++ b/.maint/CONTRIBUTORS.md
@@ -60,6 +60,7 @@ Before every release, unlisted contributors will be invited again to add their n
| Schaefer | Theo A.J. | | 0000-0003-4102-559X | Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany |
| Seeley | Saren | @sarenseeley | 0000-0002-9493-8512 | Department of Psychiatry, Icahn School of Medicine at Mount Sinai |
| Sitek | Kevin R. | | 0000-0002-2172-5786 | Speech & Hearing Bioscience & Technology Program, Harvard University |
+| Smith | David V. | @dvsmith.bsky.social | 0000-0001-5754-9633 | Department of Psychology and Neuroscience, Temple University |
| Smith | Robert E. | | 0000-0003-3636-4642 | Florey Institute of Neuroscience and Mental Health |
| Sneve | Markus H. | | 0000-0001-7644-7915 | Center for Lifespan Changes in Brain and Cognition, University of Oslo |
| Stojić | Hrvoje | | 0000-0002-9699-9052 | Max Planck UCL Centre for Computational Psychiatry and Ageing Research, University College London |
diff --git a/Dockerfile b/Dockerfile
index f2f17011d..af8730f37 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -23,6 +23,7 @@
# SOFTWARE.
ARG BASE_IMAGE=ghcr.io/nipreps/fmriprep-base:20251006
+ARG PIXI_LOCK_FLAGS=--frozen
#
# Build pixi environment
@@ -39,6 +40,7 @@ ARG BASE_IMAGE=ghcr.io/nipreps/fmriprep-base:20251006
# - ...
#
FROM ghcr.io/prefix-dev/pixi:0.53.0 AS build
+ARG PIXI_LOCK_FLAGS=--frozen
RUN apt-get update && \
apt-get install -y --no-install-recommends \
ca-certificates \
@@ -51,7 +53,7 @@ RUN pixi config set --global run-post-link-scripts insecure
RUN mkdir /app
COPY pixi.lock pyproject.toml /app
WORKDIR /app
-RUN --mount=type=cache,target=/root/.cache/rattler pixi install -e fmriprep -e test --frozen --skip fmriprep
+RUN --mount=type=cache,target=/root/.cache/rattler pixi install -e fmriprep -e test ${PIXI_LOCK_FLAGS} --skip fmriprep
RUN --mount=type=cache,target=/root/.npm pixi run --as-is -e fmriprep npm install -g svgo@^3.2.0 bids-validator@1.14.10
# Note that PATH gets hard-coded. Remove it and re-apply in final image
RUN pixi shell-hook -e fmriprep --as-is | grep -v PATH > /shell-hook.sh
@@ -59,7 +61,7 @@ RUN pixi shell-hook -e test --as-is | grep -v PATH > /test-shell-hook.sh
# Finally, install the package
COPY . /app
-RUN --mount=type=cache,target=/root/.cache/rattler pixi install -e fmriprep -e test --frozen
+RUN --mount=type=cache,target=/root/.cache/rattler pixi install -e fmriprep -e test ${PIXI_LOCK_FLAGS}
#
# Pre-fetch templates
diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py
index fbbdb51c6..47b8cabd7 100644
--- a/fmriprep/cli/parser.py
+++ b/fmriprep/cli/parser.py
@@ -103,6 +103,13 @@ def _min_one(value, parser):
raise parser.error("Argument can't be less than one.")
return value
+ def _min_zero(value, parser):
+ """Ensure an argument is not lower than 0."""
+ value = int(value)
+ if value < 0:
+ raise parser.error("Argument can't be negative.")
+ return value
+
def _to_gb(value):
scale = {'G': 1, 'T': 10**3, 'M': 1e-3, 'K': 1e-6, 'B': 1e-9}
digits = ''.join([c for c in value if c.isdigit()])
@@ -177,6 +184,7 @@ def _fallback_trt(value, parser):
PathExists = partial(_path_exists, parser=parser)
IsFile = partial(_is_file, parser=parser)
PositiveInt = partial(_min_one, parser=parser)
+ NonnegativeInt = partial(_min_zero, parser=parser)
BIDSFilter = partial(_bids_filter, parser=parser)
SliceTimeRef = partial(_slice_time_ref, parser=parser)
FallbackTRT = partial(_fallback_trt, parser=parser)
@@ -485,6 +493,27 @@ def _fallback_trt(value, parser):
'It is faster and less memory intensive, but may be less accurate.'
),
)
+ g_conf.add_argument(
+ '--me-use-warpkit',
+ action='store_true',
+ default=False,
+ help=(
+ 'Use warpkit MEDIC for susceptibility distortion correction of compatible '
+ 'multi-echo BOLD runs. Requires phase companions for each echo and a '
+ 'warpkit installation (for example, `fmriprep[warpkit]`) on Python 3.11+.'
+ ),
+ )
+ g_conf.add_argument(
+ '--me-warpkit-noise-frames',
+ action='store',
+ type=NonnegativeInt,
+ default=0,
+ help=(
+ 'Number of trailing non-imaging/noise frames to trim from each '
+ 'multi-echo magnitude and phase file before running warpkit MEDIC. '
+ 'Only applies when --me-use-warpkit is enabled.'
+ ),
+ )
g_outputs = parser.add_argument_group('Options for modulating outputs')
g_outputs.add_argument(
diff --git a/fmriprep/cli/tests/test_parser.py b/fmriprep/cli/tests/test_parser.py
index 8da654a15..09c7376e9 100644
--- a/fmriprep/cli/tests/test_parser.py
+++ b/fmriprep/cli/tests/test_parser.py
@@ -199,6 +199,27 @@ def test_slice_time_ref(tmp_path, st_ref):
_reset_config()
+def test_me_use_warpkit(tmp_path):
+ bids_path = tmp_path / 'data'
+ out_path = tmp_path / 'out'
+ args = [
+ str(bids_path),
+ str(out_path),
+ 'participant',
+ '--me-use-warpkit',
+ '--me-warpkit-noise-frames',
+ '3',
+ ]
+ bids_path.mkdir()
+
+ parser = _build_parser()
+ opts = parser.parse_args(args)
+
+ assert opts.me_use_warpkit is True
+ assert opts.me_warpkit_noise_frames == 3
+ _reset_config()
+
+
@pytest.mark.parametrize(
('args', 'expectation'),
[
diff --git a/fmriprep/config.py b/fmriprep/config.py
index c91951285..21030cef4 100644
--- a/fmriprep/config.py
+++ b/fmriprep/config.py
@@ -642,6 +642,10 @@ class workflow(_Config):
in the absence of any alternatives."""
me_t2s_fit_method = 'curvefit'
"""The method by which to estimate T2*/S0 for multi-echo data"""
+ me_use_warpkit = False
+ """Run warpkit's MEDIC workflow for multi-echo susceptibility correction."""
+ me_warpkit_noise_frames = 0
+ """Number of trailing noise frames to trim before running warpkit MEDIC."""
class loggers:
diff --git a/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib
index 50289ce1f..de3f28b41 100644
--- a/fmriprep/data/boilerplate.bib
+++ b/fmriprep/data/boilerplate.bib
@@ -337,6 +337,16 @@ @article{posse_t2s
year = 1999
}
+@article{van2023medic,
+ author = {Van, Andrew N. and Montez, David F. and Laumann, Timothy O. and Suljic, Vahdeta and Madison, Thomas and Baden, Noah J. and Ramirez-Perez, Nadeshka and Scheidter, Kristen M. and Monk, Julia S. and Whiting, Forrest I. and others},
+ title = {Framewise multi-echo distortion correction for superior functional {MRI}},
+ journal = {bioRxiv},
+ year = {2023},
+ doi = {10.1101/2023.11.28.568744},
+ url = {https://doi.org/10.1101/2023.11.28.568744},
+ publisher = {Cold Spring Harbor Laboratory}
+}
+
@article{topup,
author = {Jesper L.R. Andersson and Stefan Skare and John Ashburner},
title = {How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging},
diff --git a/fmriprep/data/reports-spec-func.yml b/fmriprep/data/reports-spec-func.yml
index dbe5b830e..fe71ddf75 100644
--- a/fmriprep/data/reports-spec-func.yml
+++ b/fmriprep/data/reports-spec-func.yml
@@ -59,17 +59,12 @@ sections:
static: false
subtitle: Alignment between the anatomical reference of the fieldmap and the BOLD reference
- bids: {datatype: figures, desc: fieldmap, suffix: bold}
- caption: Estimated fieldmap, as reconstructed on the target BOLD run space to allow
- the assessment of its alignment with the distorted data.
- The anatomical reference is the fieldmap's reference moved into the BOLD reference's grid through
- the estimated transformation.
- In other words, this plot should be equivalent to that of the
- Preprocessed estimation with varying Phase-Encoding (PE) blips shown above in the
- fieldmap section.
- Therefore, the fieldmap should be positioned relative to the anatomical reference exactly
- as it is positioned in the reportlet above.
+ caption: Estimated fieldmap shown in the corresponding run's space to allow
+ visual assessment of its alignment with the distorted BOLD reference.
+ If a dynamic fieldmap series was used, such as a framewise MEDIC estimate,
+ a representative nonzero volume from that series is shown.
static: false
- subtitle: "Reconstructed B0 map in the corresponding run's space (debug mode)"
+ subtitle: "Estimated B0 map in the corresponding run's space"
- bids: {datatype: figures, desc: sdc, suffix: bold}
caption: Results of performing susceptibility distortion correction (SDC) on the
BOLD reference image. The "distorted" image is the image that would be used to
diff --git a/fmriprep/interfaces/maths.py b/fmriprep/interfaces/maths.py
index 5de7cadb4..f5ae27293 100644
--- a/fmriprep/interfaces/maths.py
+++ b/fmriprep/interfaces/maths.py
@@ -82,3 +82,34 @@ def _run_interface(self, runtime):
self._results['out_file'] = out_file
return runtime
+
+
+class TemporalMeanInputSpec(TraitedSpec):
+ in_file = File(exists=True, mandatory=True, desc='Input 3D or 4D imaging file')
+
+
+class TemporalMeanOutputSpec(TraitedSpec):
+ out_file = File(desc='Temporal mean image')
+
+
+class TemporalMean(SimpleInterface):
+ """Collapse a time series to its temporal mean."""
+
+ input_spec = TemporalMeanInputSpec
+ output_spec = TemporalMeanOutputSpec
+
+ def _run_interface(self, runtime):
+ import nibabel as nb
+
+ img = nb.load(self.inputs.in_file)
+ data = img.get_fdata(dtype='f4')
+ if data.ndim > 3:
+ data = data.mean(axis=3, dtype='f4')
+
+ out_img = img.__class__(data, img.affine, img.header)
+ out_img.set_data_dtype(np.float32)
+ out_file = fname_presuffix(self.inputs.in_file, suffix='_mean', newpath=runtime.cwd)
+ out_img.to_filename(out_file)
+
+ self._results['out_file'] = out_file
+ return runtime
diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py
index 05df737f3..ac2eda34c 100644
--- a/fmriprep/interfaces/resampling.py
+++ b/fmriprep/interfaces/resampling.py
@@ -378,6 +378,15 @@ async def resample_series_async(
The resampled array, with shape ``coordinates.shape[1:] + (N,)``,
where N is the number of volumes in ``data``.
"""
+ fmap_is_series = fmap_hz.ndim > 3
+ if data.ndim == 3 and fmap_is_series:
+ raise ValueError('A 3D source image cannot be resampled with a 4D fieldmap series.')
+ if fmap_is_series and fmap_hz.shape[-1] != data.shape[-1]:
+ raise ValueError(
+ 'Fieldmap series and source series must have matching numbers of volumes '
+ f'(got {fmap_hz.shape[-1]} and {data.shape[-1]}).'
+ )
+
if data.ndim == 3:
return resample_vol(
data,
@@ -385,7 +394,7 @@ async def resample_series_async(
pe_info[0],
jacobian,
hmc_xfms[0] if hmc_xfms else None,
- fmap_hz,
+ fmap_hz[..., 0] if fmap_is_series else fmap_hz,
output_dtype,
order,
mode,
@@ -409,7 +418,7 @@ async def resample_series_async(
pe_info=pe_info[volid],
jacobian=jacobian,
hmc_xfm=hmc_xfms[volid] if hmc_xfms else None,
- fmap_hz=fmap_hz,
+ fmap_hz=fmap_hz[..., volid] if fmap_is_series else fmap_hz,
output=out_array[..., volid],
order=order,
mode=mode,
diff --git a/fmriprep/interfaces/tests/test_maths.py b/fmriprep/interfaces/tests/test_maths.py
index 209878cd5..9d859cd2d 100644
--- a/fmriprep/interfaces/tests/test_maths.py
+++ b/fmriprep/interfaces/tests/test_maths.py
@@ -2,7 +2,7 @@
import numpy as np
from nipype.pipeline import engine as pe
-from fmriprep.interfaces.maths import Clip
+from fmriprep.interfaces.maths import Clip, TemporalMean
def test_Clip(tmp_path):
@@ -41,3 +41,21 @@ def test_Clip(tmp_path):
assert ret.outputs.out_file == str(tmp_path / 'nonpositive/input_clipped.nii')
out_img = nb.load(ret.outputs.out_file)
assert np.allclose(out_img.get_fdata(), [[[-1.0, 0.0], [-2.0, 0.0]]])
+
+
+def test_TemporalMean(tmp_path):
+ in_file = str(tmp_path / 'timeseries.nii')
+ data = np.stack(
+ (
+ np.array([[[1.0, 2.0], [3.0, 4.0]]]),
+ np.array([[[5.0, 6.0], [7.0, 8.0]]]),
+ ),
+ axis=3,
+ )
+ nb.Nifti1Image(data, np.eye(4)).to_filename(in_file)
+
+ meaner = pe.Node(TemporalMean(in_file=in_file), name='meaner', base_dir=tmp_path)
+ ret = meaner.run()
+
+ out_img = nb.load(ret.outputs.out_file)
+ assert np.allclose(out_img.get_fdata(), [[[3.0, 4.0], [5.0, 6.0]]])
diff --git a/fmriprep/interfaces/tests/test_resampling.py b/fmriprep/interfaces/tests/test_resampling.py
new file mode 100644
index 000000000..3ed8380ff
--- /dev/null
+++ b/fmriprep/interfaces/tests/test_resampling.py
@@ -0,0 +1,30 @@
+import numpy as np
+
+from fmriprep.interfaces.resampling import resample_series
+
+
+def test_resample_series_uses_volume_specific_fieldmaps():
+ data = np.zeros((3, 1, 1, 2), dtype='f4')
+ data[:, 0, 0, 0] = [10.0, 20.0, 30.0]
+ data[:, 0, 0, 1] = [100.0, 200.0, 300.0]
+
+ coordinates = np.zeros((3, 1, 1, 1), dtype='f4')
+ pe_info = [(0, 1.0), (0, 1.0)]
+ fmap_hz = np.zeros((1, 1, 1, 2), dtype='f4')
+ fmap_hz[..., 1] = 1.0
+
+ resampled = resample_series(
+ data=data,
+ coordinates=coordinates,
+ pe_info=pe_info,
+ jacobian=False,
+ hmc_xfms=None,
+ fmap_hz=fmap_hz,
+ output_dtype='f4',
+ order=0,
+ mode='nearest',
+ prefilter=False,
+ )
+
+ assert np.allclose(resampled[0, 0, 0, 0], 10.0)
+ assert np.allclose(resampled[0, 0, 0, 1], 200.0)
diff --git a/fmriprep/interfaces/warpkit.py b/fmriprep/interfaces/warpkit.py
new file mode 100644
index 000000000..221629381
--- /dev/null
+++ b/fmriprep/interfaces/warpkit.py
@@ -0,0 +1,137 @@
+from __future__ import annotations
+
+from pathlib import Path
+
+import nibabel as nb
+import numpy as np
+from nipype.interfaces.base import File, InputMultiObject, SimpleInterface, TraitedSpec, traits
+
+
+def _nvols(img: nb.spatialimages.SpatialImage) -> int:
+ return img.shape[3] if img.ndim > 3 else 1
+
+
+def _pad_to_nvols(in_file: Path, target_nvols: int, noise_frames: int) -> Path:
+ img = nb.load(in_file)
+ nvols = _nvols(img)
+ if nvols == target_nvols:
+ return in_file
+ if nvols > target_nvols or target_nvols - nvols != noise_frames:
+ raise RuntimeError(
+ f'warpkit MEDIC returned {nvols} volume(s) for <{in_file}>, but '
+ f'{target_nvols} were expected after trimming {noise_frames} noise frame(s).'
+ )
+
+ data = np.asanyarray(img.dataobj)
+ padded = np.concatenate(
+ (data, np.repeat(data[..., -1:], target_nvols - nvols, axis=-1)),
+ axis=-1,
+ )
+ if in_file.name.endswith('.nii.gz'):
+ out_name = in_file.name[:-7] + '_padded.nii.gz'
+ elif in_file.name.endswith('.nii'):
+ out_name = in_file.name[:-4] + '_padded.nii'
+ else:
+ out_name = in_file.name + '_padded'
+ out_file = in_file.parent / out_name
+ out_img = img.__class__(
+ padded.astype(img.get_data_dtype(), copy=False),
+ img.affine,
+ img.header,
+ )
+ out_img.header.set_data_shape(padded.shape)
+ out_img.to_filename(out_file)
+ return out_file
+
+
+class WarpkitMEDICInputSpec(TraitedSpec):
+ phase = InputMultiObject(File(exists=True), mandatory=True, desc='Per-echo phase series')
+ magnitude = InputMultiObject(
+ File(exists=True), mandatory=True, desc='Per-echo magnitude series'
+ )
+ tes = traits.List(
+ traits.Float,
+ mandatory=True,
+ minlen=3,
+ desc='Echo times in milliseconds, one per echo',
+ )
+ total_readout_time = traits.Float(mandatory=True, desc='Total readout time in seconds')
+ phase_encoding_direction = traits.Enum(
+ 'i',
+ 'i-',
+ 'j',
+ 'j-',
+ 'k',
+ 'k-',
+ 'x',
+ 'x-',
+ 'y',
+ 'y-',
+ 'z',
+ 'z-',
+ mandatory=True,
+ desc='Phase encoding direction',
+ )
+ n_cpus = traits.Int(1, usedefault=True, desc='Number of worker threads for MEDIC')
+ noise_frames = traits.Int(0, usedefault=True, desc='Trailing noise frames to trim')
+ wrap_limit = traits.Bool(False, usedefault=True, desc='Disable phase unwrapping heuristics')
+ debug = traits.Bool(False, usedefault=True, desc='Enable warpkit debug mode')
+
+
+class WarpkitMEDICOutputSpec(TraitedSpec):
+ fieldmap_native = File(exists=True, desc='Framewise fieldmaps in distorted/native space')
+ displacement_map = File(exists=True, desc='Framewise displacement maps')
+ fieldmap = File(exists=True, desc='Framewise fieldmaps in undistorted space')
+
+
+class WarpkitMEDIC(SimpleInterface):
+ """Run warpkit's MEDIC workflow for a multi-echo BOLD series."""
+
+ input_spec = WarpkitMEDICInputSpec
+ output_spec = WarpkitMEDICOutputSpec
+
+ def _run_interface(self, runtime):
+ try:
+ from warpkit.api import medic
+ except ImportError as exc:
+ raise RuntimeError(
+ 'warpkit is required for MEDIC-based multi-echo SDC. '
+ 'Install fMRIPrep with the optional warpkit extra on Python 3.11+.'
+ ) from exc
+
+ out_prefix = Path(runtime.cwd) / 'warpkit_medic'
+ result = medic(
+ phase=self.inputs.phase,
+ magnitude=self.inputs.magnitude,
+ out_prefix=out_prefix,
+ tes=self.inputs.tes,
+ total_readout_time=self.inputs.total_readout_time,
+ phase_encoding_direction=self.inputs.phase_encoding_direction,
+ n_cpus=self.inputs.n_cpus,
+ noise_frames=self.inputs.noise_frames,
+ wrap_limit=self.inputs.wrap_limit,
+ debug=self.inputs.debug,
+ )
+
+ fieldmap_native = Path(result.fieldmap_native)
+ displacement_map = Path(result.displacement_map)
+ fieldmap = Path(result.fieldmap)
+
+ if self.inputs.noise_frames:
+ target_nvols = _nvols(nb.load(self.inputs.magnitude[0]))
+ fieldmap_native = _pad_to_nvols(
+ fieldmap_native,
+ target_nvols,
+ self.inputs.noise_frames,
+ )
+ displacement_map = _pad_to_nvols(
+ displacement_map,
+ target_nvols,
+ self.inputs.noise_frames,
+ )
+ fieldmap = _pad_to_nvols(fieldmap, target_nvols, self.inputs.noise_frames)
+
+ self._results['fieldmap_native'] = str(fieldmap_native)
+ self._results['displacement_map'] = str(displacement_map)
+ self._results['fieldmap'] = str(fieldmap)
+ return runtime
diff --git a/fmriprep/utils/bids.py b/fmriprep/utils/bids.py
index 708f283c5..eb409fdfb 100644
--- a/fmriprep/utils/bids.py
+++ b/fmriprep/utils/bids.py
@@ -381,6 +381,43 @@ def _unique(inlist):
return {k: _unique(v) for k, v in entities.items()}
+def select_bold_magnitude_files(file_list):
+ """Return BOLD files that should be treated as magnitude data.
+
+ This filters out explicit phase companions so pipelines that operate on
+ ``subject_data['bold']`` do not accidentally schedule ``part-phase`` series
+ as if they were standard BOLD inputs.
+ """
+ from bids.layout import parse_file_entities
+
+ return [
+ fname
+ for fname in listify(file_list)
+ if parse_file_entities(fname).get('part') in (None, 'mag')
+ ]
+
+
+def collect_bold_part_files(layout, bold_series, part):
+ """Collect BOLD part-files matching a run, sorted by echo time."""
+ magnitude_files = select_bold_magnitude_files(bold_series)
+ if not magnitude_files:
+ return []
+
+ entities = extract_entities(magnitude_files)
+ entities.pop('part', None)
+ entities.update(
+ datatype='func',
+ suffix='bold',
+ part=part,
+ extension=['.nii', '.nii.gz'],
+ )
+
+ return sorted(
+ layout.get(return_type='file', **entities),
+ key=lambda fname: layout.get_metadata(fname).get('EchoTime', 0),
+ )
+
+
def dismiss_echo(entities=None):
"""Set entities to dismiss in a DerivativesDataSink."""
if entities is None:
diff --git a/fmriprep/utils/tests/test_derivative_cache.py b/fmriprep/utils/tests/test_derivative_cache.py
index 87382f73f..5313ab1e1 100644
--- a/fmriprep/utils/tests/test_derivative_cache.py
+++ b/fmriprep/utils/tests/test_derivative_cache.py
@@ -1,6 +1,7 @@
from pathlib import Path
import pytest
+from bids.layout import BIDSLayout
from fmriprep.utils import bids
@@ -56,3 +57,31 @@ def test_transforms_found_as_str(tmp_path: Path, xfm: str):
fieldmap_id='auto_00000',
)
assert derivs == {'transforms': {xfm: str(to_find)}}
+
+
+def test_select_bold_magnitude_files():
+ files = [
+ 'sub-01/func/sub-01_task-rest_echo-1_part-mag_bold.nii.gz',
+ 'sub-01/func/sub-01_task-rest_echo-1_part-phase_bold.nii.gz',
+ 'sub-01/func/sub-01_task-rest_echo-2_part-mag_bold.nii.gz',
+ ]
+ assert bids.select_bold_magnitude_files(files) == [files[0], files[2]]
+
+
+def test_collect_bold_part_files(tmp_path: Path):
+ (tmp_path / 'dataset_description.json').write_text('{"Name": "test", "BIDSVersion": "1.8.0"}')
+ func_dir = tmp_path / 'sub-01' / 'func'
+ func_dir.mkdir(parents=True)
+
+ magnitude_files = []
+ phase_files = []
+ for echo in range(1, 4):
+ mag = func_dir / f'sub-01_task-rest_echo-{echo}_part-mag_bold.nii.gz'
+ phase = func_dir / f'sub-01_task-rest_echo-{echo}_part-phase_bold.nii.gz'
+ mag.touch()
+ phase.touch()
+ magnitude_files.append(str(mag))
+ phase_files.append(str(phase))
+
+ layout = BIDSLayout(tmp_path, validate=False)
+ assert bids.collect_bold_part_files(layout, magnitude_files, part='phase') == phase_files
diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py
index 0d8fa53d3..f6c590acc 100644
--- a/fmriprep/workflows/base.py
+++ b/fmriprep/workflows/base.py
@@ -44,7 +44,7 @@
from .. import config
from ..interfaces import DerivativesDataSink
from ..interfaces.reports import AboutSummary, SubjectSummary
-from ..utils.bids import dismiss_echo
+from ..utils.bids import dismiss_echo, select_bold_magnitude_files
def init_fmriprep_wf():
@@ -248,11 +248,19 @@ def init_single_subject_wf(
bold_runs = [
sorted(
- listify(run),
+ select_bold_magnitude_files(listify(run)),
key=lambda file: config.execution.layout.get_metadata(file).get('EchoTime', 0),
)
for run in subject_data['bold']
]
+ bold_runs = [run for run in bold_runs if run]
+
+ if not anat_only and not bold_runs:
+ task_id = config.execution.task_id or ''
+ raise RuntimeError(
+ f'No magnitude BOLD images found for participant {subject_id} and '
+ f'task {task_id}. All workflows require BOLD images.'
+ )
if subject_data['roi']:
warnings.warn(
@@ -633,11 +641,14 @@ def init_single_subject_wf(
)
fmap_cache.update(fmaps)
+ fieldmap_bold_runs = [
+ run for run in bold_runs if not (config.workflow.me_use_warpkit and len(run) > 1)
+ ]
all_estimators, estimator_map = map_fieldmap_estimation(
layout=config.execution.layout,
subject_id=subject_id,
session_id=session_id,
- bold_data=bold_runs,
+ bold_data=fieldmap_bold_runs,
ignore_fieldmaps='fieldmaps' in config.workflow.ignore,
use_syn=config.workflow.use_syn_sdc,
force_syn='syn-sdc' in config.workflow.force,
@@ -835,7 +846,8 @@ def init_single_subject_wf(
for bold_series in bold_runs:
bold_file = bold_series[0]
- fieldmap_id = estimator_map.get(bold_file)
+ use_warpkit = config.workflow.me_use_warpkit and len(bold_series) > 1
+ fieldmap_id = None if use_warpkit else estimator_map.get(bold_file)
jacobian = False
if len(bold_series) == 2:
@@ -848,7 +860,7 @@ def init_single_subject_wf(
'Please set "--echo-idx" to process one echo at a time.'
)
- if fieldmap_id:
+ if fieldmap_id and not use_warpkit:
if 'fmap-jacobian' in config.workflow.force:
jacobian = True
elif 'fmap-jacobian' not in config.workflow.ignore:
@@ -877,6 +889,7 @@ def init_single_subject_wf(
bold_series=bold_series,
precomputed=functional_cache,
fieldmap_id=fieldmap_id,
+ use_warpkit=use_warpkit,
jacobian=jacobian,
)
if bold_wf is None:
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 871be8104..ebdc85c06 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -52,11 +52,35 @@
from .t2s import init_t2s_reporting_wf
+def _build_postdesc(use_warpkit: bool) -> str:
+ """Build post-description boilerplate for functional preprocessing."""
+ postdesc = """\
+All resamplings can be performed with *a single interpolation
+step* by composing all the pertinent transformations (i.e. head-motion
+transform matrices, susceptibility distortion correction when available,
+and co-registrations to anatomical and output spaces).
+Gridded (volumetric) resamplings were performed using `nitransforms`,
+configured with cubic B-spline interpolation.
+"""
+ if use_warpkit:
+ postdesc += """\
+For compatible multi-echo runs with accompanying phase images, susceptibility
+distortion correction was performed with the Multi-Echo DIstortion Correction
+(MEDIC) procedure as implemented in *warpkit* [@van2023medic]. MEDIC estimated
+a time-varying B0 fieldmap series by fitting phase evolution across echoes
+separately for each frame of the run. The resulting fieldmap series was aligned
+to the head-motion-correction reference and applied during native-space
+resampling together with motion correction.
+"""
+ return postdesc
+
+
def init_bold_wf(
*,
bold_series: list[str],
precomputed: dict | None = None,
fieldmap_id: str | None = None,
+ use_warpkit: bool = False,
jacobian: bool = False,
) -> pe.Workflow:
"""
@@ -86,6 +110,8 @@ def init_bold_wf(
fieldmap_id
ID of the fieldmap to use to correct this BOLD series. If :obj:`None`,
no correction will be applied.
+ use_warpkit
+ Use warpkit MEDIC for compatible multi-echo BOLD runs.
Inputs
------
@@ -187,14 +213,7 @@ def init_bold_wf(
)
workflow = Workflow(name=_get_wf_name(bold_file, 'bold'))
- workflow.__postdesc__ = """\
-All resamplings can be performed with *a single interpolation
-step* by composing all the pertinent transformations (i.e. head-motion
-transform matrices, susceptibility distortion correction when available,
-and co-registrations to anatomical and output spaces).
-Gridded (volumetric) resamplings were performed using `nitransforms`,
-configured with cubic B-spline interpolation.
-"""
+ workflow.__postdesc__ = _build_postdesc(use_warpkit)
inputnode = pe.Node(
niu.IdentityInterface(
@@ -246,6 +265,7 @@ def init_bold_wf(
bold_series=bold_series,
precomputed=precomputed,
fieldmap_id=fieldmap_id,
+ use_warpkit=use_warpkit,
jacobian=jacobian,
omp_nthreads=omp_nthreads,
)
@@ -287,6 +307,7 @@ def init_bold_wf(
bold_native_wf = init_bold_native_wf(
bold_series=bold_series,
fieldmap_id=fieldmap_id,
+ use_warpkit=use_warpkit,
jacobian=jacobian,
omp_nthreads=omp_nthreads,
)
@@ -302,6 +323,7 @@ def init_bold_wf(
('outputnode.bold_mask', 'inputnode.bold_mask'),
('outputnode.motion_xfm', 'inputnode.motion_xfm'),
('outputnode.boldref2fmap_xfm', 'inputnode.boldref2fmap_xfm'),
+ ('outputnode.fieldmap', 'inputnode.fieldmap'),
('outputnode.dummy_scans', 'inputnode.dummy_scans'),
]),
]) # fmt:skip
diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py
index d929a1d93..9c54ac351 100644
--- a/fmriprep/workflows/bold/fit.py
+++ b/fmriprep/workflows/bold/fit.py
@@ -35,13 +35,15 @@
from sdcflows.workflows.apply.registration import init_coeff2epi_wf
from ... import config
+from ...interfaces.maths import TemporalMean
from ...interfaces.reports import FunctionalSummary
from ...interfaces.resampling import (
DistortionParameters,
ReconstructFieldmap,
ResampleSeries,
)
-from ...utils.bids import extract_entities
+from ...interfaces.warpkit import WarpkitMEDIC
+from ...utils.bids import collect_bold_part_files, extract_entities
from ...utils.misc import estimate_bold_mem_usage
# BOLD workflows
@@ -97,6 +99,7 @@ def init_bold_fit_wf(
bold_series: list[str],
precomputed: dict | None = None,
fieldmap_id: str | None = None,
+ use_warpkit: bool = False,
jacobian: bool = False,
omp_nthreads: int = 1,
name: str = 'bold_fit_wf',
@@ -126,6 +129,8 @@ def init_bold_fit_wf(
fieldmap_id
ID of the fieldmap to use to correct this BOLD series. If :obj:`None`,
no correction will be applied.
+ use_warpkit
+ Run warpkit MEDIC for compatible multi-echo BOLD series.
Inputs
------
@@ -179,6 +184,8 @@ def init_bold_fit_wf(
boldref2fmap_xfm
Affine transform mapping from BOLD reference space to the fieldmap
space, if applicable.
+ fieldmap
+ Direct fieldmap series in BOLD reference space, if generated in-workflow.
dummy_scans
The number of dummy scans declared or detected at the beginning of the series.
@@ -226,13 +233,50 @@ def init_bold_fit_wf(
# Get metadata from BOLD file(s)
entities = extract_entities(bold_series)
- metadata = layout.get_metadata(bold_file)
+ all_metadata = [layout.get_metadata(fname) for fname in bold_series]
+ metadata = all_metadata[0]
orientation = ''.join(nb.aff2axcodes(nb.load(bold_file).affine))
_bold_tlen, mem_gb = estimate_bold_mem_usage(bold_file)
# Boolean used to update workflow self-descriptions
multiecho = len(bold_series) > 1
+ warpkit_enabled = use_warpkit and multiecho
+ phase_files = []
+ tes_ms = []
+ if warpkit_enabled:
+ phase_files = collect_bold_part_files(layout, bold_series, part='phase')
+ if len(phase_files) != len(bold_series):
+ raise RuntimeError(
+ 'warpkit MEDIC requires one part-phase BOLD file per echo for '
+ f'<{bold_file}> (found {len(phase_files)} phase file(s) for '
+ f'{len(bold_series)} echo(es)).'
+ )
+
+ try:
+ tes_ms = [float(md['EchoTime']) * 1000 for md in all_metadata]
+ except KeyError as exc:
+ raise RuntimeError(
+ f'warpkit MEDIC requires EchoTime metadata for every echo of <{bold_file}>.'
+ ) from exc
+
+ # `bold_series` and `phase_files` are independently sorted by EchoTime
+ # (in workflows/base.py and collect_bold_part_files respectively). A
+ # missing/zero EchoTime sidecar would silently swap echo order on one
+ # side without the other, feeding warpkit a transposed pair set. Check
+ # the orderings match before handing them off.
+ if layout is None:
+ raise RuntimeError('warpkit MEDIC requires a BIDSLayout.')
+
+ phase_tes_ms = [
+ float(layout.get_metadata(p).get('EchoTime', 0.0)) * 1000 for p in phase_files
+ ]
+ if phase_tes_ms != tes_ms:
+ raise RuntimeError(
+ 'warpkit MEDIC echo-ordering mismatch between magnitude and '
+ f'phase for <{bold_file}>: magnitude TEs (ms)={tes_ms} '
+ f'phase TEs (ms)={phase_tes_ms}.'
+ )
hmc_boldref = precomputed.get('hmc_boldref')
coreg_boldref = precomputed.get('coreg_boldref')
@@ -281,6 +325,7 @@ def init_bold_fit_wf(
'motion_xfm',
'boldref2anat_xfm',
'boldref2fmap_xfm',
+ 'fieldmap',
],
),
name='outputnode',
@@ -338,6 +383,8 @@ def init_bold_fit_wf(
run_without_submitting=True,
)
summary.inputs.dummy_scans = config.workflow.dummy_scans
+ if warpkit_enabled:
+ summary.inputs.distortion_correction = 'MEDIC (warpkit)'
if config.workflow.level == 'full':
# Hack. More pain than it's worth to connect this up at a higher level.
# We can consider separating out fit and transform summaries,
@@ -349,7 +396,9 @@ def init_bold_fit_wf(
func_fit_reports_wf = init_func_fit_reports_wf(
source_file=bold_file,
- sdc_correction=fieldmap_id is not None,
+ sdc_correction=fieldmap_id is not None or warpkit_enabled,
+ fieldmap_registration=fieldmap_id is not None and not warpkit_enabled,
+ fieldmap_reportlet=warpkit_enabled,
freesurfer=config.workflow.run_reconall,
output_dir=config.execution.fmriprep_dir,
)
@@ -463,8 +512,64 @@ def init_bold_fit_wf(
else:
config.loggers.workflow.info('Found motion correction transforms - skipping Stage 2')
- # Stage 3: Register fieldmap to boldref and reconstruct in BOLD space
- if fieldmap_id:
+ # Stage 3: Prepare distortion correction in BOLD space
+ warpkit_distortion_params = None
+ if warpkit_enabled:
+ config.loggers.workflow.info('Stage 3: Adding warpkit MEDIC workflow')
+ warpkit_distortion_params = pe.Node(
+ DistortionParameters(
+ metadata=metadata,
+ in_file=bold_file,
+ fallback=config.workflow.fallback_total_readout_time,
+ ),
+ name='warpkit_distortion_params',
+ run_without_submitting=True,
+ )
+ warpkit_medic = pe.Node(
+ WarpkitMEDIC(
+ magnitude=bold_series,
+ phase=phase_files,
+ tes=tes_ms,
+ n_cpus=omp_nthreads,
+ noise_frames=config.workflow.me_warpkit_noise_frames,
+ ),
+ name='warpkit_medic',
+ # MEDIC holds all N echoes of phase + magnitude in memory as
+ # float64, plus warpkit's internal state. `mem_gb['resampled']`
+ # is the size of a single resampled BOLD volume — wildly low
+ # for high-res / many-echo datasets.
+ mem_gb=2 * len(bold_series) * mem_gb['filesize'],
+ )
+ warpkit_fieldmap = pe.Node(
+ ResampleSeries(jacobian=False),
+ name='warpkit_fieldmap',
+ n_procs=omp_nthreads,
+ mem_gb=mem_gb['resampled'],
+ )
+ warpkit_mean = pe.Node(TemporalMean(), name='warpkit_mean', mem_gb=1)
+
+ workflow.connect([
+ (warpkit_distortion_params, warpkit_medic, [
+ ('readout_time', 'total_readout_time'),
+ ('pe_direction', 'phase_encoding_direction'),
+ ]),
+ # `fieldmap` is the undistorted-space, consumer-facing output from
+ # warpkit.api.medic — it lives on the same grid as the corrected EPI.
+ # `fieldmap_native` is a debug-only field in *distorted* space; using
+ # it here causes fmriprep's pull-resampler to sample the field at
+ # the wrong physical location (off by the SDC displacement itself).
+ # See resample_vol: fmap_hz is read on the target/undistorted grid.
+ (warpkit_medic, warpkit_fieldmap, [('fieldmap', 'in_file')]),
+ (hmcref_buffer, warpkit_fieldmap, [('boldref', 'ref_file')]),
+ (hmc_buffer, warpkit_fieldmap, [('hmc_xforms', 'transforms')]),
+ (warpkit_fieldmap, warpkit_mean, [('out_file', 'in_file')]),
+ (warpkit_fieldmap, outputnode, [('out_file', 'fieldmap')]),
+ # FieldmapReportlet expects a 3D field; pass the temporal mean,
+ # not the 4D series (which would render only frame 0 or crash).
+ (warpkit_mean, func_fit_reports_wf, [('out_file', 'inputnode.fieldmap')]),
+ (fmapref_buffer, func_fit_reports_wf, [('out', 'inputnode.sdc_boldref')]),
+ ]) # fmt:skip
+ elif fieldmap_id:
config.loggers.workflow.info('Stage 3: Adding fieldmap reconstruction workflow')
fmap_select = pe.Node(
KeySelect(
@@ -586,7 +691,34 @@ def init_bold_fit_wf(
(ds_boldmask_wf, regref_buffer, [('outputnode.boldmask', 'boldmask')]),
]) # fmt:skip
- if fieldmap_id:
+ if warpkit_enabled:
+ unwarp_boldref = pe.Node(
+ ResampleSeries(jacobian=jacobian),
+ name='unwarp_boldref',
+ n_procs=omp_nthreads,
+ mem_gb=mem_gb['resampled'],
+ )
+ skullstrip_bold_wf = init_skullstrip_bold_wf()
+
+ workflow.connect([
+ (fmapref_buffer, unwarp_boldref, [('out', 'ref_file')]),
+ (enhance_boldref_wf, unwarp_boldref, [
+ ('outputnode.bias_corrected_file', 'in_file'),
+ ]),
+ (warpkit_mean, unwarp_boldref, [('out_file', 'fieldmap')]),
+ (warpkit_distortion_params, unwarp_boldref, [
+ ('readout_time', 'ro_time'),
+ ('pe_direction', 'pe_dir'),
+ ]),
+ (unwarp_boldref, ds_coreg_boldref_wf, [('out_file', 'inputnode.boldref')]),
+ (ds_coreg_boldref_wf, skullstrip_bold_wf, [
+ ('outputnode.boldref', 'inputnode.in_file'),
+ ]),
+ (skullstrip_bold_wf, ds_boldmask_wf, [
+ ('outputnode.mask_file', 'inputnode.boldmask'),
+ ]),
+ ]) # fmt:skip
+ elif fieldmap_id:
distortion_params = pe.Node(
DistortionParameters(
metadata=metadata,
@@ -717,6 +849,7 @@ def init_bold_native_wf(
*,
bold_series: list[str],
fieldmap_id: str | None = None,
+ use_warpkit: bool = False,
jacobian: bool = False,
omp_nthreads: int = 1,
name: str = 'bold_native_wf',
@@ -749,6 +882,9 @@ def init_bold_native_wf(
fieldmap_id
ID of the fieldmap to use to correct this BOLD series. If :obj:`None`,
no correction will be applied.
+ use_warpkit
+ Use a direct warpkit MEDIC fieldmap series instead of reconstructing
+ a static SDCFlows fieldmap.
Inputs
------
@@ -768,6 +904,8 @@ def init_bold_native_wf(
List of fieldmap reference files (collated with fmap_id)
fmap_coeff
List of lists of spline coefficient files (collated with fmap_id)
+ fieldmap
+ Framewise fieldmap series already aligned to ``boldref``.
Outputs
-------
@@ -810,6 +948,7 @@ def init_bold_native_wf(
all_metadata = [layout.get_metadata(bold_file) for bold_file in bold_series]
echo_times = [md.get('EchoTime') for md in all_metadata]
multiecho = len(bold_series) > 1
+ warpkit_enabled = use_warpkit and multiecho
bold_file = bold_series[0]
metadata = all_metadata[0]
@@ -846,6 +985,7 @@ def init_bold_native_wf(
'fmap_ref',
'fmap_coeff',
'fmap_id',
+ 'fieldmap',
],
),
name='inputnode',
@@ -900,13 +1040,7 @@ def init_bold_native_wf(
workflow.connect([(validate_bold, boldbuffer, [('out_file', 'bold_file')])])
# Prepare fieldmap metadata
- if fieldmap_id:
- fmap_select = pe.Node(
- KeySelect(fields=['fmap_ref', 'fmap_coeff'], key=fieldmap_id),
- name='fmap_select',
- run_without_submitting=True,
- )
-
+ if fieldmap_id or warpkit_enabled:
distortion_params = pe.Node(
DistortionParameters(
metadata=metadata,
@@ -916,16 +1050,25 @@ def init_bold_native_wf(
name='distortion_params',
run_without_submitting=True,
)
+ workflow.connect([
+ (distortion_params, boldbuffer, [
+ ('readout_time', 'ro_time'),
+ ('pe_direction', 'pe_dir'),
+ ]),
+ ]) # fmt:skip
+
+ if fieldmap_id and not warpkit_enabled:
+ fmap_select = pe.Node(
+ KeySelect(fields=['fmap_ref', 'fmap_coeff'], key=fieldmap_id),
+ name='fmap_select',
+ run_without_submitting=True,
+ )
workflow.connect([
(inputnode, fmap_select, [
('fmap_ref', 'fmap_ref'),
('fmap_coeff', 'fmap_coeff'),
('fmap_id', 'keys'),
]),
- (distortion_params, boldbuffer, [
- ('readout_time', 'ro_time'),
- ('pe_direction', 'pe_dir'),
- ]),
]) # fmt:skip
# Resample to boldref
@@ -948,7 +1091,9 @@ def init_bold_native_wf(
]),
]) # fmt:skip
- if fieldmap_id:
+ if warpkit_enabled:
+ workflow.connect([(inputnode, boldref_bold, [('fieldmap', 'fieldmap')])])
+ elif fieldmap_id:
boldref_fmap = pe.Node(ReconstructFieldmap(inverse=[True]), name='boldref_fmap', mem_gb=1)
workflow.connect([
(inputnode, boldref_fmap, [
diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py
index f8c9ebbb3..025efa9b2 100644
--- a/fmriprep/workflows/bold/outputs.py
+++ b/fmriprep/workflows/bold/outputs.py
@@ -143,6 +143,8 @@ def init_func_fit_reports_wf(
*,
source_file: str,
sdc_correction: bool,
+ fieldmap_registration: bool,
+ fieldmap_reportlet: bool,
freesurfer: bool,
output_dir: str,
name='func_fit_reports_wf',
@@ -308,7 +310,7 @@ def init_func_fit_reports_wf(
# Before: T1w brain with white matter mask
# After: Resampled boldref with white matter mask
- if sdc_correction:
+ if fieldmap_registration:
fmapref_boldref = pe.Node(
ApplyTransforms(
dimension=3,
@@ -344,7 +346,53 @@ def init_func_fit_reports_wf(
name='ds_sdcreg_report',
)
- # SDC2
+ workflow.connect([
+ (inputnode, fmapref_boldref, [
+ ('fmap_ref', 'input_image'),
+ ('coreg_boldref', 'reference_image'),
+ ('boldref2fmap_xfm', 'transforms'),
+ ]),
+ (inputnode, sdcreg_report, [
+ ('sdc_boldref', 'reference'),
+ ('fieldmap', 'fieldmap'),
+ ('bold_mask', 'mask'),
+ ]),
+ (fmapref_boldref, sdcreg_report, [('output_image', 'moving')]),
+ (sdcreg_report, ds_sdcreg_report, [('out_report', 'in_file')]),
+ ]) # fmt:skip
+
+ if fieldmap_reportlet:
+ fieldmap_report = pe.Node(
+ FieldmapReportlet(
+ reference_label='Distorted BOLD reference',
+ show=0,
+ ),
+ name='fieldmap_report',
+ mem_gb=0.1,
+ )
+
+ ds_fieldmap_report = pe.Node(
+ DerivativesDataSink(
+ source_file=source_file,
+ base_directory=output_dir,
+ desc='fieldmap',
+ suffix='bold',
+ datatype='figures',
+ dismiss_entities=dismiss_echo(),
+ ),
+ name='ds_fieldmap_report',
+ )
+
+ workflow.connect([
+ (inputnode, fieldmap_report, [
+ ('bold_mask', 'mask'),
+ ('sdc_boldref', 'reference'),
+ ('fieldmap', 'fieldmap'),
+ ]),
+ (fieldmap_report, ds_fieldmap_report, [('out_report', 'in_file')]),
+ ]) # fmt:skip
+
+ if sdc_correction:
sdc_report = pe.Node(
SimpleBeforeAfter(
before_label='Distorted',
@@ -368,18 +416,6 @@ def init_func_fit_reports_wf(
)
workflow.connect([
- (inputnode, fmapref_boldref, [
- ('fmap_ref', 'input_image'),
- ('coreg_boldref', 'reference_image'),
- ('boldref2fmap_xfm', 'transforms'),
- ]),
- (inputnode, sdcreg_report, [
- ('sdc_boldref', 'reference'),
- ('fieldmap', 'fieldmap'),
- ('bold_mask', 'mask'),
- ]),
- (fmapref_boldref, sdcreg_report, [('output_image', 'moving')]),
- (sdcreg_report, ds_sdcreg_report, [('out_report', 'in_file')]),
(inputnode, sdc_report, [
('sdc_boldref', 'before'),
('coreg_boldref', 'after'),
diff --git a/fmriprep/workflows/bold/tests/test_base.py b/fmriprep/workflows/bold/tests/test_base.py
index ff4d677f0..7ceea1719 100644
--- a/fmriprep/workflows/bold/tests/test_base.py
+++ b/fmriprep/workflows/bold/tests/test_base.py
@@ -9,7 +9,7 @@
from .... import config
from ...tests import mock_config
from ...tests.layouts import get_layout
-from ..base import init_bold_wf
+from ..base import _build_postdesc, init_bold_wf
@pytest.fixture(scope='module', autouse=True)
@@ -81,3 +81,11 @@ def test_bold_wf(
flatgraph = wf._create_flat_graph()
generate_expanded_graph(flatgraph)
+
+
+def test_bold_wf_warpkit_postdesc():
+ postdesc = _build_postdesc(True)
+
+ assert 'MEDIC' in postdesc
+ assert 'time-varying B0 fieldmap series' in postdesc
+ assert 'van2023medic' in postdesc
diff --git a/fmriprep/workflows/bold/tests/test_fit.py b/fmriprep/workflows/bold/tests/test_fit.py
index 31bdfde69..a1b55bd8c 100644
--- a/fmriprep/workflows/bold/tests/test_fit.py
+++ b/fmriprep/workflows/bold/tests/test_fit.py
@@ -3,6 +3,8 @@
import nibabel as nb
import numpy as np
import pytest
+from nipype.interfaces import utility as niu
+from nipype.pipeline import engine as pe
from nipype.pipeline.engine.utils import generate_expanded_graph
from niworkflows.utils.testing import generate_bids_skeleton
@@ -47,6 +49,17 @@ def _make_params(
)
+def _stub_enhance_and_skullstrip_bold_wf(*, omp_nthreads: int = 1):
+ workflow = pe.Workflow(name='stub_enhance_and_skullstrip_bold_wf')
+ inputnode = pe.Node(niu.IdentityInterface(fields=['in_file']), name='inputnode')
+ outputnode = pe.Node(
+ niu.IdentityInterface(fields=['bias_corrected_file', 'mask_file']),
+ name='outputnode',
+ )
+ workflow.connect([(inputnode, outputnode, [('in_file', 'bias_corrected_file')])])
+ return workflow
+
+
@pytest.mark.parametrize('task', ['rest', 'nback'])
@pytest.mark.parametrize('fieldmap_id', ['phasediff', None])
@pytest.mark.parametrize(
@@ -175,3 +188,50 @@ def test_bold_native_precomputes(
flatgraph = wf._create_flat_graph()
generate_expanded_graph(flatgraph)
+
+
+def test_bold_warpkit_graph(bids_root: Path, tmp_path: Path, monkeypatch):
+ img = nb.Nifti1Image(np.zeros((10, 10, 10, 10)), np.eye(4))
+ bold_series = [
+ str(bids_root / 'sub-01' / 'func' / f'sub-01_task-nback_echo-{i}_bold.nii.gz')
+ for i in range(1, 4)
+ ]
+ phase_series = [
+ str(bids_root / 'sub-01' / 'func' / f'sub-01_task-nback_echo-{i}_part-phase_bold.nii.gz')
+ for i in range(1, 4)
+ ]
+
+ for path in bold_series + phase_series:
+ img.to_filename(path)
+
+ monkeypatch.setattr(
+ 'fmriprep.workflows.bold.fit.init_enhance_and_skullstrip_bold_wf',
+ _stub_enhance_and_skullstrip_bold_wf,
+ )
+
+ with mock_config(bids_dir=bids_root):
+ config.workflow.bold2anat_init = 't1w'
+ config.workflow.me_warpkit_noise_frames = 3
+ fit_wf = init_bold_fit_wf(
+ bold_series=bold_series,
+ fieldmap_id=None,
+ use_warpkit=True,
+ omp_nthreads=1,
+ )
+ native_wf = init_bold_native_wf(
+ bold_series=bold_series,
+ fieldmap_id=None,
+ use_warpkit=True,
+ omp_nthreads=1,
+ )
+
+ generate_expanded_graph(fit_wf._create_flat_graph())
+ generate_expanded_graph(native_wf._create_flat_graph())
+
+ fit_nodes = fit_wf.list_node_names()
+ assert fit_wf.get_node('warpkit_medic').inputs.noise_frames == 3
+ assert any(name.endswith('fieldmap_report') for name in fit_nodes)
+ assert any(name.endswith('ds_fieldmap_report') for name in fit_nodes)
+ assert any(name.endswith('sdc_report') for name in fit_nodes)
+ assert any(name.endswith('ds_sdc_report') for name in fit_nodes)
+ assert not any(name.endswith('sdcreg_report') for name in fit_nodes)
diff --git a/pixi.lock b/pixi.lock
index dcf5aa4db..39a699fe3 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -6,8 +6,6 @@ environments:
- url: https://conda.anaconda.org/conda-forge/
indexes:
- https://pypi.org/simple
- options:
- pypi-prerelease-mode: if-necessary-or-explicit
packages:
linux-64:
- conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-4_kmp_llvm.conda
@@ -734,6 +732,7 @@ environments:
- pypi: https://files.pythonhosted.org/packages/c2/14/e2a54fabd4f08cd7af1c07030603c3356b74da07f7cc056e600436edfa17/tzlocal-5.3.1-py3-none-any.whl
- pypi: https://files.pythonhosted.org/packages/e5/4d/2e577f6db7aa0f932d19f799c18f604b2b302c65f733419b900ec07dbade/universal_pathlib-0.2.6-py3-none-any.whl
- pypi: https://files.pythonhosted.org/packages/a7/c2/fe1e52489ae3122415c51f387e221dd0773709bad6c6cdaa599e8a2c5185/urllib3-2.5.0-py3-none-any.whl
+ - pypi: https://files.pythonhosted.org/packages/39/ff/b279045420056a817049a8f745f22f68085a11f01ff221ea22bc0dafebee/warpkit-1.2.1-cp312-cp312-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl
- pypi: https://files.pythonhosted.org/packages/af/b5/123f13c975e9f27ab9c0770f514345bd406d0e8d3b7a0723af9d43f710af/wcwidth-0.2.14-py2.py3-none-any.whl
- pypi: https://files.pythonhosted.org/packages/9f/81/5d931d78d0eb732b95dc3ddaeeb71c8bb572fb01356e9133916cd729ecdd/wrapt-1.17.3-cp312-cp312-manylinux1_x86_64.manylinux_2_28_x86_64.manylinux_2_5_x86_64.whl
- pypi: https://files.pythonhosted.org/packages/d6/7d/b77455d7c7c51255b2992b429107fab811b2e36ceaf76da1e55a045dc568/xyzservices-2025.4.0-py3-none-any.whl
@@ -1095,6 +1094,7 @@ environments:
- pypi: https://files.pythonhosted.org/packages/c2/14/e2a54fabd4f08cd7af1c07030603c3356b74da07f7cc056e600436edfa17/tzlocal-5.3.1-py3-none-any.whl
- pypi: https://files.pythonhosted.org/packages/e5/4d/2e577f6db7aa0f932d19f799c18f604b2b302c65f733419b900ec07dbade/universal_pathlib-0.2.6-py3-none-any.whl
- pypi: https://files.pythonhosted.org/packages/a7/c2/fe1e52489ae3122415c51f387e221dd0773709bad6c6cdaa599e8a2c5185/urllib3-2.5.0-py3-none-any.whl
+ - pypi: https://files.pythonhosted.org/packages/39/ff/b279045420056a817049a8f745f22f68085a11f01ff221ea22bc0dafebee/warpkit-1.2.1-cp312-cp312-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl
- pypi: https://files.pythonhosted.org/packages/af/b5/123f13c975e9f27ab9c0770f514345bd406d0e8d3b7a0723af9d43f710af/wcwidth-0.2.14-py2.py3-none-any.whl
- pypi: https://files.pythonhosted.org/packages/9f/81/5d931d78d0eb732b95dc3ddaeeb71c8bb572fb01356e9133916cd729ecdd/wrapt-1.17.3-cp312-cp312-manylinux1_x86_64.manylinux_2_28_x86_64.manylinux_2_5_x86_64.whl
- pypi: https://files.pythonhosted.org/packages/d6/7d/b77455d7c7c51255b2992b429107fab811b2e36ceaf76da1e55a045dc568/xyzservices-2025.4.0-py3-none-any.whl
@@ -2134,8 +2134,8 @@ packages:
requires_python: '>=3.9'
- pypi: ./
name: fmriprep
- version: 26.0.0.dev55+g37e151d3e.d20260114
- sha256: e3e9590d09725c2cc97a110f6d30c2646d02408d0d8d056b158a175bb48c8398
+ version: 26.0.0.dev70+g43eecd295
+ sha256: f0ed814f6665164eb99485619272d0c63aefd77295c8998c30c456671fe5c29b
requires_dist:
- acres>=0.2.0
- apscheduler>=3.10
@@ -2172,6 +2172,7 @@ packages:
- sphinx-argparse!=0.5.0 ; extra == 'all'
- sphinx-rtd-theme>=0.5.2 ; extra == 'all'
- sphinx>=5 ; extra == 'all'
+ - warpkit>=1.2.1 ; python_full_version >= '3.11' and extra == 'all'
- datalad ; extra == 'container'
- datalad-osf ; extra == 'container'
- migas>=0.4.0 ; extra == 'container'
@@ -2201,6 +2202,80 @@ packages:
- pytest-env ; extra == 'tests'
- pytest-xdist>=2.5 ; extra == 'tests'
- pytest>=8.1 ; extra == 'tests'
+ - warpkit>=1.2.1 ; python_full_version >= '3.11' and extra == 'warpkit'
+ requires_python: '>=3.10'
+ editable: true
+- pypi: ./
+ name: fmriprep
+ version: 26.0.0.dev70+g43eecd295
+ sha256: f0ed814f6665164eb99485619272d0c63aefd77295c8998c30c456671fe5c29b
+ requires_dist:
+ - acres>=0.2.0
+ - apscheduler>=3.10
+ - codecarbon>=2
+ - looseversion>=1.3
+ - nibabel>=5.1.1
+ - nipype>=1.9.0
+ - nireports>=24.1.0
+ - nitime>=0.9
+ - nitransforms>=25.0.1
+ - niworkflows>=1.14.4
+ - numpy>=2.0
+ - packaging>=24
+ - pandas>=2.2
+ - psutil>=5.4
+ - pybids>=0.16
+ - requests>=2.27
+ - sdcflows>=2.15.0
+ - smriprep>=0.19.2
+ - tedana>=25.1.0
+ - templateflow>=24.2.2
+ - toml>=0.10
+ - transforms3d>=0.4.2
+ - coverage[toml]>=5.2.1 ; extra == 'all'
+ - fuzzywuzzy ; extra == 'all'
+ - migas>=0.4.0 ; extra == 'all'
+ - pydot>=1.2.3 ; extra == 'all'
+ - pytest-cov>=2.11 ; extra == 'all'
+ - pytest-env ; extra == 'all'
+ - pytest-xdist>=2.5 ; extra == 'all'
+ - pytest>=8.1 ; extra == 'all'
+ - python-levenshtein ; extra == 'all'
+ - sentry-sdk>=1.3 ; extra == 'all'
+ - sphinx-argparse!=0.5.0 ; extra == 'all'
+ - sphinx-rtd-theme>=0.5.2 ; extra == 'all'
+ - sphinx>=5 ; extra == 'all'
+ - warpkit>=1.2.1 ; python_full_version >= '3.11' and extra == 'all'
+ - datalad ; extra == 'container'
+ - datalad-osf ; extra == 'container'
+ - migas>=0.4.0 ; extra == 'container'
+ - sentry-sdk>=1.3 ; extra == 'container'
+ - pre-commit ; extra == 'dev'
+ - ruff ; extra == 'dev'
+ - pydot>=1.2.3 ; extra == 'doc'
+ - sphinx-argparse!=0.5.0 ; extra == 'doc'
+ - sphinx-rtd-theme>=0.5.2 ; extra == 'doc'
+ - sphinx>=5 ; extra == 'doc'
+ - pydot>=1.2.3 ; extra == 'docs'
+ - sphinx-argparse!=0.5.0 ; extra == 'docs'
+ - sphinx-rtd-theme>=0.5.2 ; extra == 'docs'
+ - sphinx>=5 ; extra == 'docs'
+ - duecredit ; extra == 'duecredit'
+ - fuzzywuzzy ; extra == 'maint'
+ - python-levenshtein ; extra == 'maint'
+ - migas>=0.4.0 ; extra == 'telemetry'
+ - sentry-sdk>=1.3 ; extra == 'telemetry'
+ - coverage[toml]>=5.2.1 ; extra == 'test'
+ - pytest-cov>=2.11 ; extra == 'test'
+ - pytest-env ; extra == 'test'
+ - pytest-xdist>=2.5 ; extra == 'test'
+ - pytest>=8.1 ; extra == 'test'
+ - coverage[toml]>=5.2.1 ; extra == 'tests'
+ - pytest-cov>=2.11 ; extra == 'tests'
+ - pytest-env ; extra == 'tests'
+ - pytest-xdist>=2.5 ; extra == 'tests'
+ - pytest>=8.1 ; extra == 'tests'
+ - warpkit>=1.2.1 ; python_full_version >= '3.11' and extra == 'warpkit'
requires_python: '>=3.10'
- conda: https://conda.anaconda.org/conda-forge/noarch/font-ttf-dejavu-sans-mono-2.37-hab24e00_0.tar.bz2
sha256: 58d7f40d2940dd0a8aa28651239adbf5613254df0f75789919c4e6762054403b
@@ -6696,6 +6771,18 @@ packages:
- pysocks>=1.5.6,!=1.5.7,<2.0 ; extra == 'socks'
- zstandard>=0.18.0 ; extra == 'zstd'
requires_python: '>=3.9'
+- pypi: https://files.pythonhosted.org/packages/39/ff/b279045420056a817049a8f745f22f68085a11f01ff221ea22bc0dafebee/warpkit-1.2.1-cp312-cp312-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl
+ name: warpkit
+ version: 1.2.1
+ sha256: ccbffe94bef9ebaecd507e39e8c8b59184da3745a3c3a2f506a5eb3341f028b1
+ requires_dist:
+ - nibabel>=4.0.2
+ - numpy>=1.23.3
+ - scikit-image>=0.20.0
+ - scipy>=1.8.1
+ - transforms3d>=0.4.1
+ - indexed-gzip>=1.7.0
+ requires_python: '>=3.11'
- conda: https://conda.anaconda.org/conda-forge/linux-64/wayland-1.24.0-h3e06ad9_0.conda
sha256: ba673427dcd480cfa9bbc262fd04a9b1ad2ed59a159bd8f7e750d4c52282f34c
md5: 0f2ca7906bf166247d1d760c3422cb8a
diff --git a/pyproject.toml b/pyproject.toml
index b6a9da3a8..d4a5813cf 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -79,6 +79,9 @@ telemetry = [
"migas >= 0.4.0",
"sentry-sdk >= 1.3",
]
+warpkit = [
+ "warpkit >= 1.2.1 ; python_version >= '3.11'",
+]
test = [
"coverage[toml] >= 5.2.1",
"pytest >= 8.1",
@@ -93,7 +96,7 @@ maint = [
# Aliases
docs = ["fmriprep[doc]"]
tests = ["fmriprep[test]"]
-all = ["fmriprep[doc,maint,telemetry,test]"]
+all = ["fmriprep[doc,maint,telemetry,test,warpkit]"]
[project.scripts]
fmriprep = "fmriprep.cli.run:main"
@@ -133,9 +136,9 @@ FSLDIR = "$CONDA_PREFIX"
[tool.pixi.environments]
# Local installs will use editable
default = { features = ["editable"], solve-group = "default" }
-test = { features = ["editable", "test"], solve-group = "default" }
+test = { features = ["editable", "test", "warpkit"], solve-group = "default" }
# Use "fmriprep" feature for production environment, so it is identifiable
-fmriprep = { features = ["production", "container", "telemetry"], solve-group = "production" }
+fmriprep = { features = ["production", "container", "telemetry", "warpkit"], solve-group = "production" }
[tool.pixi.feature.editable.pypi-dependencies]
fmriprep = { path = ".", editable = true }