ENH: Add M-CRIB-S -> dHCP spherical registration workflow

Author: mgxd

nibabies/workflows/anatomical/base.py

@@ -162,6 +162,9 @@ def init_infant_anat_wf(
                 "anat_aparc",
                 "anat_ribbon",
                 "template",
+                # registration sphere space is dependent on surface recon method
+                "sphere_reg",
+                "sphere_reg_fsLR",
             ]
         ),
         name="outputnode",
@@ -413,6 +416,9 @@ def init_infant_anat_wf(
             mcribs_dir=str(config.execution.mcribs_dir),  # Needed to preserve runs
         )
 
+        # M-CRIB-S to dHCP42week (32k)
+        sphere_reg_wf = init_mcribs_sphere_reg_wf()
+
         # Transformed gives
         if precomp_aseg:
             surface_recon_wf.inputs.inputnode.ants_segs = precomp_aseg
@@ -421,11 +427,18 @@ def init_infant_anat_wf(
         # fmt:off
         wf.connect([
             (inputnode, denoise_raw_t2w, [('t2w', 'input_image')]),
-            (denoise_raw_t2w, surface_recon_wf, [('output_image', 'inputnode.t2w')])
+            (denoise_raw_t2w, surface_recon_wf, [('output_image', 'inputnode.t2w')]),
         ])
         # fmt:on
+    else:
+        raise NotImplementedError
 
     if config.workflow.surface_recon_method in ('freesurfer', 'infantfs'):
+        from smriprep.workflows.surfaces import init_sphere_reg_wf
+
+        # fsaverage to fsLR
+        sphere_reg_wf = init_sphere_reg_wf()
+
         # fmt:off
         wf.connect([
             (t2w_preproc_wf, surface_recon_wf, [
@@ -472,6 +485,10 @@ def init_infant_anat_wf(
         (anat_ribbon_wf, anat_derivatives_wf, [
             ("outputnode.anat_ribbon", "inputnode.anat_ribbon"),
         ]),
+        (surface_recon_wf, sphere_reg_wf, [
+            ('outputnode.subject_id', 'inputnode.subject_id'),
+            ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
+        ]),
         (surface_recon_wf, anat_reports_wf, [
             ("outputnode.subject_id", "inputnode.subject_id"),
             ("outputnode.subjects_dir", "inputnode.subjects_dir"),
@@ -484,6 +501,12 @@ def init_infant_anat_wf(
             ("outputnode.surfaces", "inputnode.surfaces"),
             ("outputnode.morphometrics", "inputnode.morphometrics"),
         ]),
+        (sphere_reg_wf, outputnode, [
+            ('outputnode.sphere_reg', 'sphere_reg'),
+            ('outputnode.sphere_reg_fsLR', 'sphere_reg_fsLR')]),
+        (sphere_reg_wf, anat_derivatives_wf, [
+            ('outputnode.sphere_reg', 'inputnode.sphere_reg'),
+            ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR')]),
     ])
     # fmt: on
 

nibabies/workflows/anatomical/surfaces.py

@@ -12,6 +12,7 @@
 from smriprep.workflows.surfaces import init_gifti_surface_wf
 
 from ...config import DEFAULT_MEMORY_MIN_GB
+from ...data import load_resource
 
 SURFACE_INPUTS = [
     "subjects_dir",
@@ -201,6 +202,81 @@ def init_mcribs_surface_recon_wf(
     return wf
 
 
+def init_mcribs_sphere_reg_wf(*, name="mcribs_sphere_reg_wf"):
+    """
+    Generate GIFTI registration sphere files from MCRIBS template to dHCP42 (32k).
+
+    TODO: Clarify any distinction with fsLR
+    """
+    from smriprep.interfaces.surf import FixGiftiMetadata
+    from smriprep.interfaces.workbench import SurfaceSphereProjectUnproject
+
+    workflow = LiterateWorkflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(["subjects_dir", "subject_id"]),
+        name="inputnode",
+    )
+    outputnode = pe.Node(
+        niu.IdentityInterface(["sphere_reg", "sphere_reg_fsLR"]),
+        name="outputnode",
+    )
+
+    get_spheres = pe.Node(
+        niu.Function(function=_get_dhcp_spheres),
+        name='get_spheres',
+        run_without_submitting=True,
+    )
+
+    # Via FreeSurfer2CaretConvertAndRegisterNonlinear.sh#L270-L273
+    #
+    # See https://github.com/DCAN-Labs/DCAN-HCP/tree/9291324
+    sphere_gii = pe.MapNode(
+        fs.MRIsConvert(out_datatype="gii"), iterfield="in_file", name="sphere_gii"
+    )
+
+    fix_meta = pe.MapNode(FixGiftiMetadata(), iterfield="in_file", name="fix_meta")
+
+    # load template files
+    atlases = load_resource('atlases')
+
+    # SurfaceSphereProjectUnProject
+    # project to 41k dHCP atlas sphere
+    #   - sphere-in: Individual native sphere in surf directory registered to 41k atlas sphere
+    #   - sphere-to: the 41k atlas sphere, in the fsaverage directory
+    #   - sphere-unproject-from: 41k atlas sphere registered to dHCP 42wk sphere, in the fsaverage directory
+    #   - sphere-out: lh.sphere.reg2.dHCP42.native.surf.gii
+    project_unproject = pe.MapNode(
+        SurfaceSphereProjectUnproject(),
+        iterfield=["sphere_in", "sphere_project_to", "sphere_unproject_from"],
+        name="project_unproject",
+    )
+    project_unproject.inputs.sphere_project_to = [
+        atlases / 'mcribs' / 'lh.sphere.reg2.surf.gii',
+        atlases / 'mcribs' / 'rh.sphere.reg2.surf.gii',
+    ]
+    project_unproject.inputs.sphere_unproject_from = [
+        atlases / 'mcribs' / 'lh.sphere.reg.dHCP42.surf.gii',
+        atlases / 'mcribs' / 'rh.sphere.reg.dHCP42.surf.gii',
+    ]
+
+    # fmt:off
+    workflow.connect([
+        (inputnode, get_spheres, [
+            ('subjects_dir', 'subjects_dir'),
+            ('subject_id', 'subject_id'),
+        ]),
+        (get_spheres, sphere_gii, [(('out', _sorted_by_basename), 'in_file')]),
+        (sphere_gii, fix_meta, [('converted', 'in_file')]),
+        (fix_meta, project_unproject, [('out_file', 'sphere_in')]),
+        (sphere_gii, outputnode, [('converted', 'sphere_reg')]),
+        (project_unproject, outputnode, [('sphere_out', 'sphere_reg_fsLR')]),
+    ])
+    # fmt:on
+
+    return workflow
+
+
 def init_infantfs_surface_recon_wf(
     *, age_months, use_aseg=False, name="infantfs_surface_recon_wf"
 ):
@@ -483,3 +559,15 @@ def _sorted_by_basename(inlist):
     from os.path import basename
 
     return sorted(inlist, key=lambda x: str(basename(x)))
+
+
+def _get_dhcp_spheres(subject_id: str, subjects_dir: str) -> list:
+    from pathlib import Path
+
+    out = []
+    for hemi in 'lr':
+        sphere = Path(subjects_dir) / subject_id / 'surf' / f'{hemi}h.sphere.reg2'
+        if not sphere.exists():
+            raise OSError("MCRIBS spherical registration not found.")
+        out.append(str(sphere))
+    return out