Merge pull request #388 from effigies/rf/template-resampling-workflows

RF: Split template and fsLR resampling and sinking into isolated workflows

Author: effigies

smriprep/interfaces/gifti.py

@@ -5,58 +5,76 @@
 from nipype.interfaces.base import File, SimpleInterface, TraitedSpec, isdefined, traits
 
 
-class InvertShapeInputSpec(TraitedSpec):
+class MetricMathInputSpec(TraitedSpec):
     subject_id = traits.Str(desc='subject ID')
     hemisphere = traits.Enum(
         "L",
         "R",
         mandatory=True,
         desc='hemisphere',
     )
-    shape = traits.Str(desc='name of shape to invert')
-    shape_file = File(exists=True, mandatory=True, desc='input GIFTI file')
+    metric = traits.Str(desc='name of metric to invert')
+    metric_file = File(exists=True, mandatory=True, desc='input GIFTI file')
+    operation = traits.Enum(
+        "invert",
+        "abs",
+        "bin",
+        mandatory=True,
+        desc='operation to perform',
+    )
 
 
-class InvertShapeOutputSpec(TraitedSpec):
-    shape_file = File(desc='output GIFTI file')
+class MetricMathOutputSpec(TraitedSpec):
+    metric_file = File(desc='output GIFTI file')
 
 
-class InvertShape(SimpleInterface):
-    """Prepare GIFTI shape file for use in MSMSulc
+class MetricMath(SimpleInterface):
+    """Prepare GIFTI metric file for use in MSMSulc
 
     This interface mirrors the action of the following portion
     of FreeSurfer2CaretConvertAndRegisterNonlinear.sh::
 
-        wb_command -set-structure ${shape_file} CORTEX_[LEFT|RIGHT]
-        wb_command -metric-math "var * -1" ${shape_file} -var var ${shape_file}
-        wb_command -set-map-names ${shape_file} -map 1 ${subject}_[L|R]_${shape}
+        wb_command -set-structure ${metric_file} CORTEX_[LEFT|RIGHT]
+        wb_command -metric-math "var * -1" ${metric_file} -var var ${metric_file}
+        wb_command -set-map-names ${metric_file} -map 1 ${subject}_[L|R]_${metric}
+        # If abs:
+        wb_command -metric-math "abs(var)" ${metric_file} -var var ${metric_file}
 
     We do not add palette information to the output file.
     """
 
-    input_spec = InvertShapeInputSpec
-    output_spec = InvertShapeOutputSpec
+    input_spec = MetricMathInputSpec
+    output_spec = MetricMathOutputSpec
 
     def _run_interface(self, runtime):
-        subject, hemi, shape = self.inputs.subject_id, self.inputs.hemisphere, self.inputs.shape
+        subject, hemi, metric = self.inputs.subject_id, self.inputs.hemisphere, self.inputs.metric
         if not isdefined(subject):
             subject = 'sub-XYZ'
 
-        img = nb.GiftiImage.from_filename(self.inputs.shape_file)
+        img = nb.GiftiImage.from_filename(self.inputs.metric_file)
         # wb_command -set-structure
         img.meta["AnatomicalStructurePrimary"] = {'L': 'CortexLeft', 'R': 'CortexRight'}[hemi]
         darray = img.darrays[0]
         # wb_command -set-map-names
         meta = darray.meta
-        meta['Name'] = f"{subject}_{hemi}_{shape}"
-
-        # wb_command -metric-math "var * -1"
-        inv = -darray.data
+        meta['Name'] = f"{subject}_{hemi}_{metric}"
+
+        datatype = darray.datatype
+        if self.inputs.operation == "abs":
+            # wb_command -metric-math "abs(var)"
+            data = abs(darray.data)
+        elif self.inputs.operation == "invert":
+            # wb_command -metric-math "var * -1"
+            data = -darray.data
+        elif self.inputs.operation == "bin":
+            # wb_command -metric-math "var > 0"
+            data = darray.data > 0
+            datatype = 'uint8'
 
         darray = nb.gifti.GiftiDataArray(
-            inv,
+            data,
             intent=darray.intent,
-            datatype=darray.datatype,
+            datatype=datatype,
             encoding=darray.encoding,
             endian=darray.endian,
             coordsys=darray.coordsys,
@@ -65,7 +83,7 @@ def _run_interface(self, runtime):
         )
         img.darrays[0] = darray
 
-        out_filename = os.path.join(runtime.cwd, f"{subject}.{hemi}.{shape}.native.shape.gii")
+        out_filename = os.path.join(runtime.cwd, f"{subject}.{hemi}.{metric}.native.shape.gii")
         img.to_filename(out_filename)
-        self._results["shape_file"] = out_filename
+        self._results["metric_file"] = out_filename
         return runtime

smriprep/workflows/anatomical.py

@@ -55,6 +55,8 @@
 from .fit.registration import init_register_template_wf
 from .outputs import (
     init_anat_reports_wf,
+    init_ds_anat_volumes_wf,
+    init_ds_grayord_metrics_wf,
     init_ds_surface_metrics_wf,
     init_ds_surfaces_wf,
     init_ds_template_wf,
@@ -64,16 +66,20 @@
     init_ds_template_registration_wf,
     init_ds_fs_registration_wf,
     init_anat_second_derivatives_wf,
+    init_template_iterator_wf,
 )
 from .surfaces import (
     init_anat_ribbon_wf,
     init_fsLR_reg_wf,
     init_gifti_morphometrics_wf,
+    init_hcp_morphometrics_wf,
     init_gifti_surfaces_wf,
+    init_morph_grayords_wf,
     init_msm_sulc_wf,
     init_surface_derivatives_wf,
     init_surface_recon_wf,
     init_refinement_wf,
+    init_resample_midthickness_wf,
 )
 
 LOGGER = logging.getLogger("nipype.workflow")
@@ -113,6 +119,8 @@ def init_anat_preproc_wf(
 
             from niworkflows.utils.spaces import SpatialReferences, Reference
             from smriprep.workflows.anatomical import init_anat_preproc_wf
+            spaces = SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5'])
+            spaces.checkpoint()
             wf = init_anat_preproc_wf(
                 bids_root='.',
                 output_dir='.',
@@ -124,7 +132,7 @@ def init_anat_preproc_wf(
                 t2w=[],
                 skull_strip_mode='force',
                 skull_strip_template=Reference('OASIS30ANTs'),
-                spaces=SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']),
+                spaces=spaces,
                 precomputed={},
                 omp_nthreads=1,
             )
@@ -261,14 +269,13 @@ def init_anat_preproc_wf(
         omp_nthreads=omp_nthreads,
         skull_strip_fixed_seed=skull_strip_fixed_seed,
     )
-    anat_second_derivatives_wf = init_anat_second_derivatives_wf(
+    template_iterator_wf = init_template_iterator_wf(spaces=spaces)
+    ds_std_volumes_wf = init_ds_anat_volumes_wf(
         bids_root=bids_root,
-        freesurfer=freesurfer,
         output_dir=output_dir,
-        spaces=spaces,
-        cifti_output=cifti_output,
+        name="ds_std_volumes_wf",
     )
-    # fmt:off
+
     workflow.connect([
         (inputnode, anat_fit_wf, [
             ("t1w", "inputnode.t1w"),
@@ -293,19 +300,32 @@ def init_anat_preproc_wf(
             (f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}", "sphere_reg_fsLR"),
             ("outputnode.anat_ribbon", "anat_ribbon"),
         ]),
-        (anat_fit_wf, anat_second_derivatives_wf, [
-            ('outputnode.template', 'inputnode.template'),
+        (anat_fit_wf, template_iterator_wf, [
+            ("outputnode.template", "inputnode.template"),
+            ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"),
+        ]),
+        (anat_fit_wf, ds_std_volumes_wf, [
             ('outputnode.t1w_valid_list', 'inputnode.source_files'),
             ("outputnode.t1w_preproc", "inputnode.t1w_preproc"),
             ("outputnode.t1w_mask", "inputnode.t1w_mask"),
             ("outputnode.t1w_dseg", "inputnode.t1w_dseg"),
             ("outputnode.t1w_tpms", "inputnode.t1w_tpms"),
-            ('outputnode.anat2std_xfm', 'inputnode.anat2std_xfm'),
         ]),
-    ])
-    # fmt:on
-    if freesurfer:
+        (template_iterator_wf, ds_std_volumes_wf, [
+            ("outputnode.std_t1w", "inputnode.ref_file"),
+            ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"),
+            ("outputnode.space", "inputnode.space"),
+            ("outputnode.cohort", "inputnode.cohort"),
+            ("outputnode.resolution", "inputnode.resolution"),
+        ]),
+    ])  # fmt:skip
 
+    if freesurfer:
+        anat_second_derivatives_wf = init_anat_second_derivatives_wf(
+            bids_root=bids_root,
+            output_dir=output_dir,
+            cifti_output=cifti_output,
+        )
         surface_derivatives_wf = init_surface_derivatives_wf(
             cifti_output=cifti_output,
         )
@@ -316,7 +336,6 @@ def init_anat_preproc_wf(
             bids_root=bids_root, output_dir=output_dir, metrics=["curv"], name="ds_curv_wf"
         )
 
-        # fmt:off
         workflow.connect([
             (anat_fit_wf, surface_derivatives_wf, [
                 ('outputnode.t1w_preproc', 'inputnode.reference'),
@@ -336,18 +355,78 @@ def init_anat_preproc_wf(
             (surface_derivatives_wf, ds_curv_wf, [
                 ('outputnode.curv', 'inputnode.curv'),
             ]),
+            (anat_fit_wf, anat_second_derivatives_wf, [
+                ('outputnode.t1w_valid_list', 'inputnode.source_files'),
+            ]),
             (surface_derivatives_wf, anat_second_derivatives_wf, [
                 ('outputnode.out_aseg', 'inputnode.t1w_fs_aseg'),
                 ('outputnode.out_aparc', 'inputnode.t1w_fs_aparc'),
-                ('outputnode.cifti_morph', 'inputnode.cifti_morph'),
-                ('outputnode.cifti_metadata', 'inputnode.cifti_metadata'),
             ]),
             (surface_derivatives_wf, outputnode, [
                 ('outputnode.out_aseg', 't1w_aseg'),
                 ('outputnode.out_aparc', 't1w_aparc'),
             ]),
-        ])
-        # fmt:on
+        ])  # fmt:skip
+
+        if cifti_output:
+            hcp_morphometrics_wf = init_hcp_morphometrics_wf(omp_nthreads=omp_nthreads)
+            resample_midthickness_wf = init_resample_midthickness_wf(grayord_density=cifti_output)
+            morph_grayords_wf = init_morph_grayords_wf(
+                grayord_density=cifti_output, omp_nthreads=omp_nthreads
+            )
+
+            ds_grayord_metrics_wf = init_ds_grayord_metrics_wf(
+                bids_root=bids_root,
+                output_dir=output_dir,
+                metrics=['curv', 'thickness', 'sulc'],
+                cifti_output=cifti_output,
+            )
+
+            workflow.connect([
+                (anat_fit_wf, hcp_morphometrics_wf, [
+                    ('outputnode.subject_id', 'inputnode.subject_id'),
+                    ('outputnode.sulc', 'inputnode.sulc'),
+                    ('outputnode.thickness', 'inputnode.thickness'),
+                    ('outputnode.midthickness', 'inputnode.midthickness'),
+                ]),
+                (surface_derivatives_wf, hcp_morphometrics_wf, [
+                    ('outputnode.curv', 'inputnode.curv'),
+                ]),
+                (anat_fit_wf, resample_midthickness_wf, [
+                    ('outputnode.midthickness', 'inputnode.midthickness'),
+                    (
+                        f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}",
+                        "inputnode.sphere_reg_fsLR",
+                    ),
+                ]),
+                (anat_fit_wf, morph_grayords_wf, [
+                    ('outputnode.midthickness', 'inputnode.midthickness'),
+                    (
+                        f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}",
+                        "inputnode.sphere_reg_fsLR",
+                    ),
+                ]),
+                (hcp_morphometrics_wf, morph_grayords_wf, [
+                    ('outputnode.curv', 'inputnode.curv'),
+                    ('outputnode.sulc', 'inputnode.sulc'),
+                    ('outputnode.thickness', 'inputnode.thickness'),
+                    ('outputnode.roi', 'inputnode.roi'),
+                ]),
+                (resample_midthickness_wf, morph_grayords_wf, [
+                    ('outputnode.midthickness_fsLR', 'inputnode.midthickness_fsLR'),
+                ]),
+                (anat_fit_wf, ds_grayord_metrics_wf, [
+                    ('outputnode.t1w_valid_list', 'inputnode.source_files'),
+                ]),
+                (morph_grayords_wf, ds_grayord_metrics_wf, [
+                    ('outputnode.curv_fsLR', 'inputnode.curv'),
+                    ('outputnode.curv_metadata', 'inputnode.curv_metadata'),
+                    ('outputnode.thickness_fsLR', 'inputnode.thickness'),
+                    ('outputnode.thickness_metadata', 'inputnode.thickness_metadata'),
+                    ('outputnode.sulc_fsLR', 'inputnode.sulc'),
+                    ('outputnode.sulc_metadata', 'inputnode.sulc_metadata'),
+                ]),
+            ])  # fmt:skip
 
     return workflow
 
@@ -1264,7 +1343,7 @@ def init_anat_fit_wf(
         # fmt:on
     elif msm_sulc:
         LOGGER.info("ANAT Stage 10: Found pre-computed MSM-Sulc registration sphere")
-        fsLR_buffer.inputs.sphere_reg_msm = sorted(precomputed["sphere_reg_msm"])
+        msm_buffer.inputs.sphere_reg_msm = sorted(precomputed["sphere_reg_msm"])
     else:
         LOGGER.info("ANAT Stage 10: MSM-Sulc disabled")
 

smriprep/workflows/base.py

@@ -81,6 +81,8 @@ def init_smriprep_wf(
             os.environ['FREESURFER_HOME'] = os.getcwd()
             from smriprep.workflows.base import init_smriprep_wf
             from niworkflows.utils.spaces import SpatialReferences, Reference
+            spaces = SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5'])
+            spaces.checkpoint()
             wf = init_smriprep_wf(
                 sloppy=False,
                 debug=False,
@@ -98,7 +100,7 @@ def init_smriprep_wf(
                 skull_strip_fixed_seed=False,
                 skull_strip_mode='force',
                 skull_strip_template=Reference('OASIS30ANTs'),
-                spaces=SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']),
+                spaces=spaces,
                 subject_list=['smripreptest'],
                 work_dir='.',
                 bids_filters=None,
@@ -250,6 +252,8 @@ def init_single_subject_wf(
             from niworkflows.utils.spaces import SpatialReferences, Reference
             from smriprep.workflows.base import init_single_subject_wf
             BIDSLayout = namedtuple('BIDSLayout', ['root'])
+            spaces = SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5'])
+            spaces.checkpoint()
             wf = init_single_subject_wf(
                 sloppy=False,
                 debug=False,
@@ -266,7 +270,7 @@ def init_single_subject_wf(
                 skull_strip_fixed_seed=False,
                 skull_strip_mode='force',
                 skull_strip_template=Reference('OASIS30ANTs'),
-                spaces=SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']),
+                spaces=spaces,
                 subject_id='test',
                 bids_filters=None,
                 cifti_output=None,