[ENH] allow topup+drbuddi for hbcd (#667)

.circleci/HBCD_preproc.sh

@@ -0,0 +1,44 @@
+#!/bin/bash
+
+cat << DOC
+
+Test paired DWI series with DRBUDDI
+===================================
+
+This tests the following features:
+ - Eddy is run on a CPU
+ - DRBUDDI is run with two DWI series
+
+DOC
+
+set +e
+source ./get_data.sh
+TESTDIR=${PWD}
+TESTNAME=HBCD
+get_config_data ${TESTDIR}
+get_bids_data ${TESTDIR} HBCD
+CFG=${TESTDIR}/data/nipype.cfg
+EDDY_CFG=${TESTDIR}/data/eddy_config.json
+
+# For the run
+setup_dir ${TESTDIR}/${TESTNAME}
+TEMPDIR=${TESTDIR}/${TESTNAME}/work
+OUTPUT_DIR=${TESTDIR}/${TESTNAME}/derivatives
+BIDS_INPUT_DIR=${TESTDIR}/data/hbcd_sim
+export FS_LICENSE=${TESTDIR}/data/license.txt
+QSIPREP_CMD=$(run_qsiprep_cmd ${BIDS_INPUT_DIR} ${OUTPUT_DIR})
+
+# Do the HBCD-style run
+${QSIPREP_CMD} \
+	-w ${TEMPDIR} \
+	--sloppy \
+	--anat-modality T2w \
+	--infant \
+	--denoise-method dwidenoise \
+	--b1_biascorrect_stage none \
+	--pepolar-method DRBUDDI \
+	--eddy_config ${EDDY_CFG} \
+	--output-resolution 5 \
+    -vv --stop-on-first-crash
+
+

.circleci/get_data.sh

@@ -306,6 +306,15 @@ get_bids_data() {
     mkdir -p ${WORKDIR}/data
     cd ${WORKDIR}/data
 
+    # Phantom hbcd-like data
+    if [[ ${DS} = HBCD ]]; then
+      ${WGET} \
+        -O HBCD.tar.xz \
+	    "https://upenn.box.com/shared/static/gn1ec8x7mtk1f07l97d0th9idn4qv3yx.xz"
+      tar xvfJ HBCD.tar.xz -C ${WORKDIR}/data/
+      rm HBCD.tar.xz
+    fi
+
     # Down-sampled compressed sensing DSI
     if [[ ${DS} = DSCSDSI ]]; then
       ${WGET} \

qsiprep/cli/run.py

@@ -459,8 +459,8 @@ def get_parser():
         '--pepolar-method', '--pepolar_method',
         action='store',
         default='TOPUP',
-        choices=['TOPUP', 'DRBUDDI'],
-        help='select which SDC method to use for PEPOLAR fieldmaps (default: OASIS)')
+        choices=['TOPUP', 'DRBUDDI', 'TOPUP+DRBUDDI'],
+        help='select which SDC method to use for PEPOLAR fieldmaps (default: TOPUP)')
 
     g_fmap.add_argument(
         '--denoised_image_sdc', '--denoised_image_sdc',

qsiprep/interfaces/eddy.py

@@ -260,6 +260,7 @@ def _run_interface(self, runtime):
 def boilerplate_from_eddy_config(eddy_config, fieldmap_type, pepolar_method):
     """Write boilerplate text based on an eddy config dict.
     """
+    doing_2stage = "drbuddi" in pepolar_method.lower()
     ext_eddy = ExtendedEddy(**eddy_config)
     desc = [
         "FSL (version %s)'s eddy was used for head motion correction and "
@@ -359,8 +360,8 @@ def boilerplate_from_eddy_config(eddy_config, fieldmap_type, pepolar_method):
                     % (ext_eddy.inputs.mporder, niter, s2v_interp, lam))
 
     # distortion correction
-    if pepolar_method.lower() == 'topup':
-        desc.append(topup_boilerplate(fieldmap_type))
+    if "topup" in pepolar_method.lower():
+        desc.append(topup_boilerplate(fieldmap_type, pepolar_method))
     # DRBUDDI is described in its own workflow
 
     # move by susceptibility
@@ -379,29 +380,42 @@ def boilerplate_from_eddy_config(eddy_config, fieldmap_type, pepolar_method):
 
     # Format the interpolation
     lsr_ref = ' [@fsllsr]' if ext_eddy.inputs.method == 'lsr' else ''
-    desc.append("Final interpolation was performed using the `%s` method%s.\n\n" % (
-        ext_eddy.inputs.method, lsr_ref))
-
+    if doing_2stage:
+        desc.append("Interpolation after head motion and initial susceptibility "
+                    "distortion correction")
+    else:
+        desc.append("Final interpolation")
+    desc.append("was performed using the `%s` method%s." % (
+                ext_eddy.inputs.method, lsr_ref))
+    if not doing_2stage:
+        desc.append("\n\n")
     return " ".join(desc)
 
 
-def topup_boilerplate(fieldmap_type):
+def topup_boilerplate(fieldmap_type, pepolar_method):
     """Write boilerplate text based on fieldmaps
     """
+    if fieldmap_type not in ("rpe_series", "epi"):
+        return ""
     desc = []
-    if fieldmap_type in ("rpe_series", "epi"):
-        desc.append("Data was collected with reversed phase-encode blips, resulting "
-                    "in pairs of images with distortions going in opposite directions.")
-        if fieldmap_type == "epi":
-            desc.append("Here, b=0 reference images with reversed "
-                        "phase encoding directions were used "
-                        "along with an equal number of b=0 images extracted "
-                        "from the DWI scans.")
-        else:
-            desc.append("Here, multiple DWI series were acquired with opposite phase encoding "
-                        "directions, so b=0 images were extracted from each.")
-        desc.append("From these pairs the susceptibility-induced off-resonance field was "
-                    "estimated using a method similar to that described in [@topup]. "
-                    "The fieldmaps were ultimately incorporated into the "
-                    "Eddy current and head motion correction interpolation.")
+    desc.append("\n\nData was collected with reversed phase-encode blips, resulting "
+                "in pairs of images with distortions going in opposite directions.")
+
+    if 'drbuddi' in pepolar_method.lower():
+        desc.append(
+            "Distortion correction was performed in two stages. In the first stage, "
+            "FSL's TOPUP [@topup]")
+    else:
+        desc.append("FSL's TOPUP [@topup]")
+
+    desc.append("was used to estimate a susceptibility-induced off-resonance field based on")
+    if fieldmap_type == "epi":
+        desc.append("b=0 reference images with reversed "
+                    "phase encoding directions.")
+    else:
+        desc.append("b=0 images extracted from multiple DWI series  "
+                    "with reversed phase encoding directions.")
+    desc.append("The TOPUP-estimated fieldmap was incorporated into the "
+                "Eddy current and head motion correction interpolation.")
+
     return " ".join(desc)

qsiprep/interfaces/tortoise.py

@@ -635,3 +635,36 @@ def make_bmat_file(bvals, bvecs):
         ["FSLBVecsToTORTOISEBmatrix", op.abspath(bvals), op.abspath(bvecs)])
     print(pout)
     return bvals.replace("bval", "bmtxt")
+
+
+def generate_drbuddi_boilerplate(fieldmap_type, t2w_sdc, with_topup=False):
+    """Generate boilerplate that describes how DRBUDDI is being used."""
+
+    desc = ["\n\nDRBUDDI [@drbuddi], part of the TORTOISE [@tortoisev3] software package,"]
+    if not with_topup:
+        # Until now there will have been no description of the SDC procedure.
+        # Add extra details about the input data.
+        desc.append(
+            "was used to perform susceptibility distortion correction. "
+            "Data was collected with reversed phase-encode blips, resulting "
+            "in pairs of images with distortions going in opposite directions.")
+    else:
+        desc += ["was used to perform a second stage of distortion correction."]
+
+    # Describe what's going on
+    if fieldmap_type == "epi":
+        desc.append("DRBUDDI used b=0 reference images with reversed "
+                    "phase encoding directions to estimate")
+    else:
+        desc.append("DRBUDDI used multiple motion-corrected DWI series acquired "
+                    "with opposite phase encoding "
+                    "directions. A b=0 image **and** the Fractional Anisotropy "
+                    "images from both phase encoding diesctions were used together in "
+                    "a multi-modal registration to estimate")
+    desc.append("the susceptibility-induced off-resonance field.")
+
+    if t2w_sdc:
+        desc.append("A T2-weighted image was included in the multimodal registration.")
+    desc.append("Signal intensity was adjusted "
+                "in the final interpolated images using a method similar to LSR.\n\n" )
+    return " ".join(desc)

qsiprep/viz/config.json

@@ -163,15 +163,22 @@
             {
                 "name": "epi/unwarp",
                 "file_pattern": "dwi/.*sdc_b0.*\\.",
-                "title": "Susceptibility distortion correction",
-                "description": "Results of performing susceptibility distortion correction (SDC) using b=0 images",
+                "title": "Susceptibility distortion correction (TOPUP)",
+                "description": "Results of performing susceptibility distortion correction (SDC) using b=0 images in TOPUP",
+                "imgtype": "svg+xml"
+            },
+            {
+                "name": "epi/unwarp",
+                "file_pattern": "dwi/.*sdcdrbuddi_b0.*\\.",
+                "title": "Susceptibility distortion correction (DRBUDDI)",
+                "description": "Results of performing susceptibility distortion correction (SDC) using DRBUDDI. b=0 images are shown",
                 "imgtype": "svg+xml"
             },
             {
                 "name": "epi/unwarp",
                 "file_pattern": "dwi/.*sdc_b0t2w.*\\.",
-                "title": "Susceptibility distortion correction",
-                "description": "Results of performing susceptibility distortion correction (SDC) using b=0 images. The overlay shown is an ad-hoc segmentation of a t2w image and is only for display purposes.",
+                "title": "Susceptibility distortion correction (DRBUDDI+T2w)",
+                "description": "Results of performing susceptibility distortion correction (SDC) using DRBUDDI. The overlay shown is an ad-hoc segmentation of a t2w image and is only for display purposes.",
                 "imgtype": "svg+xml"
             },
             {
@@ -191,7 +198,7 @@
             {
                 "name": "epi_mean_t1_registration/b0_coreg",
                 "file_pattern": "dwi/.*_coreg\\.",
-                "title": "b=0 to T1 registration",
+                "title": "b=0 to anatomical reference registration",
                 "description": "<code>antsRegistration</code> was used to generate transformations from the b=0 reference image to the T1w-image.",
                 "imgtype": "svg+xml"
             },

qsiprep/workflows/base.py

@@ -166,7 +166,7 @@ def init_qsiprep_wf(
             How to combine multiple distortion groups after correction. 'concat', 'average' or
             'none'
         pepolar_method : str
-            Either 'DRBUDDI' or 'TOPUP'. The method for SDC when EPI fieldmaps are used.
+            Either 'DRBUDDI', 'TOPUP' or 'TOPUP+DRBUDDI'. The method for SDC when EPI fieldmaps are used.
         omp_nthreads : int
             Maximum number of threads an individual process may use
         skull_strip_template : str

qsiprep/workflows/dwi/base.py

@@ -142,7 +142,7 @@ def init_dwi_preproc_wf(dwi_only,
         unringing_method : str
             algorithm to use for removing Gibbs ringing. Options: none, mrdegibbs
         pepolar_method : str
-            Either 'DRBUDDI' or 'TOPUP'. The method for SDC when EPI fieldmaps are used.
+            Either 'DRBUDDI', 'TOPUP' or 'TOPUP+DRBUDDI'. The method for SDC when EPI fieldmaps are used.
         b1_biascorrect_stage : str
             'final', 'none' or 'legacy'
         no_b0_harmonization : bool
@@ -405,7 +405,7 @@ def init_dwi_preproc_wf(dwi_only,
 
     # Fieldmap reports should vary depending on which type of correction is performed
     # PEPOLAR (epi, rpe series) will produce potentially much more detailed reports
-    doing_topup = fieldmap_type in ("epi", "rpe_series") and pepolar_method.lower() == "topup"
+    doing_topup = fieldmap_type in ("epi", "rpe_series") and "topup" in pepolar_method.lower()
     if fieldmap_type not in ("epi", "rpe_series", None) or doing_topup:
         fmap_unwarp_report_wf = init_fmap_unwarp_report_wf()
         ds_report_sdc = pe.Node(
@@ -427,7 +427,8 @@ def init_dwi_preproc_wf(dwi_only,
             (fmap_unwarp_report_wf, ds_report_sdc, [('outputnode.report', 'in_file')])
         ])
 
-    elif fieldmap_type in ("epi", "rpe_series"):
+    # DRBUDDI has some extra reports that we want to save. Make sure we get them!
+    if fieldmap_type in ("epi", "rpe_series") and "drbuddi" in pepolar_method.lower():
 
         if os.path.exists(t2w_sdc):
             extended_pepolar_report_wf = init_extended_pepolar_report_wf(
@@ -446,7 +447,7 @@ def init_dwi_preproc_wf(dwi_only,
 
         ds_report_b0_sdc = pe.Node(
             DerivativesMaybeDataSink(
-                desc="sdc",
+                desc="sdcdrbuddi",
                 suffix='b0' if not t2w_sdc else 'b0t2w',
                 source_file=source_file),
             name='ds_report_b0_sdc',

qsiprep/workflows/dwi/finalize.py

@@ -90,7 +90,7 @@ def init_dwi_finalize_wf(scan_groups,
         output_resolution : float
             Output voxel resolution in mm
         pepolar_method : str
-            Either 'DRBUDDI' or 'TOPUP'. The method for SDC when EPI fieldmaps are used.
+            Either 'DRBUDDI', 'TOPUP' or 'TOPUP+DRBUDDI'. The method for SDC when EPI fieldmaps are used.
         omp_nthreads : int
             Maximum number of threads an individual process may use
         low_mem : bool

qsiprep/workflows/dwi/fsl.py

@@ -71,7 +71,7 @@ def init_fsl_hmc_wf(scan_groups,
             threshold for a slice to be replaced with imputed values. Overrides the
             parameter in ``eddy_config`` if set to a number > 0.
         pepolar_method : str
-            Either 'DRBUDDI' or 'TOPUP'. The method for SDC when EPI fieldmaps are used.
+            Either 'DRBUDDI', 'TOPUP' or 'DRBUDDI+TOPUP'. The method for SDC when EPI fieldmaps are used.
         eddy_config: str
             Path to a JSON file containing settings for the call to ``eddy``.
 
@@ -100,7 +100,7 @@ def init_fsl_hmc_wf(scan_groups,
     fsl_check = os.environ.get('FSL_BUILD')
     if fsl_check=="no_fsl":
         raise Exception(
-            """Container in use does not have FSL. To use this workflow, 
+            """Container in use does not have FSL. To use this workflow,
             please download the qsiprep container with FSL installed.""")
     inputnode = pe.Node(
         niu.IdentityInterface(
@@ -220,7 +220,9 @@ def init_fsl_hmc_wf(scan_groups,
     fieldmap_type = scan_groups['fieldmap_info']['suffix'] or ''
     workflow.__desc__ = boilerplate_from_eddy_config(eddy_args, fieldmap_type,
                                                      pepolar_method)
-    if fieldmap_type in ('epi', 'rpe_series') and pepolar_method.lower() == "topup":
+
+    # Are we running TOPUP?
+    if fieldmap_type in ('epi', 'rpe_series') and "topup" in pepolar_method.lower():
         # If there are EPI fieldmaps in fmaps/, make sure they get to TOPUP. It will always use
         # b=0 images from the DWI series regardless
         gather_inputs.inputs.topup_requested = True
@@ -247,8 +249,6 @@ def init_fsl_hmc_wf(scan_groups,
         topup_to_eddy_reg = pe.Node(fsl.FLIRT(dof=6, output_type="NIFTI_GZ"),
                                     name="topup_to_eddy_reg")
         workflow.connect([
-            # There will be no SDC warps, they are applied by eddy
-            (gather_inputs, outputnode, [('forward_warps', 'to_dwi_ref_warps')]),
             (gather_inputs, topup, [
                 ('topup_datain', 'encoding_file'),
                 ('topup_imain', 'in_file'),
@@ -263,27 +263,31 @@ def init_fsl_hmc_wf(scan_groups,
             # Use corrected images from TOPUP to make a mask for eddy
             (topup, unwarped_mean, [('out_corrected', 'in_files')]),
             (unwarped_mean, pre_eddy_b0_ref_wf, [('out_avg', 'inputnode.b0_template')]),
-            (b0_ref_to_lps, outputnode, [('dwi_file', 'b0_template')]),
             # Save reports
             (gather_inputs, topup_summary, [('topup_report', 'summary')]),
             (topup_summary, ds_report_topupsummary, [('out_report', 'in_file')])
-
         ])
 
-        return workflow
-
-    # All the following distortion correction methods operate on images
-    # *already processed by eddy*. Therefore we need to make a mask for
-    # eddy that contains both distorted brain shapes
-    distorted_merge = pe.Node(
-        IntraModalMerge(hmc=True, to_lps=False), name='distorted_merge')
-    workflow.connect([
+        if "drbuddi" not in pepolar_method.lower():
+            LOGGER.info("Using single-stage SDC, TOPUP-only")
+            workflow.connect([
+                # There will be no SDC warps, they are applied by eddy
+                (gather_inputs, outputnode, [('forward_warps', 'to_dwi_ref_warps')]),
+                (b0_ref_to_lps, outputnode, [('dwi_file', 'b0_template')]),
+            ])
+    else:
+        # If we're not using TOPUP we need to make a mask for eddy based on the distorted brain shapes
+        distorted_merge = pe.Node(
+            IntraModalMerge(hmc=True, to_lps=False), name='distorted_merge')
         # Use the distorted mask for eddy
-        (gather_inputs, distorted_merge, [('topup_imain', 'in_files')]),
-        (distorted_merge, pre_eddy_b0_ref_wf, [('out_avg', 'inputnode.b0_template')])])
+        workflow.connect([
+            (gather_inputs, distorted_merge, [('topup_imain', 'in_files')]),
+            (distorted_merge, pre_eddy_b0_ref_wf, [('out_avg', 'inputnode.b0_template')])])
+
 
-    if fieldmap_type in ('epi', 'rpe_series') and pepolar_method.lower() == "drbuddi":
+    if fieldmap_type in ('epi', 'rpe_series') and 'drbuddi' in pepolar_method.lower():
         outputnode.inputs.sdc_method = "DRBUDDI"
+        LOGGER.info("Running DRBUDDI for SDC")
 
         # Let gather_inputs know we're doing pepolar, even though it's not topup
         gather_inputs.inputs.topup_requested = True
@@ -292,6 +296,7 @@ def init_fsl_hmc_wf(scan_groups,
 
         drbuddi_wf = init_drbuddi_wf(scan_groups=scan_groups, omp_nthreads=omp_nthreads,
                                      t2w_sdc=t2w_sdc, b0_threshold=b0_threshold,
+                                     pepolar_method=pepolar_method,
                                      raw_image_sdc=raw_image_sdc,
                                      sloppy=sloppy)
 
@@ -323,10 +328,8 @@ def init_fsl_hmc_wf(scan_groups,
 
         return workflow
 
-
-
     if fieldmap_type in ('fieldmap', 'syn') or fieldmap_type.startswith("phase"):
-
+        LOGGER.info(f"Computing fieldmap directly from {fieldmap_type}")
         outputnode.inputs.sdc_method = fieldmap_type
         b0_sdc_wf = init_sdc_wf(
             scan_groups['fieldmap_info'], dwi_metadata, omp_nthreads=omp_nthreads,
@@ -348,7 +351,7 @@ def init_fsl_hmc_wf(scan_groups,
                 ('outputnode.method', 'sdc_method'),
                 ('outputnode.b0_ref', 'b0_template')])])
 
-    else:
+    if not fieldmap_type:
         outputnode.inputs.sdc_method = "None"
         workflow.connect([
             (b0_ref_to_lps, outputnode, [