diff --git a/.circleci/HBCD_preproc.sh b/.circleci/HBCD_preproc.sh new file mode 100644 index 00000000..a6728fee --- /dev/null +++ b/.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 + + diff --git a/.circleci/get_data.sh b/.circleci/get_data.sh index 1395fa8b..eca70f4d 100644 --- a/.circleci/get_data.sh +++ b/.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} \ diff --git a/qsiprep/cli/run.py b/qsiprep/cli/run.py index 8b10e6b0..d8137a0a 100644 --- a/qsiprep/cli/run.py +++ b/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', diff --git a/qsiprep/data/mni_1mm_t2w_lps.nii.gz b/qsiprep/data/mni_1mm_t2w_lps.nii.gz new file mode 100644 index 00000000..e3e64191 Binary files /dev/null and b/qsiprep/data/mni_1mm_t2w_lps.nii.gz differ diff --git a/qsiprep/interfaces/eddy.py b/qsiprep/interfaces/eddy.py index 00a56700..fbe23a1f 100644 --- a/qsiprep/interfaces/eddy.py +++ b/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) \ No newline at end of file diff --git a/qsiprep/interfaces/tortoise.py b/qsiprep/interfaces/tortoise.py index 00fef2f0..29075fdd 100644 --- a/qsiprep/interfaces/tortoise.py +++ b/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) \ No newline at end of file diff --git a/qsiprep/viz/config.json b/qsiprep/viz/config.json index 56c0abd5..1664c47c 100644 --- a/qsiprep/viz/config.json +++ b/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": "antsRegistration was used to generate transformations from the b=0 reference image to the T1w-image.", "imgtype": "svg+xml" }, diff --git a/qsiprep/workflows/base.py b/qsiprep/workflows/base.py index 2a8d7cfe..1ed87feb 100644 --- a/qsiprep/workflows/base.py +++ b/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 diff --git a/qsiprep/workflows/dwi/base.py b/qsiprep/workflows/dwi/base.py index 1a05f6b5..5ffe07ad 100644 --- a/qsiprep/workflows/dwi/base.py +++ b/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', diff --git a/qsiprep/workflows/dwi/finalize.py b/qsiprep/workflows/dwi/finalize.py index f7a43b52..8268d5cf 100644 --- a/qsiprep/workflows/dwi/finalize.py +++ b/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 diff --git a/qsiprep/workflows/dwi/fsl.py b/qsiprep/workflows/dwi/fsl.py index c520c9ae..c223d433 100644 --- a/qsiprep/workflows/dwi/fsl.py +++ b/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, [ diff --git a/qsiprep/workflows/dwi/hmc_sdc.py b/qsiprep/workflows/dwi/hmc_sdc.py index 3dff035e..5124e57a 100644 --- a/qsiprep/workflows/dwi/hmc_sdc.py +++ b/qsiprep/workflows/dwi/hmc_sdc.py @@ -163,11 +163,12 @@ def init_qsiprep_hmcsdc_wf(scan_groups, fieldmap_type = fieldmap_info['suffix'] if fieldmap_type in ('epi', 'rpe_series'): - if pepolar_method.lower() == "topup": + if 'topup' in pepolar_method.lower(): raise Exception("TOPUP is not supported with SHORELine ") drbuddi_wf = init_drbuddi_wf(scan_groups=scan_groups, omp_nthreads=omp_nthreads, b0_threshold=b0_threshold, raw_image_sdc=raw_image_sdc, + pepolar_method=pepolar_method, sloppy=sloppy, t2w_sdc=t2w_sdc) # apply the head motion correction transforms diff --git a/qsiprep/workflows/fieldmap/drbuddi.py b/qsiprep/workflows/fieldmap/drbuddi.py index 3f50f8aa..fc44e02b 100644 --- a/qsiprep/workflows/fieldmap/drbuddi.py +++ b/qsiprep/workflows/fieldmap/drbuddi.py @@ -28,14 +28,15 @@ from ...engine import Workflow from ...interfaces.fmap import get_distortion_grouping from ...interfaces.tortoise import ( - GatherDRBUDDIInputs, DRBUDDI, DRBUDDIAggregateOutputs) + GatherDRBUDDIInputs, DRBUDDI, DRBUDDIAggregateOutputs, + generate_drbuddi_boilerplate) from ...interfaces.reports import TopupSummary LOGGER = logging.getLogger('nipype.workflow') DEFAULT_MEMORY_MIN_GB = 0.01 -def init_drbuddi_wf(scan_groups, b0_threshold, raw_image_sdc, t2w_sdc, omp_nthreads=1, +def init_drbuddi_wf(scan_groups, b0_threshold, pepolar_method, raw_image_sdc, t2w_sdc, omp_nthreads=1, name="drbuddi_sdc_wf", sloppy=False): """ This workflow implements the heuristics to choose a @@ -114,29 +115,15 @@ def init_drbuddi_wf(scan_groups, b0_threshold, raw_image_sdc, t2w_sdc, omp_nthre "down_fa_corrected_image", "t2w_image" ]), name='outputnode') - desc = ["Data was collected with reversed phase-encode blips, resulting " - "in pairs of images with distortions going in opposite directions."] fieldmap_info = scan_groups['fieldmap_info'] if fieldmap_info['suffix'] not in ('epi', 'rpe_series', 'dwi'): raise Exception("DRBUDDI workflow requires epi, rpe_series or dwi fieldmaps") - # Describe what's going on - if fieldmap_info["suffix"] == "epi": - desc.append("Here, b=0 reference images with reversed " - "phase encoding directions were used to estimate ") - else: - desc.append("Here, multiple DWI series were 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("An updated version of DRBUDDI [@drbuddi], part of the TORTOISE [@tortoisev3] " - "software package was used to estimate distortion. Signal intensity was adjusted" - "in the final interpolated images using a method similar to LSR.\n" ) - workflow.__desc__ = " ".join(desc) + workflow.__desc__ = generate_drbuddi_boilerplate( + fieldmap_type=fieldmap_info['suffix'], + t2w_sdc=t2w_sdc, + with_topup="topup" in pepolar_method.lower()) outputnode.inputs.method = \ 'PEB/PEPOLAR (phase-encoding based / PE-POLARity): %s' % fieldmap_info['suffix'] @@ -211,4 +198,4 @@ def init_drbuddi_wf(scan_groups, b0_threshold, raw_image_sdc, t2w_sdc, omp_nthre ("b0_ref", "b0_ref")]) ]) - return workflow + return workflow \ No newline at end of file