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.

.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/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,52 @@ 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 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,
+                iterfield=["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
+
+        raise RuntimeError(
+            f"Unknown fieldmap estimation strategy <{estimator}>."
+        )
 
     return workflow