From a3660c860d21781611d419e2b06b74a7c40b17ba Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Mon, 4 May 2026 06:45:53 -0400 Subject: [PATCH 01/14] Add optional warpkit MEDIC support for multiecho BOLD --- Dockerfile | 5 +- fmriprep/cli/parser.py | 10 ++ fmriprep/cli/tests/test_parser.py | 13 ++ fmriprep/config.py | 2 + fmriprep/interfaces/maths.py | 31 ++++ fmriprep/interfaces/resampling.py | 13 +- fmriprep/interfaces/tests/test_maths.py | 20 ++- fmriprep/utils/bids.py | 37 +++++ fmriprep/utils/tests/test_derivative_cache.py | 29 ++++ fmriprep/workflows/base.py | 23 ++- fmriprep/workflows/bold/base.py | 6 + fmriprep/workflows/bold/fit.py | 150 +++++++++++++++--- fmriprep/workflows/bold/tests/test_fit.py | 32 ++++ pyproject.toml | 9 +- 14 files changed, 347 insertions(+), 33 deletions(-) diff --git a/Dockerfile b/Dockerfile index f2f17011d..0e484f773 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 @@ -51,7 +52,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 +60,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..e7c36f870 100644 --- a/fmriprep/cli/parser.py +++ b/fmriprep/cli/parser.py @@ -485,6 +485,16 @@ 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_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..f13552b8c 100644 --- a/fmriprep/cli/tests/test_parser.py +++ b/fmriprep/cli/tests/test_parser.py @@ -199,6 +199,19 @@ 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'] + bids_path.mkdir() + + parser = _build_parser() + opts = parser.parse_args(args) + + assert opts.me_use_warpkit is True + _reset_config() + + @pytest.mark.parametrize( ('args', 'expectation'), [ diff --git a/fmriprep/config.py b/fmriprep/config.py index c91951285..9ccabbe2d 100644 --- a/fmriprep/config.py +++ b/fmriprep/config.py @@ -642,6 +642,8 @@ 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.""" class loggers: 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/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..367bb49c0 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -57,6 +57,7 @@ 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 +87,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 ------ @@ -246,6 +249,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 +291,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 +307,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..4ca8d25d5 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,32 @@ 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 hmc_boldref = precomputed.get('hmc_boldref') coreg_boldref = precomputed.get('coreg_boldref') @@ -281,6 +307,7 @@ def init_bold_fit_wf( 'motion_xfm', 'boldref2anat_xfm', 'boldref2fmap_xfm', + 'fieldmap', ], ), name='outputnode', @@ -314,7 +341,7 @@ def init_bold_fit_wf( if coreg_boldref: regref_buffer.inputs.boldref = coreg_boldref config.loggers.workflow.debug(f'Reusing coregistration reference: {coreg_boldref}') - fmapref_buffer.inputs.sbref_files = sbref_files + fmapref_buffer.inputs.sbref_files = [] if warpkit_enabled else sbref_files summary = pe.Node( FunctionalSummary( @@ -338,6 +365,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 +378,7 @@ 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 and not warpkit_enabled, freesurfer=config.workflow.run_reconall, output_dir=config.execution.fmriprep_dir, ) @@ -463,8 +492,49 @@ 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, + ), + name='warpkit_medic', + mem_gb=mem_gb['resampled'], + ) + 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'), + ]), + (warpkit_medic, warpkit_fieldmap, [('fieldmap_native', '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')]), + ]) # fmt:skip + elif fieldmap_id: config.loggers.workflow.info('Stage 3: Adding fieldmap reconstruction workflow') fmap_select = pe.Node( KeySelect( @@ -550,7 +620,7 @@ def init_bold_fit_wf( config.loggers.workflow.info('Stage 4: Adding coregistration boldref workflow') # If sbref files are available, add them to the list of sources - if sbref_files and nb.load(sbref_files[0]).ndim > 3: + if sbref_files and not warpkit_enabled and nb.load(sbref_files[0]).ndim > 3: raw_sbref_wf = init_raw_boldref_wf( name='raw_sbref_wf', bold_file=sbref_files[0], @@ -586,7 +656,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 +814,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 +847,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 +869,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 +913,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 +950,7 @@ def init_bold_native_wf( 'fmap_ref', 'fmap_coeff', 'fmap_id', + 'fieldmap', ], ), name='inputnode', @@ -900,13 +1005,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 +1015,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 +1056,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/tests/test_fit.py b/fmriprep/workflows/bold/tests/test_fit.py index 31bdfde69..6f2ea6bdb 100644 --- a/fmriprep/workflows/bold/tests/test_fit.py +++ b/fmriprep/workflows/bold/tests/test_fit.py @@ -175,3 +175,35 @@ 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): + 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) + + with mock_config(bids_dir=bids_root): + 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()) diff --git a/pyproject.toml b/pyproject.toml index 1c446a2f2..5ebe9be04 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 } From 2e2ad294ddec43f7ef3c9de60efaf0bf47adb42f Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Mon, 4 May 2026 08:36:42 -0400 Subject: [PATCH 02/14] Fix PIXI_LOCK_FLAGS scope in Docker build --- Dockerfile | 1 + 1 file changed, 1 insertion(+) diff --git a/Dockerfile b/Dockerfile index 0e484f773..af8730f37 100644 --- a/Dockerfile +++ b/Dockerfile @@ -40,6 +40,7 @@ ARG PIXI_LOCK_FLAGS=--frozen # - ... # 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 \ From 3933446cb3720beada9b268cb6eefbd6a35b9448 Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Mon, 4 May 2026 11:06:52 -0400 Subject: [PATCH 03/14] Add missing warpkit interface files --- fmriprep/interfaces/tests/test_resampling.py | 30 ++++++++ fmriprep/interfaces/warpkit.py | 78 ++++++++++++++++++++ 2 files changed, 108 insertions(+) create mode 100644 fmriprep/interfaces/tests/test_resampling.py create mode 100644 fmriprep/interfaces/warpkit.py 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..2c634b08e --- /dev/null +++ b/fmriprep/interfaces/warpkit.py @@ -0,0 +1,78 @@ +from __future__ import annotations + +from pathlib import Path + +from nipype.interfaces.base import File, InputMultiObject, SimpleInterface, TraitedSpec, traits + + +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') + 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, + wrap_limit=self.inputs.wrap_limit, + debug=self.inputs.debug, + ) + + self._results['fieldmap_native'] = str(result.fieldmap_native) + self._results['displacement_map'] = str(result.displacement_map) + self._results['fieldmap'] = str(result.fieldmap) + return runtime From 41769d6883b5c11de5f746df3142af8044247d25 Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Mon, 4 May 2026 19:28:36 -0400 Subject: [PATCH 04/14] Add warpkit MEDIC reportlets and citation boilerplate --- fmriprep/data/boilerplate.bib | 10 +++ fmriprep/data/reports-spec-func.yml | 15 ++--- fmriprep/workflows/bold/base.py | 32 +++++++--- fmriprep/workflows/bold/fit.py | 6 +- fmriprep/workflows/bold/outputs.py | 73 ++++++++++++++++------ fmriprep/workflows/bold/tests/test_base.py | 9 ++- fmriprep/workflows/bold/tests/test_fit.py | 28 ++++++++- 7 files changed, 133 insertions(+), 40 deletions(-) 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/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 367bb49c0..56eed1c6d 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -52,6 +52,29 @@ 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 estimates +a framewise B0 fieldmap from the phase evolution across echoes at each +repetition time, and 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], @@ -190,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( diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 4ca8d25d5..0eb427263 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -378,7 +378,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 and not warpkit_enabled, + 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, ) @@ -533,6 +535,8 @@ def init_bold_fit_wf( (hmc_buffer, warpkit_fieldmap, [('hmc_xforms', 'transforms')]), (warpkit_fieldmap, warpkit_mean, [('out_file', 'in_file')]), (warpkit_fieldmap, outputnode, [('out_file', 'fieldmap')]), + (warpkit_fieldmap, 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') diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index df8de999c..58c10c321 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', @@ -300,7 +302,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, @@ -348,7 +350,57 @@ 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, fmap_boldref, [ + ('fieldmap', 'input_image'), + ('coreg_boldref', 'reference_image'), + ('boldref2fmap_xfm', 'transforms'), + ]), + (inputnode, sdcreg_report, [ + ('sdc_boldref', 'reference'), + ('bold_mask', 'mask'), + ]), + (fmapref_boldref, sdcreg_report, [('output_image', 'moving')]), + (fmap_boldref, sdcreg_report, [('output_image', 'fieldmap')]), + (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, [ + ('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', @@ -372,23 +424,6 @@ def init_func_fit_reports_wf( ) workflow.connect([ - (inputnode, fmapref_boldref, [ - ('fmap_ref', 'input_image'), - ('coreg_boldref', 'reference_image'), - ('boldref2fmap_xfm', 'transforms'), - ]), - (inputnode, fmap_boldref, [ - ('fieldmap', 'input_image'), - ('coreg_boldref', 'reference_image'), - ('boldref2fmap_xfm', 'transforms'), - ]), - (inputnode, sdcreg_report, [ - ('sdc_boldref', 'reference'), - ('bold_mask', 'mask'), - ]), - (fmapref_boldref, sdcreg_report, [('output_image', 'moving')]), - (fmap_boldref, sdcreg_report, [('output_image', 'fieldmap')]), - (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..2e30b7ae1 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,10 @@ 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 'van2023medic' in postdesc diff --git a/fmriprep/workflows/bold/tests/test_fit.py b/fmriprep/workflows/bold/tests/test_fit.py index 6f2ea6bdb..9609a8d0f 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( @@ -177,7 +190,7 @@ def test_bold_native_precomputes( generate_expanded_graph(flatgraph) -def test_bold_warpkit_graph(bids_root: Path, tmp_path: Path): +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') @@ -191,7 +204,13 @@ def test_bold_warpkit_graph(bids_root: Path, tmp_path: Path): 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' fit_wf = init_bold_fit_wf( bold_series=bold_series, fieldmap_id=None, @@ -207,3 +226,10 @@ def test_bold_warpkit_graph(bids_root: Path, tmp_path: Path): generate_expanded_graph(fit_wf._create_flat_graph()) generate_expanded_graph(native_wf._create_flat_graph()) + + fit_nodes = fit_wf.list_node_names() + 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) From eff7d66c82820796ac4491fcb25ddebc1d0e3efc Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Thu, 7 May 2026 11:05:02 -0400 Subject: [PATCH 05/14] Polish warpkit MEDIC reporting and boilerplate --- fmriprep/workflows/bold/base.py | 10 +++++----- fmriprep/workflows/bold/outputs.py | 1 + fmriprep/workflows/bold/tests/test_base.py | 1 + 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 56eed1c6d..ebdc85c06 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -66,11 +66,11 @@ def _build_postdesc(use_warpkit: bool) -> str: 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 estimates -a framewise B0 fieldmap from the phase evolution across echoes at each -repetition time, and the resulting fieldmap series was aligned to the -head-motion-correction reference and applied during native-space resampling -together with motion 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 diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index 58c10c321..b1dc876f9 100644 --- a/fmriprep/workflows/bold/outputs.py +++ b/fmriprep/workflows/bold/outputs.py @@ -394,6 +394,7 @@ def init_func_fit_reports_wf( workflow.connect([ (inputnode, fieldmap_report, [ + ('bold_mask', 'mask'), ('sdc_boldref', 'reference'), ('fieldmap', 'fieldmap'), ]), diff --git a/fmriprep/workflows/bold/tests/test_base.py b/fmriprep/workflows/bold/tests/test_base.py index 2e30b7ae1..7ceea1719 100644 --- a/fmriprep/workflows/bold/tests/test_base.py +++ b/fmriprep/workflows/bold/tests/test_base.py @@ -87,4 +87,5 @@ 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 From 43eecd2959f94a436c10e7b28dff2e3acf77784a Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Thu, 7 May 2026 20:02:29 -0400 Subject: [PATCH 06/14] MAINT: Add David V. Smith to contributors --- .maint/CONTRIBUTORS.md | 1 + 1 file changed, 1 insertion(+) 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 | From 66baf820bb31aefd05cc85a5b2ee4711cca569cc Mon Sep 17 00:00:00 2001 From: "nipreps[bot]" Date: Fri, 8 May 2026 01:35:16 +0000 Subject: [PATCH 07/14] [DATALAD RUNCMD] bash -c '! pixi lock --check || git chec... === Do not change lines below === { "chain": [], "cmd": "bash -c '! pixi lock --check || git checkout .'", "exit": 0, "extra_inputs": [], "inputs": [ "pixi.lock", "pyproject.toml" ], "outputs": [ "pixi.lock" ], "pwd": "." } ^^^ Do not change lines above ^^^ --- pixi.lock | 95 ++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 91 insertions(+), 4 deletions(-) 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 From 4eb308f034ad8986923a0d9244fa2ec7cb261d51 Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Sun, 10 May 2026 22:35:06 -0400 Subject: [PATCH 08/14] ENH: Expose warpkit MEDIC noise-frame trimming --- fmriprep/cli/parser.py | 19 +++++++ fmriprep/cli/tests/test_parser.py | 10 +++- fmriprep/config.py | 2 + fmriprep/interfaces/warpkit.py | 65 +++++++++++++++++++++-- fmriprep/workflows/bold/fit.py | 1 + fmriprep/workflows/bold/tests/test_fit.py | 2 + 6 files changed, 95 insertions(+), 4 deletions(-) diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py index e7c36f870..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) @@ -495,6 +503,17 @@ def _fallback_trt(value, parser): '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 f13552b8c..09c7376e9 100644 --- a/fmriprep/cli/tests/test_parser.py +++ b/fmriprep/cli/tests/test_parser.py @@ -202,13 +202,21 @@ def test_slice_time_ref(tmp_path, st_ref): 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'] + 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() diff --git a/fmriprep/config.py b/fmriprep/config.py index 9ccabbe2d..21030cef4 100644 --- a/fmriprep/config.py +++ b/fmriprep/config.py @@ -644,6 +644,8 @@ class workflow(_Config): """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/interfaces/warpkit.py b/fmriprep/interfaces/warpkit.py index 2c634b08e..221629381 100644 --- a/fmriprep/interfaces/warpkit.py +++ b/fmriprep/interfaces/warpkit.py @@ -2,9 +2,48 @@ 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( @@ -34,6 +73,7 @@ class WarpkitMEDICInputSpec(TraitedSpec): 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') @@ -68,11 +108,30 @@ def _run_interface(self, runtime): 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, ) - self._results['fieldmap_native'] = str(result.fieldmap_native) - self._results['displacement_map'] = str(result.displacement_map) - self._results['fieldmap'] = str(result.fieldmap) + 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/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 0eb427263..cf0a66e28 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -513,6 +513,7 @@ def init_bold_fit_wf( phase=phase_files, tes=tes_ms, n_cpus=omp_nthreads, + noise_frames=config.workflow.me_warpkit_noise_frames, ), name='warpkit_medic', mem_gb=mem_gb['resampled'], diff --git a/fmriprep/workflows/bold/tests/test_fit.py b/fmriprep/workflows/bold/tests/test_fit.py index 9609a8d0f..a1b55bd8c 100644 --- a/fmriprep/workflows/bold/tests/test_fit.py +++ b/fmriprep/workflows/bold/tests/test_fit.py @@ -211,6 +211,7 @@ def test_bold_warpkit_graph(bids_root: Path, tmp_path: Path, monkeypatch): 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, @@ -228,6 +229,7 @@ def test_bold_warpkit_graph(bids_root: Path, tmp_path: Path, monkeypatch): 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) From 89f431d4086f36b38f0730c4b14d50ebad7c479e Mon Sep 17 00:00:00 2001 From: Andrew Van Date: Sun, 10 May 2026 22:46:13 -0400 Subject: [PATCH 09/14] FIX: feed undistorted MEDIC fieldmap (`fmap`) to the resampler MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `warpkit.api.medic` exposes three outputs: `fieldmap` (undistorted-space, the consumer-facing output), `fieldmap_native` (a distorted-space debug output), and `displacement_map`. The Stage-3 wiring used `fieldmap_native`, which makes fmriprep's pull-resampler look up the field at the wrong physical location at each target voxel (it's off by the SDC displacement itself). `resample_vol` documents fmap_hz as "sampled to the target space, in Hz" (interfaces/resampling.py:272-273) — i.e., the undistorted/output grid. That's `fieldmap`, not `fieldmap_native`. Verified empirically by running `resample_series_async` against both fields on two datasets: * warpkit's 3.5mm fixture: RMS Δ ≈ 358, max|Δ| ≈ 17.5k vs mean signal ~7050 * ds007637 sub-04 at 2mm: RMS Δ ≈ 211, max|Δ| ≈ 21k vs mean signal ~4175 Relative error grows with resolution, as expected — a fixed physical displacement spans more voxels in a finer grid. --- fmriprep/workflows/bold/fit.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index cf0a66e28..55af2a24e 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -531,7 +531,13 @@ def init_bold_fit_wf( ('readout_time', 'total_readout_time'), ('pe_direction', 'phase_encoding_direction'), ]), - (warpkit_medic, warpkit_fieldmap, [('fieldmap_native', 'in_file')]), + # `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')]), From 6c590002eb9f1ee782ff55b8b56471c6c36ecace Mon Sep 17 00:00:00 2001 From: Andrew Van Date: Sun, 10 May 2026 22:46:41 -0400 Subject: [PATCH 10/14] FIX: feed 3D temporal-mean field to FieldmapReportlet `func_fit_reports_wf.inputnode.fieldmap` flows into FieldmapReportlet (workflows/bold/outputs.py:373-401), which expects a 3D static field. The Stage-3 warpkit branch was passing the 4D HMC-aligned per-frame field directly, so the reportlet either crashes on the 4D NIfTI or silently renders frame 0 only. Pipe the temporal mean (already computed by `warpkit_mean` on line 526) instead. The 4D series remains the outputnode.fieldmap that the bold native workflow consumes for per-frame correction. --- fmriprep/workflows/bold/fit.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 55af2a24e..867227352 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -542,7 +542,9 @@ def init_bold_fit_wf( (hmc_buffer, warpkit_fieldmap, [('hmc_xforms', 'transforms')]), (warpkit_fieldmap, warpkit_mean, [('out_file', 'in_file')]), (warpkit_fieldmap, outputnode, [('out_file', 'fieldmap')]), - (warpkit_fieldmap, func_fit_reports_wf, [('out_file', 'inputnode.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: From 81ca35b9c618363207e594aa13c655554017b87b Mon Sep 17 00:00:00 2001 From: Andrew Van Date: Sun, 10 May 2026 22:47:19 -0400 Subject: [PATCH 11/14] FIX: scale warpkit_medic mem_gb by echo count MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `mem_gb['resampled']` is the memory footprint of one resampled BOLD volume (`filesize * 4` in `estimate_bold_mem_usage`). MEDIC actually holds all N echoes of phase + magnitude in memory as float64, plus warpkit's internal state (unwrapped phases, masks, per-frame fields). For typical multi-echo data this underestimates real usage by an order of magnitude. ds007637 sub-04 is 5 echoes × ~320 MB compressed — close to 3 GB resident before warpkit's working set. The nipype scheduler will happily oversubscribe and OOM. Use `2 * n_echoes * filesize` instead, which roughly doubles the input footprint to cover outputs and warpkit's internal arrays. --- fmriprep/workflows/bold/fit.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 867227352..027dd681f 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -516,7 +516,11 @@ def init_bold_fit_wf( noise_frames=config.workflow.me_warpkit_noise_frames, ), name='warpkit_medic', - mem_gb=mem_gb['resampled'], + # 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), From 197bc44c16c67419043fe2d81891c35af24de10c Mon Sep 17 00:00:00 2001 From: Andrew Van Date: Sun, 10 May 2026 22:47:56 -0400 Subject: [PATCH 12/14] FIX: stop silently disabling SBRef when MEDIC is enabled MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Stage-3 warpkit branch cleared `fmapref_buffer.sbref_files` and gated `raw_sbref_wf` with `not warpkit_enabled`. SBRef is acquired separately from the BOLD run and is orthogonal to MEDIC — dropping it silently is a regression for sites that acquire it. If there's a concrete incompatibility between SBRef and the MEDIC field-application path, it should be a stated limitation with a clear error, not an invisible behavior change. Removing the gate so SBRef flows through the same way it does on every other distortion-correction path. --- fmriprep/workflows/bold/fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 027dd681f..e4d1b5160 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -341,7 +341,7 @@ def init_bold_fit_wf( if coreg_boldref: regref_buffer.inputs.boldref = coreg_boldref config.loggers.workflow.debug(f'Reusing coregistration reference: {coreg_boldref}') - fmapref_buffer.inputs.sbref_files = [] if warpkit_enabled else sbref_files + fmapref_buffer.inputs.sbref_files = sbref_files summary = pe.Node( FunctionalSummary( @@ -637,7 +637,7 @@ def init_bold_fit_wf( config.loggers.workflow.info('Stage 4: Adding coregistration boldref workflow') # If sbref files are available, add them to the list of sources - if sbref_files and not warpkit_enabled and nb.load(sbref_files[0]).ndim > 3: + if sbref_files and nb.load(sbref_files[0]).ndim > 3: raw_sbref_wf = init_raw_boldref_wf( name='raw_sbref_wf', bold_file=sbref_files[0], From 53d8fc7593de41ac00059b28101c0a5ee336b4f9 Mon Sep 17 00:00:00 2001 From: Andrew Van Date: Sun, 10 May 2026 22:49:06 -0400 Subject: [PATCH 13/14] FIX: assert echo ordering matches between magnitude and phase `bold_series` is sorted by EchoTime in `workflows/base.py:250` and `phase_files` is independently sorted by EchoTime in `collect_bold_part_files`. They *should* line up, but a missing or zero-default EchoTime sidecar on one side and not the other silently transposes the pair set fed to warpkit. Add an explicit cross-check so that misalignment fails fast with a readable message instead of producing garbage fieldmaps. --- fmriprep/workflows/bold/fit.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index e4d1b5160..9468ba795 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -260,6 +260,22 @@ def init_bold_fit_wf( 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. Assert + # the orderings match before handing them off. + assert layout is not None # narrowed for type-checker; warpkit needs BIDS + 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') # Can contain From a57f163716244a2e49d85fb6c0f5e3fee7bb1ac1 Mon Sep 17 00:00:00 2001 From: "David V. Smith" Date: Mon, 18 May 2026 08:08:45 -0400 Subject: [PATCH 14/14] FIX: Avoid assert in MEDIC echo-order check --- fmriprep/workflows/bold/fit.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 9468ba795..9c54ac351 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -263,9 +263,11 @@ def init_bold_fit_wf( # `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. Assert + # side without the other, feeding warpkit a transposed pair set. Check # the orderings match before handing them off. - assert layout is not None # narrowed for type-checker; warpkit needs BIDS + 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 ]