FIX: Goodvoxel projection

.circleci/bcp_anat_outputs.txt

@@ -23,16 +23,20 @@ sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_from-MNIInfant+1_to-T1w_mode-image_xf
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_from-T1w_to-fsnative_mode-image_xfm.txt
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_from-T1w_to-MNIInfant+1_mode-image_xfm.h5
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_curv.shape.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_inflated.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_midthickness.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_pial.surf.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_space-fsLR_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_sulc.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_thickness.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_white.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_curv.shape.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_inflated.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_midthickness.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_pial.surf.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_space-fsLR_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_sulc.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_thickness.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_white.surf.gii

.circleci/bcp_full_outputs.txt

@@ -23,16 +23,20 @@ sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_from-MNIInfant+1_to-T1w_mode-image_xf
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_from-T1w_to-fsnative_mode-image_xfm.txt
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_from-T1w_to-MNIInfant+1_mode-image_xfm.h5
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_curv.shape.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_inflated.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_midthickness.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_pial.surf.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_space-fsLR_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_sulc.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_thickness.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-L_white.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_curv.shape.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_inflated.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_midthickness.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_pial.surf.gii
+sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_space-fsLR_desc-reg_sphere.surf.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_sulc.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_thickness.shape.gii
 sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_hemi-R_white.surf.gii

nibabies/data/__init__.py

@@ -0,0 +1,25 @@
+import atexit
+from contextlib import ExitStack
+from pathlib import Path
+
+try:
+    from functools import cache
+except ImportError:  # PY38
+    from functools import lru_cache as cache
+
+try:  # Prefer backport to leave consistency to dependency spec
+    from importlib_resources import as_file, files
+except ImportError:
+    from importlib.resources import as_file, files
+
+__all__ = ["load_resource"]
+
+exit_stack = ExitStack()
+atexit.register(exit_stack.close)
+
+path = files(__package__)
+
+
+@cache
+def load_resource(fname: str) -> Path:
+    return exit_stack.enter_context(as_file(path.joinpath(fname)))

nibabies/interfaces/utils.py

@@ -0,0 +1,64 @@
+import os
+import re
+
+from nipype.interfaces.base import (
+    BaseInterfaceInputSpec,
+    File,
+    InputMultiObject,
+    OutputMultiObject,
+    SimpleInterface,
+    TraitedSpec,
+    traits,
+)
+
+
+class CiftiSelectInputSpec(BaseInterfaceInputSpec):
+    hemi = traits.Enum("L", "R", desc="Hemisphere")
+    surfaces = InputMultiObject(File(exists=True), desc="Surfaces")
+    morphometrics = InputMultiObject(File(exists=True), desc="Surface morphometrics")
+    spherical_registrations = InputMultiObject(
+        File(exists=True), desc="Spherical registration to fsLR"
+    )
+    template_spheres = InputMultiObject(File(exists=True), desc="fsLR sphere")
+    template_surfaces = InputMultiObject(File(exists=True), desc="fsLR midthickness")
+    template_rois = InputMultiObject(File(exists=True), desc="fsLR ROIs")
+
+
+class CiftiSelectOutputSpec(TraitedSpec):
+    white = OutputMultiObject(File, desc="white surface")
+    pial = OutputMultiObject(File, desc="pial surface")
+    midthickness = OutputMultiObject(File, desc="midthickness surface")
+    thickness = OutputMultiObject(File, desc="thickness surface")
+    sphere_reg = OutputMultiObject(File, desc="fsLR spherical regisration")
+    template_sphere = OutputMultiObject(File, desc="fsLR sphere")
+    template_surface = OutputMultiObject(File, desc="fsLR surface (midthickness)")
+    template_roi = OutputMultiObject(File, desc="fsLR ROIs")
+
+
+class CiftiSelect(SimpleInterface):
+    input_spec = CiftiSelectInputSpec
+    output_spec = CiftiSelectOutputSpec
+
+    def _run_interface(self, runtime):
+        idx = 0 if self.inputs.hemi == "L" else 1
+        all_surfaces = (self.inputs.surfaces or []) + (self.inputs.morphometrics or [])
+        container = {
+            'white': [],
+            'pial': [],
+            'midthickness': [],
+            'thickness': [],
+            'sphere_reg': self.inputs.spherical_registrations or [],
+            'template_sphere': self.inputs.template_spheres or [],
+            'template_surface': self.inputs.template_surfaces or [],
+            'template_roi': self.inputs.template_rois or [],
+        }
+        find_name = re.compile(r'(?:^|[^d])(?P<name>white|pial|midthickness|thickness)')
+        for surface in all_surfaces:
+            match = find_name.search(os.path.basename(surface))
+            if match:
+                container[match.group('name')].append(surface)
+
+        for name, vals in container.items():
+            if vals:
+                self._results[name] = sorted(vals, key=os.path.basename)[idx]
+        return runtime

nibabies/interfaces/workbench.py

@@ -1570,3 +1570,113 @@ class CreateSignedDistanceVolume(WBCommand):
     input_spec = CreateSignedDistanceVolumeInputSpec
     output_spec = CreateSignedDistanceVolumeOutputSpec
     _cmd = "wb_command -create-signed-distance-volume"
+
+
+class SurfaceAverageInputSpec(CommandLineInputSpec):
+    out_file = File(
+        name_template="averaged.surf.gii",
+        position=0,
+        desc="output file",
+    )
+    surfaces = InputMultiObject(
+        traits.Either(
+            File(exists=True),
+            traits.Tuple(File(exists=True), traits.Float),
+        ),
+        argstr="%s",
+        position=3,
+        desc="Surface, or surface and weighted average tuple, to include in the average",
+    )
+
+
+class SurfaceAverageOutputSpec(TraitedSpec):
+    out_file = File(desc="The output averaged surface")
+    # stddev_metric = File(desc="The output metric for 3D sample standard deviation")
+    # uncert_metric = File(desc="The output metric for uncertainty")
+
+
+class SurfaceAverage(WBCommand):
+    """
+    AVERAGE SURFACE FILES TOGETHER
+    wb_command -surface-average
+      <surface-out> - output - the output averaged surface
+
+      [-stddev] - compute 3D sample standard deviation
+         <stddev-metric-out> - output - the output metric for 3D sample
+            standard deviation
+
+      [-uncertainty] - compute caret5 'uncertainty'
+         <uncert-metric-out> - output - the output metric for uncertainty
+
+      [-surf] - repeatable - specify a surface to include in the average
+         <surface> - a surface file to average
+
+         [-weight] - specify a weighted average
+            <weight> - the weight to use (default 1)
+
+      The 3D sample standard deviation is computed as
+      'sqrt(sum(squaredlength(xyz - mean(xyz)))/(n - 1))'.
+
+      Uncertainty is a legacy measure used in caret5, and is computed as
+      'sum(length(xyz - mean(xyz)))/n'.
+
+      When weights are used, the 3D sample standard deviation treats them as
+      reliability weights.
+    """
+
+    input_spec = SurfaceAverageInputSpec
+    output_spec = SurfaceAverageOutputSpec
+    _cmd = "wb_command -surface-average"
+
+    def _format_arg(self, name, trait_spec, value):
+        if name == 'surfaces':
+            cmd = []
+            for val in value:
+                if len(val) == 2:
+                    cmd.append(f"{val[0]} -weight {val[-1]}")
+                else:
+                    cmd.append(val)
+            return '-surf ' + ' -surf '.join(cmd)
+        return super()._format_arg(name, trait_spec, value)
+
+    def _list_output(self):
+        outputs = self.output_spec().get()
+        outputs["out_file"] = os.path.abspath(self.inputs.out_file)
+        return outputs
+
+
+class SurfaceVertexAreasInputSpec(CommandLineInputSpec):
+    in_file = File(
+        exists=True,
+        mandatory=True,
+        position=0,
+        argstr="%s",
+        desc="Input surface",
+    )
+    out_file = File(
+        name_template="%s.shape.gii",
+        name_source="in_file",
+        position=1,
+        argstr="%s",
+        desc="Output vertex areas",
+    )
+
+
+class SurfaceVertexAreasOutputSpec(TraitedSpec):
+    out_file = File(desc="Output vertex areas")
+
+
+class SurfaceVertexAreas(WBCommand):
+    """
+    MEASURE SURFACE AREA EACH VERTEX IS RESPONSIBLE FOR
+    wb_command -surface-vertex-areas
+      <surface> - the surface to measure
+      <metric> - output - the output metric
+
+      Each vertex gets one third of the area of each triangle it is a part of.
+      Units are mm^2.
+    """
+
+    input_spec = SurfaceVertexAreasInputSpec
+    output_spec = SurfaceVertexAreasOutputSpec
+    _cmd = "wb_command -surface-vertex-areas"

nibabies/tests/test_config.py

@@ -100,7 +100,6 @@ def test_config_spaces():
 
     assert [str(s) for s in spaces.get_standard(full_spec=True)] == [
         'MNIInfant:cohort-1:res-native',  # Default output space
-        'fsaverage:den-164k',
         'MNI152NLin6Asym:res-2',
     ]
 

nibabies/workflows/anatomical/base.py

@@ -162,6 +162,10 @@ def init_infant_anat_wf(
                 "anat_aparc",
                 "anat_ribbon",
                 "template",
+                # registration sphere space is dependent on surface recon method
+                "sphere_reg",
+                "sphere_reg_fsLR",
+                "midthickness_fsLR",
             ]
         ),
         name="outputnode",
@@ -399,7 +403,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(
@@ -413,6 +417,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 +428,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 +486,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,24 +502,39 @@ 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
 
     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, [
+            (sphere_reg_wf, anat_fsLR_resampling_wf, [
+                ('outputnode.sphere_reg', 'inputnode.sphere_reg'),
+                ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR')]),
+            (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.surfaces', 'inputnode.surfaces'),
+                ('outputnode.morphometrics', 'inputnode.morphometrics')]),
+            (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.midthickness_fsLR", "midthickness_fsLR")])
         ])
         # fmt:on