ENH: Calculate T1/T2 ratio

src/smriprep/interfaces/calc.py

@@ -0,0 +1,90 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2024 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Image calculation interfaces."""
+
+import os
+from pathlib import Path
+
+import nibabel as nb
+import numpy as np
+from nipype.interfaces.base import (
+    File,
+    SimpleInterface,
+    TraitedSpec,
+)
+
+
+class T1T2RatioInputSpec(TraitedSpec):
+    t1w_file = File(exists=True, mandatory=True, desc='T1-weighted image')
+    t2w_file = File(exists=True, mandatory=True, desc='T2-weighted image')
+    mask_file = File(exists=True, desc='Brain mask')
+
+
+class T1T2RatioOutputSpec(TraitedSpec):
+    t1t2_file = File(exists=True, desc='T1/T2 ratio image')
+
+
+class T1T2Ratio(SimpleInterface):
+    input_spec = T1T2RatioInputSpec
+    output_spec = T1T2RatioOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['t1t2_file'] = make_t1t2_ratio(
+            self.inputs.t1w_file, self.inputs.t2w_file, self.inputs.mask_file, newpath=runtime.cwd
+        )
+        return runtime
+
+
+def make_t1t2_ratio(
+    t1w_file: str,
+    t2w_file: str,
+    mask_file: str | None = None,
+    newpath: str | None = None,
+) -> str:
+    t1w = nb.load(t1w_file)
+    t2w = nb.load(t2w_file)
+    if mask_file is not None:
+        mask = np.asanyarray(nb.load(mask_file).dataobj) != 0
+    else:
+        mask = np.ones(t1w.shape, dtype=bool)
+
+    t1w_data = t1w.get_fdata(dtype=np.float32)
+    t2w_data = t2w.get_fdata(dtype=np.float32)
+
+    t1t2_data = np.zeros_like(t1w_data)
+
+    ratio = t1w_data[mask] / t2w_data[mask]
+    ratio[~np.isfinite(ratio)] = 0
+    minval = ratio.min()
+    maxval = ratio.max()
+
+    t1t2_data[mask] = (ratio - minval) / (maxval - minval) * 100
+
+    t1t2 = nb.Nifti1Image(t1t2_data, t1w.affine, t1w.header)
+    t1t2.header.set_data_dtype(np.float32)
+
+    t1t2_path = Path(newpath or os.getcwd()) / 't1t2_ratio.nii.gz'
+
+    t1t2.to_filename(t1t2_path)
+
+    return str(t1t2_path)

src/smriprep/interfaces/tests/test_calc.py

@@ -0,0 +1,24 @@
+import nibabel as nb
+from nipype.pipeline import engine as pe
+from templateflow import api as tf
+
+from ..calc import T1T2Ratio
+
+
+def test_T1T2Ratio(tmp_path):
+    t1w = tf.get('MNI152NLin2009cAsym', desc=None, resolution=1, suffix='T1w')
+    t2w = tf.get('MNI152NLin2009cAsym', desc=None, resolution=1, suffix='T2w')
+    mask = tf.get('MNI152NLin2009cAsym', desc='brain', resolution=1, suffix='mask')
+
+    t1t2 = pe.Node(
+        T1T2Ratio(t1w_file=t1w, t2w_file=t2w, mask_file=mask),
+        name='t1t2',
+        base_dir=tmp_path,
+    )
+
+    result = t1t2.run()
+
+    t1t2ratio = nb.load(result.outputs.t1t2_file)
+    assert t1t2ratio.shape == (193, 229, 193)
+    assert t1t2ratio.get_fdata().min() == 0.0
+    assert t1t2ratio.get_fdata().max() == 100.0

src/smriprep/workflows/anatomical.py

@@ -57,6 +57,7 @@
 import smriprep
 
 from ..interfaces import DerivativesDataSink
+from ..interfaces.calc import T1T2Ratio
 from ..utils.misc import apply_lut as _apply_bids_lut
 from ..utils.misc import fs_isRunning as _fs_isRunning
 from .fit.registration import init_register_template_wf
@@ -83,6 +84,7 @@
     init_hcp_morphometrics_wf,
     init_morph_grayords_wf,
     init_msm_sulc_wf,
+    init_myelinmap_fsLR_wf,
     init_refinement_wf,
     init_resample_surfaces_wf,
     init_surface_derivatives_wf,
@@ -459,6 +461,52 @@ def init_anat_preproc_wf(
                 ]),
             ])  # fmt:skip
 
+            if t2w:
+                myelinmap_sources = pe.Node(niu.Merge(2), name='myelinmap_sources')
+                t1t2_ratio = pe.Node(T1T2Ratio(), name='t1t2_ratio')
+                myelinmap_fsLR_wf = init_myelinmap_fsLR_wf(
+                    grayord_density=cifti_output, omp_nthreads=omp_nthreads, mem_gb=1
+                )
+                ds_grayord_myelinmap_wf = init_ds_grayord_metrics_wf(
+                    bids_root=bids_root,
+                    output_dir=output_dir,
+                    metrics=['myelinmap'],
+                    cifti_output=cifti_output,
+                )
+
+                workflow.connect([
+                    (anat_fit_wf, t1t2_ratio, [
+                        ('outputnode.t1w_preproc', 't1w_file'),
+                        ('outputnode.t2w_preproc', 't2w_file'),
+                    ]),
+                    (t1t2_ratio, myelinmap_fsLR_wf, [('out_file', 'inputnode.in_file')]),
+                    (anat_fit_wf, myelinmap_fsLR_wf, [
+                        ('outputnode.midthickness', 'inputnode.midthickness'),
+                        (
+                            f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}",
+                            'inputnode.sphere_reg_fsLR',
+                        ),
+                    ]),
+                    (hcp_morphometrics_wf, myelinmap_fsLR_wf, [
+                        ('outputnode.thickness', 'inputnode.thickness'),
+                        ('outputnode.roi', 'inputnode.cortex_mask'),
+                    ]),
+                    (resample_surfaces_wf, myelinmap_fsLR_wf, [
+                        ('outputnode.midthickness_fsLR', 'inputnode.midthickness_fsLR'),
+                    ]),
+                    (anat_fit_wf, myelinmap_sources, [
+                        ('outputnode.t1w_preproc', 'in1'),
+                        ('outputnode.t2w_preproc', 'in2'),
+                    ]),
+                    (myelinmap_sources, ds_grayord_myelinmap_wf, [
+                        ('out', 'inputnode.source_files'),
+                    ]),
+                    (myelinmap_fsLR_wf, ds_grayord_myelinmap_wf, [
+                        ('outputnode.out_file', 'inputnode.myelinmap'),
+                        ('outputnode.out_metadata', 'inputnode.myelinmap_metadata'),
+                    ]),
+                ])  # fmt:skip
+
     return workflow
 
 

src/smriprep/workflows/surfaces.py

@@ -30,6 +30,7 @@
 
 import typing as ty
 
+from nibabel.processing import fwhm2sigma
 from nipype.interfaces import freesurfer as fs
 from nipype.interfaces import io as nio
 from nipype.interfaces import utility as niu
@@ -54,6 +55,7 @@
 )
 
 import smriprep
+from smriprep.interfaces.cifti import GenerateDScalar
 from smriprep.interfaces.surf import MakeRibbon
 from smriprep.interfaces.workbench import SurfaceResample
 
@@ -1476,19 +1478,14 @@ def init_morph_grayords_wf(
     thickness_fsLR
         HCP-style thickness file in CIFTI-2 format, resampled to fsLR
     """
-    import templateflow.api as tf
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
-    from smriprep.interfaces.cifti import GenerateDScalar
-
     workflow = Workflow(name=name)
     workflow.__desc__ = f"""\
 *Grayordinate* "dscalar" files containing {grayord_density} samples were
 resampled onto fsLR using the Connectome Workbench [@hcppipelines].
 """
 
-    fslr_density = '32k' if grayord_density == '91k' else '59k'
-
     inputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
@@ -1527,18 +1524,289 @@ def init_morph_grayords_wf(
         name='outputnode',
     )
 
-    atlases = smriprep.load_data('atlases')
+    metrics = ['curv', 'sulc', 'thickness']
+    select_surfaces = pe.Node(
+        KeySelect(fields=metrics, keys=['L', 'R']),
+        name='select_surfaces',
+        run_without_submitting=True,
+    )
+
+    workflow.connect([
+        (hemisource, select_surfaces, [('hemi', 'key')]),
+    ]) # fmt:skip
+
+    for metric in metrics:
+        resample_and_mask_wf = init_resample_and_mask_wf(
+            grayord_density=grayord_density,
+            omp_nthreads=omp_nthreads,
+            mem_gb=1,
+            name=f'resample_and_mask_{metric}_wf',
+        )
+        cifti_metric = pe.JoinNode(
+            GenerateDScalar(grayordinates=grayord_density, scalar_name=metric),
+            name=f'cifti_{metric}',
+            joinfield=['scalar_surfs'],
+            joinsource='hemisource',
+        )
+
+        workflow.connect([
+            (inputnode, select_surfaces, [(metric, metric)]),
+            (inputnode, resample_and_mask_wf, [
+                ('midthickness', 'inputnode.midthickness'),
+                ('midthickness_fsLR', 'inputnode.midthickness_fsLR'),
+                ('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
+                ('roi', 'inputnode.cortex_mask'),
+            ]),
+            (hemisource, resample_and_mask_wf, [('hemi', 'inputnode.hemi')]),
+            (select_surfaces, resample_and_mask_wf, [(metric, 'inputnode.in_file')]),
+            (resample_and_mask_wf, cifti_metric, [('outputnode.out_file', 'scalar_surfs')]),
+            (cifti_metric, outputnode, [
+                ('out_file', f'{metric}_fsLR'),
+                ('out_metadata', f'{metric}_metadata'),
+            ]),
+        ])  # fmt:skip
+
+    return workflow
+
+
+def init_myelinmap_fsLR_wf(
+    grayord_density: ty.Literal['91k', '170k'],
+    omp_nthreads: int,
+    mem_gb: float,
+    name: str = 'myelinmap_fsLR_wf',
+):
+    """Resample myelinmap volume to fsLR surface.
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: colored
+            :simple_form: yes
+
+            from smriprep.workflows.surfaces import init_myelinmap_fsLR_wf
+            wf = init_myelinmap_fsLR_wf(grayord_density='91k', omp_nthreads=1, mem_gb=1)
+
+    Parameters
+    ----------
+    grayord_density : :class:`str`
+        Either ``"91k"`` or ``"170k"``, representing the total *grayordinates*.
+    omp_nthreads : :class:`int`
+        Maximum number of threads an individual process may use
+    mem_gb : :class:`float`
+        Size of BOLD file in GB
+    name : :class:`str`
+        Name of workflow (default: ``"myelinmap_fsLR_wf"``)
+
+    Inputs
+    ------
+    in_file : :class:`str`
+        Path to the myelin map in subject volume space
+    thickness : :class:`list` of :class:`str`
+        Path to left and right hemisphere thickness GIFTI shape files
+    midthickness : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surface files
+    midthickness_fsLR : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surface files in fsLR space
+    sphere_reg_fsLR : :class:`list` of :class:`str`
+        Path to left and right hemisphere sphere.reg GIFTI surface files,
+        mapping from subject to fsLR
+    cortex_mask : :class:`list` of :class:`str`
+        Path to left and right hemisphere cortex mask GIFTI files
+
+    Outputs
+    -------
+    out_fsLR : :class:`str`
+        Path to the resampled myelin map in fsLR space
+
+    """
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+    from niworkflows.interfaces.utility import KeySelect
+    from niworkflows.interfaces.workbench import VolumeToSurfaceMapping
+
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'in_file',
+                'thickness',
+                'midthickness',
+                'midthickness_fsLR',
+                'sphere_reg_fsLR',
+                'cortex_mask',
+                'volume_roi',
+            ]
+        ),
+        name='inputnode',
+    )
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=['out_file', 'out_metadata']),
+        name='outputnode',
+    )
+
+    hemisource = pe.Node(
+        niu.IdentityInterface(fields=['hemi']),
+        name='hemisource',
+        iterables=[('hemi', ['L', 'R'])],
+    )
+
     select_surfaces = pe.Node(
         KeySelect(
             fields=[
-                'curv',
-                'sulc',
                 'thickness',
-                'roi',
+                'midthickness',
+            ],
+            keys=['L', 'R'],
+        ),
+        name='select_surfaces',
+        run_without_submitting=True,
+    )
+
+    volume_to_surface = pe.Node(
+        VolumeToSurfaceMapping(method='myelin-style', sigma=fwhm2sigma(5)),
+        name='volume_to_surface',
+        mem_gb=mem_gb * 3,
+        n_procs=omp_nthreads,
+    )
+    # For SmoothedMyelinMap
+    # smooth = pe.Node(
+    #     MetricSmooth(sigma=fwhm2sigma(4), nearest=True),
+    #     name='metric_dilate',
+    #     mem_gb=1,
+    #     n_procs=omp_nthreads,
+    # )
+    resample_and_mask_wf = init_resample_and_mask_wf(
+        grayord_density=grayord_density,
+        omp_nthreads=omp_nthreads,
+        mem_gb=mem_gb,
+    )
+    cifti_myelinmap = pe.JoinNode(
+        GenerateDScalar(grayordinates=grayord_density, scalar_name='MyelinMap'),
+        name='cifti_myelinmap',
+        joinfield=['scalar_surfs'],
+        joinsource='hemisource',
+    )
+
+    workflow.connect([
+        (inputnode, select_surfaces, [
+            ('midthickness', 'midthickness'),
+            ('thickness', 'thickness'),
+        ]),
+        (hemisource, select_surfaces, [('hemi', 'key')]),
+        # Resample volume to native surface
+        (inputnode, volume_to_surface, [
+            ('in_file', 'volume_file'),
+            ('ribbon_file', 'ribbon_roi'),
+        ]),
+        (select_surfaces, volume_to_surface, [
+            ('midthickness', 'surface_file'),
+            ('thickness', 'thickness'),
+        ]),
+        (inputnode, resample_and_mask_wf, [
+            ('midthickness', 'inputnode.midthickness'),
+            ('midthickness_fsLR', 'inputnode.midthickness_fsLR'),
+            ('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
+            ('cortex_mask', 'inputnode.cortex_mask'),
+        ]),
+        (hemisource, resample_and_mask_wf, [('hemi', 'inputnode.hemi')]),
+        (volume_to_surface, resample_and_mask_wf, [('out_file', 'inputnode.in_file')]),
+        (resample_and_mask_wf, cifti_myelinmap, [('outputnode.out_file', 'scalar_surfs')]),
+        (cifti_myelinmap, outputnode, [
+            ('out_file', 'out_file'),
+            ('out_metadata', 'out_metadata'),
+        ]),
+    ])  # fmt:skip
+
+    return workflow
+
+
+def init_resample_and_mask_wf(
+    grayord_density: ty.Literal['91k', '170k'],
+    omp_nthreads: int,
+    mem_gb: float,
+    name: str = 'resample_and_mask_wf',
+):
+    """Resample GIFTI surfaces to fsLR space and mask with fsLR ROI.
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: colored
+            :simple_form: yes
+
+            from smriprep.workflows.surfaces import init_resample_and_mask_wf
+            wf = init_resample_and_mask_wf(
+                grayord_density='91k',
+                omp_nthreads=1,
+                mem_gb=1,
+            )
+
+    Parameters
+    ----------
+    grayord_density : :class:`str`
+        Either ``"91k"`` or ``"170k"``, representing the total *grayordinates*.
+    omp_nthreads : :class:`int`
+        Maximum number of threads an individual process may use
+    mem_gb : :class:`float`
+        Size of BOLD file in GB
+    name : :class:`str`
+        Name of workflow (default: ``"resample_and_mask_wf"``)
+
+    Inputs
+    ------
+    in_file : :class:`str`
+        Path to metric file in subject space
+    hemi : :class:`str`
+        Hemisphere identifier (``"L"`` or ``"R"``)
+    midthickness : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surfaces.
+    midthickness_fsLR : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surfaces in fsLR space.
+    sphere_reg_fsLR : :class:`list` of :class:`str`
+        Path to left and right hemisphere sphere.reg GIFTI surfaces, mapping from subject to fsLR
+    cortex_mask : :class:`list` of :class:`str`
+        Path to left and right hemisphere cortex mask GIFTI files
+
+    Outputs
+    -------
+    metric_fsLR : :class:`str`
+        Path to metric resampled as GIFTI file in fsLR space
+
+    """
+    import templateflow.api as tf
+    from nipype.pipeline import engine as pe
+    from niworkflows.interfaces.utility import KeySelect
+
+    fslr_density = '32k' if grayord_density == '91k' else '59k'
+
+    workflow = pe.Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'in_file',
+                'hemi',
+                'midthickness',
+                'midthickness_fsLR',
+                'sphere_reg_fsLR',
+                'cortex_mask',
+            ]
+        ),
+        name='inputnode',
+    )
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=['out_file']),
+        name='outputnode',
+    )
+
+    select_surfaces = pe.Node(
+        KeySelect(
+            fields=[
                 'midthickness',
                 'midthickness_fsLR',
                 'sphere_reg_fsLR',
                 'template_sphere',
+                'cortex_mask',
                 'template_roi',
             ],
             keys=['L', 'R'],
@@ -1559,60 +1827,41 @@ def init_morph_grayords_wf(
         )
         for hemi in ['L', 'R']
     ]
+    atlases = smriprep.load_data('atlases')
     select_surfaces.inputs.template_roi = [
         str(atlases / f'L.atlasroi.{fslr_density}_fs_LR.shape.gii'),
         str(atlases / f'R.atlasroi.{fslr_density}_fs_LR.shape.gii'),
     ]
 
+    resample_to_fsLR = pe.Node(
+        MetricResample(method='ADAP_BARY_AREA', area_surfs=True),
+        name='resample_to_fsLR',
+        mem_gb=1,
+        n_procs=omp_nthreads,
+    )
+    mask_fsLR = pe.Node(MetricMask(), name='mask_fsLR')
+
     workflow.connect([
         (inputnode, select_surfaces, [
-            ('curv', 'curv'),
-            ('sulc', 'sulc'),
-            ('thickness', 'thickness'),
-            ('roi', 'roi'),
             ('midthickness', 'midthickness'),
             ('midthickness_fsLR', 'midthickness_fsLR'),
             ('sphere_reg_fsLR', 'sphere_reg_fsLR'),
+            ('cortex_mask', 'cortex_mask'),
+            ('hemi', 'key'),
         ]),
-        (hemisource, select_surfaces, [('hemi', 'key')]),
+        (inputnode, resample_to_fsLR, [('in_file', 'in_file')]),
+        (select_surfaces, resample_to_fsLR, [
+            ('sphere_reg_fsLR', 'current_sphere'),
+            ('template_sphere', 'new_sphere'),
+            ('midthickness', 'current_area'),
+            ('midthickness_fsLR', 'new_area'),
+            ('cortex_mask', 'roi_metric'),
+        ]),
+        (select_surfaces, mask_fsLR, [('template_roi', 'mask')]),
+        (resample_to_fsLR, mask_fsLR, [('out_file', 'in_file')]),
+        (mask_fsLR, outputnode, [('out_file', 'out_file')]),
     ])  # fmt:skip
 
-    for metric in ('curv', 'sulc', 'thickness'):
-        resampler = pe.Node(
-            MetricResample(method='ADAP_BARY_AREA', area_surfs=True),
-            name=f'resample_{metric}',
-            n_procs=omp_nthreads,
-        )
-        mask_fsLR = pe.Node(
-            MetricMask(),
-            name=f'mask_{metric}',
-            n_procs=omp_nthreads,
-        )
-        cifti_metric = pe.JoinNode(
-            GenerateDScalar(grayordinates=grayord_density, scalar_name=metric),
-            name=f'cifti_{metric}',
-            joinfield=['scalar_surfs'],
-            joinsource='hemisource',
-        )
-
-        workflow.connect([
-            (select_surfaces, resampler, [
-                (metric, 'in_file'),
-                ('sphere_reg_fsLR', 'current_sphere'),
-                ('template_sphere', 'new_sphere'),
-                ('midthickness', 'current_area'),
-                ('midthickness_fsLR', 'new_area'),
-                ('roi', 'roi_metric'),
-            ]),
-            (select_surfaces, mask_fsLR, [('template_roi', 'mask')]),
-            (resampler, mask_fsLR, [('out_file', 'in_file')]),
-            (mask_fsLR, cifti_metric, [('out_file', 'scalar_surfs')]),
-            (cifti_metric, outputnode, [
-                ('out_file', f'{metric}_fsLR'),
-                ('out_metadata', f'{metric}_metadata'),
-            ]),
-        ])  # fmt:skip
-
     return workflow