Use standard-space anatomical brain mask (#1204)

docs/usage.rst

@@ -165,11 +165,14 @@ We recommend NOT setting the datatype, suffix, or file extension in the filter f
 ``"t1w"`` selects a native T1w-space, preprocessed T1w file.
 
 ``"t2w"`` selects a native T1w-space, preprocessed T2w file.
+If not T1w file is available, this file will be in T2w space.
 
-``"anat_dseg"`` selects a native T1w-space segmentation file.
-This file is primarily used for figures.
+``"anat_dseg"`` selects a segmentation file in the same space as the BOLD data.
+This file is not used for anything.
 
-``"anat_brainmask"`` selects a native T1w-space brain mask.
+``"anat_brainmask"`` selects an anatomically-derived brain mask in the same space as the BOLD data.
+This file is used (1) to estimate head radius for FD calculation and
+(2) to calculate coregistration quality metrics.
 
 ``"anat_to_template_xfm"`` selects a transform from T1w (or T2w, if no T1w image is available)
 space to standard space.

xcp_d/interfaces/plotting.py

@@ -247,10 +247,14 @@ class _QCPlotsInputSpec(BaseInterfaceInputSpec):
 
     # Inputs used only for nifti data
     seg_file = File(exists=True, mandatory=False, desc="Seg file for nifti")
-    anat_brainmask = File(exists=True, mandatory=False, desc="Mask in T1W")
+    anat_brainmask = File(
+        exists=True,
+        mandatory=False,
+        desc="Anatomically-derived brain mask in anatomical space.",
+    )
     template_mask = File(exists=True, mandatory=False, desc="Template mask")
-    bold2T1w_mask = File(exists=True, mandatory=False, desc="Bold mask in MNI")
-    bold2temp_mask = File(exists=True, mandatory=False, desc="Bold mask in T1W")
+    bold2T1w_mask = File(exists=True, mandatory=False, desc="BOLD mask in MNI space")
+    bold2temp_mask = File(exists=True, mandatory=False, desc="BOLD mask in anatomical space")
 
 
 class _QCPlotsOutputSpec(TraitedSpec):

xcp_d/utils/bids.py

@@ -199,10 +199,9 @@ def collect_data(
             "to": ["T1w", "T2w"],
             "suffix": "xfm",
         },
-        # native T1w-space brain mask
+        # brain mask in same standard space as BOLD data
         "anat_brainmask": {
             "datatype": "anat",
-            "space": None,
             "desc": "brain",
             "suffix": "mask",
             "extension": ".nii.gz",
@@ -289,10 +288,12 @@ def collect_data(
 
         queries["anat_to_template_xfm"]["to"] = volspace
         queries["template_to_anat_xfm"]["from"] = volspace
+        queries["anat_brainmask"]["space"] = volspace
     else:
         # use the BOLD file's space if the BOLD file is a nifti.
         queries["anat_to_template_xfm"]["to"] = queries["bold"]["space"]
         queries["template_to_anat_xfm"]["from"] = queries["bold"]["space"]
+        queries["anat_brainmask"]["space"] = queries["bold"]["space"]
 
     # Grab the first (and presumably best) density and resolution if there are multiple.
     # This probably works well for resolution (1 typically means 1x1x1,

xcp_d/utils/qcmetrics.py

@@ -26,9 +26,9 @@ def compute_registration_qc(bold2t1w_mask, anat_brainmask, bold2template_mask, t
     Parameters
     ----------
     bold2t1w_mask : :obj:`str`
-        Path to the BOLD mask in T1w space.
+        Path to the BOLD mask in anatomical (T1w or T2w) space.
     anat_brainmask : :obj:`str`
-        Path to the T1w mask.
+        Path to the anatomically-derived brain mask in anatomical space.
     bold2template_mask : :obj:`str`
         Path to the BOLD mask in template space.
     template_mask : :obj:`str`

xcp_d/workflows/base.py

@@ -20,6 +20,7 @@
 
 from xcp_d import config
 from xcp_d.__about__ import __version__
+from xcp_d.interfaces.ants import ApplyTransforms
 from xcp_d.interfaces.bids import DerivativesDataSink
 from xcp_d.interfaces.report import AboutSummary, SubjectSummary
 from xcp_d.utils.bids import (
@@ -132,6 +133,7 @@ def init_single_subject_wf(subject_id: str):
     )
     t1w_available = subj_data["t1w"] is not None
     t2w_available = subj_data["t2w"] is not None
+    anat_mod = "t1w" if t1w_available else "t2w"
 
     mesh_available, standard_space_mesh, mesh_files = collect_mesh_data(
         layout=config.execution.layout,
@@ -156,7 +158,7 @@ def init_single_subject_wf(subject_id: str):
                 "subj_data",  # not currently used, but will be in future
                 "t1w",
                 "t2w",  # optional
-                "anat_brainmask",  # not used by cifti workflow
+                "anat_brainmask",  # used to estimate head radius and for QC metrics
                 "anat_dseg",
                 "template_to_anat_xfm",  # not used by cifti workflow
                 "anat_to_template_xfm",
@@ -359,8 +361,24 @@ def init_single_subject_wf(subject_id: str):
                 ])  # fmt:skip
 
     # Estimate head radius, if necessary
+    # Need to warp the standard-space brain mask to the anatomical space to estimate head radius
+    warp_brainmask = ApplyTransforms(
+        input_image=subj_data["anat_brainmask"],
+        transforms=[subj_data["template_to_anat_xfm"]],
+        reference_image=subj_data[anat_mod],
+        num_threads=2,
+        interpolation="GenericLabel",
+        input_image_type=3,
+        dimension=3,
+    )
+    os.makedirs(config.execution.work_dir / workflow.fullname, exist_ok=True)
+    warp_brainmask_results = warp_brainmask.run(
+        cwd=(config.execution.work_dir / workflow.fullname),
+    )
+    anat_brainmask_in_anat_space = warp_brainmask_results.outputs.output_image
+
     head_radius = estimate_brain_radius(
-        mask_file=subj_data["anat_brainmask"],
+        mask_file=anat_brainmask_in_anat_space,
         head_radius=config.workflow.head_radius,
     )
 
@@ -473,6 +491,11 @@ def init_single_subject_wf(subject_id: str):
                     ]),
                 ])  # fmt:skip
 
+                # The post-processing workflow needs a native anatomical-space image as a reference
+                workflow.connect([
+                    (inputnode, postprocess_bold_wf, [(anat_mod, "inputnode.anat_native")]),
+                ])  # fmt:skip
+
             if config.workflow.combine_runs and (n_task_runs > 1):
                 for io_name, node in merge_dict.items():
                     workflow.connect([
@@ -490,6 +513,7 @@ def init_single_subject_wf(subject_id: str):
                 (inputnode, concatenate_data_wf, [
                     ("anat_brainmask", "inputnode.anat_brainmask"),
                     ("template_to_anat_xfm", "inputnode.template_to_anat_xfm"),
+                    (anat_mod, "inputnode.anat_native"),
                 ]),
             ])  # fmt:skip
 

xcp_d/workflows/bold.py

@@ -111,7 +111,7 @@ def init_postprocess_nifti_wf(
         Preprocessed T2w image, warped to standard space.
         Fed from the subject workflow.
     anat_brainmask
-        T1w brain mask, used for transforms in the QC report workflow.
+        T1w brain mask in standard space, used for transforms in the QC report workflow.
         Fed from the subject workflow.
     %(fmriprep_confounds_file)s
     fmriprep_confounds_json
@@ -159,6 +159,7 @@ def init_postprocess_nifti_wf(
                 "template_to_anat_xfm",
                 "t1w",
                 "t2w",
+                "anat_native",
                 "anat_brainmask",
                 "fmriprep_confounds_file",
                 "fmriprep_confounds_json",
@@ -325,6 +326,7 @@ def init_postprocess_nifti_wf(
             ("bold_file", "inputnode.name_source"),
             ("boldref", "inputnode.boldref"),
             ("bold_mask", "inputnode.bold_mask"),
+            ("anat_native", "inputnode.anat"),
             ("anat_brainmask", "inputnode.anat_brainmask"),
             ("template_to_anat_xfm", "inputnode.template_to_anat_xfm"),
         ]),

xcp_d/workflows/concatenation.py

@@ -64,6 +64,7 @@ def init_concatenate_data_wf(TR, head_radius, name="concatenate_data_wf"):
         Brain mask files for each of the BOLD runs.
         This will be a list of paths for NIFTI inputs, or a list of Undefineds for CIFTI ones.
     anat_brainmask : :obj:`str`
+        Anatomically-derived brain mask in the same standard space as the BOLD mask.
     %(template_to_anat_xfm)s
     %(boldref)s
     %(timeseries)s
@@ -98,6 +99,7 @@ def init_concatenate_data_wf(TR, head_radius, name="concatenate_data_wf"):
                 "smoothed_denoised_bold",
                 "bold_mask",  # only for niftis, from postproc workflows
                 "boldref",  # only for niftis, from postproc workflows
+                "anat_native",  # only for niftis, from data collection
                 "anat_brainmask",  # only for niftis, from data collection
                 "template_to_anat_xfm",  # only for niftis, from data collection
                 "timeseries",
@@ -166,6 +168,7 @@ def init_concatenate_data_wf(TR, head_radius, name="concatenate_data_wf"):
     workflow.connect([
         (inputnode, qc_report_wf, [
             ("template_to_anat_xfm", "inputnode.template_to_anat_xfm"),
+            ("anat_native", "inputnode.anat"),
             ("anat_brainmask", "inputnode.anat_brainmask"),
         ]),
         (clean_name_source, qc_report_wf, [("name_source", "inputnode.name_source")]),

xcp_d/workflows/plotting.py

@@ -61,8 +61,10 @@ def init_qc_report_wf(
     %(boldref)s
         Only used with non-CIFTI data.
     bold_mask
+        Path to the BOLD run's brain mask in the same space as ``preprocessed_bold``.
         Only used with non-CIFTI data.
     anat_brainmask
+        Path to the anatomical brain mask in the same standard space as ``bold_mask``.
         Only used with non-CIFTI data.
     %(template_to_anat_xfm)s
         Only used with non-CIFTI data.
@@ -94,6 +96,7 @@ def init_qc_report_wf(
                 "run_index",  # will only be set for concatenated data
                 # nifti-only inputs
                 "bold_mask",
+                "anat",  # T1w/T2w image in anatomical space
                 "anat_brainmask",
                 "boldref",
                 "template_to_anat_xfm",
@@ -158,7 +161,7 @@ def init_qc_report_wf(
         workflow.connect([
             (inputnode, warp_boldmask_to_t1w, [
                 ("bold_mask", "input_image"),
-                ("anat_brainmask", "reference_image"),
+                ("anat", "reference_image"),
             ]),
             (get_native2space_transforms, warp_boldmask_to_t1w, [
                 ("bold_to_t1w_xfms", "transforms"),
@@ -185,6 +188,27 @@ def init_qc_report_wf(
             ]),
         ])  # fmt:skip
 
+        # Warp the standard-space anatomical brain mask to the anatomical space
+        warp_anatmask_to_t1w = pe.Node(
+            ApplyTransforms(
+                dimension=3,
+                interpolation="NearestNeighbor",
+            ),
+            name="warp_anatmask_to_t1w",
+            n_procs=omp_nthreads,
+            mem_gb=1,
+        )
+        workflow.connect([
+            (inputnode, warp_anatmask_to_t1w, [
+                ("bold_mask", "input_image"),
+                ("anat", "reference_image"),
+            ]),
+            (get_native2space_transforms, warp_anatmask_to_t1w, [
+                ("bold_to_t1w_xfms", "transforms"),
+                ("bold_to_t1w_xfms_invert", "invert_transform_flags"),
+            ]),
+        ])  # fmt:skip
+
         # NIFTI files require a tissue-type segmentation in the same space as the BOLD data.
         # Get the set of transforms from MNI152NLin6Asym (the dseg) to the BOLD space.
         # Given that xcp-d doesn't process native-space data, this transform will never be used.
@@ -276,13 +300,11 @@ def init_qc_report_wf(
 
         if config.workflow.file_format == "nifti":
             workflow.connect([
-                (inputnode, qc_report, [
-                    ("anat_brainmask", "anat_brainmask"),
-                    ("bold_mask", "mask_file"),
-                ]),
+                (inputnode, qc_report, [("bold_mask", "mask_file")]),
                 (warp_dseg_to_bold, qc_report, [("output_image", "seg_file")]),
                 (warp_boldmask_to_t1w, qc_report, [("output_image", "bold2T1w_mask")]),
                 (warp_boldmask_to_mni, qc_report, [("output_image", "bold2temp_mask")]),
+                (warp_anatmask_to_t1w, qc_report, [("output_image", "anat_brainmask")]),
             ])  # fmt:skip
         else:
             qc_report.inputs.mask_file = None