ENH: Enable resampling morphometrics to fsLR CIFTI-2 files

setup.cfg

@@ -21,6 +21,7 @@ python_requires = >=3.8
 install_requires =
     indexed_gzip >= 0.8.8
     lockfile
+    looseversion
     matplotlib >= 2.2.0
     nibabel >= 4.0.1
     nipype >= 1.7.0

smriprep/conftest.py

@@ -1,3 +1,22 @@
 import os
+from pathlib import Path
+from tempfile import TemporaryDirectory
+import pytest
 
 os.environ['NO_ET'] = '1'
+
+
+@pytest.fixture(autouse=True, scope="session")
+def default_cwd():
+    cwd = os.getcwd()
+    with TemporaryDirectory(prefix="smriprepTest") as tmpdir:
+        try:
+            os.chdir(tmpdir)
+            yield Path(tmpdir)
+        finally:
+            os.chdir(cwd)
+
+
+@pytest.fixture(autouse=True, scope="session")
+def populate_default_cwd(default_cwd):
+    Path.write_bytes(default_cwd / 'lh.white', b'')

smriprep/interfaces/cifti.py

@@ -0,0 +1,196 @@
+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 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(f"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 = ci.Cifti2MetaData(metadata)
+
+    img = ci.Cifti2Image(dataobj=np.atleast_2d(np.concatenate(arrays)), header=hdr)
+    img.nifti_header.set_intent("NIFTI_INTENT_CONNECTIVITY_DENSE_SCALARS")
+
+    stem = Path(scalar_surfs[0]).name.split(".")[0]
+    cifti_stem = "_".join(ent for ent in stem.split("_") if not ent.startswith("hemi-"))
+    out_file = Path.cwd() / f"{cifti_stem}.dscalar.nii"
+    img.to_filename(out_file)
+    return out_file

smriprep/interfaces/freesurfer.py

@@ -22,10 +22,12 @@
 #
 """Nipype's recon-all replacement."""
 import os
+from looseversion import LooseVersion
 from nipype import logging
 from nipype.utils.filemanip import check_depends
 from nipype.interfaces.base import traits, InputMultiObject, isdefined, File
 from nipype.interfaces import freesurfer as fs
+from niworkflows.interfaces import freesurfer as nwfs
 
 iflogger = logging.getLogger("nipype.interface")
 
@@ -174,6 +176,11 @@ def cmdline(self):
                 no_run = False
                 continue
 
+            # FreeSurfer changed the meaning and order of -apas2aseg without
+            # updating the recon table on the wiki. Hack it until fixed in nipype.
+            if step == 'apas2aseg' and fs.Info.looseversion() >= LooseVersion("7.3.0"):
+                infiles = []
+
             subj_dir = os.path.join(subjects_dir, self.inputs.subject_id)
             if check_depends(
                 [os.path.join(subj_dir, f) for f in outfiles],
@@ -283,3 +290,34 @@ def _gen_filename(self, name):
             return self.inputs.in_file
 
         return super()._gen_filename(name)
+
+
+class MakeMidthickness(nwfs.MakeMidthickness):
+    """Patched MakeMidthickness interface
+
+    Ensures output filenames are specified with hemisphere labels, when appropriate.
+    This may not cover all use-cases in MRIsExpand, but we're just making midthickness
+    files.
+
+    >>> from smriprep.interfaces.freesurfer import MakeMidthickness
+    >>> mris_expand = MakeMidthickness(thickness=True, distance=0.5)
+    >>> mris_expand.inputs.in_file = 'lh.white'
+    >>> mris_expand.cmdline
+    'mris_expand -thickness lh.white 0.5 lh.expanded'
+    >>> mris_expand.inputs.out_name = 'graymid'
+    >>> mris_expand.cmdline
+    'mris_expand -thickness lh.white 0.5 lh.graymid'
+
+        Explicit hemisphere labels should still be respected:
+
+    >>> mris_expand.inputs.out_name = 'rh.graymid'
+    >>> mris_expand.cmdline
+    'mris_expand -thickness lh.white 0.5 rh.graymid'
+    """
+
+    def _format_arg(self, name, trait_spec, value):
+        # FreeSurfer at some point changed whether it would add the hemi label onto the
+        # surface. Therefore we'll do it ourselves.
+        if name == "out_name":
+            value = self._associated_file(self.inputs.in_file, value)
+        return super()._format_arg(name, trait_spec, value)

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
 
 

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