diff --git a/.circleci/config.yml b/.circleci/config.yml index 714f60774..6e28eded5 100755 --- a/.circleci/config.yml +++ b/.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 diff --git a/dmriprep/cli/parser.py b/dmriprep/cli/parser.py index 65fb03230..a1e7e4542 100644 --- a/dmriprep/cli/parser.py +++ b/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", diff --git a/dmriprep/interfaces/vectors.py b/dmriprep/interfaces/vectors.py index 1329bd629..6e7010d10 100644 --- a/dmriprep/interfaces/vectors.py +++ b/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() diff --git a/dmriprep/workflows/base.py b/dmriprep/workflows/base.py index 9378186b3..80eb7cb00 100755 --- a/dmriprep/workflows/base.py +++ b/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 diff --git a/setup.cfg b/setup.cfg index 526caeba0..d29df45b1 100644 --- a/setup.cfg +++ b/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