ENH: Add optional warpkit MEDIC support for multi-echo BOLD distortion correction

.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                                            |

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,15 +53,15 @@ 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
 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

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(

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'),
     [

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:

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},

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 
-      <em>Preprocessed estimation with varying Phase-Encoding (PE) blips</em> 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 <em>B<sub>0</sub></em> map in the corresponding run's space (debug mode)"
+    subtitle: "Estimated <em>B<sub>0</sub></em> 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

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

fmriprep/interfaces/resampling.py

@@ -378,14 +378,23 @@ 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,
             coordinates,
             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,

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]]])

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)