ENH: Output thickness, curvature, and sulcal depth files

.circleci/ds005_outputs.txt

@@ -19,14 +19,20 @@ smriprep/sub-01/anat/sub-01_desc-preproc_T1w.nii.gz
 smriprep/sub-01/anat/sub-01_dseg.nii.gz
 smriprep/sub-01/anat/sub-01_from-fsnative_to-T1w_mode-image_xfm.txt
 smriprep/sub-01/anat/sub-01_from-T1w_to-fsnative_mode-image_xfm.txt
+smriprep/sub-01/anat/sub-01_hemi-L_curv.shape.gii
 smriprep/sub-01/anat/sub-01_hemi-L_inflated.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_midthickness.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_pial.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_smoothwm.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-L_sulc.shape.gii
+smriprep/sub-01/anat/sub-01_hemi-L_thickness.shape.gii
+smriprep/sub-01/anat/sub-01_hemi-R_curv.shape.gii
 smriprep/sub-01/anat/sub-01_hemi-R_inflated.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_midthickness.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_pial.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_smoothwm.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-R_sulc.shape.gii
+smriprep/sub-01/anat/sub-01_hemi-R_thickness.shape.gii
 smriprep/sub-01/anat/sub-01_label-CSF_probseg.nii.gz
 smriprep/sub-01/anat/sub-01_label-GM_probseg.nii.gz
 smriprep/sub-01/anat/sub-01_label-WM_probseg.nii.gz

setup.cfg

@@ -23,7 +23,7 @@ install_requires =
     matplotlib >= 2.2.0
     nibabel >= 4.0.1
     nipype >= 1.7.0
-    niworkflows ~= 1.6.0
+    niworkflows @ git+https://github.com/nipreps/niworkflows.git@enh/path2bids-regex
     numpy
     packaging
     pybids >= 0.11.1

smriprep/data/io_spec.json

@@ -71,6 +71,12 @@
         "hemi": ["L", "R"],
         "extension": "surf.gii",
         "suffix": [ "inflated", "midthickness", "pial", "smoothwm"]
+      },
+      "morphometrics": {
+        "datatype": "anat",
+        "hemi": ["L", "R"],
+        "extension": "shape.gii",
+        "suffix": ["thickness", "sulc", "curv"]
       }
     }
   },

smriprep/interfaces/freesurfer.py

@@ -24,7 +24,7 @@
 import os
 from nipype import logging
 from nipype.utils.filemanip import check_depends
-from nipype.interfaces.base import traits, InputMultiObject, isdefined
+from nipype.interfaces.base import traits, InputMultiObject, isdefined, File
 from nipype.interfaces import freesurfer as fs
 
 iflogger = logging.getLogger("nipype.interface")
@@ -199,3 +199,94 @@ def _format_arg(self, name, trait_spec, value):
         if name == "hemi":
             return trait_spec.argstr % value
         return super()._format_arg(name, trait_spec, value)
+
+
+class _MRIsConvertDataInputSpec(fs.utils.MRIsConvertInputSpec):
+    in_file = File(
+        exists=True,
+        position=-2,
+        genfile=True,
+        argstr="%s",
+        desc="Surface file",
+    )
+    _xor = ('annot_file', 'parcstats_file', 'label_file', 'scalarcurv_file', 'functional_file')
+    annot_file = File(
+        exists=True,
+        argstr="--annot %s",
+        mandatory=True,
+        xor=_xor,
+        desc="input is annotation or gifti label data",
+    )
+
+    parcstats_file = File(
+        exists=True,
+        argstr="--parcstats %s",
+        mandatory=True,
+        xor=_xor,
+        desc="infile is name of text file containing label/val pairs",
+    )
+
+    label_file = File(
+        exists=True,
+        argstr="--label %s",
+        mandatory=True,
+        xor=_xor,
+        desc="infile is .label file, label is name of this label",
+    )
+
+    scalarcurv_file = File(
+        exists=True,
+        argstr="-c %s",
+        mandatory=True,
+        xor=_xor,
+        desc="input is scalar curv overlay file (must still specify surface)",
+    )
+
+    functional_file = File(
+        exists=True,
+        argstr="-f %s",
+        mandatory=True,
+        xor=_xor,
+        desc="input is functional time-series or other multi-frame data (must specify surface)",
+    )
+    target_surface = traits.Enum(
+        "white",
+        "smoothwm",
+        usedefault=True,
+        mandatory=True,
+        desc="Target surface type to return if deriving",
+    )
+
+
+class MRIsConvertData(fs.utils.MRIsConvert):
+    """Convert surface data files (label, curvature, functional, etc)
+    Wraps mris_convert to automatically select the correct ?h.white surface if
+    passed a file from the subject's surf/ directory
+    """
+    input_spec = _MRIsConvertDataInputSpec
+
+    def _gen_filename(self, name):
+        if name == "in_file":
+            if isdefined(self.inputs.in_file):
+                return self.inputs.in_file
+
+            # Find file we're trying to convert
+            fname = None
+            for opt in ('annot', 'parcstats', 'label', 'scalarcurv', 'functional'):
+                input_file = getattr(self.inputs, f"{opt}_file")
+                if isdefined(input_file):
+                    fname = input_file
+                    break
+
+            if fname is None:
+                raise ValueError("Missing file to derive filename from.")
+
+            # $SUB/lh.curv -> $SUB/lh.white, etc
+            dirname, basename = os.path.split(fname)
+            hemi = basename.split('.', 1)[0]
+            if hemi not in ('lh', 'rh'):
+                return None
+            self.inputs.in_file = os.path.join(dirname, f"{hemi}.{self.inputs.target_surface}")
+            return self.inputs.in_file
+
+        return super()._gen_filename(name)

smriprep/utils/bids.py

@@ -40,7 +40,7 @@ def get_outputnode_spec():
     'anat2std_xfm', 'std2anat_xfm',
     't1w_aseg', 't1w_aparc',
     't1w2fsnative_xfm', 'fsnative2t1w_xfm',
-    'surfaces']
+    'surfaces', 'morphometrics']
 
     """
     spec = loads(Path(pkgrf("smriprep", "data/io_spec.json")).read_text())["queries"]

smriprep/workflows/anatomical.py

@@ -198,6 +198,8 @@ def init_anat_preproc_wf(
         subject space to T1w
     surfaces
         GIFTI surfaces (gray/white boundary, midthickness, pial, inflated)
+    morphometrics
+        GIFTIs of cortical thickness, curvature, and sulcal depth
 
     See Also
     --------
@@ -545,6 +547,7 @@ def _check_img(img):
             ('outputnode.t1w2fsnative_xfm', 't1w2fsnative_xfm'),
             ('outputnode.fsnative2t1w_xfm', 'fsnative2t1w_xfm'),
             ('outputnode.surfaces', 'surfaces'),
+            ('outputnode.morphometrics', 'morphometrics'),
             ('outputnode.out_aseg', 't1w_aseg'),
             ('outputnode.out_aparc', 't1w_aparc')]),
         (applyrefined, buffernode, [('out_file', 't1w_brain')]),
@@ -561,6 +564,7 @@ def _check_img(img):
             ('t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'),
             ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'),
             ('surfaces', 'inputnode.surfaces'),
+            ('morphometrics', 'inputnode.morphometrics'),
         ]),
     ])
     # fmt:on

smriprep/workflows/outputs.py

@@ -262,6 +262,8 @@ def init_anat_derivatives_wf(
         subject space to T1w
     surfaces
         GIFTI surfaces (gray/white boundary, midthickness, pial, inflated)
+    morphometrics
+        GIFTIs of cortical thickness, curvature, and sulcal depth
     t1w_fs_aseg
         FreeSurfer's aseg segmentation, in native T1w space
     t1w_fs_aparc
@@ -287,6 +289,7 @@ def init_anat_derivatives_wf(
                 "t1w2fsnative_xfm",
                 "fsnative2t1w_xfm",
                 "surfaces",
+                "morphometrics",
                 "t1w_fs_aseg",
                 "t1w_fs_aparc",
             ]
@@ -599,6 +602,16 @@ def init_anat_derivatives_wf(
         name="ds_surfs",
         run_without_submitting=True,
     )
+    # Morphometrics
+    name_morphs = pe.MapNode(
+        Path2BIDS(), iterfield="in_file", name="name_morphs", run_without_submitting=True,
+    )
+    ds_morphs = pe.MapNode(
+        DerivativesDataSink(base_directory=output_dir, extension=".shape.gii"),
+        iterfield=["in_file", "hemi", "suffix"],
+        name="ds_morphs",
+        run_without_submitting=True,
+    )
     # Parcellations
     ds_t1w_fsaseg = pe.Node(
         DerivativesDataSink(
@@ -628,6 +641,11 @@ def init_anat_derivatives_wf(
                                ('source_files', 'source_file')]),
         (name_surfs, ds_surfs, [('hemi', 'hemi'),
                                 ('suffix', 'suffix')]),
+        (inputnode, name_morphs, [('morphometrics', 'in_file')]),
+        (inputnode, ds_morphs, [('morphometrics', 'in_file'),
+                                ('source_files', 'source_file')]),
+        (name_morphs, ds_morphs, [('hemi', 'hemi'),
+                                  ('suffix', 'suffix')]),
         (inputnode, ds_t1w_fsaseg, [('t1w_fs_aseg', 'in_file'),
                                     ('source_files', 'source_file')]),
         (inputnode, ds_t1w_fsparc, [('t1w_fs_aparc', 'in_file'),

smriprep/workflows/surfaces.py

@@ -36,7 +36,6 @@
 )
 
 from ..interfaces.freesurfer import ReconAll
-from ..interfaces.surf import NormalizeSurf
 
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 from niworkflows.interfaces.freesurfer import (
@@ -159,6 +158,8 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"):
         FreeSurfer's aseg segmentation, in native T1w space
     out_aparc
         FreeSurfer's aparc+aseg segmentation, in native T1w space
+    morphometrics
+        GIFTIs of cortical thickness, curvature, and sulcal depth
 
     See also
     --------
@@ -203,6 +204,7 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"):
                 "out_brainmask",
                 "out_aseg",
                 "out_aparc",
+                "morphometrics",
             ]
         ),
         name="outputnode",
@@ -270,8 +272,6 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"):
         # reoriented image
         (inputnode, fsnative2t1w_xfm, [('t1w', 'target_file')]),
         (autorecon1, fsnative2t1w_xfm, [('T1', 'source_file')]),
-        (fsnative2t1w_xfm, gifti_surface_wf, [
-            ('out_reg_file', 'inputnode.fsnative2t1w_xfm')]),
         (fsnative2t1w_xfm, t1w2fsnative_xfm, [('out_reg_file', 'in_lta')]),
         # Refine ANTs mask, deriving new mask from FS' aseg
         (inputnode, refine, [('corrected_t1', 'in_anat'),
@@ -291,7 +291,8 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"):
         # Output
         (autorecon_resume_wf, outputnode, [('outputnode.subjects_dir', 'subjects_dir'),
                                            ('outputnode.subject_id', 'subject_id')]),
-        (gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces')]),
+        (gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces'),
+                                        ('outputnode.morphometrics', 'morphometrics')]),
         (t1w2fsnative_xfm, outputnode, [('out_lta', 't1w2fsnative_xfm')]),
         (fsnative2t1w_xfm, outputnode, [('out_reg_file', 'fsnative2t1w_xfm')]),
         (refine, outputnode, [('out_file', 'out_brainmask')]),
@@ -504,15 +505,20 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"):
     surfaces
         GIFTI surfaces for gray/white matter boundary, pial surface,
         midthickness (or graymid) surface, and inflated surfaces
+    morphometrics
+        GIFTIs of cortical thickness, curvature, and sulcal depth
 
     """
+    from ..interfaces.freesurfer import MRIsConvertData
+    from ..interfaces.surf import NormalizeSurf
+
     workflow = Workflow(name=name)
 
     inputnode = pe.Node(
         niu.IdentityInterface(["subjects_dir", "subject_id", "fsnative2t1w_xfm"]),
         name="inputnode",
     )
-    outputnode = pe.Node(niu.IdentityInterface(["surfaces"]), name="outputnode")
+    outputnode = pe.Node(niu.IdentityInterface(["surfaces", "morphometrics"]), name="outputnode")
 
     get_surfaces = pe.Node(nio.FreeSurferSource(), name="get_surfaces")
 
@@ -532,9 +538,18 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"):
         run_without_submitting=True,
     )
     fs2gii = pe.MapNode(
-        fs.MRIsConvert(out_datatype="gii", to_scanner=True), iterfield="in_file", name="fs2gii"
+        fs.MRIsConvert(out_datatype="gii", to_scanner=True), iterfield="in_file", name="fs2gii",
     )
     fix_surfs = pe.MapNode(NormalizeSurf(), iterfield="in_file", name="fix_surfs")
+    surfmorph_list = pe.Node(
+        niu.Merge(3, ravel_inputs=True),
+        name="surfmorph_list",
+        run_without_submitting=True,
+    )
+    fscurv2funcgii = pe.MapNode(
+        MRIsConvertData(out_datatype="gii", to_scanner=True, target_surface="smoothwm"),
+        iterfield="scalarcurv_file", name="fscurv2funcgii",
+    )
 
     # fmt:off
     workflow.connect([
@@ -555,6 +570,11 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"):
         (fs2gii, fix_surfs, [('converted', 'in_file')]),
         (inputnode, fix_surfs, [('fsnative2t1w_xfm', 'transform_file')]),
         (fix_surfs, outputnode, [('out_file', 'surfaces')]),
+        (get_surfaces, surfmorph_list, [('thickness', 'in1'),
+                                        ('sulc', 'in2'),
+                                        ('curv', 'in3')]),
+        (surfmorph_list, fscurv2funcgii, [('out', 'scalarcurv_file')]),
+        (fscurv2funcgii, outputnode, [('converted', 'morphometrics')]),
     ])
     # fmt:on
     return workflow