ENH: *fieldmap-less* integration from *SDCFlows* 2.0

.circleci/config.yml

@@ -315,6 +315,7 @@ jobs:
                 --user $(id -u):$(id -g) \
                 nipreps/dmriprep:latest /data /out participant -vv $eddy \
                 --fs-subjects-dir /data/derivatives/freesurfer-6.0.1 --sloppy \
+                --output-spaces MNI152NLin2009cAsym --use-syn-sdc \
                 --notrack --skip-bids-validation -w /work --omp-nthreads 2 --nprocs 2
       - store_artifacts:
           path: /tmp/ds000206/derivatives/dmriprep

dmriprep/cli/parser.py

@@ -225,7 +225,7 @@ def _bids_filter(value):
         default="register",
         choices=["register", "header"],
         help='Either "register" (the default) to initialize volumes at center or "header"'
-        " to use the header information when coregistering DWI to T1w images.",
+        " to use the header information when co-registering DWI to T1w images.",
     )
 
     #  ANTs options
@@ -243,35 +243,27 @@ def _bids_filter(value):
         "run-to-run replicability when used with --omp-nthreads 1",
     )
 
-    # Fieldmap options
-    g_fmap = parser.add_argument_group("Specific options for handling fieldmaps")
-    g_fmap.add_argument(
-        "--fmap-bspline",
-        action="store_true",
-        default=False,
-        help="fit a B-Spline field using least-squares (experimental)",
-    )
-    g_fmap.add_argument(
-        "--fmap-no-demean",
-        action="store_false",
-        default=True,
-        help="do not remove median (within mask) from fieldmap",
-    )
-
     # SyN-unwarp options
     g_syn = parser.add_argument_group("Specific options for SyN distortion correction")
     g_syn.add_argument(
         "--use-syn-sdc",
         action="store_true",
+        dest="use_syn",
         default=False,
-        help="EXPERIMENTAL: Use fieldmap-free distortion correction",
+        help="""\
+Attempt to set-up fieldmap-less estimation of fieldmaps via nonlinear registration with ANTs \
+if no other fieldmap estimation method is available. Fieldmap-less estimation will not be used \
+when sufficient fieldmap information (B0 mapping with SEI or GRE, or PEPOLAR estimation with \
+EPIs) is retrieved from the BIDS structure for a given subject.
+""",
     )
     g_syn.add_argument(
         "--force-syn",
         action="store_true",
         default=False,
-        help="EXPERIMENTAL/TEMPORARY: Use SyN correction in addition to "
-        "fieldmap correction, if available",
+        help="""\
+Force estimation of fieldmaps with the fieldmap-less approach. The use of this feature \
+is discouraged.""",
     )
 
     # FreeSurfer options
@@ -334,7 +326,7 @@ def _bids_filter(value):
         action="store_true",
         default=False,
         help="only generate reports, don't run workflows. This will only rerun report "
-        "aggregation, not reportlet generation for specific nodes.",
+        "aggregation, not 'reportlet' generation for specific nodes.",
     )
     g_other.add_argument(
         "--run-uuid",

dmriprep/config/reports-spec.yml

@@ -26,7 +26,7 @@ sections:
     caption: Surfaces (white and pial) reconstructed with FreeSurfer (<code>recon-all</code>)
       overlaid on the participant's T1w template.
     subtitle: Surface reconstruction
-- name: Fieldmaps
+- name: <em>B<sub>0</sub></em> field mapping
   ordering: session,run,fmapid
   reportlets:
   - bids: {datatype: figures, desc: mapped, suffix: fieldmap}
@@ -38,7 +38,7 @@ sections:
     description: Hover over the panels with the mouse pointer to also visualize the intensity of the
       field inhomogeneity in Hertz.
     static: false
-    subtitle: "Susceptibility-derived Distortion Correction (SDC): field inhomogeneity estimation"
+    subtitle: "Preprocessed <em>B<sub>0</sub></em> mapping acquisition"
   - bids: {datatype: figures, desc: phasediff, suffix: fieldmap}
     caption: Inhomogeneities of the <em>B<sub>0</sub></em> field introduce (oftentimes severe) spatial distortions
       along the phase-encoding direction of the image. A Gradient-Recalled Echo (GRE) scheme for the 
@@ -48,7 +48,7 @@ sections:
     description: Hover over the panels with the mouse pointer to also visualize the intensity of the
       field inhomogeneity in Hertz.
     static: false
-    subtitle: "Susceptibility-derived Distortion Correction (SDC): field inhomogeneity estimation"
+    subtitle: "Preprocessed mapping of phase-difference acquisition"
   - bids: {datatype: figures, desc: pepolar, suffix: fieldmap}
     caption: Inhomogeneities of the <em>B<sub>0</sub></em> field introduce (oftentimes severe) spatial distortions
       along the phase-encoding direction of the image. Utilizing two or more images with different
@@ -58,7 +58,17 @@ sections:
     description: Hover on the panels with the mouse pointer to also visualize the intensity of the
       inhomogeneity of the field in Hertz.
     static: false
-    subtitle: "Susceptibility-derived Distortion Correction (SDC): field inhomogeneity estimation"
+    subtitle: "Preprocessed estimation with varying Phase-Endocing (PE) blips"
+  - bids: {datatype: figures, desc: anat, suffix: fieldmap}
+    caption: Inhomogeneities of the <em>B<sub>0</sub></em> field introduce (oftentimes severe) spatial distortions
+      along the phase-encoding direction of the image. Utilizing an <em>anatomically-correct</em> acquisition
+      (for instance, T1w or T2w), it is possible to estimate the inhomogeneity of the field by means of nonlinear
+      registration. The plot below shows a reference EPI (echo-planar imaging) volume generated
+      using two or more EPI images with the same PE encoding, after alignment to the anatomical scan.
+    description: Hover on the panels with the mouse pointer to also visualize the intensity of the
+      inhomogeneity of the field in Hertz.
+    static: false
+    subtitle: "Preprocessed estimation by nonlinear registration to an anatomical scan (&ldquo;<em>fieldmap-less</em>&rdquo;)"
 - name: Diffusion
   ordering: session,acquisition,run
   reportlets:

dmriprep/interfaces/vectors.py

@@ -30,6 +30,7 @@ class _CheckGradientTableOutputSpec(TraitedSpec):
     full_sphere = traits.Bool()
     pole = traits.Tuple(traits.Float, traits.Float, traits.Float)
     b0_ixs = traits.List(traits.Int)
+    b0_mask = traits.List(traits.Bool)
 
 
 class CheckGradientTable(SimpleInterface):
@@ -81,6 +82,7 @@ def _run_interface(self, runtime):
         pole = table.pole
         self._results["pole"] = tuple(pole)
         self._results["full_sphere"] = np.all(pole == 0.0)
+        self._results["b0_mask"] = table.b0mask.tolist()
         self._results["b0_ixs"] = np.where(table.b0mask)[0].tolist()
 
         cwd = Path(runtime.cwd).absolute()

dmriprep/workflows/base.py

@@ -313,13 +313,19 @@ def init_single_subject_wf(subject_id):
         fmap_estimators = find_estimators(
             layout=config.execution.layout,
             subject=subject_id,
-            fmapless=False,
+            fmapless=config.workflow.use_syn,
+            force_fmapless=config.workflow.force_syn,
         )
 
-        # Add fieldmap-less estimators
-        if not fmap_estimators and config.workflow.use_syn:
-            # estimators = [fm.FieldmapEstimation()]
-            raise NotImplementedError
+        if (
+            any(f.method == fm.EstimatorType.ANAT for f in fmap_estimators)
+            and "MNI152NLin2009cAsym" not in spaces.get_spaces(nonstandard=False, dim=(3,))
+        ):
+            # Although this check would go better within parser, allow datasets with fieldmaps
+            # not to require spatial standardization of the T1w.
+            raise RuntimeError("""\
+Argument '--use-sdc-syn' requires having 'MNI152NLin2009cAsym' as one output standard space. \
+Please add the 'MNI152NLin2009cAsym' keyword to the '--output-spaces' argument""")
 
     # Nuts and bolts: initialize individual run's pipeline
     dwi_preproc_list = []
@@ -354,11 +360,15 @@ def init_single_subject_wf(subject_id):
         dwi_preproc_list.append(dwi_preproc_wf)
 
     if not fmap_estimators:
+        config.loggers.workflow.warning(
+            "Data for fieldmap estimation not present. Please note that these data "
+            "will not be corrected for susceptibility distortions."
+        )
         return workflow
 
     # SDC Step 2: Manually add further estimators (e.g., fieldmap-less)
     fmap_wf = init_fmap_preproc_wf(
-        debug=config.execution.debug,
+        debug=config.execution.debug is True,
         estimators=fmap_estimators,
         omp_nthreads=config.nipype.omp_nthreads,
         output_dir=str(output_dir),
@@ -371,7 +381,6 @@ def init_single_subject_wf(subject_id):
 BIDS structure for this particular subject.
 """
 
-    # TODO: Requires nipreps/sdcflows#147
     for dwi_preproc_wf in dwi_preproc_list:
         # fmt: off
         workflow.connect([
@@ -392,38 +401,70 @@ def init_single_subject_wf(subject_id):
 
     # Step 3: Manually connect PEPOLAR
     for estimator in fmap_estimators:
-        if estimator.method != fm.EstimatorType.PEPOLAR:
+        config.loggers.workflow.info(f"""\
+Setting-up fieldmap "{estimator.bids_id}" ({estimator.method}) with \
+<{', '.join(s.path.name for s in estimator.sources)}>""")
+        if estimator.method in (fm.EstimatorType.MAPPED, fm.EstimatorType.PHASEDIFF):
             continue
 
         suffices = set(s.suffix for s in estimator.sources)
 
-        if sorted(suffices) == ["epi"]:
+        if estimator.method == fm.EstimatorType.PEPOLAR and sorted(suffices) == ["epi"]:
             getattr(fmap_wf.inputs, f"in_{estimator.bids_id}").in_data = [
                 str(s.path) for s in estimator.sources
             ]
             getattr(fmap_wf.inputs, f"in_{estimator.bids_id}").metadata = [
                 s.metadata for s in estimator.sources
             ]
-        else:
-            raise NotImplementedError
-            # from niworkflows.interfaces.utility import KeySelect
-            # est_id = estimator.bids_id
-            # estim_select = pe.MapNode(
-            #     KeySelect(fields=["metadata", "dwi_reference", "dwi_mask", "gradients_rasb",]),
-            #     name=f"fmap_select_{est_id}",
-            #     run_without_submitting=True,
-            #     iterfields=["key"]
-            # )
-            # estim_select.inputs.key = [
-            #     str(s.path) for s in estimator.sources if s.suffix in ("epi", "dwi", "sbref")
-            # ]
-            # # fmt:off
-            # workflow.connect([
-            #     (referencenode, estim_select, [("dwi_file", "keys"),
-            #                                   ("metadata", "metadata"),
-            #                                   ("dwi_reference", "dwi_reference"),
-            #                                   ("gradients_rasb", "gradients_rasb")]),
-            # ])
-            # # fmt:on
+            continue
+
+        if estimator.method == fm.EstimatorType.PEPOLAR:
+            raise NotImplementedError(
+                "Sophisticated PEPOLAR schemes (e.g., using DWI+EPI) are unsupported."
+            )
+
+        if estimator.method == fm.EstimatorType.ANAT:
+            from sdcflows.workflows.fit.syn import init_syn_preprocessing_wf
+            from ..interfaces.vectors import CheckGradientTable
+
+            sources = [
+                str(s.path) for s in estimator.sources
+                if s.suffix in ("dwi",)
+            ]
+            layout = config.execution.layout
+            syn_preprocessing_wf = init_syn_preprocessing_wf(
+                omp_nthreads=config.nipype.omp_nthreads,
+                debug=config.execution.debug is True,
+                auto_bold_nss=False,
+                t1w_inversion=True,
+                name=f"syn_preprocessing_{estimator.bids_id}",
+            )
+            syn_preprocessing_wf.inputs.inputnode.in_epis = sources
+            syn_preprocessing_wf.inputs.inputnode.in_meta = [
+                layout.get_metadata(s) for s in sources
+            ]
+            b0_masks = pe.MapNode(CheckGradientTable(), name=f"b0_masks_{estimator.bids_id}",
+                                  iterfield=("dwi_file", "in_bvec", "in_bval"))
+            b0_masks.inputs.dwi_file = sources
+            b0_masks.inputs.in_bvec = [str(layout.get_bvec(s)) for s in sources]
+            b0_masks.inputs.in_bval = [str(layout.get_bval(s)) for s in sources]
+
+            # fmt:off
+            workflow.connect([
+                (anat_preproc_wf, syn_preprocessing_wf, [
+                    ("outputnode.t1w_preproc", "inputnode.in_anat"),
+                    ("outputnode.t1w_mask", "inputnode.mask_anat"),
+                    ("outputnode.std2anat_xfm", "inputnode.std2anat_xfm"),
+                ]),
+                (b0_masks, syn_preprocessing_wf, [("b0_mask", "inputnode.t_masks")]),
+                (syn_preprocessing_wf, fmap_wf, [
+                    ("outputnode.epi_ref", f"in_{estimator.bids_id}.epi_ref"),
+                    ("outputnode.epi_mask", f"in_{estimator.bids_id}.epi_mask"),
+                    ("outputnode.anat_ref", f"in_{estimator.bids_id}.anat_ref"),
+                    ("outputnode.anat_mask", f"in_{estimator.bids_id}.anat_mask"),
+                    ("outputnode.sd_prior", f"in_{estimator.bids_id}.sd_prior"),
+                ]),
+            ])
+            # fmt:on
 
     return workflow

docs/requirements.txt

@@ -6,11 +6,11 @@ niworkflows >=1.4.0rc5,<1.5
 packaging
 pydot>=1.2.3
 pydotplus
-sdcflows ~= 2.0.0
+sdcflows ~= 2.0.1
 smriprep >= 0.8.0rc2
 sphinx >=2.1.2,<3.0
 sphinx-argparse
 sphinx_rtd_theme
 sphinxcontrib-apidoc ~= 0.3.0
 templateflow
-toml
+toml

setup.cfg

@@ -29,7 +29,7 @@ install_requires =
     numpy
     pybids >= 0.11.1
     pyyaml
-    sdcflows ~= 2.0.0
+    sdcflows ~= 2.0.1
     smriprep >= 0.8.0rc2
     svgutils != 0.3.2
     templateflow ~= 0.6