Add SynB0 DISCO

qsiprep/cli/parser.py

@@ -576,6 +576,14 @@
         help='EXPERIMENTAL/TEMPORARY: Use SyN correction in addition to '
         'fieldmap correction, if available',
     )
+    g_syn.add_argument(
+        "--syn-method",
+        "--syn_method",
+        choices=["synb0", "legacy"],
+        action="store",
+        default="synb0",
+        help="Mehtod for SYN unwarping.  Options are synb0 (default) or legacy",
+    )
 
     g_other = parser.add_argument_group('Other options')
     g_other.add_argument('--version', action='version', version=verstr)

qsiprep/config.py

@@ -616,6 +616,9 @@ class workflow(_Config):
     use_syn_sdc = None
     """Run *fieldmap-less* susceptibility-derived distortions estimation
     in the absence of any alternatives."""
+    syn_method = None
+    """Which method to use for synthetic distortion correction. Must be
+    synb0 or legacy."""
 
 
 class loggers:

qsiprep/interfaces/synb0.py

@@ -0,0 +1,24 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+Interfaces for using SynB0-DISCO
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+"""
+import os
+import os.path as op
+
+
+def get_synb0_atlas(masked=True, res="low"):
+
+    atlas_dir = os.getenv("SYNB0_ATLASES")
+    if not atlas_dir:
+        raise Exception("Unable to locate SynB0 atlases. Define a SYNB0_ATLASES variable.")
+
+    res_str = "" if res == "high" else "_2_5"
+    mask_str = "_mask" if masked else ""
+
+    atlas_file = f"mni_icbm152_t1_tal_nlin_asym_09c{mask_str}{res_str}.nii.gz"
+    return op.join(atlas_dir, atlas_file)

qsiprep/tests/test_workflows.py

@@ -0,0 +1,11 @@
+from pathlib import Path
+from qsiprep.workflows.anatomical.volume import init_synb0_anat_wf
+
+input_nii = "/data/qsiprep_testing_data/anat_preprocessing/" \
+    "rigid_acpc_resample_brain_t1w/masked_brain_trans.nii"
+
+wf = init_synb0_anat_wf()
+wf.inputs.inputnode.t1w_brain_acpc = input_nii
+wf.base_dir = "/data/qsiprep_testing/synb0"
+wf.run()
+

qsiprep/workflows/anatomical/volume.py

@@ -32,7 +32,7 @@
 
 """
 
-from nipype.interfaces import afni, ants, mrtrix3
+from nipype.interfaces import afni, ants, freesurfer, mrtrix3
 from nipype.interfaces import utility as niu
 from nipype.interfaces.ants import BrainExtraction, N4BiasFieldCorrection
 from nipype.interfaces.base import traits
@@ -54,6 +54,7 @@
 )
 from ...interfaces.itk import AffineToRigid, DisassembleTransform
 from ...interfaces.niworkflows import RobustMNINormalizationRPT
+from ...interfaces.synb0 import get_synb0_atlas
 from ...utils.misc import fix_multi_source_name
 
 ANTS_VERSION = BrainExtraction().version or '<ver>'
@@ -137,6 +138,8 @@ def init_anat_preproc_wf(
         ANTs-compatible affine-and-warp transform file (inverse)
     t1_resampling_grid
         Image of the preprocessed t1 to be used as the reference output for dwis
+    t1_for_synb0
+        T1w image to be used for synb0-DISCO (if selected)
 
     """
 
@@ -881,7 +884,81 @@ def init_anat_normalization_wf(anatomical_template, has_rois=False) -> Workflow:
     return workflow
 
 
-def init_dl_prep_wf(name='dl_prep_wf') -> Workflow:
+def init_synb0_anat_wf(name="synb0_anat_wf") -> Workflow:
+    """Creates a T1w image that has been prepared/scaled for synb0 disco.
+
+    The steps come from the normalize_T1.sh script in v3.1 of Synb0-DISCO
+
+    Inputs
+    ------
+    t1w_brain_acpc
+        T1-weighted structural image to skull-strip
+
+    Outputs
+    -------
+    t1w_brain_acpc_nu
+        Full resolution, scaled brain image
+    t1w_brain_acpc_nu_2_5
+        The ``t1w_brain_acpc_nu`` image resampled into SynB0-DISCO's
+        input 2.5mm grid.
+
+    """
+    workflow = Workflow(name=name)
+    inputnode = pe.Node(niu.IdentityInterface(fields=["t1w_brain_acpc"]), name="inputnode")
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=["t1w_brain_acpc_nu", "t1w_brain_acpc_nu_2_5"]),
+        name="outputnode",
+    )
+
+    # mri_convert $JOB_PATH/T1.nii.gz $JOB_PATH/T1.mgz
+    convert_to_mgz = pe.Node(
+        freesurfer.MRIConvert(out_type="mgz"),
+        name="convert_to_mgz",
+    )
+
+    # mri_nu_correct.mni --i $JOB_PATH/T1.mgz --o $JOB_PATH/T1_N3.mgz --n 2
+    nu_mni = pe.Node(
+        freesurfer.MNIBiasCorrection(iterations=2),
+        name="nu_mni",
+    )
+
+    # mri_normalize -g 1 -mprage $JOB_PATH/T1_N3.mgz $JOB_PATH/T1_norm.mgz
+    mri_normalize = pe.Node(
+        freesurfer.Normalize(gradient=1, args="-mprage", out_file="T1_norm.mgz"),
+        name="mri_normalize",
+    )
+
+    # mri_convert $JOB_PATH/T1_norm.mgz $T1_NORM_PATH
+    convert_to_nii = pe.Node(
+        freesurfer.MRIConvert(out_type="niigz"),
+        name="convert_to_nii",
+    )
+
+    # From prepare_input.sh:
+    # "antsApplyTransforms -d 3 -i $T1_NORM_PATH -r $T1_ATLAS_2_5_PATH \
+    #   -n BSpline -t "$ANTS_OUT"0GenericAffine.mat -o $T1_NORM_LIN_ATLAS_2_5_PATH"
+    resample_to_model_grid = pe.Node(
+        ants.ApplyTransforms(
+            reference_image=get_synb0_atlas(masked=True, res="low"),
+            dimension=3,
+            interpolation="BSpline",
+            transforms="identity",
+        ),
+        name="resample_to_model_grid",
+    )
+    workflow.connect([
+        (inputnode, convert_to_mgz, [('t1w_brain_acpc', 'in_file')]),
+        (convert_to_mgz, nu_mni, [('out_file', 'in_file')]),
+        (nu_mni, mri_normalize, [('out_file', 'in_file')]),
+        (mri_normalize, convert_to_nii, [('out_file', 'in_file')]),
+        (convert_to_nii, outputnode, [('out_file', 't1w_brain_acpc_nu')]),
+        (convert_to_nii, resample_to_model_grid, [('out_file', 'input_image')]),
+        (resample_to_model_grid, outputnode, [('output_image', 't1w_brain_acpc_nu_2_5')])
+    ])  # fmt:skip
+    return workflow
+
+
+def init_dl_prep_wf(name="dl_prep_wf") -> Workflow:
     """Prepare images for use in the FreeSurfer deep learning functions"""
     workflow = Workflow(name=name)
     inputnode = pe.Node(niu.IdentityInterface(fields=['image']), name='inputnode')