Collect templates from TemplateFlow instead of package data (#49)

* Rename QSIRecon to QSIPrep as necessary.

* Remove unused template files.

* Load templates from TemplateFlow.

* Remove remaining unused templates.

* Whoops!

* Apply suggestions from code review

Co-authored-by: Matt Cieslak <[email protected]>

* Use MNIInfant cohort-2, not cohort-1.

* Mask the template file.

* Fix orientation.

* Update anatomical.py

---------

Co-authored-by: Matt Cieslak <[email protected]>

qsirecon/interfaces/anatomical.py

@@ -25,15 +25,14 @@
     traits,
 )
 from nipype.utils.filemanip import fname_presuffix
-from pkg_resources import resource_filename as pkgrf
 
 from ..utils.ingress import ukb_dirname_to_bids
 from .images import to_lps
 
 LOGGER = logging.getLogger("nipype.interface")
 
 
-class QSIReconAnatomicalIngressInputSpec(BaseInterfaceInputSpec):
+class QSIPrepAnatomicalIngressInputSpec(BaseInterfaceInputSpec):
     recon_input_dir = traits.Directory(
         exists=True, mandatory=True, help="directory containing subject results directories"
     )
@@ -42,7 +41,7 @@ class QSIReconAnatomicalIngressInputSpec(BaseInterfaceInputSpec):
     infant_mode = traits.Bool(mandatory=True)
 
 
-class QSIReconAnatomicalIngressOutputSpec(TraitedSpec):
+class QSIPrepAnatomicalIngressOutputSpec(TraitedSpec):
     # sub-1_desc-aparcaseg_dseg.nii.gz
     t1_aparc = File()
     # sub-1_dseg.nii.gz
@@ -65,19 +64,17 @@ class QSIReconAnatomicalIngressOutputSpec(TraitedSpec):
     t1_2_mni_reverse_transform = File()
     # sub-1_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5
     t1_2_mni_forward_transform = File()
-    # Generic: what template was used?
-    template_image = File()
 
 
-class QSIReconAnatomicalIngress(SimpleInterface):
-    """Get only the useful files from a QSIRecon anatomical output.
+class QSIPrepAnatomicalIngress(SimpleInterface):
+    """Get only the useful files from a QSIPrep anatomical output.
 
     Many of the preprocessed outputs aren't useful for reconstruction
     (mainly anything that has been mapped forward into template space).
     """
 
-    input_spec = QSIReconAnatomicalIngressInputSpec
-    output_spec = QSIReconAnatomicalIngressOutputSpec
+    input_spec = QSIPrepAnatomicalIngressInputSpec
+    output_spec = QSIPrepAnatomicalIngressOutputSpec
 
     def _run_interface(self, runtime):
         # The path to the output from the qsirecon run
@@ -137,15 +134,6 @@ def _run_interface(self, runtime):
             "t1_2_mni_forward_transform",
             "%s/sub-%s*_from-T*w_to-MNI152NLin2009cAsym_mode-image_xfm.h5" % (anat_root, sub),
         )
-        if not self.inputs.infant_mode:
-            self._results["template_image"] = pkgrf(
-                "qsirecon",
-                "data/mni_1mm_t1w_lps_brain.nii.gz",
-            )
-        else:
-            self._results["template_image"] = pkgrf(
-                "qsirecon", "data/mni_1mm_t1w_lps_brain_infant.nii.gz"
-            )
 
         return runtime
 
@@ -163,13 +151,13 @@ def _get_if_exists(self, name, pattern, excludes=None):
             self._results[name] = files[0]
 
 
-class UKBAnatomicalIngressInputSpec(QSIReconAnatomicalIngressInputSpec):
+class UKBAnatomicalIngressInputSpec(QSIPrepAnatomicalIngressInputSpec):
     recon_input_dir = traits.Directory(
         exists=True, mandatory=True, help="directory containing a single subject's results"
     )
 
 
-class UKBAnatomicalIngress(QSIReconAnatomicalIngress):
+class UKBAnatomicalIngress(QSIPrepAnatomicalIngress):
     input_spec = UKBAnatomicalIngressInputSpec
 
     def _run_interface(self, runtime):
@@ -344,3 +332,50 @@ def _run_interface(self, runtime):
 
         self._results["voxel_size"] = voxel_size
         return runtime
+
+
+class _GetTemplateInputSpec(BaseInterfaceInputSpec):
+    template_name = traits.Enum(
+        "MNI152NLin2009cAsym",
+        "MNIInfant",
+        mandatory=True,
+    )
+
+
+class _GetTemplateOutputSpec(BaseInterfaceInputSpec):
+    template_file = File(exists=True)
+    mask_file = File(exists=True)
+
+
+class GetTemplate(SimpleInterface):
+    input_spec = _GetTemplateInputSpec
+    output_spec = _GetTemplateOutputSpec
+
+    def _run_interface(self, runtime):
+        from templateflow.api import get as get_template
+
+        template_file = str(
+            get_template(
+                self.inputs.template_name,
+                cohort=[None, "2"],
+                resolution="1",
+                desc=None,
+                suffix="T1w",
+                extension=".nii.gz",
+            ),
+        )
+        mask_file = str(
+            get_template(
+                self.inputs.template_name,
+                cohort=[None, "2"],
+                resolution="1",
+                desc="brain",
+                suffix="mask",
+                extension=".nii.gz",
+            ),
+        )
+
+        self._results["template_file"] = template_file
+        self._results["mask_file"] = mask_file
+
+        return runtime

qsirecon/interfaces/ingress.py

@@ -43,7 +43,7 @@
 LOGGER = logging.getLogger("nipype.interface")
 
 
-class QsiReconDWIIngressInputSpec(BaseInterfaceInputSpec):
+class QSIPrepDWIIngressInputSpec(BaseInterfaceInputSpec):
     # DWI files
     dwi_file = File(exists=True)
     bval_file = File(exists=True)
@@ -52,7 +52,7 @@ class QsiReconDWIIngressInputSpec(BaseInterfaceInputSpec):
     atlas_names = traits.List()
 
 
-class QsiReconDWIIngressOutputSpec(TraitedSpec):
+class QSIPrepDWIIngressOutputSpec(TraitedSpec):
     subject_id = traits.Str()
     session_id = traits.Str()
     space_id = traits.Str()
@@ -73,9 +73,9 @@ class QsiReconDWIIngressOutputSpec(TraitedSpec):
     slice_qc_file = File(exists=True)
 
 
-class QsiReconDWIIngress(SimpleInterface):
-    input_spec = QsiReconDWIIngressInputSpec
-    output_spec = QsiReconDWIIngressOutputSpec
+class QSIPrepDWIIngress(SimpleInterface):
+    input_spec = QSIPrepDWIIngressInputSpec
+    output_spec = QSIPrepDWIIngressOutputSpec
 
     def _run_interface(self, runtime):
         params = get_bids_params(self.inputs.dwi_file)
@@ -121,7 +121,7 @@ def _get_qc_filename(self, out_root, params, desc, suffix):
         return out_root + "/" + fname + "_desc-%s_dwi.%s" % (desc, suffix)
 
 
-class _UKBioBankDWIIngressInputSpec(QsiReconDWIIngressInputSpec):
+class _UKBioBankDWIIngressInputSpec(QSIPrepDWIIngressInputSpec):
     dwi_file = File(exists=False, help="The name of what a BIDS dwi file may have been")
     data_dir = traits.Directory(
         exists=True, help="The UKB data directory for a subject. Must contain DTI/ and T1/"
@@ -130,7 +130,7 @@ class _UKBioBankDWIIngressInputSpec(QsiReconDWIIngressInputSpec):
 
 class UKBioBankDWIIngress(SimpleInterface):
     input_spec = _UKBioBankDWIIngressInputSpec
-    output_spec = QsiReconDWIIngressOutputSpec
+    output_spec = QSIPrepDWIIngressOutputSpec
 
     def _run_interface(self, runtime):
         runpath = Path(runtime.cwd)

qsirecon/interfaces/interchange.py

@@ -5,8 +5,8 @@
     traits,
 )
 
-from qsirecon.interfaces.anatomical import QSIReconAnatomicalIngress
-from qsirecon.interfaces.ingress import QsiReconDWIIngress
+from qsirecon.interfaces.anatomical import QSIPrepAnatomicalIngress
+from qsirecon.interfaces.ingress import QSIPrepDWIIngress
 
 # Anatomical (t1w/t2w) slots
 FS_FILES_TO_REGISTER = [
@@ -23,9 +23,9 @@
     "fs_to_qsiprep_transform_mrtrix",
 ]
 
-# These come directly from QSIRecon outputs. They're aligned to the DWIs in AC-PC
+# These come directly from QSIPrep outputs. They're aligned to the DWIs in AC-PC
 qsiprep_highres_anatomical_ingressed_fields = (
-    QSIReconAnatomicalIngress.output_spec.class_editable_traits()
+    QSIPrepAnatomicalIngress.output_spec.class_editable_traits()
 )
 
 # The init_recon_anatomical anatomical workflow can create additional
@@ -38,7 +38,7 @@
 )
 
 # These are read directly from QSIRecon's dwi results.
-qsiprep_output_names = QsiReconDWIIngress().output_spec.class_editable_traits()
+qsiprep_output_names = QSIPrepDWIIngress().output_spec.class_editable_traits()
 
 # dMRI + registered anatomical fields
 recon_workflow_anatomical_input_fields = anatomical_workflow_outputs + [

qsirecon/workflows/base.py

@@ -76,7 +76,7 @@ def init_single_subject_recon_wf(subject_id):
     subject_id : str
         Single subject label
     """
-    from ..interfaces.ingress import QsiReconDWIIngress, UKBioBankDWIIngress
+    from ..interfaces.ingress import QSIPrepDWIIngress, UKBioBankDWIIngress
     from ..interfaces.interchange import (
         ReconWorkflowInputs,
         anatomical_workflow_outputs,
@@ -152,7 +152,7 @@ def init_single_subject_recon_wf(subject_id):
         # Get the preprocessed DWI and all the related preprocessed images
         if config.workflow.recon_input_pipeline == "qsiprep":
             dwi_ingress_nodes[dwi_file] = pe.Node(
-                QsiReconDWIIngress(dwi_file=dwi_file), name=wf_name + "_ingressed_dwi_data"
+                QSIPrepDWIIngress(dwi_file=dwi_file), name=wf_name + "_ingressed_dwi_data"
             )
 
         elif config.workflow.recon_input_pipeline == "ukb":

qsirecon/workflows/recon/anatomical.py

@@ -11,12 +11,12 @@
 from nipype.interfaces.base import traits
 from nipype.pipeline import engine as pe
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
-from niworkflows.interfaces.reportlets.registration import SpatialNormalizationRPT
 from pkg_resources import resource_filename as pkgrf
 
 from ... import config
 from ...interfaces.anatomical import (
-    QSIReconAnatomicalIngress,
+    GetTemplate,
+    QSIPrepAnatomicalIngress,
     UKBAnatomicalIngress,
     VoxelSizeChooser,
 )
@@ -272,18 +272,18 @@ def gather_qsiprep_anatomical_data(subject_id):
         "has_freesurfer": False,
     }
     recon_input_dir = config.execution.bids_dir
-    # Check to see if we have a T1w preprocessed by QSIRecon
+    # Check to see if we have a T1w preprocessed by QSIPrep
     missing_qsiprep_anats = check_qsiprep_anatomical_outputs(recon_input_dir, subject_id, "T1w")
     has_qsiprep_t1w = not missing_qsiprep_anats
     status["has_qsiprep_t1w"] = has_qsiprep_t1w
     if missing_qsiprep_anats:
         config.loggers.workflow.info(
-            "Missing T1w QSIRecon outputs found: %s", " ".join(missing_qsiprep_anats)
+            "Missing T1w QSIPrep outputs found: %s", " ".join(missing_qsiprep_anats)
         )
     else:
-        config.loggers.workflow.info("Found usable QSIRecon-preprocessed T1w image and mask.")
+        config.loggers.workflow.info("Found usable QSIPrep-preprocessed T1w image and mask.")
     anat_ingress = pe.Node(
-        QSIReconAnatomicalIngress(subject_id=subject_id, recon_input_dir=recon_input_dir),
+        QSIPrepAnatomicalIngress(subject_id=subject_id, recon_input_dir=recon_input_dir),
         name="qsiprep_anat_ingress",
     )
 
@@ -296,7 +296,7 @@ def gather_qsiprep_anatomical_data(subject_id):
 
     if missing_qsiprep_transforms:
         config.loggers.workflow.info(
-            "Missing T1w QSIRecon outputs: %s", " ".join(missing_qsiprep_transforms)
+            "Missing T1w QSIPrep outputs: %s", " ".join(missing_qsiprep_transforms)
         )
 
     return anat_ingress, status
@@ -523,13 +523,32 @@ def _get_status():
             "has_qsiprep_t1w_transforms": has_qsiprep_t1w_transforms,
         }
 
+    # XXX: This is a temporary solution until QSIRecon supports flexible output spaces.
+    get_template = pe.Node(
+        GetTemplate(
+            template_name="MNI152NLin2009cAsym" if not config.workflow.infant else "MNIInfant",
+        ),
+        name="get_template",
+    )
+    mask_template = pe.Node(
+        afni.Calc(expr="a*b", outputtype="NIFTI_GZ"),
+        name="mask_template",
+    )
+    reorient_to_lps = pe.Node(
+        afni.Resample(orientation="RAI", outputtype="NIFTI_GZ"),
+        name="reorient_to_lps",
+    )
+
     reference_grid_wf = init_output_grid_wf()
     workflow.connect([
-        (inputnode, reference_grid_wf, [
-            ('template_image', 'inputnode.template_image'),
-            ('dwi_ref', 'inputnode.input_image')]),
-        (reference_grid_wf, buffernode, [
-            ('outputnode.grid_image', 'resampling_template')])
+        (get_template, mask_template, [
+            ('template_file', 'in_file_a'),
+            ('mask_file', 'in_file_b'),
+        ]),
+        (mask_template, reorient_to_lps, [('out_file', 'in_file')]),
+        (inputnode, reference_grid_wf, [('dwi_ref', 'inputnode.input_image')]),
+        (reorient_to_lps, reference_grid_wf, [('out_file', 'inputnode.template_image')]),
+        (reference_grid_wf, buffernode, [('outputnode.grid_image', 'resampling_template')]),
     ])  # fmt:skip
 
     # Missing Freesurfer AND QSIRecon T1ws, or the user wants a DWI-based mask
@@ -818,46 +837,6 @@ def _get_resampled(atlas_configs, atlas_name, to_retrieve):
     return atlas_configs[atlas_name][to_retrieve]
 
 
-def get_t1w_registration_node(infant_mode, sloppy, omp_nthreads):
-
-    # Gets an ants interface for t1w-based normalization
-    if sloppy:
-        config.loggers.workflow.info("Using QuickSyN")
-        # Requires a warp file: make an inaccurate one
-        settings = pkgrf("qsirecon", "data/quick_syn.json")
-        t1_2_mni = pe.Node(
-            SpatialNormalizationRPT(
-                float=True,
-                generate_report=True,
-                settings=[settings],
-            ),
-            name="t1_2_mni",
-            n_procs=omp_nthreads,
-            mem_gb=2,
-        )
-    else:
-        t1_2_mni = pe.Node(
-            SpatialNormalizationRPT(
-                float=True,
-                generate_report=True,
-                flavor="precise",
-            ),
-            name="t1_2_mni",
-            n_procs=omp_nthreads,
-            mem_gb=2,
-        )
-        # Get the template image
-    if not infant_mode:
-        ref_img_brain = pkgrf("qsirecon", "data/mni_1mm_t1w_lps_brain.nii.gz")
-    else:
-        ref_img_brain = pkgrf("qsirecon", "data/mni_1mm_t1w_lps_brain_infant.nii.gz")
-
-    t1_2_mni.inputs.template = "MNI152NLin2009cAsym"
-    t1_2_mni.inputs.reference_image = ref_img_brain
-    t1_2_mni.inputs.orientation = "LPS"
-    return t1_2_mni
-
-
 def init_output_grid_wf() -> Workflow:
     """Generate a non-oblique, uniform voxel-size grid around a brain."""
     workflow = Workflow(name="output_grid_wf")

qsirecon/workflows/reports.py

@@ -12,7 +12,7 @@
 
 from .. import config
 from ..interfaces import DerivativesDataSink
-from ..interfaces.ingress import QsiReconDWIIngress
+from ..interfaces.ingress import QSIPrepDWIIngress
 from ..interfaces.interchange import qsiprep_output_names, recon_workflow_input_fields
 from ..interfaces.reports import InteractiveReport
 
@@ -86,7 +86,7 @@ def init_single_subject_json_report_wf(subject_id, name):
         niu.IdentityInterface(fields=recon_workflow_input_fields), name="inputnode"
     )
     qsirecon_preprocessed_dwi_data = pe.Node(
-        QsiReconDWIIngress(), name="qsirecon_preprocessed_dwi_data"
+        QSIPrepDWIIngress(), name="qsirecon_preprocessed_dwi_data"
     )
 
     # For doctests