ENH: Enable resampling morphometrics to fsLR CIFTI-2 files

effigies · effigies · commit 73fa73bf819d · 2023-02-28T12:54:15.000-05:00
diff --git a/smriprep/interfaces/cifti.py b/smriprep/interfaces/cifti.py
@@ -0,0 +1,195 @@
+import json
+from pathlib import Path
+import typing as ty
+
+import nibabel as nb
+import numpy as np
+
+from nibabel import cifti2 as ci
+from nipype.interfaces.base import TraitedSpec, traits, File, SimpleInterface
+from nipype.utils.filemanip import split_filename
+from templateflow import api as tf
+
+
+class _GenerateDScalarInputSpec(TraitedSpec):
+    surface_target = traits.Enum(
+        "fsLR",
+        usedefault=True,
+        desc="CIFTI surface target space",
+    )
+    grayordinates = traits.Enum(
+        "91k", "170k", usedefault=True, desc="Final CIFTI grayordinates"
+    )
+    scalar_surfs = traits.List(
+        File(exists=True),
+        mandatory=True,
+        desc="list of surface BOLD GIFTI files (length 2 with order [L,R])",
+    )
+    scalar_name = traits.Str(mandatory=True, desc="Name of scalar")
+
+
+class _GenerateDScalarOutputSpec(TraitedSpec):
+    out_file = File(desc="generated CIFTI file")
+    out_metadata = File(desc="CIFTI metadata JSON")
+
+
+class GenerateDScalar(SimpleInterface):
+    """
+    Generate a HCP-style CIFTI-2 image from scalar surface files.
+    """
+    input_spec = _GenerateDScalarInputSpec
+    output_spec = _GenerateDScalarOutputSpec
+
+    def _run_interface(self, runtime):
+
+        surface_labels, metadata = _prepare_cifti(self.inputs.grayordinates)
+        self._results["out_file"] = _create_cifti_image(
+            self.inputs.scalar_surfs,
+            surface_labels,
+            self.inputs.scalar_name,
+            metadata,
+        )
+        metadata_file = Path("dscalar.json").absolute()
+        metadata_file.write_text(json.dumps(metadata, indent=2))
+        self._results["out_metadata"] = str(metadata_file)
+        return runtime
+
+
+def _prepare_cifti(grayordinates: str) -> ty.Tuple[list, dict]:
+    """
+    Fetch the required templates needed for CIFTI-2 generation, based on input surface density.
+
+    Parameters
+    ----------
+    grayordinates :
+        Total CIFTI grayordinates (91k, 170k)
+
+    Returns
+    -------
+    surface_labels
+        Surface label files for vertex inclusion/exclusion.
+    metadata
+        Dictionary with BIDS metadata.
+
+    Examples
+    --------
+    >>> surface_labels, metadata = _prepare_cifti('91k')
+    >>> surface_labels  # doctest: +ELLIPSIS
+    ['.../tpl-fsLR_hemi-L_den-32k_desc-nomedialwall_dparc.label.gii', \
+     '.../tpl-fsLR_hemi-R_den-32k_desc-nomedialwall_dparc.label.gii']
+    >>> metadata # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
+    {'Density': '91,282 grayordinates corresponding to all of the grey matter sampled at a \
+2mm average vertex spacing...', 'SpatialReference': {'CIFTI_STRUCTURE_CORTEX_LEFT': ...}}
+
+    """
+
+    grayord_key = {
+        "91k": {
+            "surface-den": "32k",
+            "tf-res": "02",
+            "grayords": "91,282",
+            "res-mm": "2mm"
+        },
+        "170k": {
+            "surface-den": "59k",
+            "tf-res": "06",
+            "grayords": "170,494",
+            "res-mm": "1.6mm"
+        }
+    }
+    if grayordinates not in grayord_key:
+        raise NotImplementedError("Grayordinates {grayordinates} is not supported.")
+
+    total_grayords = grayord_key[grayordinates]['grayords']
+    res_mm = grayord_key[grayordinates]['res-mm']
+    surface_density = grayord_key[grayordinates]['surface-den']
+    # Fetch templates
+    surface_labels = [
+        str(
+            tf.get(
+                "fsLR",
+                density=surface_density,
+                hemi=hemi,
+                desc="nomedialwall",
+                suffix="dparc",
+                raise_empty=True,
+            )
+        )
+        for hemi in ("L", "R")
+    ]
+
+    tf_url = "https://templateflow.s3.amazonaws.com"
+    surfaces_url = (  # midthickness is the default, but varying levels of inflation are all valid
+        f"{tf_url}/tpl-fsLR/tpl-fsLR_den-{surface_density}_hemi-%s_midthickness.surf.gii"
+    )
+    metadata = {
+        "Density": (
+            f"{total_grayords} grayordinates corresponding to all of the grey matter sampled at a "
+            f"{res_mm} average vertex spacing on the surface"
+        ),
+        "SpatialReference": {
+            "CIFTI_STRUCTURE_CORTEX_LEFT": surfaces_url % "L",
+            "CIFTI_STRUCTURE_CORTEX_RIGHT": surfaces_url % "R",
+        }
+    }
+    return surface_labels, metadata
+
+
+def _create_cifti_image(
+    scalar_surfs: ty.Tuple[str, str],
+    surface_labels: ty.Tuple[str, str],
+    scalar_name: str,
+    metadata: ty.Optional[dict] = None,
+):
+    """
+    Generate CIFTI image in target space.
+
+    Parameters
+    ----------
+    scalar_surfs
+        Surface scalar files (L,R)
+    surface_labels
+        Surface label files used to remove medial wall (L,R)
+    metadata
+        Metadata to include in CIFTI header
+    scalar_name
+        Name to apply to scalar map
+
+    Returns
+    -------
+    out :
+        BOLD data saved as CIFTI dtseries
+    """
+    brainmodels = []
+    arrays = []
+
+    for idx, hemi in enumerate(('left', 'right')):
+        labels = nb.load(surface_labels[idx])
+        mask = np.bool_[labels.darrays[0].data]
+
+        struct = f'cortex_{hemi}'
+        brainmodels.append(
+            ci.BrainModelAxis(struct, vertex=np.nonzero(mask)[0], nvertices={struct: len(mask)})
+        )
+
+        morph_scalar = nb.load(scalar_surfs[idx])
+        arrays.append(morph_scalar.darrays[0].data[mask].astype("float32"))
+
+    # provide some metadata to CIFTI matrix
+    if not metadata:
+        metadata = {
+            "surface": "fsLR",
+        }
+
+    # generate and save CIFTI image
+    hdr = ci.Cifti2Header.from_axes(
+        (ci.ScalarAxis([scalar_name]), brainmodels[0] + brainmodels[1])
+    )
+    hdr.matrix.metadata.update(metadata)
+
+    img = ci.Cifti2Image(dataobj=np.atleast_2d(np.concatenate(arrays)), header=hdr)
+    img.nifti_header.set_intent("NIFTI_INTENT_CONNECTIVITY_DENSE_SCALARS")
+
+    out_file = "{}.dscalar.nii".format(split_filename(scalar_surfs[0])[1])
+    img.to_filename(out_file)
+    return Path.cwd() / out_file
diff --git a/smriprep/workflows/anatomical.py b/smriprep/workflows/anatomical.py
@@ -51,7 +51,7 @@
 from ..utils.misc import apply_lut as _apply_bids_lut, fs_isRunning as _fs_isRunning
 from .norm import init_anat_norm_wf
 from .outputs import init_anat_reports_wf, init_anat_derivatives_wf
-from .surfaces import init_anat_ribbon_wf, init_surface_recon_wf
+from .surfaces import init_anat_ribbon_wf, init_surface_recon_wf, init_morph_grayords_wf
 
 LOGGER = logging.getLogger("nipype.workflow")
 
@@ -69,6 +69,7 @@ def init_anat_preproc_wf(
     skull_strip_mode,
     skull_strip_template,
     spaces,
+    cifti_output=False,
     debug=False,
     existing_derivatives=None,
     name="anat_preproc_wf",
@@ -464,6 +465,7 @@ def _check_img(img):
         t2w=t2w,
         output_dir=output_dir,
         spaces=spaces,
+        cifti_output=cifti_output,
     )
 
     # fmt:off
@@ -639,6 +641,22 @@ def _check_img(img):
     ])
     # fmt:on
 
+    if cifti_output:
+        morph_grayords_wf = init_morph_grayords_wf(grayord_density=cifti_output)
+        anat_derivatives_wf.get_node('inputnode').inputs.cifti_density = cifti_output
+        # fmt:off
+        workflow.connect([
+            (surface_recon_wf, morph_grayords_wf, [
+                ('outputnode.subject_id', 'inputnode.subject_id'),
+                ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
+            ]),
+            (morph_grayords_wf, anat_derivatives_wf, [
+                ("outputnode.cifti_morph", "inputnode.cifti_morph"),
+                ("outputnode.cifti_metadata", "inputnode.cifti_metadata"),
+            ]),
+        ])
+        # fmt:on
+
     return workflow
 
 
diff --git a/smriprep/workflows/base.py b/smriprep/workflows/base.py
@@ -414,6 +414,7 @@ def init_single_subject_wf(
         skull_strip_mode=skull_strip_mode,
         skull_strip_template=skull_strip_template,
         spaces=spaces,
+        cifti_output=False,  # Enabling this needs a CLI flag
     )
 
     # fmt:off
diff --git a/smriprep/workflows/outputs.py b/smriprep/workflows/outputs.py
@@ -204,6 +204,7 @@ def init_anat_derivatives_wf(
     t2w,
     output_dir,
     spaces,
+    cifti_output,
     name="anat_derivatives_wf",
     tpm_labels=BIDS_TISSUE_ORDER,
 ):
@@ -224,6 +225,8 @@ def init_anat_derivatives_wf(
         Workflow name (default: anat_derivatives_wf)
     tpm_labels : :obj:`tuple`
         Tissue probability maps in order
+    cifti_output : :obj:`bool`
+        Whether the ``--cifti-output`` flag was set.
 
     Inputs
     ------
@@ -273,6 +276,12 @@ def init_anat_derivatives_wf(
         FreeSurfer's aseg segmentation, in native T1w space
     t1w_fs_aparc
         FreeSurfer's aparc+aseg segmentation, in native T1w space
+    cifti_morph
+        Morphometric CIFTI-2 dscalar files
+    cifti_density
+        Grayordinate density
+    cifti_metadata
+        JSON files containing metadata dictionaries
 
     """
     from niworkflows.interfaces.utility import KeySelect
@@ -299,6 +308,9 @@ def init_anat_derivatives_wf(
                 "anat_ribbon",
                 "t1w_fs_aseg",
                 "t1w_fs_aparc",
+                'cifti_metadata',
+                'cifti_density',
+                'cifti_morph',
             ]
         ),
         name="inputnode",
@@ -696,6 +708,27 @@ def init_anat_derivatives_wf(
 
     ])
     # fmt:on
+
+    if cifti_output:
+        ds_cifti_morph = pe.MapNode(
+            DerivativesDataSink(
+                base_directory=output_dir,
+                suffix=['curv', 'sulc', 'thickness'],
+                compress=False,
+                space='fsLR',
+            ),
+            name='ds_cifti_morph',
+            run_without_submitting=True,
+            iterfield=["in_file", "meta_dict", "suffix"],
+        )
+        # fmt:off
+        workflow.connect([
+            (inputnode, ds_cifti_morph, [('cifti_morph', 'in_file'),
+                                         ('source_file', 'source_file'),
+                                         ('cifti_density', 'density'),
+                                         (('cifti_metadata', _read_json), 'meta_dict')])
+        ])
+        # fmt:on
     return workflow
 
 
@@ -796,3 +829,10 @@ def _combine_cohort(in_template):
             return template
         return f"{template}+{in_template.split('cohort-')[-1].split(':')[0]}"
     return [_combine_cohort(v) for v in in_template]
+
+
+def _read_json(in_file):
+    from json import loads
+    from pathlib import Path
+
+    return loads(Path(in_file).read_text())
diff --git a/smriprep/workflows/surfaces.py b/smriprep/workflows/surfaces.py

Original file line number	Diff line number	Diff line change
`@@ -414,6 +414,7 @@ def init_single_subject_wf(`
`414`	`414`	`skull_strip_mode=skull_strip_mode,`
`415`	`415`	`skull_strip_template=skull_strip_template,`
`416`	`416`	`spaces=spaces,`
	`417`	`+ cifti_output=False, # Enabling this needs a CLI flag`
`417`	`418`	`)`
`418`	`419`
`419`	`420`	`# fmt:off`