ENH: Integrate morphometric fsLR resampling

Author: mgxd

nibabies/workflows/anatomical/base.py

@@ -402,7 +402,7 @@ def init_infant_anat_wf(
     elif config.workflow.surface_recon_method == 'mcribs':
         from nipype.interfaces.ants import DenoiseImage
 
-        from .surfaces import init_mcribs_surface_recon_wf
+        from .surfaces import init_mcribs_sphere_reg_wf, init_mcribs_surface_recon_wf
 
         # Denoise raw T2w, since using the template / preproc resulted in intersection errors
         denoise_raw_t2w = pe.Node(
@@ -511,20 +511,24 @@ def init_infant_anat_wf(
     # fmt: on
 
     if cifti_output:
-        from smriprep.workflows.surfaces import init_morph_grayords_wf
+        from nibabies.workflows.anatomical.resampling import (
+            init_anat_fsLR_resampling_wf,
+        )
 
-        morph_grayords_wf = init_morph_grayords_wf(grayord_density=cifti_output)
+        is_mcribs = config.workflow.surface_recon_method == "mcribs"
+        # handles morph_grayords_wf
+        anat_fsLR_resampling_wf = init_anat_fsLR_resampling_wf(cifti_output, mcribs=is_mcribs)
         anat_derivatives_wf.get_node('inputnode').inputs.cifti_density = cifti_output
         # fmt:off
         wf.connect([
-            (surface_recon_wf, morph_grayords_wf, [
+            (surface_recon_wf, anat_fsLR_resampling_wf, [
                 ('outputnode.subject_id', 'inputnode.subject_id'),
-                ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
-            ]),
-            (morph_grayords_wf, anat_derivatives_wf, [
+                ('outputnode.subjects_dir', 'inputnode.subjects_dir')]),
+            (anat_fsLR_resampling_wf, anat_derivatives_wf, [
                 ("outputnode.cifti_morph", "inputnode.cifti_morph"),
-                ("outputnode.cifti_metadata", "inputnode.cifti_metadata"),
-            ]),
+                ("outputnode.cifti_metadata", "inputnode.cifti_metadata")]),
+            (anat_fsLR_resampling_wf, outputnode, [
+                ("outputnode.fsLR_midthickness", "fsLR_midthickness")])
         ])
         # fmt:on
 

nibabies/workflows/anatomical/resampling.py

@@ -47,7 +47,16 @@ def init_anat_fsLR_resampling_wf(
         iterables=[('hemi', ['L', 'R'])],
     )
 
-    outputnode = pe.Node(niu.IdentityInterface(fields=['fsLR_midthickness']), name='outputnode')
+    outputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'fsLR_midthickness',
+                'cifti_morph',
+                'cifti_metadata',
+            ]
+        ),
+        name='outputnode',
+    )
 
     # select white, midthickness and pial surfaces based on hemi
     select_surfaces = pe.Node(CiftiSelect(), name='select_surfaces')
@@ -96,50 +105,33 @@ def init_anat_fsLR_resampling_wf(
     else:
         morph_grayords_wf = init_morph_grayords_wf(grayord_density)
 
-    workflow.connect(
-        [
-            (
-                inputnode,
-                select_surfaces,
-                [("surfaces", "surfaces"), ("sphere_reg_fsLR", "spherical_registrations")],
-            ),
-            (itersource, select_surfaces, [("hemi", "hemi")]),
-            # Downsample midthickness to fsLR density
-            (
-                select_surfaces,
-                downsampled_midthickness,
-                [
-                    ("midthickness", "surface_in"),
-                    ("sphere_reg", "current_sphere"),
-                    ("template_sphere", "new_sphere"),
-                ],
-            ),
-            (downsampled_midthickness, joinnode, [("surface_out", "fsLR_midthickness")]),
-            (joinnode, outputnode, [("surface_out", "fsLR_midthickness")]),
-            # resample surfaces
-            (
-                inputnode,
-                morph_grayords_wf,
-                [
-                    ("subject_id", "inputnode.subject_id"),
-                    ("subjects_dir", "inputnode.subjects_dir"),
-                ],
-            ),
-            (
-                morph_grayords_wf,
-                outputnode,
-                [
-                    ("outputnode.cifti_morph", "cifti_morph"),
-                    ("outputnode.cifti_metadata", "cifti_metadata"),
-                ],
-            ),
-        ]
-    )
+    # fmt:off
+    workflow.connect([
+        (inputnode, select_surfaces, [
+            ("surfaces", "surfaces"),
+            ("sphere_reg_fsLR", "spherical_registrations")]),
+        (itersource, select_surfaces, [("hemi", "hemi")]),
+        # Downsample midthickness to fsLR density
+        (select_surfaces, downsampled_midthickness, [
+            ("midthickness", "surface_in"),
+            ("sphere_reg", "current_sphere"),
+            ("template_sphere", "new_sphere")]),
+        (downsampled_midthickness, joinnode, [("surface_out", "fsLR_midthickness")]),
+        (joinnode, outputnode, [("surface_out", "fsLR_midthickness")]),
+        # resample morphometrics to fsLR 32k
+        (inputnode, morph_grayords_wf, [
+            ("subject_id", "inputnode.subject_id"),
+            ("subjects_dir", "inputnode.subjects_dir")]),
+        (morph_grayords_wf, outputnode, [
+            ("outputnode.cifti_morph", "cifti_morph"),
+            ("outputnode.cifti_metadata", "cifti_metadata")]),
+    ])
+    # fmt:on
     return workflow
 
 
 def init_mcribs_morph_grayords_wf(
-    grayord_density: ty.Literal['91k', '170k'],
+    grayord_density: ty.Literal['91k'],  # Only 91k supported ATM
     name: str = "morph_grayords_wf",
 ):
     """
@@ -180,21 +172,17 @@ def init_mcribs_morph_grayords_wf(
         Paths to JSON files containing metadata corresponding to ``cifti_morph``
 
     """
-    import templateflow.api as tf
-    from nipype.interfaces.io import FreeSurferSource
     from nipype.interfaces.workbench import MetricResample
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from smriprep.interfaces.cifti import GenerateDScalar
 
     workflow = Workflow(name=name)
     workflow.__desc__ = f"""\
 *Grayordinate* "dscalar" files [@hcppipelines] containing {grayord_density} samples were
-also generated using the highest-resolution ``fsaverage`` as an intermediate standardized
+also generated using `M-CRIB-S` as an intermediate standardized
 surface space.
 """
 
-    fslr_density = "32k" if grayord_density == "91k" else "59k"
-
     inputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
@@ -213,17 +201,12 @@ def init_mcribs_morph_grayords_wf(
         name="outputnode",
     )
 
-    get_surfaces = pe.Node(FreeSurferSource(), name="get_surfaces")
-
     surfmorph_list = pe.Node(
         niu.Merge(3, ravel_inputs=True),
         name="surfmorph_list",
         run_without_submitting=True,
     )
 
-    # TODO: Extract L/R midthickness from surfaces
-    # TODO: Coerce morphometrics into curv-L, curv-R, sulc-L, sulc-R, thickness-L, thickness-R
-
     # Setup Workbench command. LR ordering for hemi can be assumed, as it is imposed
     # by the iterfield of the MapNode in the surface sampling workflow above.
     resample = pe.MapNode(
@@ -244,7 +227,6 @@ def init_mcribs_morph_grayords_wf(
         str(atlases / 'mcribs' / 'lh.sphere.reg.dHCP42.surf.gii'),
         str(atlases / 'mcribs' / 'rh.sphere.reg.dHCP42.surf.gii'),
     ] * 3
-    # current area: FreeSurfer (M-CRIB-S) midthickness
     resample.inputs.new_sphere = [
         str(atlases / 'dHCP' / 'dHCP.week42.L.sphere.surf.gii'),
         str(atlases / 'dHCP' / 'dHCP.week42.R.sphere.surf.gii'),
@@ -267,23 +249,23 @@ def init_mcribs_morph_grayords_wf(
 
     # fmt: off
     workflow.connect([
-        (inputnode, get_surfaces, [
-            ('subject_id', 'subject_id'),
-            ('subjects_dir', 'subjects_dir'),
+        (inputnode, resample, [
+            ("fsLR_midthickness", "new_area"),
+            (('surfaces', _get_surf, "midthickness", 3), "current_area")]),
+        (inputnode, surfmorph_list, [
+            (('morphometrics', _get_surf, "curv"), "in1"),
+            (('morphometrics', _get_surf, "sulc"), "in2"),
+            (('morphometrics', _get_surf, "thickness"), "in3"),
         ]),
-        (get_surfaces, surfmorph_list, [
-            (('curv', _sorted_by_basename), 'in1'),
-            (('sulc', _sorted_by_basename), 'in2'),
-            (('thickness', _sorted_by_basename), 'in3'),
-        ]),
-        # (surfmorph_list, surf2surf, [('out', 'source_file')]),
-        # (surf2surf, resample, [('out_file', 'in_file')]),
-        (inputnode, resample, [("fsLR_midthickness", "new_area")]),
         (resample, gen_cifti, [
             (("out_file", _collate), "scalar_surfs")]),
         (gen_cifti, outputnode, [("out_file", "cifti_morph"),
                                  ("out_metadata", "cifti_metadata")]),
     ])
     # fmt: on
-
     return workflow
+
+
+def _get_surf(surfaces, name, mult=1):
+    "Select a specific surface by name, and optionally multiple it."
+    return [surf for surf in _sorted_by_basename(surfaces) if name in surf] * mult