diff --git a/smriprep/interfaces/calc.py b/smriprep/interfaces/calc.py new file mode 100644 index 0000000000..248ceae358 --- /dev/null +++ b/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 +# +# 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) diff --git a/smriprep/interfaces/tests/test_calc.py b/smriprep/interfaces/tests/test_calc.py new file mode 100644 index 0000000000..17d068b1f0 --- /dev/null +++ b/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 diff --git a/smriprep/workflows/anatomical.py b/smriprep/workflows/anatomical.py index 974d2c57a4..4c23da329a 100644 --- a/smriprep/workflows/anatomical.py +++ b/smriprep/workflows/anatomical.py @@ -57,6 +57,8 @@ 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 @@ -645,6 +647,7 @@ def init_anat_fit_wf( 'sphere_reg_fsLR', 'sphere_reg_msm', 'anat_ribbon', + 't1t2_ratio', # Reverse transform; not computable from forward transform 'std2anat_xfm', # Metadata @@ -1308,6 +1311,16 @@ def init_anat_fit_wf( LOGGER.info('ANAT Stage 8a: Found pre-computed cortical ribbon mask') outputnode.inputs.anat_ribbon = precomputed['anat_ribbon'] + if t2w and 't1t2ratio' not in precomputed: + LOGGER.info('ANAT Stage 8b: Creating T1w/T2w ratio map') + # Commented out to pacify linter. + # t1t2_ratio = pe.Node(T1T2Ratio(), name='t1t2_ratio') + + elif not t2w: + LOGGER.info('ANAT No T2w images provided - skipping Stage 8b') + else: + LOGGER.info('ANAT Found precomputed T1w/T2w ratio map - skipping Stage 8b') + # Stage 9: Baseline fsLR registration if len(precomputed.get('sphere_reg_fsLR', [])) < 2: LOGGER.info('ANAT Stage 9: Creating fsLR registration sphere') diff --git a/smriprep/workflows/surfaces.py b/smriprep/workflows/surfaces.py index 48895377a3..2424f4593a 100644 --- a/smriprep/workflows/surfaces.py +++ b/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 @@ -1447,19 +1449,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=[ @@ -1498,18 +1495,285 @@ 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, + ) + + 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)]), + (hemisource, select_surfaces, [('hemi', 'key')]), + (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, + ) + # 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'], @@ -1530,60 +1794,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