Merge pull request #488 from nipreps/fix/correctly-handle-priors

FIX: Re-enable priors respecting ``sd_priors`` argument

Author: effigies

sdcflows/workflows/fit/syn.py

@@ -26,7 +26,6 @@
 
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu
-from nipype.interfaces.base import Undefined
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
 from ... import data
@@ -44,13 +43,11 @@
 
 def init_syn_sdc_wf(
     *,
-    atlas_threshold=3,
     sloppy=False,
     debug=False,
     name="syn_sdc_wf",
     omp_nthreads=1,
     laplacian_weight=None,
-    sd_prior=True,
     **kwargs,
 ):
     """
@@ -59,10 +56,6 @@ def init_syn_sdc_wf(
     SyN deformation is restricted to the phase-encoding (PE) direction.
     If no PE direction is specified, anterior-posterior PE is assumed.
 
-    SyN deformation is also restricted to regions that are expected to have a
-    >3mm (approximately 1 voxel) warp, based on the fieldmap atlas.
-
-
     Workflow Graph
         .. workflow ::
             :graph2use: orig
@@ -73,9 +66,6 @@ def init_syn_sdc_wf(
 
     Parameters
     ----------
-    atlas_threshold : :obj:`float`
-        Exclude from the registration metric computation areas with average distortions
-        below this threshold (in mm).
     sloppy : :obj:`bool`
         Whether a fast (less accurate) configuration of the workflow should be applied.
     debug : :obj:`bool`
@@ -102,9 +92,6 @@ def init_syn_sdc_wf(
         A preprocessed, skull-stripped anatomical (T1w or T2w) image resampled in EPI space.
     anat_mask : :obj:`str`
         Path to the brain mask corresponding to ``anat_ref`` in EPI space.
-    sd_prior : :obj:`str`
-        A template map of areas with strong susceptibility distortions (SD) to regularize
-        the cost function of SyN
 
     Outputs
     -------
@@ -150,16 +137,15 @@ def init_syn_sdc_wf(
     workflow = Workflow(name=name)
     workflow.__desc__ = f"""\
 A deformation field to correct for susceptibility distortions was estimated
-based on *fMRIPrep*'s *fieldmap-less* approach.
+based on *SDCFlows*' *fieldmap-less* approach.
 The deformation field is that resulting from co-registering the EPI reference
-to the same-subject T1w-reference with its intensity inverted [@fieldmapless1;
-@fieldmapless2].
+to the same-subject's T1w-reference [@fieldmapless1; @fieldmapless2].
 Registration is performed with `antsRegistration`
 (ANTs {ants_version or "-- version unknown"}), and
 the process regularized by constraining deformation to be nonzero only
-along the phase-encoding direction, and modulated with an average fieldmap
-template [@fieldmapless3].
+along the phase-encoding direction.
 """
+
     inputnode = pe.Node(niu.IdentityInterface(INPUT_FIELDS), name="inputnode")
     outputnode = pe.Node(
         niu.IdentityInterface(
@@ -211,10 +197,11 @@ def init_syn_sdc_wf(
 
     epi_umask = pe.Node(Union(), name="epi_umask")
     moving_masks = pe.Node(
-        niu.Merge(3),
+        niu.Merge(2),
         name="moving_masks",
         run_without_submitting=True,
     )
+    moving_masks.inputs.in1 = "NULL"
 
     fixed_masks = pe.Node(
         niu.Merge(2),
@@ -287,7 +274,7 @@ def init_syn_sdc_wf(
         (inputnode, amask2epi, [("epi_mask", "reference_image")]),
         (inputnode, zooms_bmask, [("anat_mask", "input_image")]),
         (inputnode, fixed_masks, [("anat_mask", "in1"),
-                                  ("anat_mask", "in2")]),
+                                  ("sd_prior", "in2")]),
         (inputnode, anat_dilmsk, [("anat_mask", "in_file")]),
         (inputnode, warp_dir, [("anat_ref", "fixed_image")]),
         (inputnode, anat_merge, [("anat_ref", "in1")]),
@@ -298,9 +285,7 @@ def init_syn_sdc_wf(
         (inputnode, epi_umask, [("epi_mask", "in1")]),
         (lap_anat, lap_anat_norm, [("output_image", "in_file")]),
         (lap_anat_norm, anat_merge, [("out", "in2")]),
-        (epi_umask, moving_masks, [("out_file", "in1"),
-                                   ("out_file", "in2"),
-                                   ("out_file", "in3")]),
+        (epi_umask, moving_masks, [("out_file", "in2")]),
         (clip_epi, epi_merge, [("out_file", "in1")]),
         (clip_epi, lap_epi, [("out_file", "op1")]),
         (clip_epi, zooms_epi, [("out_file", "in_file")]),
@@ -339,6 +324,7 @@ def init_syn_sdc_wf(
 
 def init_syn_preprocessing_wf(
     *,
+    atlas_threshold=3,
     debug=False,
     name="syn_preprocessing_wf",
     omp_nthreads=1,
@@ -359,6 +345,9 @@ def init_syn_preprocessing_wf(
 
     Parameters
     ----------
+    atlas_threshold : :obj:`float`
+        Mask excluding areas with average distortions below this threshold (in mm)
+        on the prior.
     debug : :obj:`bool`
         Whether a fast (less accurate) configuration of the workflow should be applied.
     name : :obj:`str`
@@ -499,6 +488,8 @@ def _remove_first_mask(in_file):
     sampling_ref = pe.Node(GenerateSamplingReference(), name="sampling_ref")
 
     if sd_prior:
+        from niworkflows.interfaces.nibabel import Binarize
+
         # Mapping & preparing prior knowledge
         # Concatenate transform files:
         # 1) MNI -> anat; 2) ATLAS -> MNI
@@ -521,19 +512,20 @@ def _remove_first_mask(in_file):
             mem_gb=0.3,
         )
 
+        prior_msk = pe.Node(Binarize(thresh_low=atlas_threshold), name="prior_msk")
+
         workflow.connect([
             (inputnode, transform_list, [("std2anat_xfm", "in2")]),
             (epi2anat, transform_list, [("forward_transforms", "in1")]),
             (transform_list, prior2epi, [("out", "transforms")]),
             (sampling_ref, prior2epi, [("out_file", "reference_image")]),
-            (prior2epi, outputnode, [("output_image", "sd_prior")]),
+            (prior2epi, prior_msk, [("output_image", "in_file")]),
+            (prior_msk, outputnode, [("out_mask", "sd_prior")]),
         ])  # fmt:skip
 
     else:
-        # no prior to be used
-        # MG: Future goal is to allow using alternative mappings
-        # i.e. in the case of infants, where priors change depending on development
-        outputnode.inputs.sd_prior = Undefined
+        # no prior to be used -> set anatomical mask as prior
+        workflow.connect(mask_dtype, "out", outputnode, "sd_prior")
 
     workflow.connect([
         (inputnode, epi_reference_wf, [("in_epis", "inputnode.in_files")]),

sdcflows/workflows/fit/tests/test_syn.py

@@ -31,8 +31,7 @@
 
 @pytest.mark.veryslow
 @pytest.mark.slow
-@pytest.mark.parametrize("sd_prior", [True, False])
-def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
+def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode):
     """Build and run an SDC-SyN workflow."""
     derivs_path = datadir / "ds000054" / "derivatives"
     smriprep = derivs_path / "smriprep-0.6" / "sub-100185" / "anat"
@@ -44,7 +43,6 @@ def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
         debug=sloppy_mode,
         auto_bold_nss=True,
         t1w_inversion=True,
-        sd_prior=sd_prior,
     )
     prep_wf.inputs.inputnode.in_epis = [
         str(
@@ -79,7 +77,6 @@ def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
         debug=sloppy_mode,
         sloppy=sloppy_mode,
         omp_nthreads=4,
-        sd_prior=sd_prior,
     )
 
     # fmt: off
@@ -170,6 +167,24 @@ def test_syn_wf_inputs(sloppy, laplacian_weight):
     assert wf.inputs.syn.metric_weight == metric_weight
 
 
+@pytest.mark.parametrize("sd_prior", [True, False])
+def test_syn_preprocessing_wf_inputs(sd_prior):
+    """Test appropriate instantiation of the SDC-SyN preprocessing workflow."""
+
+    prep_wf = init_syn_preprocessing_wf(
+        omp_nthreads=4,
+        sd_prior=sd_prior,
+        auto_bold_nss=True,
+        t1w_inversion=True,
+    )
+
+    if not sd_prior:
+        with pytest.raises(AttributeError):
+            prep_wf.inputs.prior_msk.in_file
+    else:
+        assert prep_wf.inputs.prior_msk.thresh_low
+
+
 @pytest.mark.parametrize("ants_version", ["2.2.0", "2.1.0", None])
 def test_syn_wf_version(monkeypatch, ants_version):
     """Ensure errors are triggered with ANTs < 2.2."""