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

Currently, this is blocked by the absence of a general solution to the
reference generation issue (see nipreps/niworkflows#601).

Once adequate DWI references are connected to the estimation workflow,
this feature should work out.
We can test the new feature with ds000206.

Resolves: #111.

Author: oesteban

.circleci/config.yml

@@ -314,7 +314,7 @@ jobs:
                 -v /tmp/ds000206/work:/work \
                 --user $(id -u):$(id -g) \
                 nipreps/dmriprep:latest /data /out participant -vv $eddy \
-                --fs-subjects-dir /data/derivatives/freesurfer-6.0.1 --sloppy \
+                --fs-subjects-dir /data/derivatives/freesurfer-6.0.1 --sloppy --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/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,6 +360,10 @@ 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)
@@ -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,75 @@ 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.interfaces.brainmask import BrainExtraction
+            from sdcflows.workflows.fit.syn import init_syn_preprocessing_wf
+            from ...interfaces.vectors import CheckGradientTable
+
+            sources = [str(s.path) for s in estimator.sources]
+            layout = config.execution.layout
+            syn_preprocessing_wf = init_syn_preprocessing_wf(
+                omp_nthreads=config.nipype.omp_nthreads,
+                auto_bold_nss=False,
+            )
+            syn_preprocessing_wf.inputs.in_epis = sources
+            syn_preprocessing_wf.inputs.in_meta = [
+                layout.get_metadata(s) for s in sources
+            ]
+
+            b0_masks = pe.MapNode(CheckGradientTable(), name="b0_masks",
+                                  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]
+
+            epi_brain = pe.Node(BrainExtraction(), name="epi_brain")
+
+            # fmt:off
+            workflow.connect([
+                (anat_preproc_wf, syn_preprocessing_wf, [
+                    ("outputnode.t1w_preproc", "inputnode.in_anat"),
+                    ("outputnode.t1w_mask", "inputnode.mask_anat"),
+                ]),
+                (b0_masks, syn_preprocessing_wf, [("b0_mask", "inputnode.t_masks")]),
+                (syn_preprocessing_wf, fmap_wf, [
+                    ("outputnode.epi_ref", "inputnode.epi_ref"),
+                    ("outputnode.anat_ref", "inputnode.anat_ref"),
+                    ("outputnode.anat_mask", "inputnode.anat_mask"),
+                    ("outputnode.anat2epi_xfm", "inputnode.anat2epi_xfm"),
+                ]),
+                (syn_preprocessing_wf, epi_brain, [("outputnode.epi_ref", "in_file")]),
+                (anat_preproc_wf, fmap_wf, [
+                    ("outputnode.std2anat_xfm", "inputnode.std2anat_xfm"),
+                ]),
+                (epi_brain, fmap_wf, [("out_mask", "inputnode.epi_mask")]),
+            ])
+            # fmt:on
+
+        raise RuntimeError(
+            f"Unknown fieldmap estimation strategy <{estimator}>."
+        )
 
     return workflow

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