Output T2w images if available (#648)

docs/output.rst

@@ -42,6 +42,7 @@ The  ``xcp_d`` outputs are written out in BIDS format and consist of three main
           sub-<label>/[ses-<label>/]
                anat/
                     <source_entities>_space-MNI152NLin6Asym_desc-preproc_T1w.nii.gz
+                    <source_entities>_space-MNI152NLin6Asym_desc-preproc_T2w.nii.gz
                     <source_entities>_space-MNI152NLin6Asym_dseg.nii.gz
 
    If the ``--warp-surfaces-native2std`` option is selected, and reconstructed surfaces are available in the preprocessed dataset,

xcp_d/tests/conftest.py

@@ -87,6 +87,8 @@ def fmriprep_with_freesurfer_data(datasets):
         func_dir,
         "sub-01_task-rest_space-MNI152NLin2009cAsym_res-2_boldref.nii.gz",
     )
+    files["t1w"] = os.path.join(anat_dir, "sub-01_desc-preproc_T1w.nii.gz")
+    files["t1seg"] = os.path.join(anat_dir, "sub-01_desc-aseg_dseg.nii.gz")
 
     return files
 

xcp_d/tests/test_workflows_anatomical.py

@@ -0,0 +1,51 @@
+"""Tests for the xcp_d.workflows.anatomical module."""
+import os
+import shutil
+
+from xcp_d.workflows import anatomical
+
+
+def test_warp_anats_to_template_wf(fmriprep_with_freesurfer_data, tmp_path_factory):
+    """Test xcp_d.workflows.anatomical.init_warp_anats_to_template_wf."""
+    tmpdir = tmp_path_factory.mktemp("test_nifti_conn")
+
+    t1w_to_template_xform = fmriprep_with_freesurfer_data["t1w_to_template_xform"]
+    t1w = fmriprep_with_freesurfer_data["t1w"]
+    t1seg = fmriprep_with_freesurfer_data["t1seg"]
+    t2w = os.path.join(tmpdir, "sub-01_desc-preproc_T2w.nii.gz")  # pretend t1w is t2w
+    shutil.copyfile(t1w, t2w)
+
+    wf = anatomical.init_warp_anats_to_template_wf(
+        output_dir=tmpdir,
+        input_type="fmriprep",
+        t2w_available=True,
+        target_space="MNI152NLin2009cAsym",
+        omp_nthreads=1,
+        mem_gb=0.1,
+        name="warp_anats_to_template_wf",
+    )
+    wf.inputs.inputnode.t1w_to_template = t1w_to_template_xform
+    wf.inputs.inputnode.t1w = t1w
+    wf.inputs.inputnode.t1seg = t1seg
+    wf.inputs.inputnode.t2w = t2w
+    wf.base_dir = tmpdir
+    wf.run()
+
+    out_anat_dir = os.path.join(tmpdir, "xcp_d", "sub-01", "anat")
+    out_t1w = os.path.join(
+        out_anat_dir,
+        "sub-01_space-MNI152NLin2009cAsym_desc-preproc_T1w.nii.gz",
+    )
+    assert os.path.isfile(out_t1w), os.listdir(out_anat_dir)
+
+    out_t2w = os.path.join(
+        out_anat_dir,
+        "sub-01_space-MNI152NLin2009cAsym_desc-preproc_T2w.nii.gz",
+    )
+    assert os.path.isfile(out_t2w), os.listdir(out_anat_dir)
+
+    out_t1seg = os.path.join(
+        out_anat_dir,
+        "sub-01_space-MNI152NLin2009cAsym_desc-aseg_dseg.nii.gz",
+    )
+    assert os.path.isfile(out_t1seg), os.listdir(out_anat_dir)

xcp_d/utils/bids.py

@@ -190,7 +190,21 @@ def collect_data(
         # all preprocessed BOLD files in the right space/resolution/density
         "bold": {"datatype": "func", "suffix": "bold", "desc": ["preproc", None]},
         # native T1w-space, preprocessed T1w file
-        "t1w": {"datatype": "anat", "space": None, "suffix": "T1w", "extension": ".nii.gz"},
+        "t1w": {
+            "datatype": "anat",
+            "space": None,
+            "desc": "preproc",
+            "suffix": "T1w",
+            "extension": ".nii.gz",
+        },
+        # native T2w-space, preprocessed T1w file
+        "t2w": {
+            "datatype": "anat",
+            "space": [None, "T1w"],
+            "desc": "preproc",
+            "suffix": "T2w",
+            "extension": ".nii.gz",
+        },
         # native T1w-space dseg file, but not aseg or aparcaseg
         "t1w_seg": {
             "datatype": "anat",
@@ -203,7 +217,7 @@ def collect_data(
         # from entity will be set later
         "template_to_t1w_xform": {
             "datatype": "anat",
-            "to": "T1w",
+            "to": ["T1w", "T2w"],
             "suffix": "xfm",
         },
         # native T1w-space brain mask
@@ -218,7 +232,7 @@ def collect_data(
         # to entity will be set later
         "t1w_to_template_xform": {
             "datatype": "anat",
-            "from": "T1w",
+            "from": ["T1w", "T2w"],
             "suffix": "xfm",
         },
     }
@@ -227,6 +241,9 @@ def collect_data(
     else:
         queries["bold"]["extension"] = ".nii.gz"
 
+    # List any fields that aren't required
+    OPTIONAL_FIELDS = ["t2w"]
+
     # Apply filters. These may override anything.
     bids_filters = bids_filters or {}
     for acq, entities in bids_filters.items():
@@ -302,10 +319,17 @@ def collect_data(
     for field, filenames in subj_data.items():
         # All fields except the BOLD data should have a single file
         if field != "bold" and isinstance(filenames, list):
-            if not filenames:
+            if field not in OPTIONAL_FIELDS and not filenames:
                 raise FileNotFoundError(f"No {field} found with query: {queries[field]}")
 
-            subj_data[field] = filenames[0]
+            if len(filenames) == 1:
+                subj_data[field] = filenames[0]
+            elif len(filenames) > 1:
+                filenames_str = "\n\t".join(filenames)
+                LOGGER.warning(f"Multiple files found for query '{field}':\n\t{filenames_str}")
+                subj_data[field] = filenames[0]
+            else:
+                subj_data[field] = None
 
     LOGGER.debug(f"Collected data:\n{yaml.dump(subj_data, default_flow_style=False, indent=4)}")
 

xcp_d/workflows/anatomical.py

@@ -39,12 +39,13 @@
 def init_warp_anats_to_template_wf(
     output_dir,
     input_type,
+    t2w_available,
     target_space,
     omp_nthreads,
     mem_gb,
     name="warp_anats_to_template_wf",
 ):
-    """Copy T1w and segmentation to the derivative directory.
+    """Copy T1w, segmentation, and, optionally, T2w to the derivative directory.
 
     If necessary, this workflow will also warp the images to standard space.
 
@@ -57,6 +58,7 @@ def init_warp_anats_to_template_wf(
             wf = init_warp_anats_to_template_wf(
                 output_dir=".",
                 input_type="fmriprep",
+                t2w_available=True,
                 target_space="MNI152NLin6Asym",
                 omp_nthreads=1,
                 mem_gb=0.1,
@@ -67,6 +69,8 @@ def init_warp_anats_to_template_wf(
     ----------
     %(output_dir)s
     %(input_type)s
+    t2w_available : bool
+        True if a preprocessed T2w is available, False if not.
     target_space : str
         Target NIFTI template for T1w.
     %(omp_nthreads)s
@@ -77,7 +81,11 @@ def init_warp_anats_to_template_wf(
     Inputs
     ------
     t1w : str
-        Path to the T1w file.
+        Path to the preprocessed T1w file.
+        This file may be in standard space or native T1w space.
+    t2w : str
+        Path to the preprocessed T2w file.
+        This file may be in standard space or native T1w space.
     t1seg : str
         Path to the T1w segmentation file.
     %(t1w_to_template)s
@@ -86,7 +94,7 @@ def init_warp_anats_to_template_wf(
     workflow = Workflow(name=name)
 
     inputnode = pe.Node(
-        niu.IdentityInterface(fields=["t1w", "t1seg", "t1w_to_template"]),
+        niu.IdentityInterface(fields=["t1w", "t2w", "t1seg", "t1w_to_template"]),
         name="inputnode",
     )
 
@@ -105,37 +113,56 @@ def init_warp_anats_to_template_wf(
     )
 
     if input_type in ("dcan", "hcp"):
-        # Assume that the T1w and T1w segmentation files are in standard space,
+        # Assume that the T1w, T1w segmentation, and T2w files are in standard space,
         # but don't have the "space" entity, for the "dcan" and "hcp" derivatives.
         # This is a bug, and the converted filenames are inaccurate, so we have this
         # workaround in place.
-        ds_t1w = pe.Node(
+        ds_t1w_std = pe.Node(
             DerivativesDataSink(
                 base_directory=output_dir,
                 extension=".nii.gz",
             ),
-            name="ds_t1w",
+            name="ds_t1w_std",
             run_without_submitting=False,
         )
 
-        ds_t1seg = pe.Node(
+        ds_t1seg_std = pe.Node(
             DerivativesDataSink(
                 base_directory=output_dir,
                 extension=".nii.gz",
             ),
-            name="ds_t1seg",
+            name="ds_t1seg_std",
             run_without_submitting=False,
         )
 
         # fmt:off
         workflow.connect([
-            (inputnode, ds_t1w, [("t1w", "in_file")]),
-            (inputnode, ds_t1seg, [("t1seg", "in_file")]),
+            (inputnode, ds_t1w_std, [("t1w", "in_file")]),
+            (inputnode, ds_t1seg_std, [("t1seg", "in_file")]),
         ])
         # fmt:on
 
+        if t2w_available:
+            ds_t2w_std = pe.Node(
+                DerivativesDataSink(
+                    base_directory=output_dir,
+                    extension=".nii.gz",
+                ),
+                name="ds_t2w_std",
+                run_without_submitting=False,
+            )
+
+            # fmt:off
+            workflow.connect([
+                (inputnode, ds_t2w_std, [
+                    ("t2w", "in_file"),
+                    ("t2w", "source_file"),
+                ]),
+            ])
+            # fmt:on
+
     else:
-        # Warp the native T1w-space T1w and T1w segmentation files to the selected standard space.
+        # Warp the native T1w-space T1w, T1w segmentation, and T2w files to standard space.
         warp_t1w_to_template = pe.Node(
             ApplyTransforms(
                 num_threads=2,
@@ -180,42 +207,74 @@ def init_warp_anats_to_template_wf(
         ])
         # fmt:on
 
-        ds_t1w = pe.Node(
+        if t2w_available:
+            t2w_transform = pe.Node(
+                ApplyTransforms(
+                    num_threads=2,
+                    reference_image=template_file,
+                    interpolation="LanczosWindowedSinc",
+                    input_image_type=3,
+                    dimension=3,
+                ),
+                name="t2w_transform",
+                mem_gb=mem_gb,
+                n_procs=omp_nthreads,
+            )
+
+            ds_t2w_std = pe.Node(
+                DerivativesDataSink(
+                    base_directory=output_dir,
+                    space=target_space,
+                    cohort=cohort,
+                    extension=".nii.gz",
+                ),
+                name="ds_t2w_std",
+                run_without_submitting=False,
+            )
+
+            # fmt:off
+            workflow.connect([
+                (inputnode, t2w_transform, [
+                    ("t2w", "input_image"),
+                    ("t1w_to_template", "transforms"),
+                ]),
+                (t2w_transform, ds_t2w_std, [("output_image", "in_file")]),
+                (inputnode, ds_t2w_std, [("t2w", "source_file")]),
+            ])
+            # fmt:on
+
+        ds_t1w_std = pe.Node(
             DerivativesDataSink(
                 base_directory=output_dir,
                 space=target_space,
                 cohort=cohort,
                 extension=".nii.gz",
             ),
-            name="ds_t1w",
+            name="ds_t1w_std",
             run_without_submitting=False,
         )
 
-        # fmt:off
-        workflow.connect([(warp_t1w_to_template, ds_t1w, [("output_image", "in_file")])])
-        # fmt:on
+        workflow.connect([(warp_t1w_to_template, ds_t1w_std, [("output_image", "in_file")])])
 
-        ds_t1seg = pe.Node(
+        ds_t1seg_std = pe.Node(
             DerivativesDataSink(
                 base_directory=output_dir,
                 space=target_space,
                 cohort=cohort,
                 extension=".nii.gz",
             ),
-            name="ds_t1seg",
+            name="ds_t1seg_std",
             run_without_submitting=False,
         )
 
-        # fmt:off
-        workflow.connect([(warp_t1seg_to_template, ds_t1seg, [("output_image", "in_file")])])
-        # fmt:on
+        workflow.connect([(warp_t1seg_to_template, ds_t1seg_std, [("output_image", "in_file")])])
 
     # fmt:off
     workflow.connect([
-        (inputnode, ds_t1w, [("t1w", "source_file")]),
-        (inputnode, ds_t1seg, [("t1seg", "source_file")]),
-        (ds_t1w, outputnode, [("out_file", "t1w")]),
-        (ds_t1seg, outputnode, [("out_file", "t1seg")]),
+        (inputnode, ds_t1w_std, [("t1w", "source_file")]),
+        (inputnode, ds_t1seg_std, [("t1seg", "source_file")]),
+        (ds_t1w_std, outputnode, [("out_file", "t1w")]),
+        (ds_t1seg_std, outputnode, [("out_file", "t1seg")]),
     ])
     # fmt:on
 

xcp_d/workflows/base.py

@@ -373,6 +373,7 @@ def init_subject_wf(
             fields=[
                 "subj_data",  # not currently used, but will be in future
                 "t1w",
+                "t2w",  # optional
                 "t1w_mask",  # not used by cifti workflow
                 "t1w_seg",
                 "template_to_t1w_xform",  # not used by cifti workflow
@@ -395,6 +396,7 @@ def init_subject_wf(
     )
     inputnode.inputs.subj_data = subj_data
     inputnode.inputs.t1w = subj_data["t1w"]
+    inputnode.inputs.t2w = subj_data["t2w"]
     inputnode.inputs.t1w_mask = subj_data["t1w_mask"]
     inputnode.inputs.t1w_seg = subj_data["t1w_seg"]
     inputnode.inputs.template_to_t1w_xform = subj_data["template_to_t1w_xform"]
@@ -484,6 +486,7 @@ def init_subject_wf(
     warp_anats_to_template_wf = init_warp_anats_to_template_wf(
         output_dir=output_dir,
         input_type=input_type,
+        t2w_available=subj_data["t2w"] is not None,
         target_space=target_space,
         omp_nthreads=omp_nthreads,
         mem_gb=5,  # RF: need to change memory size
@@ -493,6 +496,7 @@ def init_subject_wf(
     workflow.connect([
         (inputnode, warp_anats_to_template_wf, [
             ("t1w", "inputnode.t1w"),
+            ("t2w", "inputnode.t2w"),
             ("t1w_seg", "inputnode.t1seg"),
             ("t1w_to_template_xform", "inputnode.t1w_to_template"),
         ]),