From 42481ea028f61902703a58e7fb710c407e132f33 Mon Sep 17 00:00:00 2001 From: madisoth <56737848+madisoth@users.noreply.github.com> Date: Sat, 3 Sep 2022 19:38:32 -0700 Subject: [PATCH 01/11] add options "--project-goodvoxels", "--surface-sampler", plus supporting Workbench interfaces --- nibabies/cli/parser.py | 18 + nibabies/config.py | 5 + nibabies/data/tests/config.toml | 2 + nibabies/interfaces/metric.py | 305 ++++++++++++ nibabies/interfaces/volume.py | 164 +++++++ nibabies/interfaces/wbvoltosurf.py | 298 ++++++++++++ nibabies/workflows/base.py | 2 + nibabies/workflows/bold/base.py | 17 +- nibabies/workflows/bold/resampling.py | 673 +++++++++++++++++++++++--- 9 files changed, 1424 insertions(+), 60 deletions(-) create mode 100644 nibabies/interfaces/metric.py create mode 100644 nibabies/interfaces/volume.py create mode 100644 nibabies/interfaces/wbvoltosurf.py diff --git a/nibabies/cli/parser.py b/nibabies/cli/parser.py index 25b847ed..6be65545 100644 --- a/nibabies/cli/parser.py +++ b/nibabies/cli/parser.py @@ -333,6 +333,24 @@ def _slice_time_ref(value, parser): help="Replace medial wall values with NaNs on functional GIFTI files. Only " "performed for GIFTI files mapped to a freesurfer subject (fsaverage or fsnative).", ) + g_conf.add_argument( + "--project-goodvoxels", + required=False, + action="store_true", + default=False, + help="Exclude voxels whose timeseries have locally high coefficient of variation " + "from surface resampling. Only performed for GIFTI files mapped to a freesurfer subject " + "(fsaverage or fsnative).", + ) + g_conf.add_argument( + "--surface-sampler", + action="store", + choices=("fs", "wb"), + default="fs", + help="select sampling method used for projection of BOLD timeseries to surface. " + "'fs' to use FreeSurfer mri_vol2surf (default). " + "'wb' to use Workbench wb_command -volume-to-surface-mapping.", + ) g_conf.add_argument( "--slice-time-ref", required=False, diff --git a/nibabies/config.py b/nibabies/config.py index 0b3847d2..333f0d23 100644 --- a/nibabies/config.py +++ b/nibabies/config.py @@ -553,6 +553,11 @@ class workflow(_Config): """Run FreeSurfer ``recon-all`` with the ``-logitudinal`` flag.""" medial_surface_nan = None """Fill medial surface with :abbr:`NaNs (not-a-number)` when sampling.""" + project_goodvoxels = False + """Exclude voxels with locally high coefficient of variation from sampling.""" + surface_sampler = "fs" + """``fs`` (default) or ``wb`` to select FreeSurfer-based or Workbench-based + method for sampling BOLD to surface.""" regressors_all_comps = None """Return all CompCor components.""" regressors_dvars_th = None diff --git a/nibabies/data/tests/config.toml b/nibabies/data/tests/config.toml index b4b398f8..6cb874b1 100644 --- a/nibabies/data/tests/config.toml +++ b/nibabies/data/tests/config.toml @@ -52,6 +52,8 @@ hires = true ignore = [ "slicetiming",] longitudinal = false medial_surface_nan = false +project_goodvoxels = false +surface_sampler = "fs" regressors_all_comps = false regressors_dvars_th = 1.5 regressors_fd_th = 0.5 diff --git a/nibabies/interfaces/metric.py b/nibabies/interfaces/metric.py new file mode 100644 index 00000000..023fda7e --- /dev/null +++ b/nibabies/interfaces/metric.py @@ -0,0 +1,305 @@ +# -*- coding: utf-8 -*- +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""This module provides interfaces for workbench surface commands""" +import os + +from nipype.interfaces.base import TraitedSpec, File, traits, CommandLineInputSpec +from nipype.interfaces.workbench.base import WBCommand +from nipype import logging + +iflogger = logging.getLogger("nipype.interface") + + +class MetricDilateInputSpec(CommandLineInputSpec): + + in_file = File( + exists=True, + mandatory=True, + argstr="%s ", + position=0, + desc="the metric to dilate", + ) + + surface = File( + exists=True, + mandatory=True, + argstr="%s ", + position=1, + desc="the surface to compute on", + ) + + distance = traits.Float( + mandatory=True, + argstr="%f ", + position=2, + desc="distance in mm to dilate", + ) + + out_file = File( + name_source=["in_file"], + name_template="%s.func.gii", + keep_extension=False, + argstr="%s ", + position=3, + desc="output - the output metric", + ) + + bad_vertex_roi = File( + argstr="-bad-vertex-roi %s ", + position=4, + desc="metric file, positive values denote vertices to have their values replaced", + ) + + data_roi = File( + argstr="-data-roi %s ", + position=5, + desc="metric file, positive values denote vertices that have data", + ) + + column = traits.Int( + position=6, + argstr="-column %d ", + desc="the column number", + ) + + nearest = traits.Bool( + position=7, + argstr="-nearest ", + desc="use the nearest good value instead of a weighted average", + ) + + linear = traits.Bool( + position=8, + argstr="-linear ", + desc="fill in values with linear interpolation along strongest gradient", + ) + + exponent = traits.Float( + argstr="-exponent %f ", + position=9, + default=6.0, + desc="exponent n to use in (area / (distance ^ n)) as the " + "weighting function (default 6)", + ) + + corrected_areas = File( + argstr="-corrected-areas %s ", + position=10, + desc="vertex areas to use instead of computing them from the surface", + ) + + legacy_cutoff = traits.Bool( + position=11, + argstr="-legacy-cutoff ", + desc="use the v1.3.2 method of choosing how many vertices to " + "use when calulating the dilated value with weighted method", + ) + + +class MetricDilateOutputSpec(TraitedSpec): + out_file = File(exists=True, desc="output file") + + +class MetricDilate(WBCommand): + """ + For all data values designated as bad, if they neighbor a good value or + are within the specified distance of a good value in the same kind of + model, replace the value with a distance weighted average of nearby good + values, otherwise set the value to zero. If -nearest is specified, it + will use the value from the closest good value within range instead of a + weighted average. When the input file contains label data, nearest + dilation is used on the surface, and weighted popularity is used in the + volume. + + The -corrected-areas options are intended for dilating on group average + surfaces, but it is only an approximate correction for the reduction of + structure in a group average surface. + + If -bad-vertex-roi is specified, all values, including those with + value zero, are good, except for locations with a positive value in the + ROI. If it is not specified, only values equal to zero are bad. + + >>> from nipype.interfaces.workbench import MetricDilate + >>> metdil = MetricDilate() + >>> metdil.inputs.in_file = 'sub-01_task-rest_bold_space-fsaverage5.L.func.gii' + >>> metdil.inputs.surface = 'sub-01.L.midthickness.surf.gii' + >>> metdil.inputs.distance = 10 + >>> metdil.inputs.nearest = True + >>> metdil.cmdline + 'wb_command -metric-dilate \ + sub-01_task-rest_bold_space-fsaverage5.L.func.gii \ + sub-01.L.midthickness.surf.gii \ + 10 \ + sub-01_task-rest_bold_space-fsaverage5.L.func.gii \ + -nearest' + """ + input_spec = MetricDilateInputSpec + output_spec = MetricDilateOutputSpec + _cmd = "wb_command -metric-dilate " + + +class MetricResampleInputSpec(CommandLineInputSpec): + in_file = File( + exists=True, + mandatory=True, + argstr="%s", + position=0, + desc="The metric file to resample", + ) + current_sphere = File( + exists=True, + mandatory=True, + argstr="%s", + position=1, + desc="A sphere surface with the mesh that the metric is currently on", + ) + new_sphere = File( + exists=True, + mandatory=True, + argstr="%s", + position=2, + desc="A sphere surface that is in register with and" + " has the desired output mesh", + ) + method = traits.Enum( + "ADAP_BARY_AREA", + "BARYCENTRIC", + argstr="%s", + mandatory=True, + position=3, + desc="The method name - ADAP_BARY_AREA method is recommended for" + " ordinary metric data, because it should use all data while" + " downsampling, unlike BARYCENTRIC. If ADAP_BARY_AREA is used," + " exactly one of area_surfs or area_metrics must be specified", + ) + out_file = File( + name_source=["new_sphere"], + name_template="%s.out", + keep_extension=True, + argstr="%s", + position=4, + desc="The output metric", + ) + area_surfs = traits.Bool( + position=5, + argstr="-area-surfs", + xor=["area_metrics"], + desc="Specify surfaces to do vertex area correction based on", + ) + area_metrics = traits.Bool( + position=5, + argstr="-area-metrics", + xor=["area_surfs"], + desc="Specify vertex area metrics to do area correction based on", + ) + current_area = File( + exists=True, + position=6, + argstr="%s", + desc="A relevant anatomical surface with mesh OR" + " a metric file with vertex areas for mesh", + ) + new_area = File( + exists=True, + position=7, + argstr="%s", + desc="A relevant anatomical surface with mesh OR" + " a metric file with vertex areas for mesh", + ) + roi_metric = File( + exists=True, + position=8, + argstr="-current-roi %s", + desc="Input roi on the current mesh used to exclude non-data vertices", + ) + valid_roi_out = traits.Bool( + position=9, + argstr="-valid-roi-out", + desc="Output the ROI of vertices that got data from valid source vertices", + ) + largest = traits.Bool( + position=10, + argstr="-largest", + desc="Use only the value of the vertex with the largest weight", + ) + + +class MetricResampleOutputSpec(TraitedSpec): + out_file = File(exists=True, desc="the output metric") + roi_file = File(desc="ROI of vertices that got data from valid source vertices") + + +class MetricResample(WBCommand): + """ + Resample a metric file to a different mesh + + Resamples a metric file, given two spherical surfaces that are in + register. If ``ADAP_BARY_AREA`` is used, exactly one of -area-surfs or + ``-area-metrics`` must be specified. + + The ``ADAP_BARY_AREA`` method is recommended for ordinary metric data, + because it should use all data while downsampling, unlike ``BARYCENTRIC``. + The recommended areas option for most data is individual midthicknesses + for individual data, and averaged vertex area metrics from individual + midthicknesses for group average data. + + The ``-current-roi`` option only masks the input, the output may be slightly + dilated in comparison, consider using ``-metric-mask`` on the output when + using ``-current-roi``. + + The ``-largest option`` results in nearest vertex behavior when used with + ``BARYCENTRIC``. When resampling a binary metric, consider thresholding at + 0.5 after resampling rather than using ``-largest``. + + >>> from nipype.interfaces.workbench import MetricResample + >>> metres = MetricResample() + >>> metres.inputs.in_file = 'sub-01_task-rest_bold_space-fsaverage5.L.func.gii' + >>> metres.inputs.method = 'ADAP_BARY_AREA' + >>> metres.inputs.current_sphere = 'fsaverage5_std_sphere.L.10k_fsavg_L.surf.gii' + >>> metres.inputs.new_sphere = 'fs_LR-deformed_to-fsaverage.L.sphere.32k_fs_LR.surf.gii' + >>> metres.inputs.area_metrics = True + >>> metres.inputs.current_area = 'fsaverage5.L.midthickness_va_avg.10k_fsavg_L.shape.gii' + >>> metres.inputs.new_area = 'fs_LR.L.midthickness_va_avg.32k_fs_LR.shape.gii' + >>> metres.cmdline + 'wb_command -metric-resample sub-01_task-rest_bold_space-fsaverage5.L.func.gii \ + fsaverage5_std_sphere.L.10k_fsavg_L.surf.gii \ + fs_LR-deformed_to-fsaverage.L.sphere.32k_fs_LR.surf.gii \ + ADAP_BARY_AREA fs_LR-deformed_to-fsaverage.L.sphere.32k_fs_LR.surf.out \ + -area-metrics fsaverage5.L.midthickness_va_avg.10k_fsavg_L.shape.gii \ + fs_LR.L.midthickness_va_avg.32k_fs_LR.shape.gii' + """ + + input_spec = MetricResampleInputSpec + output_spec = MetricResampleOutputSpec + _cmd = "wb_command -metric-resample" + + def _format_arg(self, opt, spec, val): + if opt in ["current_area", "new_area"]: + if not self.inputs.area_surfs and not self.inputs.area_metrics: + raise ValueError( + "{} was set but neither area_surfs or" + " area_metrics were set".format(opt) + ) + if opt == "method": + if ( + val == "ADAP_BARY_AREA" + and not self.inputs.area_surfs + and not self.inputs.area_metrics + ): + raise ValueError( + "Exactly one of area_surfs or area_metrics" " must be specified" + ) + if opt == "valid_roi_out" and val: + # generate a filename and add it to argstr + roi_out = self._gen_filename(self.inputs.in_file, suffix="_roi") + iflogger.info("Setting roi output file as", roi_out) + spec.argstr += " " + roi_out + return super(MetricResample, self)._format_arg(opt, spec, val) + + def _list_outputs(self): + outputs = super(MetricResample, self)._list_outputs() + if self.inputs.valid_roi_out: + roi_file = self._gen_filename(self.inputs.in_file, suffix="_roi") + outputs["roi_file"] = os.path.abspath(roi_file) + return outputs diff --git a/nibabies/interfaces/volume.py b/nibabies/interfaces/volume.py new file mode 100644 index 00000000..8866f990 --- /dev/null +++ b/nibabies/interfaces/volume.py @@ -0,0 +1,164 @@ +# -*- coding: utf-8 -*- +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""This module provides interfaces for workbench volume commands""" +import os + +from nipype.interfaces.base import TraitedSpec, File, traits, CommandLineInputSpec +from nipype.interfaces.workbench.base import WBCommand +from nipype import logging + +iflogger = logging.getLogger("nipype.interface") + + +class CreateSignedDistanceVolumeInputSpec(CommandLineInputSpec): + surface = File( + exists=True, + mandatory=True, + argstr="%s ", + position=0, + desc="the input surface", + ) + ref_space = File( + exists=True, + mandatory=True, + argstr="%s ", + position=1, + desc="a volume in the desired output space (dims, spacing, origin)", + ) + out_vol = File( + name_source=["surface"], + name_template="%s.distvol.nii.gz", + argstr="%s ", + position=2, + desc="output - the output volume", + ) + roi_out = File( + name_source=["surface"], + name_template="%s.distvolroi.nii.gz", + argstr="-roi-out %s ", + position=3, + desc="output roi volume of where the output has a computed value", + ) + fill_value = traits.Float( + mandatory=False, + argstr="-fill-value %f ", + position=4, + desc="specify a value to put in all voxels that don't get assigned a distance, default 0", + ) + exact_limit = traits.Float( + mandatory=False, + argstr="-exact-limit %f ", + position=5, + desc="specify distance for exact output in mm, default 5", + ) + approx_limit = traits.Float( + mandatory=False, + argstr="-approx-limit %f ", + position=6, + desc="specify distance for approximate output in mm, default 20", + ) + approx_neighborhood = traits.Int( + mandatory=False, + argstr="-approx-neighborhood %d ", + position=7, + desc="size of neighborhood cube measured from center to face, default 2 = 5x5x5", + ) + winding_method = traits.Enum( + "EVEN_ODD", + "NEGATIVE", + "NONZERO", + "NORMALS", + argstr="-winding %s ", + usedefault=True, + position=8, + desc="winding method for point inside surface test, choices: " + "EVEN_ODD (default) " + "NEGATIVE " + "NONZERO " + "NORMALS " + ) + +class CreateSignedDistanceVolumeOutputSpec(TraitedSpec): + out_vol = File(exists=True, desc="output - the output volume") + roi_out = File(desc="output roi volume of where the output has a computed value") + + +class CreateSignedDistanceVolume(WBCommand): + """ +CREATE SIGNED DISTANCE VOLUME FROM SURFACE + wb_command -create-signed-distance-volume + - the input surface + - a volume in the desired output space (dims, spacing, origin) + - output - the output volume + + [-roi-out] - output an roi volume of where the output has a computed + value + - output - the output roi volume + + [-fill-value] - specify a value to put in all voxels that don't get + assigned a distance + - value to fill with (default 0) + + [-exact-limit] - specify distance for exact output + - distance in mm (default 5) + + [-approx-limit] - specify distance for approximate output + - distance in mm (default 20) + + [-approx-neighborhood] - voxel neighborhood for approximate calculation + - size of neighborhood cube measured from center to face, in + voxels (default 2 = 5x5x5) + + [-winding] - winding method for point inside surface test + - name of the method (default EVEN_ODD) + + Computes the signed distance function of the surface. Exact distance is + calculated by finding the closest point on any surface triangle to the + center of the voxel. Approximate distance is calculated starting with + these distances, using dijkstra's method with a neighborhood of voxels. + Specifying too small of an exact distance may produce unexpected results. + Valid specifiers for winding methods are as follows: + + EVEN_ODD (default) + NEGATIVE + NONZERO + NORMALS + + The NORMALS method uses the normals of triangles and edges, or the + closest triangle hit by a ray from the point. This method may be + slightly faster, but is only reliable for a closed surface that does not + cross through itself. All other methods count entry (positive) and exit + (negative) crossings of a vertical ray from the point, then counts as + inside if the total is odd, negative, or nonzero, respectively. + + >>> from nipype.interfaces.workbench import CreateSignedDistanceVolume + >>> distvol = CreateSignedDistanceVolume() + >>> distvol.inputs.surface = 'sub-01.L.pial.native.surf.gii' + >>> distvol.inputs.refspace = 'sub-01_T1w.nii.gz' + >>> distvol.inputs.out_vol = 'sub-01.L.pial.native.surf.distvol.nii.gz' + >>> distvol.inputs.roi_out = 'sub-01.L.pial.native.surf.distvolroi.nii.gz' + >>> distvol.inputs.fill_value = 0 + >>> distvol.inputs.exact_limit = 5 + >>> distvol.inputs.approx_limit = 20 + >>> distvol.inputs.approx_neighborhood = 2 + >>> distvol.inputs.winding_method = 'EVEN_ODD' + >>> distvol.cmdline + 'wb_command -create-signed-distance-volume sub-01.L.pial.native.surf.gii \ + sub-01_T1w.nii.gz \ + sub-01.L.pial.native.surf.distvol.nii.gz \ + -roi-out sub-01.L.pial.native.surf.distvolroi.nii.gz \ + -fill-value 0 \ + -exact-limit 5 \ + -approx-limit 20 \ + -approx-neighborhood 2 \ + -winding EVEN_ODD' + """ + + input_spec = CreateSignedDistanceVolumeInputSpec + output_spec = CreateSignedDistanceVolumeOutputSpec + _cmd = "wb_command -create-signed-distance-volume" + + def _list_outputs(self): + outputs = super(CreateSignedDistanceVolume, self)._list_outputs() + return outputs diff --git a/nibabies/interfaces/wbvoltosurf.py b/nibabies/interfaces/wbvoltosurf.py new file mode 100644 index 00000000..12409b74 --- /dev/null +++ b/nibabies/interfaces/wbvoltosurf.py @@ -0,0 +1,298 @@ +# -*- coding: utf-8 -*- +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""This module provides interfaces for workbench surface commands""" +import os + +from nipype.interfaces.base import TraitedSpec, File, traits, CommandLineInputSpec +from nipype.interfaces.workbench.base import WBCommand +from nipype import logging + +iflogger = logging.getLogger("nipype.interface") + + +class VolumeToSurfaceMappingInputSpec(CommandLineInputSpec): + + in_file = File( + exists=True, + mandatory=True, + argstr="%s ", + position=0, + desc="the volume to map data from", + ) + + surface = File( + exists=True, + mandatory=True, + argstr="%s ", + position=1, + desc="the surface to map the data onto", + ) + + out_file = File( + name_source=["in_file"], + name_template="%s.func.gii", + keep_extension=False, + argstr="%s ", + position=2, + desc="output - the output metric file", + ) + + mapping_method = traits.Enum( + "trilinear", + "enclosing", + "cubic", + "ribbon-constrained", + "myelin-style", + position=3, + argstr="-%s ", + desc="choose mapping method: trilinear, enclosing voxel, cubic spline, " + "ribbon-constrained, or myelin-style", + ) + + inner_surf = File( + argstr="%s ", + position=4, + desc="the inner surface of the ribbon", + ) + + outer_surf = File( + argstr="%s ", + position=5, + desc="the outer surface of the ribbon", + ) + + roi_volume = File( + argstr="-volume-roi %s ", + position=6, + desc="the roi volume file to use", + ) + + weighted = traits.Bool( + position=7, + argstr="-weighted ", + desc="treat the roi values as weightings rather than binary", + ) + + subdiv_num = traits.Int( + position=8, + argstr="-voxel-subdiv %d ", + desc="number of subdivisions while estimating voxel weights, " + "default 3 if this option is not specified", + ) + + thin_columns = traits.Bool( + position=9, + argstr="-thin-columns ", + desc="use non-overlapping polyhedra", + ) + + gaussian_scale = traits.Float( + argstr="-gaussian %f ", + position=10, + desc="value to multiply the local thickness by, to get the gaussian sigma. " + "reduce weight to voxels that aren't near ", + ) + + interpolate_method = traits.Enum( + "CUBIC", + "ENCLOSING_VOXEL", + "TRILINEAR", + position=11, + argstr="-interpolate %s ", + desc="must be CUBIC, ENCLOSING_VOXEL, or TRILINEAR. " + "instead of a weighted average of voxels, " + "interpolate at subpoints inside the ribbon", + ) + + bad_vertices_out = File( + argstr="-bad-vertices-out %s ", + position=12, + desc="output an ROI of which vertices didn't intersect any valid voxels", + ) + + output_weights = traits.Bool( + position=13, + argstr="-output-weights ", + desc="the column number", + ) + + output_weights_vertex = traits.Int( + position=14, + argstr="%d ", + desc="the column number", + ) + + output_weights_weights_out = File( + argstr="%s ", + position=15, + desc="output - volume to write the weights to", + ) + + output_weights_text_out = File( + argstr="%s ", + position=16, + desc="text file name. " + "write the voxel weights for all vertices to a text file", + ) + + myelin_style_ribbon_roi = File( + argstr="%s ", + position=17, + desc="roi volume of the cortical ribbon for this hemisphere", + ) + + myelin_style_thickness = File( + argstr="%s ", + position=18, + desc="a metric file of cortical thickness", + ) + + myelin_style_sigma = traits.Float( + argstr="%f ", + position=19, + desc="gaussian kernel in mm for weighting voxels within range", + ) + + legacy_bug = traits.Bool( + position=20, + argstr="-legacy-bug ", + desc="emulate old Workbench v1.2.3 and earlier code that didn't follow a cylinder cutoff", + ) + + subvol_select = File( + argstr="-subvol-select %s ", + position=21, + desc="the subvolume number or name", + ) + + +class VolumeToSurfaceMappingOutputSpec(TraitedSpec): + out_file = File(exists=True, desc="the output metric file") + bad_vertices_out = File(desc="ROI of which vertices didn't intersect any valid voxels") + output_weights_weights_out = File(desc="volume to which the specified vertex's weights are written") + output_weights_text_out = File(desc="text file of the voxel weights for all vertices") + + +class VolumeToSurfaceMapping(WBCommand): + """ + You must specify exactly one mapping method. Enclosing voxel uses the + value from the voxel the vertex lies inside, while trilinear does a 3D + linear interpolation based on the voxels immediately on each side of the + vertex's position. + + The ribbon mapping method constructs a polyhedron from the vertex's + neighbors on each surface, and estimates the amount of this polyhedron's + volume that falls inside any nearby voxels, to use as the weights for + sampling. If -thin-columns is specified, the polyhedron uses the edge + midpoints and triangle centroids, so that neighboring vertices do not + have overlapping polyhedra. This may require increasing -voxel-subdiv to + get enough samples in each voxel to reliably land inside these smaller + polyhedra. The volume ROI is useful to exclude partial volume effects of + voxels the surfaces pass through, and will cause the mapping to ignore + voxels that don't have a positive value in the mask. The subdivision + number specifies how it approximates the amount of the volume the + polyhedron intersects, by splitting each voxel into NxNxN pieces, and + checking whether the center of each piece is inside the polyhedron. If + you have very large voxels, consider increasing this if you get zeros in + your output. The -gaussian option makes it act more like the myelin + method, where the distance of a voxel from is used to + downweight the voxel. The -interpolate suboption, instead of doing a + weighted average of voxels, interpolates from the volume at the + subdivided points inside the ribbon. If using both -interpolate and the + -weighted suboption to -volume-roi, the roi volume weights are linearly + interpolated, unless the -interpolate method is ENCLOSING_VOXEL, in which + case ENCLOSING_VOXEL is also used for sampling the roi volume weights. + + The myelin style method uses part of the caret5 myelin mapping command to + do the mapping: for each surface vertex, take all voxels that are in a + cylinder with radius and height equal to cortical thickness, centered on + the vertex and aligned with the surface normal, and that are also within + the ribbon ROI, and apply a gaussian kernel with the specified sigma to + them to get the weights to use. The -legacy-bug flag reverts to the + unintended behavior present from the initial implementation up to and + including v1.2.3, which had only the tangential cutoff and a bounding box + intended to be larger than where the cylinder cutoff should have been. + + >>> from nipype.interfaces.workbench import VolumeToSurfaceMapping + >>> vtsm = VolumeToSurfaceMapping() + >>> vtsm.inputs.in_file = 'sub-01_task-rest_run-1_bold.nii.gz' + >>> vtsm.inputs.surface = 'sub-01.L.midthickness.surf.gii' + >>> vtsm.inputs.ribbon_constrained = True + >>> vtsm.inputs.inner_surf = 'sub-01.L.white.surf.gii' + >>> vtsm.inputs.outer_surf = 'sub-01.L.pial.surf.gii' + >>> vtsm.inputs.roi_volume = 'sub-01_task-rest_run-1_goodvoxels.nii.gz' + >>> vtsm.cmdline + 'wb_command -volume-to-surface-mapping \ + sub-01_task-rest_run-1_bold.nii.gz \ + sub-01.L.midthickness.surf.gii \ + sub-01_task-rest_run-1_bold.func.gii \ + -ribbon-constrained sub-01.L.white.surf.gii sub-01.L.pial.surf.gii \ + -volume-roi sub-01_task-rest_run-1_goodvoxels.nii.gz + """ + input_spec = VolumeToSurfaceMappingInputSpec + output_spec = VolumeToSurfaceMappingOutputSpec + _cmd = "wb_command -volume-to-surface-mapping " + + def _format_arg(self, opt, spec, val): + if opt in [ + "inner_surf", + "outer_surf", + "roi_volume", + "subdiv_num", + "thin_columns", + "gaussian_scale", + "interpolate_method", + "bad_vertices_out", + "output_weights", + "output_weights_text_out", + ]: + if self.inputs.mapping_method.val is not "ribbon-constrained": + raise ValueError( + "{} was set, but ribbon-constrained mapping method was not".format(opt) + ) + + if opt == "weighted": + if (val == True) and (self.inputs.mapping_method.val is not "ribbon-constrained"): + raise ValueError( + "weighted was set, but the required roi_volume was not" + ) + + if ((self.inputs.mapping_method.val == "ribbon-constrained") + and not (self.inputs.inner_surf and self.inputs.outer_surf) + ): + raise ValueError( + "mapping method is ribbon-constrained but at least one of " + "inner_surf, outer_surf was not set" + ) + + if opt == "output_weights": + if ((val == True) and not (self.inputs.output_weights_vertex + and self.inputs.output_weights_weights_out) + ): + raise ValueError( + "output_weights was set but at least one of " + "output_weights_vertex, output_weights_weights_out was not set" + ) + + if ((self.inputs.mapping_method.val == "myelin-style") + and not (self.inputs.myelin_style_ribbon_roi + and self.inputs.myelin_style_thickness + and self.inputs.myelin_style_sigma) + ): + raise ValueError( + "mapping method is myelin-style but at least one of myelin_style_ribbon_roi, " + "myelin_style_thickness, myelin_style_sigma was not set" + ) + + if opt == "legacy_bug": + if (val == True) and (not self.inputs.mapping_method.val is not "myelin-style"): + raise ValueError( + "legacy_bug was set, but the mapping method is not myelin-style" + ) + + return super(VolumeToSurfaceMapping, self)._format_arg(opt, spec, val) + + def _list_outputs(self): + outputs = super(VolumeToSurfaceMapping, self)._list_outputs() + return outputs \ No newline at end of file diff --git a/nibabies/workflows/base.py b/nibabies/workflows/base.py index 104f1c35..a4d16864 100644 --- a/nibabies/workflows/base.py +++ b/nibabies/workflows/base.py @@ -198,6 +198,7 @@ def init_single_subject_wf(subject_id, session_id=None): derivatives = config.execution.derivatives or {} anat_modality = "t1w" if subject_data["t1w"] else "t2w" spaces = config.workflow.spaces + surface_sampler = config.workflow.surface_sampler # Make sure we always go through these two checks if not anat_only and not subject_data["bold"]: task_id = config.execution.task_id @@ -478,6 +479,7 @@ def init_single_subject_wf(subject_id, session_id=None): ('outputnode.subject_id', 'inputnode.subject_id'), ('outputnode.t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'), ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ('outputnode.surfaces', 'inputnode.anat_giftis'), ]), ]) # fmt: on diff --git a/nibabies/workflows/bold/base.py b/nibabies/workflows/bold/base.py index 4dd801b6..67d5e88c 100644 --- a/nibabies/workflows/bold/base.py +++ b/nibabies/workflows/bold/base.py @@ -302,9 +302,13 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non Gridded (volumetric) resamplings were performed using `antsApplyTransforms` (ANTs), configured with Lanczos interpolation to minimize the smoothing effects of other kernels [@lanczos]. -Non-gridded (surface) resamplings were performed using `mri_vol2surf` -(FreeSurfer). -""" +Non-gridded (surface) resamplings were performed using {method}.. +""".format( + method={"fs": "mri_vol2surf (FreeSurfer)", + "wb": "wb_command -volume-to-surface-mapping (Workbench)" + }[config.workflow.surface_sampler] + ) + inputnode = pe.Node( niu.IdentityInterface( @@ -321,6 +325,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non "anat2std_xfm", "std2anat_xfm", "template", + "anat_giftis", # from bold reference workflow "bold_ref", "bold_ref_xfm", @@ -851,6 +856,8 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non mem_gb=mem_gb["resampled"], surface_spaces=freesurfer_spaces, medial_surface_nan=config.workflow.medial_surface_nan, + project_goodvoxels=config.workflow.project_goodvoxels, + surface_sampler=config.workflow.surface_sampler, name="bold_surf_wf", ) # fmt:off @@ -858,7 +865,9 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non (inputnode, bold_surf_wf, [ ('subjects_dir', 'inputnode.subjects_dir'), ('subject_id', 'inputnode.subject_id'), - ('t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm')]), + ('t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'), + ("anat_giftis", "inputnode.anat_giftis"), + ("anat_mask", "inputnode.t1w_mask")]), (bold_t1_trans_wf, bold_surf_wf, [('outputnode.bold_t1', 'inputnode.source_file')]), (bold_surf_wf, outputnode, [('outputnode.surfaces', 'surfaces')]), (bold_surf_wf, func_derivatives_wf, [ diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index ec38ab69..eae9c2f5 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -11,18 +11,41 @@ """ import nipype.interfaces.workbench as wb from nipype.interfaces import freesurfer as fs +from nipype.interfaces import fsl from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe +from nipype import Function +from niworkflows.interfaces.freesurfer import MakeMidthickness +from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms +from niworkflows.interfaces.freesurfer import MedialNaNs +from ...interfaces.volume import CreateSignedDistanceVolume +from ...interfaces.metric import MetricDilate, MetricResample +from ...interfaces.wbvoltosurf import VolumeToSurfaceMapping + + from ...config import DEFAULT_MEMORY_MIN_GB -def init_bold_surf_wf(mem_gb, surface_spaces, medial_surface_nan, name="bold_surf_wf"): +def init_bold_surf_wf(mem_gb, + surface_spaces, + medial_surface_nan, + project_goodvoxels, + surface_sampler, + name="bold_surf_wf"): """ Sample functional images to FreeSurfer surfaces. For each vertex, the cortical ribbon is sampled at six points (spaced 20% of thickness apart) and averaged. + + If --surface-sampler wb is used, Workbench's wb_command -volume-to-surface-mapping + with -ribbon-constrained option is used instead of the default FreeSurfer mri_vol2surf. + Note that unlike HCP, no additional spatial smoothing is applied to the surface-projected + data. + + If --project-goodvoxels is used, a "goodvoxels" BOLD mask, as described in [@hcppipelines], + is generated and applied to the functional image before sampling to surface. Outputs are in GIFTI format. Workflow Graph @@ -33,7 +56,9 @@ def init_bold_surf_wf(mem_gb, surface_spaces, medial_surface_nan, name="bold_sur from fmriprep.workflows.bold import init_bold_surf_wf wf = init_bold_surf_wf(mem_gb=0.1, surface_spaces=['fsnative', 'fsaverage5'], - medial_surface_nan=False) + medial_surface_nan=False, + project_goodvoxels=False, + surface_sampler="fs") Parameters ---------- @@ -44,19 +69,29 @@ def init_bold_surf_wf(mem_gb, surface_spaces, medial_surface_nan, name="bold_sur native surface. medial_surface_nan : :obj:`bool` Replace medial wall values with NaNs on functional GIFTI files + project_goodvoxels : :obj:`bool` + Exclude voxels with locally high coefficient of variation, or that lie outside the + cortical surfaces, from the surface projection. + surface_sampler : :obj:`str` + 'fs' (default) or 'wb' to specify FreeSurfer-based or Workbench-based + volume to surface mapping Inputs ------ source_file - Motion-corrected BOLD series in T1 space - t1w_preproc - Bias-corrected structural template image + Motion-corrected BOLD series in T1w space subjects_dir FreeSurfer SUBJECTS_DIR subject_id FreeSurfer subject ID t1w2fsnative_xfm LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space + itk_bold_to_t1 + Affine transform from ``ref_bold_brain`` to T1w space (ITK format) + anat_giftis + GIFTI anatomical surfaces in T1w space + t1w_mask + Mask of the skull-stripped T1w image Outputs ------- @@ -67,26 +102,44 @@ def init_bold_surf_wf(mem_gb, surface_spaces, medial_surface_nan, name="bold_sur from nipype.interfaces.io import FreeSurferSource from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.surf import GiftiSetAnatomicalStructure + import templateflow as tf + workflow = Workflow(name=name) workflow.__desc__ = """\ The BOLD time-series were resampled onto the following surfaces (FreeSurfer reconstruction nomenclature): -{out_spaces}. +{out_spaces} with {sampling_method} """.format( out_spaces=", ".join(["*%s*" % s for s in surface_spaces]) ) + if project_goodvoxels: + workflow.__desc__ += """\ +Before resampling, a "goodvoxels" mask [@hcppipelines] was applied, +excluding voxels whose time-series have a locally high coetfficient of +variation, or that lie outside the cortical surfaces, from the +surface projection. +""" + inputnode = pe.Node( niu.IdentityInterface( - fields=["source_file", "subject_id", "subjects_dir", "t1w2fsnative_xfm"] + fields=["source_file", + "subject_id", + "subjects_dir", + "t1w2fsnative_xfm", + "anat_giftis", + "t1w_mask"] ), name="inputnode", ) + itersource = pe.Node(niu.IdentityInterface(fields=["target"]), name="itersource") itersource.iterables = [("target", surface_spaces)] - get_fsnative = pe.Node(FreeSurferSource(), name="get_fsnative", run_without_submitting=True) + get_fsnative = pe.Node( + FreeSurferSource(), name="get_fsnative", run_without_submitting=True + ) def select_target(subject_id, space): """Get the target subject ID, given a source subject ID and a target space.""" @@ -106,7 +159,11 @@ def select_target(subject_id, space): run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB, ) - itk2lta = pe.Node(niu.Function(function=_itk2lta), name="itk2lta", run_without_submitting=True) + + itk2lta = pe.Node( + niu.Function(function=_itk2lta), name="itk2lta", run_without_submitting=True + ) + sampler = pe.MapNode( fs.SampleToSurface( cortex_mask=True, @@ -117,11 +174,20 @@ def select_target(subject_id, space): sampling_range=(0, 1, 0.2), sampling_units="frac", ), + name_source=['source_file'], + keep_extension=False, + name_template='%s.func.gii', iterfield=["hemi"], name="sampler", mem_gb=mem_gb * 3, ) sampler.inputs.hemi = ["lh", "rh"] + + # Refine if medial vertices should be NaNs + medial_nans = pe.MapNode( + MedialNaNs(), iterfield=["in_file"], name="medial_nans", mem_gb=DEFAULT_MEMORY_MIN_GB + ) + update_metadata = pe.MapNode( GiftiSetAnatomicalStructure(), iterfield=["in_file"], @@ -134,49 +200,550 @@ def select_target(subject_id, space): joinsource="itersource", name="outputnode", ) + + if not project_goodvoxels: + # fmt: off + workflow.connect([ + (inputnode, get_fsnative, [('subject_id', 'subject_id'), + ('subjects_dir', 'subjects_dir')]), + (inputnode, targets, [('subject_id', 'subject_id')]), + (inputnode, rename_src, [('source_file', 'in_file')]), + (inputnode, itk2lta, [('source_file', 'src_file'), + ('t1w2fsnative_xfm', 'in_file')]), + (get_fsnative, itk2lta, [('brain', 'dst_file')]), # InfantFS: Use brain instead of T1 + (inputnode, sampler, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id')]), + (itersource, targets, [('target', 'space')]), + (itersource, rename_src, [('target', 'subject')]), + (itk2lta, sampler, [('out', 'reg_file')]), + (targets, sampler, [('out', 'target_subject')]), + (rename_src, sampler, [('out_file', 'source_file')]), + (update_metadata, outputnode, [('out_file', 'surfaces')]), + (itersource, outputnode, [('target', 'target')]), + ]) + # fmt: on - # fmt: off - workflow.connect([ - (inputnode, get_fsnative, [('subject_id', 'subject_id'), - ('subjects_dir', 'subjects_dir')]), - (inputnode, targets, [('subject_id', 'subject_id')]), - (inputnode, rename_src, [('source_file', 'in_file')]), - (inputnode, itk2lta, [('source_file', 'src_file'), - ('t1w2fsnative_xfm', 'in_file')]), - (get_fsnative, itk2lta, [('brain', 'dst_file')]), # InfantFS: Use brain instead of T1 - (inputnode, sampler, [('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id')]), - (itersource, targets, [('target', 'space')]), - (itersource, rename_src, [('target', 'subject')]), - (itk2lta, sampler, [('out', 'reg_file')]), - (targets, sampler, [('out', 'target_subject')]), - (rename_src, sampler, [('out_file', 'source_file')]), - (update_metadata, outputnode, [('out_file', 'surfaces')]), - (itersource, outputnode, [('target', 'target')]), - ]) - # fmt: on + if not medial_surface_nan: + workflow.connect(sampler, "out_file", update_metadata, "in_file") + return workflow - if not medial_surface_nan: - workflow.connect(sampler, "out_file", update_metadata, "in_file") + # fmt: off + workflow.connect([ + (inputnode, medial_nans, [('subjects_dir', 'subjects_dir')]), + (sampler, medial_nans, [('out_file', 'in_file')]), + (medial_nans, update_metadata, [('out_file', 'in_file')]), + ]) + # fmt: on return workflow - from niworkflows.interfaces.freesurfer import MedialNaNs + # 0, 1 = wm; 2, 3 = pial; 6, 7 = mid + # note that order of lh / rh within each surf type is not guaranteed due to use + # of unsorted glob by FreeSurferSource prior, but we can do a sort + # to ensure consistent ordering + select_wm = pe.Node( + niu.Select(index=[0, 1]), + name="select_wm_pial", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) - # Refine if medial vertices should be NaNs - medial_nans = pe.MapNode( - MedialNaNs(), iterfield=["in_file"], name="medial_nans", mem_gb=DEFAULT_MEMORY_MIN_GB + select_pial = pe.Node( + niu.Select(index=[2, 3],), + name="select_pial", + mem_gb=DEFAULT_MEMORY_MIN_GB, ) - # fmt: off + select_midthick = pe.Node( + niu.Select(index=[6, 7]), + name="select_midthick", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + create_wm_distvol = pe.MapNode( + CreateSignedDistanceVolume(), + iterfield=["surface"], + name="create_wm_distvol", + mem_gb=mem_gb, + ) + + create_pial_distvol = pe.MapNode( + CreateSignedDistanceVolume(), + iterfield=["surface"], + name="create_pial_distvol", + mem_gb=mem_gb, + ) + + thresh_wm_distvol = pe.MapNode( + fsl.maths.MathsCommand(args="-thr 0 -bin -mul 255"), + iterfield=["in_file"], + name="thresh_wm_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + uthresh_pial_distvol = pe.MapNode( + fsl.maths.MathsCommand(args="-uthr 0 -abs -bin -mul 255"), + iterfield=["in_file"], + name="uthresh_pial_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + bin_wm_distvol = pe.MapNode( + fsl.maths.UnaryMaths(operation="bin"), + iterfield=["in_file"], + name="bin_wm_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + bin_pial_distvol = pe.MapNode( + fsl.maths.UnaryMaths(operation="bin"), + iterfield=["in_file"], + name="bin_pial_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + split_wm_distvol = pe.Node( + niu.Split(splits=[1, 1]), + name="split_wm_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + merge_wm_distvol_no_flatten = pe.Node( + niu.Merge(2), + no_flatten=True, + name="merge_wm_distvol_no_flatten", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + make_ribbon_vol = pe.MapNode( + fsl.maths.MultiImageMaths(op_string="-mas %s -mul 255 "), + iterfield=["in_file", "operand_files"], + name="make_ribbon_vol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + bin_ribbon_vol = pe.MapNode( + fsl.maths.UnaryMaths(operation="bin"), + iterfield=["in_file"], + name="bin_ribbon_vol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + split_squeeze_ribbon_vol = pe.Node( + niu.Split(splits=[1, 1], squeeze=True), + name="split_squeeze_ribbon", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + combine_ribbon_vol_hemis = pe.Node( + fsl.maths.BinaryMaths(operation="add"), + name="combine_ribbon_vol_hemis", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + ribbon_boldsrc_xfm = pe.Node( + ApplyTransforms(interpolation='MultiLabel', + transforms='identity'), + name="ribbon_boldsrc_xfm", + mem_gb=mem_gb, + ) + + stdev_volume = pe.Node( + fsl.maths.StdImage(dimension='T'), + name="stdev_volume", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + mean_volume = pe.Node( + fsl.maths.MeanImage(dimension='T'), + name="mean_volume", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_volume = pe.Node( + fsl.maths.BinaryMaths(operation='div'), + name="cov_volume", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_ribbon = pe.Node( + fsl.ApplyMask(), + name="cov_ribbon", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_ribbon_mean = pe.Node( + fsl.ImageStats(op_string='-M '), + name="cov_ribbon_mean", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_ribbon_std = pe.Node( + fsl.ImageStats(op_string='-S '), + name="cov_ribbon_std", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_ribbon_norm = pe.Node( + fsl.maths.BinaryMaths(operation='div'), + name="cov_ribbon_norm", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + smooth_norm = pe.Node( + fsl.maths.MathsCommand(args="-bin -s 5 "), + name="smooth_norm", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + merge_smooth_norm = pe.Node( + niu.Merge(1), + name="merge_smooth_norm", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_ribbon_norm_smooth = pe.Node( + fsl.maths.MultiImageMaths(op_string='-s 5 -div %s -dilD '), + name="cov_ribbon_norm_smooth", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_norm = pe.Node( + fsl.maths.BinaryMaths(operation='div'), + name="cov_norm", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_norm_modulate = pe.Node( + fsl.maths.BinaryMaths(operation='div'), + name="cov_norm_modulate", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + cov_norm_modulate_ribbon = pe.Node( + fsl.ApplyMask(), + name="cov_norm_modulate_ribbon", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + def _calc_upper_thr(in_stats): + return in_stats[0] + (in_stats[1] * 0.5) + + upper_thr_val = pe.Node( + Function( + input_names=["in_stats"], + output_names=["upper_thresh"], + function=_calc_upper_thr + ), + name="upper_thr_val", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + def _calc_lower_thr(in_stats): + return in_stats[1] - (in_stats[0] * 0.5) + + lower_thr_val = pe.Node( + Function( + input_names=["in_stats"], + output_names=["lower_thresh"], + function=_calc_lower_thr + ), + name="lower_thr_val", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + mod_ribbon_mean = pe.Node( + fsl.ImageStats(op_string='-M '), + name="mod_ribbon_mean", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + mod_ribbon_std = pe.Node( + fsl.ImageStats(op_string='-S '), + name="mod_ribbon_std", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + merge_mod_ribbon_stats = pe.Node( + niu.Merge(2), + name="merge_mod_ribbon_stats", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + bin_mean_volume = pe.Node( + fsl.maths.UnaryMaths(operation="bin"), + name="bin_mean_volume", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + merge_goodvoxels_operands = pe.Node( + niu.Merge(2), + name="merge_goodvoxels_operands", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + goodvoxels_thr = pe.Node( + fsl.maths.Threshold(), + name="goodvoxels_thr", + mem_gb=mem_gb, + ) + + goodvoxels_mask = pe.Node( + fsl.maths.MultiImageMaths(op_string='-bin -sub %s -mul -1 '), + name="goodvoxels_mask", + mem_gb=mem_gb, + ) + + goodvoxels_ribbon_mask = pe.Node( + fsl.ApplyMask(), + name_source=['in_file'], + keep_extension=True, + name_template='%s', + name="goodvoxels_ribbon_mask", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + apply_goodvoxels_ribbon_mask = pe.Node( + fsl.ApplyMask(), + name_source=['in_file'], + keep_extension=True, + name_template='%s', + name="apply_goodvoxels_ribbon_mask", + mem_gb=mem_gb * 3, + ) + + get_target_wm = pe.MapNode( + FreeSurferSource(), + iterfield=["in_file", "surface"], + name="get_target_wm", + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + get_target_wm.inputs.hemi=["lh","rh"] + + make_target_midthick = pe.MapNode( + MakeMidthickness(thickness=True, distance=0.5), + iterfield=["in_file"], + name="make_target_midthick", + run_without_submitting=True, + mem_gb=mem_gb * 3, + ) + + target_midthick_gifti = pe.MapNode( + fs.MRIsConvert(out_datatype="gii"), + iterfield=["in_file"], + name="target_midthick_gifti", + run_without_submitting=True, + mem_gb=mem_gb, + ) + + wbsampler = pe.MapNode( + VolumeToSurfaceMapping( + mapping_method="ribbon-constrained", + ), + iterfield=["in_file"], + name="wbsampler", + mem_gb=mem_gb * 3, + ) + + metric_dilate = pe.MapNode( + MetricDilate( + distance=10, + nearest=True, + ), + iterfield=["in_file", "surface"], + name="metric_dilate", + mem_gb=mem_gb * 3, + ) + + native_to_target = pe.MapNode( + MetricResample( + method="ADAP_BARY_AREA", + area_metrics=True, + ), + iterfield=[ + "in_file", + "out_file", + "new_sphere", + "new_area", + "current_sphere", + "current_area", + ], + name="native_to_target", + ) + native_to_target.inputs.new_sphere = [ + str( + tf.api.get("fsaverage", hemi=hemi, density="164k", desc="std", suffix="sphere") + ) + for hemi in "LR" + ] + + # make FS midthick if target doesn't have them already workflow.connect([ - (inputnode, medial_nans, [('subjects_dir', 'subjects_dir')]), - (sampler, medial_nans, [('out_file', 'in_file')]), - (medial_nans, update_metadata, [('out_file', 'in_file')]), + (inputnode, get_target_wm, [('subjects_dir', 'subjects_dir')]), + (targets, get_target_wm, [('out', 'subject_id')]), + (get_target_wm, make_target_midthick, [("white", "in_file")]), + (make_target_midthick, target_midthick_gifti, [("out_file", "in_file")]), ]) - # fmt: on - return workflow + + # make HCP-style ribbon volume in T1w space + workflow.connect([ + (inputnode, select_wm, [("anat_giftis", "inlist")]), + (inputnode, select_pial, [("anat_giftis", "inlist")]), + (inputnode, select_midthick, [("anat_giftis", "inlist")]), + (select_wm, create_wm_distvol, [("out", "surface")]), + (inputnode, create_wm_distvol, [("t1w_mask", "ref_space")]), + (select_pial, create_pial_distvol, [("out", "surface")]), + (inputnode, create_pial_distvol, [("t1w_mask", "ref_space")]), + (create_wm_distvol, thresh_wm_distvol, [("out_vol", "in_file")]), + (create_pial_distvol, uthresh_pial_distvol, [("out_vol", "in_file")]), + (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), + (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), + (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), + (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), + (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), + (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), + (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), + ]) + + # make HCP-style "goodvoxels" mask in t1w space for filtering outlier voxels + # in bold timeseries, based on modulated normalized covariance + workflow.connect([ + (combine_ribbon_vol_hemis, ribbon_boldsrc_xfm, [("out_file", 'input_image')]), + (rename_src, stdev_volume, [("out_file", "in_file")]), + (rename_src, mean_volume, [("out_file", "in_file")]), + (mean_volume, ribbon_boldsrc_xfm, [('out_file', 'reference_image')]), + (stdev_volume, cov_volume, [("out_file", "in_file")]), + (mean_volume, cov_volume, [("out_file", "operand_file")]), + (cov_volume, cov_ribbon, [("out_file", "in_file")]), + (ribbon_boldsrc_xfm, cov_ribbon, [("output_image", "mask_file")]), + (cov_ribbon, cov_ribbon_mean, [("out_file", "in_file")]), + (cov_ribbon, cov_ribbon_std, [("out_file", "in_file")]), + (cov_ribbon, cov_ribbon_norm, [("out_file", "in_file")]), + (cov_ribbon_mean, cov_ribbon_norm, [("out_stat", "operand_value")]), + (cov_ribbon_norm, smooth_norm, [("out_file", "in_file")]), + (smooth_norm, merge_smooth_norm, [("out_file", "in1")]), + (cov_ribbon_norm, cov_ribbon_norm_smooth, [("out_file", "in_file")]), + (merge_smooth_norm, cov_ribbon_norm_smooth, [("out", "operand_files")]), + (cov_ribbon_mean, cov_norm, [("out_stat", "operand_value")]), + (cov_volume, cov_norm, [("out_file", "in_file")]), + (cov_norm, cov_norm_modulate, [("out_file", "in_file")]), + (cov_ribbon_norm_smooth, cov_norm_modulate, [("out_file", "operand_file")]), + (cov_norm_modulate, cov_norm_modulate_ribbon, [("out_file", "in_file")]), + (ribbon_boldsrc_xfm, cov_norm_modulate_ribbon, [("output_image", "mask_file")]), + (cov_norm_modulate_ribbon, mod_ribbon_mean, [("out_file", "in_file")]), + (cov_norm_modulate_ribbon, mod_ribbon_std, [("out_file", "in_file")]), + (mod_ribbon_mean, merge_mod_ribbon_stats, [("out_stat", "in1")]), + (mod_ribbon_std, merge_mod_ribbon_stats, [("out_stat", "in2")]), + (merge_mod_ribbon_stats, upper_thr_val, [("out", "in_stats")]), + (merge_mod_ribbon_stats, lower_thr_val, [("out", "in_stats")]), + (mean_volume, bin_mean_volume, [("out_file", "in_file")]), + (upper_thr_val, goodvoxels_thr, [("upper_thresh", "thresh")]), + (cov_norm_modulate, goodvoxels_thr, [("out_file", "in_file")]), + (bin_mean_volume, merge_goodvoxels_operands, [("out_file", "in1")]), + (goodvoxels_thr, goodvoxels_mask, [("out_file", "in_file")]), + (merge_goodvoxels_operands, goodvoxels_mask, [("out", "operand_files")]), + ]) + + if surface_sampler is "fs": + # apply goodvoxels ribbon mask to bold + workflow.connect([ + (goodvoxels_mask, goodvoxels_ribbon_mask, [("out_file", "in_file")]), + (ribbon_boldsrc_xfm, goodvoxels_ribbon_mask, [("output_image", "mask_file")]), + (goodvoxels_ribbon_mask, apply_goodvoxels_ribbon_mask, [("out_file", "mask_file")]), + (rename_src, apply_goodvoxels_ribbon_mask, [("out_file", "in_file")]), + ]) + + # project masked bold to target surfs + workflow.connect([ + (inputnode, get_fsnative, [("subject_id", "subject_id"), + ("subjects_dir", "subjects_dir")]), + (inputnode, targets, [("subject_id", "subject_id")]), + (inputnode, rename_src, [("source_file", "in_file")]), + (inputnode, itk2lta, [("source_file", "src_file"), + ("t1w2fsnative_xfm", "in_file")]), + (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 + (inputnode, sampler, [("subjects_dir", "subjects_dir"), + ("subject_id", "subject_id")]), + (itersource, targets, [("target", "space")]), + (itersource, rename_src, [("target", "subject")]), + (itk2lta, sampler, [("out", "reg_file")]), + (targets, sampler, [("out", "target_subject")]), + (apply_goodvoxels_ribbon_mask, sampler, [("out_file", "source_file")]), + (sampler, metric_dilate, [("out_file", "in_file")]), + (target_midthick_gifti, metric_dilate, [("converted", "surface")]), + (update_metadata, outputnode, [("out_file", "surfaces")]), + (itersource, outputnode, [("target", "target")]), + ]) + + # fmt:on + if not medial_surface_nan: + # fmt:off + workflow.connect([ + (metric_dilate, update_metadata, [("out_file", "in_file")]), + ]) + # fmt:on + return workflow + + # fmt:off + workflow.connect([ + (metric_dilate, medial_nans, [("out_file", "in_file")]), + (medial_nans, update_metadata, [("out_file", "in_file")]), + ]) + # fmt:on + return workflow + # wb method first projects to native surfs, then resamples to fsaverage + # (fsaverage{3,4,5,6} not supported as yet) + if surface_sampler is "wb": + workflow.connect([ + (inputnode, get_fsnative, [("subject_id", "subject_id"), + ("subjects_dir", "subjects_dir")]), + (inputnode, targets, [("subject_id", "subject_id")]), + (inputnode, rename_src, [("source_file", "in_file")]), + (inputnode, itk2lta, [("source_file", "src_file"), + ("t1w2fsnative_xfm", "in_file")]), + (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 + (itersource, targets, [("target", "space")]), + (itersource, rename_src, [("target", "subject")]), + (rename_src, wbsampler, [("out_file", "in_file")]), + (select_midthick, wbsampler, [(("out_file", _sorted), "surface")]), + (select_wm, wbsampler, [(("out_file", _sorted), "inner_surf")]), + (select_pial, wbsampler, [(("out_file", _sorted), "outer_surf")]), + (goodvoxels_mask, wbsampler, [("out_file", "roi_volume")]), + (wbsampler, metric_dilate, [("out_file", "in_file")]), + (metric_dilate, native_to_target, [("out_file", "in_file")]), + (select_midthick, native_to_target, [(("out_file", _sorted), "current_area")]), + (get_fsnative, native_to_target, [(("sphere_reg", _sorted), "current_sphere")]), + (target_midthick_gifti, native_to_target, [("converted", "new_area")]), + (update_metadata, outputnode, [("out_file", "surfaces")]), + (itersource, outputnode, [("target", "target")]), + ]) + + # fmt:on + if not medial_surface_nan: + # fmt:off + workflow.connect([ + (sampler, metric_dilate, [("out_file", "in_file")]), + (target_midthick_gifti, metric_dilate, [("converted", "surface")]), + (metric_dilate, update_metadata, [("out_file", "in_file")]), + ]) + # fmt:on + return workflow + + # fmt:off + workflow.connect([ + (sampler, metric_dilate, [("out_file", "in_file")]), + (target_midthick_gifti, metric_dilate, [("converted", "surface")]), + (metric_dilate, medial_nans, [("out_file", "in_file")]), + (medial_nans, update_metadata, [("out_file", "in_file")]), + ]) + # fmt:on + return workflow + def init_bold_std_trans_wf( freesurfer, mem_gb, @@ -620,7 +1187,7 @@ def init_bold_grayords_wf(grayord_density, mem_gb, repetition_time, name="bold_g subcortical_labels : :obj:`str` Volume file containing all subcortical labels surf_files : :obj:`str` - List of BOLD files resampled on the fsaverage (ico7) surfaces. + List of BOLD files resampled on the fsaverage (ico7) surfaces surf_refs : List of unique identifiers corresponding to the BOLD surface-conversions. @@ -640,7 +1207,6 @@ def init_bold_grayords_wf(grayord_density, mem_gb, repetition_time, name="bold_g from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.utility import KeySelect - from ...interfaces.nibabel import ReorientImage from ...interfaces.workbench import CiftiCreateDenseTimeseries workflow = Workflow(name=name) @@ -732,16 +1298,9 @@ def init_bold_grayords_wf(grayord_density, mem_gb, repetition_time, name="bold_g niu.Function(function=_split_surfaces, output_names=["left_surface", "right_surface"]), name="split_surfaces", ) - - reorient_data = pe.Node(ReorientImage(target_orientation="LAS"), name="reorient_data") - reorient_labels = reorient_data.clone(name="reorient_labels") - gen_cifti = pe.Node(CiftiCreateDenseTimeseries(timestep=repetition_time), name="gen_cifti") - gen_cifti.inputs.roi_left = tf.api.get( - "fsLR", density=fslr_density, hemi="L", desc="nomedialwall", suffix="dparc" - ) - gen_cifti.inputs.roi_right = tf.api.get( - "fsLR", density=fslr_density, hemi="R", desc="nomedialwall", suffix="dparc" + gen_cifti.inputs.volume_structure_labels = str( + tf.api.get("MNI152NLin6Asym", resolution=2, atlas="HCP", suffix="dseg") ) gen_cifti_metadata = pe.Node( niu.Function(function=_gen_metadata, output_names=["out_metadata", "variant", "density"]), @@ -751,10 +1310,9 @@ def init_bold_grayords_wf(grayord_density, mem_gb, repetition_time, name="bold_g # fmt: off workflow.connect([ - (inputnode, reorient_data, [("subcortical_volume", "in_file")]), - (inputnode, reorient_labels, [("subcortical_labels", "in_file")]), - (reorient_data, gen_cifti, [("out_file", "volume_data")]), - (reorient_labels, gen_cifti, [("out_file", "volume_structure_labels")]), + (inputnode, gen_cifti, [ + ('subcortical_volume', 'volume_data'), + ('subcortical_labels', 'volume_structure_labels')]), (inputnode, select_fs_surf, [('surf_files', 'surf_files'), ('surf_refs', 'keys')]), (select_fs_surf, resample, [('surf_files', 'in_file')]), @@ -862,3 +1420,6 @@ def _itk2lta(in_file, src_file, dst_file): in_file, fmt="fs" if in_file.endswith(".lta") else "itk", reference=src_file ).to_filename(out_file, moving=dst_file, fmt="fs") return str(out_file) + +def _sorted(inlist): + return sorted(inlist) From b14cf08cab866d256709a82ba1086c124ca49fc6 Mon Sep 17 00:00:00 2001 From: madisoth <56737848+madisoth@users.noreply.github.com> Date: Wed, 14 Sep 2022 13:50:21 -0700 Subject: [PATCH 02/11] additional work on --project-goodvoxels, revert work on --surface-sampler option (to be continued in a separate branch) --- nibabies/cli/parser.py | 9 - nibabies/config.py | 3 - nibabies/data/tests/config.toml | 1 - nibabies/interfaces/metric.py | 305 ------------------------- nibabies/interfaces/wbvoltosurf.py | 298 ------------------------ nibabies/workflows/base.py | 1 - nibabies/workflows/bold/base.py | 9 +- nibabies/workflows/bold/resampling.py | 317 ++++++++------------------ 8 files changed, 98 insertions(+), 845 deletions(-) delete mode 100644 nibabies/interfaces/metric.py delete mode 100644 nibabies/interfaces/wbvoltosurf.py diff --git a/nibabies/cli/parser.py b/nibabies/cli/parser.py index 6be65545..02ecfea3 100644 --- a/nibabies/cli/parser.py +++ b/nibabies/cli/parser.py @@ -342,15 +342,6 @@ def _slice_time_ref(value, parser): "from surface resampling. Only performed for GIFTI files mapped to a freesurfer subject " "(fsaverage or fsnative).", ) - g_conf.add_argument( - "--surface-sampler", - action="store", - choices=("fs", "wb"), - default="fs", - help="select sampling method used for projection of BOLD timeseries to surface. " - "'fs' to use FreeSurfer mri_vol2surf (default). " - "'wb' to use Workbench wb_command -volume-to-surface-mapping.", - ) g_conf.add_argument( "--slice-time-ref", required=False, diff --git a/nibabies/config.py b/nibabies/config.py index 333f0d23..9f459bf5 100644 --- a/nibabies/config.py +++ b/nibabies/config.py @@ -555,9 +555,6 @@ class workflow(_Config): """Fill medial surface with :abbr:`NaNs (not-a-number)` when sampling.""" project_goodvoxels = False """Exclude voxels with locally high coefficient of variation from sampling.""" - surface_sampler = "fs" - """``fs`` (default) or ``wb`` to select FreeSurfer-based or Workbench-based - method for sampling BOLD to surface.""" regressors_all_comps = None """Return all CompCor components.""" regressors_dvars_th = None diff --git a/nibabies/data/tests/config.toml b/nibabies/data/tests/config.toml index 6cb874b1..b9bd3e16 100644 --- a/nibabies/data/tests/config.toml +++ b/nibabies/data/tests/config.toml @@ -53,7 +53,6 @@ ignore = [ "slicetiming",] longitudinal = false medial_surface_nan = false project_goodvoxels = false -surface_sampler = "fs" regressors_all_comps = false regressors_dvars_th = 1.5 regressors_fd_th = 0.5 diff --git a/nibabies/interfaces/metric.py b/nibabies/interfaces/metric.py deleted file mode 100644 index 023fda7e..00000000 --- a/nibabies/interfaces/metric.py +++ /dev/null @@ -1,305 +0,0 @@ -# -*- coding: utf-8 -*- -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -"""This module provides interfaces for workbench surface commands""" -import os - -from nipype.interfaces.base import TraitedSpec, File, traits, CommandLineInputSpec -from nipype.interfaces.workbench.base import WBCommand -from nipype import logging - -iflogger = logging.getLogger("nipype.interface") - - -class MetricDilateInputSpec(CommandLineInputSpec): - - in_file = File( - exists=True, - mandatory=True, - argstr="%s ", - position=0, - desc="the metric to dilate", - ) - - surface = File( - exists=True, - mandatory=True, - argstr="%s ", - position=1, - desc="the surface to compute on", - ) - - distance = traits.Float( - mandatory=True, - argstr="%f ", - position=2, - desc="distance in mm to dilate", - ) - - out_file = File( - name_source=["in_file"], - name_template="%s.func.gii", - keep_extension=False, - argstr="%s ", - position=3, - desc="output - the output metric", - ) - - bad_vertex_roi = File( - argstr="-bad-vertex-roi %s ", - position=4, - desc="metric file, positive values denote vertices to have their values replaced", - ) - - data_roi = File( - argstr="-data-roi %s ", - position=5, - desc="metric file, positive values denote vertices that have data", - ) - - column = traits.Int( - position=6, - argstr="-column %d ", - desc="the column number", - ) - - nearest = traits.Bool( - position=7, - argstr="-nearest ", - desc="use the nearest good value instead of a weighted average", - ) - - linear = traits.Bool( - position=8, - argstr="-linear ", - desc="fill in values with linear interpolation along strongest gradient", - ) - - exponent = traits.Float( - argstr="-exponent %f ", - position=9, - default=6.0, - desc="exponent n to use in (area / (distance ^ n)) as the " - "weighting function (default 6)", - ) - - corrected_areas = File( - argstr="-corrected-areas %s ", - position=10, - desc="vertex areas to use instead of computing them from the surface", - ) - - legacy_cutoff = traits.Bool( - position=11, - argstr="-legacy-cutoff ", - desc="use the v1.3.2 method of choosing how many vertices to " - "use when calulating the dilated value with weighted method", - ) - - -class MetricDilateOutputSpec(TraitedSpec): - out_file = File(exists=True, desc="output file") - - -class MetricDilate(WBCommand): - """ - For all data values designated as bad, if they neighbor a good value or - are within the specified distance of a good value in the same kind of - model, replace the value with a distance weighted average of nearby good - values, otherwise set the value to zero. If -nearest is specified, it - will use the value from the closest good value within range instead of a - weighted average. When the input file contains label data, nearest - dilation is used on the surface, and weighted popularity is used in the - volume. - - The -corrected-areas options are intended for dilating on group average - surfaces, but it is only an approximate correction for the reduction of - structure in a group average surface. - - If -bad-vertex-roi is specified, all values, including those with - value zero, are good, except for locations with a positive value in the - ROI. If it is not specified, only values equal to zero are bad. - - >>> from nipype.interfaces.workbench import MetricDilate - >>> metdil = MetricDilate() - >>> metdil.inputs.in_file = 'sub-01_task-rest_bold_space-fsaverage5.L.func.gii' - >>> metdil.inputs.surface = 'sub-01.L.midthickness.surf.gii' - >>> metdil.inputs.distance = 10 - >>> metdil.inputs.nearest = True - >>> metdil.cmdline - 'wb_command -metric-dilate \ - sub-01_task-rest_bold_space-fsaverage5.L.func.gii \ - sub-01.L.midthickness.surf.gii \ - 10 \ - sub-01_task-rest_bold_space-fsaverage5.L.func.gii \ - -nearest' - """ - input_spec = MetricDilateInputSpec - output_spec = MetricDilateOutputSpec - _cmd = "wb_command -metric-dilate " - - -class MetricResampleInputSpec(CommandLineInputSpec): - in_file = File( - exists=True, - mandatory=True, - argstr="%s", - position=0, - desc="The metric file to resample", - ) - current_sphere = File( - exists=True, - mandatory=True, - argstr="%s", - position=1, - desc="A sphere surface with the mesh that the metric is currently on", - ) - new_sphere = File( - exists=True, - mandatory=True, - argstr="%s", - position=2, - desc="A sphere surface that is in register with and" - " has the desired output mesh", - ) - method = traits.Enum( - "ADAP_BARY_AREA", - "BARYCENTRIC", - argstr="%s", - mandatory=True, - position=3, - desc="The method name - ADAP_BARY_AREA method is recommended for" - " ordinary metric data, because it should use all data while" - " downsampling, unlike BARYCENTRIC. If ADAP_BARY_AREA is used," - " exactly one of area_surfs or area_metrics must be specified", - ) - out_file = File( - name_source=["new_sphere"], - name_template="%s.out", - keep_extension=True, - argstr="%s", - position=4, - desc="The output metric", - ) - area_surfs = traits.Bool( - position=5, - argstr="-area-surfs", - xor=["area_metrics"], - desc="Specify surfaces to do vertex area correction based on", - ) - area_metrics = traits.Bool( - position=5, - argstr="-area-metrics", - xor=["area_surfs"], - desc="Specify vertex area metrics to do area correction based on", - ) - current_area = File( - exists=True, - position=6, - argstr="%s", - desc="A relevant anatomical surface with mesh OR" - " a metric file with vertex areas for mesh", - ) - new_area = File( - exists=True, - position=7, - argstr="%s", - desc="A relevant anatomical surface with mesh OR" - " a metric file with vertex areas for mesh", - ) - roi_metric = File( - exists=True, - position=8, - argstr="-current-roi %s", - desc="Input roi on the current mesh used to exclude non-data vertices", - ) - valid_roi_out = traits.Bool( - position=9, - argstr="-valid-roi-out", - desc="Output the ROI of vertices that got data from valid source vertices", - ) - largest = traits.Bool( - position=10, - argstr="-largest", - desc="Use only the value of the vertex with the largest weight", - ) - - -class MetricResampleOutputSpec(TraitedSpec): - out_file = File(exists=True, desc="the output metric") - roi_file = File(desc="ROI of vertices that got data from valid source vertices") - - -class MetricResample(WBCommand): - """ - Resample a metric file to a different mesh - - Resamples a metric file, given two spherical surfaces that are in - register. If ``ADAP_BARY_AREA`` is used, exactly one of -area-surfs or - ``-area-metrics`` must be specified. - - The ``ADAP_BARY_AREA`` method is recommended for ordinary metric data, - because it should use all data while downsampling, unlike ``BARYCENTRIC``. - The recommended areas option for most data is individual midthicknesses - for individual data, and averaged vertex area metrics from individual - midthicknesses for group average data. - - The ``-current-roi`` option only masks the input, the output may be slightly - dilated in comparison, consider using ``-metric-mask`` on the output when - using ``-current-roi``. - - The ``-largest option`` results in nearest vertex behavior when used with - ``BARYCENTRIC``. When resampling a binary metric, consider thresholding at - 0.5 after resampling rather than using ``-largest``. - - >>> from nipype.interfaces.workbench import MetricResample - >>> metres = MetricResample() - >>> metres.inputs.in_file = 'sub-01_task-rest_bold_space-fsaverage5.L.func.gii' - >>> metres.inputs.method = 'ADAP_BARY_AREA' - >>> metres.inputs.current_sphere = 'fsaverage5_std_sphere.L.10k_fsavg_L.surf.gii' - >>> metres.inputs.new_sphere = 'fs_LR-deformed_to-fsaverage.L.sphere.32k_fs_LR.surf.gii' - >>> metres.inputs.area_metrics = True - >>> metres.inputs.current_area = 'fsaverage5.L.midthickness_va_avg.10k_fsavg_L.shape.gii' - >>> metres.inputs.new_area = 'fs_LR.L.midthickness_va_avg.32k_fs_LR.shape.gii' - >>> metres.cmdline - 'wb_command -metric-resample sub-01_task-rest_bold_space-fsaverage5.L.func.gii \ - fsaverage5_std_sphere.L.10k_fsavg_L.surf.gii \ - fs_LR-deformed_to-fsaverage.L.sphere.32k_fs_LR.surf.gii \ - ADAP_BARY_AREA fs_LR-deformed_to-fsaverage.L.sphere.32k_fs_LR.surf.out \ - -area-metrics fsaverage5.L.midthickness_va_avg.10k_fsavg_L.shape.gii \ - fs_LR.L.midthickness_va_avg.32k_fs_LR.shape.gii' - """ - - input_spec = MetricResampleInputSpec - output_spec = MetricResampleOutputSpec - _cmd = "wb_command -metric-resample" - - def _format_arg(self, opt, spec, val): - if opt in ["current_area", "new_area"]: - if not self.inputs.area_surfs and not self.inputs.area_metrics: - raise ValueError( - "{} was set but neither area_surfs or" - " area_metrics were set".format(opt) - ) - if opt == "method": - if ( - val == "ADAP_BARY_AREA" - and not self.inputs.area_surfs - and not self.inputs.area_metrics - ): - raise ValueError( - "Exactly one of area_surfs or area_metrics" " must be specified" - ) - if opt == "valid_roi_out" and val: - # generate a filename and add it to argstr - roi_out = self._gen_filename(self.inputs.in_file, suffix="_roi") - iflogger.info("Setting roi output file as", roi_out) - spec.argstr += " " + roi_out - return super(MetricResample, self)._format_arg(opt, spec, val) - - def _list_outputs(self): - outputs = super(MetricResample, self)._list_outputs() - if self.inputs.valid_roi_out: - roi_file = self._gen_filename(self.inputs.in_file, suffix="_roi") - outputs["roi_file"] = os.path.abspath(roi_file) - return outputs diff --git a/nibabies/interfaces/wbvoltosurf.py b/nibabies/interfaces/wbvoltosurf.py deleted file mode 100644 index 12409b74..00000000 --- a/nibabies/interfaces/wbvoltosurf.py +++ /dev/null @@ -1,298 +0,0 @@ -# -*- coding: utf-8 -*- -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -"""This module provides interfaces for workbench surface commands""" -import os - -from nipype.interfaces.base import TraitedSpec, File, traits, CommandLineInputSpec -from nipype.interfaces.workbench.base import WBCommand -from nipype import logging - -iflogger = logging.getLogger("nipype.interface") - - -class VolumeToSurfaceMappingInputSpec(CommandLineInputSpec): - - in_file = File( - exists=True, - mandatory=True, - argstr="%s ", - position=0, - desc="the volume to map data from", - ) - - surface = File( - exists=True, - mandatory=True, - argstr="%s ", - position=1, - desc="the surface to map the data onto", - ) - - out_file = File( - name_source=["in_file"], - name_template="%s.func.gii", - keep_extension=False, - argstr="%s ", - position=2, - desc="output - the output metric file", - ) - - mapping_method = traits.Enum( - "trilinear", - "enclosing", - "cubic", - "ribbon-constrained", - "myelin-style", - position=3, - argstr="-%s ", - desc="choose mapping method: trilinear, enclosing voxel, cubic spline, " - "ribbon-constrained, or myelin-style", - ) - - inner_surf = File( - argstr="%s ", - position=4, - desc="the inner surface of the ribbon", - ) - - outer_surf = File( - argstr="%s ", - position=5, - desc="the outer surface of the ribbon", - ) - - roi_volume = File( - argstr="-volume-roi %s ", - position=6, - desc="the roi volume file to use", - ) - - weighted = traits.Bool( - position=7, - argstr="-weighted ", - desc="treat the roi values as weightings rather than binary", - ) - - subdiv_num = traits.Int( - position=8, - argstr="-voxel-subdiv %d ", - desc="number of subdivisions while estimating voxel weights, " - "default 3 if this option is not specified", - ) - - thin_columns = traits.Bool( - position=9, - argstr="-thin-columns ", - desc="use non-overlapping polyhedra", - ) - - gaussian_scale = traits.Float( - argstr="-gaussian %f ", - position=10, - desc="value to multiply the local thickness by, to get the gaussian sigma. " - "reduce weight to voxels that aren't near ", - ) - - interpolate_method = traits.Enum( - "CUBIC", - "ENCLOSING_VOXEL", - "TRILINEAR", - position=11, - argstr="-interpolate %s ", - desc="must be CUBIC, ENCLOSING_VOXEL, or TRILINEAR. " - "instead of a weighted average of voxels, " - "interpolate at subpoints inside the ribbon", - ) - - bad_vertices_out = File( - argstr="-bad-vertices-out %s ", - position=12, - desc="output an ROI of which vertices didn't intersect any valid voxels", - ) - - output_weights = traits.Bool( - position=13, - argstr="-output-weights ", - desc="the column number", - ) - - output_weights_vertex = traits.Int( - position=14, - argstr="%d ", - desc="the column number", - ) - - output_weights_weights_out = File( - argstr="%s ", - position=15, - desc="output - volume to write the weights to", - ) - - output_weights_text_out = File( - argstr="%s ", - position=16, - desc="text file name. " - "write the voxel weights for all vertices to a text file", - ) - - myelin_style_ribbon_roi = File( - argstr="%s ", - position=17, - desc="roi volume of the cortical ribbon for this hemisphere", - ) - - myelin_style_thickness = File( - argstr="%s ", - position=18, - desc="a metric file of cortical thickness", - ) - - myelin_style_sigma = traits.Float( - argstr="%f ", - position=19, - desc="gaussian kernel in mm for weighting voxels within range", - ) - - legacy_bug = traits.Bool( - position=20, - argstr="-legacy-bug ", - desc="emulate old Workbench v1.2.3 and earlier code that didn't follow a cylinder cutoff", - ) - - subvol_select = File( - argstr="-subvol-select %s ", - position=21, - desc="the subvolume number or name", - ) - - -class VolumeToSurfaceMappingOutputSpec(TraitedSpec): - out_file = File(exists=True, desc="the output metric file") - bad_vertices_out = File(desc="ROI of which vertices didn't intersect any valid voxels") - output_weights_weights_out = File(desc="volume to which the specified vertex's weights are written") - output_weights_text_out = File(desc="text file of the voxel weights for all vertices") - - -class VolumeToSurfaceMapping(WBCommand): - """ - You must specify exactly one mapping method. Enclosing voxel uses the - value from the voxel the vertex lies inside, while trilinear does a 3D - linear interpolation based on the voxels immediately on each side of the - vertex's position. - - The ribbon mapping method constructs a polyhedron from the vertex's - neighbors on each surface, and estimates the amount of this polyhedron's - volume that falls inside any nearby voxels, to use as the weights for - sampling. If -thin-columns is specified, the polyhedron uses the edge - midpoints and triangle centroids, so that neighboring vertices do not - have overlapping polyhedra. This may require increasing -voxel-subdiv to - get enough samples in each voxel to reliably land inside these smaller - polyhedra. The volume ROI is useful to exclude partial volume effects of - voxels the surfaces pass through, and will cause the mapping to ignore - voxels that don't have a positive value in the mask. The subdivision - number specifies how it approximates the amount of the volume the - polyhedron intersects, by splitting each voxel into NxNxN pieces, and - checking whether the center of each piece is inside the polyhedron. If - you have very large voxels, consider increasing this if you get zeros in - your output. The -gaussian option makes it act more like the myelin - method, where the distance of a voxel from is used to - downweight the voxel. The -interpolate suboption, instead of doing a - weighted average of voxels, interpolates from the volume at the - subdivided points inside the ribbon. If using both -interpolate and the - -weighted suboption to -volume-roi, the roi volume weights are linearly - interpolated, unless the -interpolate method is ENCLOSING_VOXEL, in which - case ENCLOSING_VOXEL is also used for sampling the roi volume weights. - - The myelin style method uses part of the caret5 myelin mapping command to - do the mapping: for each surface vertex, take all voxels that are in a - cylinder with radius and height equal to cortical thickness, centered on - the vertex and aligned with the surface normal, and that are also within - the ribbon ROI, and apply a gaussian kernel with the specified sigma to - them to get the weights to use. The -legacy-bug flag reverts to the - unintended behavior present from the initial implementation up to and - including v1.2.3, which had only the tangential cutoff and a bounding box - intended to be larger than where the cylinder cutoff should have been. - - >>> from nipype.interfaces.workbench import VolumeToSurfaceMapping - >>> vtsm = VolumeToSurfaceMapping() - >>> vtsm.inputs.in_file = 'sub-01_task-rest_run-1_bold.nii.gz' - >>> vtsm.inputs.surface = 'sub-01.L.midthickness.surf.gii' - >>> vtsm.inputs.ribbon_constrained = True - >>> vtsm.inputs.inner_surf = 'sub-01.L.white.surf.gii' - >>> vtsm.inputs.outer_surf = 'sub-01.L.pial.surf.gii' - >>> vtsm.inputs.roi_volume = 'sub-01_task-rest_run-1_goodvoxels.nii.gz' - >>> vtsm.cmdline - 'wb_command -volume-to-surface-mapping \ - sub-01_task-rest_run-1_bold.nii.gz \ - sub-01.L.midthickness.surf.gii \ - sub-01_task-rest_run-1_bold.func.gii \ - -ribbon-constrained sub-01.L.white.surf.gii sub-01.L.pial.surf.gii \ - -volume-roi sub-01_task-rest_run-1_goodvoxels.nii.gz - """ - input_spec = VolumeToSurfaceMappingInputSpec - output_spec = VolumeToSurfaceMappingOutputSpec - _cmd = "wb_command -volume-to-surface-mapping " - - def _format_arg(self, opt, spec, val): - if opt in [ - "inner_surf", - "outer_surf", - "roi_volume", - "subdiv_num", - "thin_columns", - "gaussian_scale", - "interpolate_method", - "bad_vertices_out", - "output_weights", - "output_weights_text_out", - ]: - if self.inputs.mapping_method.val is not "ribbon-constrained": - raise ValueError( - "{} was set, but ribbon-constrained mapping method was not".format(opt) - ) - - if opt == "weighted": - if (val == True) and (self.inputs.mapping_method.val is not "ribbon-constrained"): - raise ValueError( - "weighted was set, but the required roi_volume was not" - ) - - if ((self.inputs.mapping_method.val == "ribbon-constrained") - and not (self.inputs.inner_surf and self.inputs.outer_surf) - ): - raise ValueError( - "mapping method is ribbon-constrained but at least one of " - "inner_surf, outer_surf was not set" - ) - - if opt == "output_weights": - if ((val == True) and not (self.inputs.output_weights_vertex - and self.inputs.output_weights_weights_out) - ): - raise ValueError( - "output_weights was set but at least one of " - "output_weights_vertex, output_weights_weights_out was not set" - ) - - if ((self.inputs.mapping_method.val == "myelin-style") - and not (self.inputs.myelin_style_ribbon_roi - and self.inputs.myelin_style_thickness - and self.inputs.myelin_style_sigma) - ): - raise ValueError( - "mapping method is myelin-style but at least one of myelin_style_ribbon_roi, " - "myelin_style_thickness, myelin_style_sigma was not set" - ) - - if opt == "legacy_bug": - if (val == True) and (not self.inputs.mapping_method.val is not "myelin-style"): - raise ValueError( - "legacy_bug was set, but the mapping method is not myelin-style" - ) - - return super(VolumeToSurfaceMapping, self)._format_arg(opt, spec, val) - - def _list_outputs(self): - outputs = super(VolumeToSurfaceMapping, self)._list_outputs() - return outputs \ No newline at end of file diff --git a/nibabies/workflows/base.py b/nibabies/workflows/base.py index a4d16864..5ac91064 100644 --- a/nibabies/workflows/base.py +++ b/nibabies/workflows/base.py @@ -198,7 +198,6 @@ def init_single_subject_wf(subject_id, session_id=None): derivatives = config.execution.derivatives or {} anat_modality = "t1w" if subject_data["t1w"] else "t2w" spaces = config.workflow.spaces - surface_sampler = config.workflow.surface_sampler # Make sure we always go through these two checks if not anat_only and not subject_data["bold"]: task_id = config.execution.task_id diff --git a/nibabies/workflows/bold/base.py b/nibabies/workflows/bold/base.py index 67d5e88c..661c89e1 100644 --- a/nibabies/workflows/bold/base.py +++ b/nibabies/workflows/bold/base.py @@ -302,12 +302,8 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non Gridded (volumetric) resamplings were performed using `antsApplyTransforms` (ANTs), configured with Lanczos interpolation to minimize the smoothing effects of other kernels [@lanczos]. -Non-gridded (surface) resamplings were performed using {method}.. -""".format( - method={"fs": "mri_vol2surf (FreeSurfer)", - "wb": "wb_command -volume-to-surface-mapping (Workbench)" - }[config.workflow.surface_sampler] - ) +Non-gridded (surface) resamplings were performed using `mri_vol2surf` (FreeSurfer) +""" inputnode = pe.Node( @@ -857,7 +853,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non surface_spaces=freesurfer_spaces, medial_surface_nan=config.workflow.medial_surface_nan, project_goodvoxels=config.workflow.project_goodvoxels, - surface_sampler=config.workflow.surface_sampler, name="bold_surf_wf", ) # fmt:off diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index eae9c2f5..d957deb5 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -9,18 +9,19 @@ .. autofunction:: init_bold_preproc_trans_wf """ + + import nipype.interfaces.workbench as wb from nipype.interfaces import freesurfer as fs from nipype.interfaces import fsl from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe from nipype import Function -from niworkflows.interfaces.freesurfer import MakeMidthickness from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms from niworkflows.interfaces.freesurfer import MedialNaNs from ...interfaces.volume import CreateSignedDistanceVolume -from ...interfaces.metric import MetricDilate, MetricResample -from ...interfaces.wbvoltosurf import VolumeToSurfaceMapping +from os.path import basename + @@ -31,7 +32,6 @@ def init_bold_surf_wf(mem_gb, surface_spaces, medial_surface_nan, project_goodvoxels, - surface_sampler, name="bold_surf_wf"): """ Sample functional images to FreeSurfer surfaces. @@ -39,11 +39,6 @@ def init_bold_surf_wf(mem_gb, For each vertex, the cortical ribbon is sampled at six points (spaced 20% of thickness apart) and averaged. - If --surface-sampler wb is used, Workbench's wb_command -volume-to-surface-mapping - with -ribbon-constrained option is used instead of the default FreeSurfer mri_vol2surf. - Note that unlike HCP, no additional spatial smoothing is applied to the surface-projected - data. - If --project-goodvoxels is used, a "goodvoxels" BOLD mask, as described in [@hcppipelines], is generated and applied to the functional image before sampling to surface. Outputs are in GIFTI format. @@ -57,8 +52,7 @@ def init_bold_surf_wf(mem_gb, wf = init_bold_surf_wf(mem_gb=0.1, surface_spaces=['fsnative', 'fsaverage5'], medial_surface_nan=False, - project_goodvoxels=False, - surface_sampler="fs") + project_goodvoxels=False) Parameters ---------- @@ -72,9 +66,6 @@ def init_bold_surf_wf(mem_gb, project_goodvoxels : :obj:`bool` Exclude voxels with locally high coefficient of variation, or that lie outside the cortical surfaces, from the surface projection. - surface_sampler : :obj:`str` - 'fs' (default) or 'wb' to specify FreeSurfer-based or Workbench-based - volume to surface mapping Inputs ------ @@ -109,10 +100,11 @@ def init_bold_surf_wf(mem_gb, workflow.__desc__ = """\ The BOLD time-series were resampled onto the following surfaces (FreeSurfer reconstruction nomenclature): -{out_spaces} with {sampling_method} +{out_spaces} """.format( - out_spaces=", ".join(["*%s*" % s for s in surface_spaces]) - ) + out_spaces=", ".join(["*%s*" % s for s in surface_spaces]), +) + if project_goodvoxels: workflow.__desc__ += """\ @@ -176,7 +168,6 @@ def select_target(subject_id, space): ), name_source=['source_file'], keep_extension=False, - name_template='%s.func.gii', iterfield=["hemi"], name="sampler", mem_gb=mem_gb * 3, @@ -188,6 +179,15 @@ def select_target(subject_id, space): MedialNaNs(), iterfield=["in_file"], name="medial_nans", mem_gb=DEFAULT_MEMORY_MIN_GB ) + # Rename the source file to the output space to simplify naming later + prepend_hemi = pe.MapNode( + niu.Rename(format_string="%(hemi)s.target.func.gii"), + iterfield=["in_file", "hemi"], + name="prepend_hemi", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + prepend_hemi.inputs.hemi=["lh", "rh"] + update_metadata = pe.MapNode( GiftiSetAnatomicalStructure(), iterfield=["in_file"], @@ -200,41 +200,6 @@ def select_target(subject_id, space): joinsource="itersource", name="outputnode", ) - - if not project_goodvoxels: - # fmt: off - workflow.connect([ - (inputnode, get_fsnative, [('subject_id', 'subject_id'), - ('subjects_dir', 'subjects_dir')]), - (inputnode, targets, [('subject_id', 'subject_id')]), - (inputnode, rename_src, [('source_file', 'in_file')]), - (inputnode, itk2lta, [('source_file', 'src_file'), - ('t1w2fsnative_xfm', 'in_file')]), - (get_fsnative, itk2lta, [('brain', 'dst_file')]), # InfantFS: Use brain instead of T1 - (inputnode, sampler, [('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id')]), - (itersource, targets, [('target', 'space')]), - (itersource, rename_src, [('target', 'subject')]), - (itk2lta, sampler, [('out', 'reg_file')]), - (targets, sampler, [('out', 'target_subject')]), - (rename_src, sampler, [('out_file', 'source_file')]), - (update_metadata, outputnode, [('out_file', 'surfaces')]), - (itersource, outputnode, [('target', 'target')]), - ]) - # fmt: on - - if not medial_surface_nan: - workflow.connect(sampler, "out_file", update_metadata, "in_file") - return workflow - - # fmt: off - workflow.connect([ - (inputnode, medial_nans, [('subjects_dir', 'subjects_dir')]), - (sampler, medial_nans, [('out_file', 'in_file')]), - (medial_nans, update_metadata, [('out_file', 'in_file')]), - ]) - # fmt: on - return workflow # 0, 1 = wm; 2, 3 = pial; 6, 7 = mid # note that order of lh / rh within each surf type is not guaranteed due to use @@ -242,12 +207,12 @@ def select_target(subject_id, space): # to ensure consistent ordering select_wm = pe.Node( niu.Select(index=[0, 1]), - name="select_wm_pial", + name="select_wm", mem_gb=DEFAULT_MEMORY_MIN_GB, ) select_pial = pe.Node( - niu.Select(index=[2, 3],), + niu.Select(index=[2, 3]), name="select_pial", mem_gb=DEFAULT_MEMORY_MIN_GB, ) @@ -496,7 +461,6 @@ def _calc_lower_thr(in_stats): fsl.ApplyMask(), name_source=['in_file'], keep_extension=True, - name_template='%s', name="goodvoxels_ribbon_mask", mem_gb=DEFAULT_MEMORY_MIN_GB, ) @@ -505,93 +469,53 @@ def _calc_lower_thr(in_stats): fsl.ApplyMask(), name_source=['in_file'], keep_extension=True, - name_template='%s', name="apply_goodvoxels_ribbon_mask", mem_gb=mem_gb * 3, ) - - get_target_wm = pe.MapNode( - FreeSurferSource(), - iterfield=["in_file", "surface"], - name="get_target_wm", - run_without_submitting=True, - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - get_target_wm.inputs.hemi=["lh","rh"] - make_target_midthick = pe.MapNode( - MakeMidthickness(thickness=True, distance=0.5), - iterfield=["in_file"], - name="make_target_midthick", - run_without_submitting=True, - mem_gb=mem_gb * 3, - ) + if not project_goodvoxels: + # fmt: off + workflow.connect([ + (inputnode, get_fsnative, [('subject_id', 'subject_id'), + ('subjects_dir', 'subjects_dir')]), + (inputnode, targets, [('subject_id', 'subject_id')]), + (inputnode, rename_src, [('source_file', 'in_file')]), + (inputnode, itk2lta, [('source_file', 'src_file'), + ('t1w2fsnative_xfm', 'in_file')]), + (get_fsnative, itk2lta, [('brain', 'dst_file')]), # InfantFS: Use brain instead of T1 + (inputnode, sampler, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id')]), + (itersource, targets, [('target', 'space')]), + (itersource, rename_src, [('target', 'subject')]), + (itk2lta, sampler, [('out', 'reg_file')]), + (targets, sampler, [('out', 'target_subject')]), + (rename_src, sampler, [('out_file', 'source_file')]), + (update_metadata, outputnode, [('out_file', 'surfaces')]), + (itersource, outputnode, [('target', 'target')]), + ]) + # fmt: on - target_midthick_gifti = pe.MapNode( - fs.MRIsConvert(out_datatype="gii"), - iterfield=["in_file"], - name="target_midthick_gifti", - run_without_submitting=True, - mem_gb=mem_gb, - ) - - wbsampler = pe.MapNode( - VolumeToSurfaceMapping( - mapping_method="ribbon-constrained", - ), - iterfield=["in_file"], - name="wbsampler", - mem_gb=mem_gb * 3, - ) + if not medial_surface_nan: + workflow.connect(sampler, "out_file", update_metadata, "in_file") + return workflow - metric_dilate = pe.MapNode( - MetricDilate( - distance=10, - nearest=True, - ), - iterfield=["in_file", "surface"], - name="metric_dilate", - mem_gb=mem_gb * 3, - ) - - native_to_target = pe.MapNode( - MetricResample( - method="ADAP_BARY_AREA", - area_metrics=True, - ), - iterfield=[ - "in_file", - "out_file", - "new_sphere", - "new_area", - "current_sphere", - "current_area", - ], - name="native_to_target", - ) - native_to_target.inputs.new_sphere = [ - str( - tf.api.get("fsaverage", hemi=hemi, density="164k", desc="std", suffix="sphere") - ) - for hemi in "LR" - ] - - # make FS midthick if target doesn't have them already - workflow.connect([ - (inputnode, get_target_wm, [('subjects_dir', 'subjects_dir')]), - (targets, get_target_wm, [('out', 'subject_id')]), - (get_target_wm, make_target_midthick, [("white", "in_file")]), - (make_target_midthick, target_midthick_gifti, [("out_file", "in_file")]), - ]) + # fmt: off + workflow.connect([ + (inputnode, medial_nans, [('subjects_dir', 'subjects_dir')]), + (sampler, medial_nans, [('out_file', 'in_file')]), + (medial_nans, update_metadata, [('out_file', 'in_file')]), + ]) + # fmt: on + return workflow # make HCP-style ribbon volume in T1w space workflow.connect([ (inputnode, select_wm, [("anat_giftis", "inlist")]), (inputnode, select_pial, [("anat_giftis", "inlist")]), (inputnode, select_midthick, [("anat_giftis", "inlist")]), - (select_wm, create_wm_distvol, [("out", "surface")]), + (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surface")]), (inputnode, create_wm_distvol, [("t1w_mask", "ref_space")]), - (select_pial, create_pial_distvol, [("out", "surface")]), + (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surface")]), (inputnode, create_pial_distvol, [("t1w_mask", "ref_space")]), (create_wm_distvol, thresh_wm_distvol, [("out_vol", "in_file")]), (create_pial_distvol, uthresh_pial_distvol, [("out_vol", "in_file")]), @@ -647,102 +571,51 @@ def _calc_lower_thr(in_stats): (merge_goodvoxels_operands, goodvoxels_mask, [("out", "operand_files")]), ]) - if surface_sampler is "fs": - # apply goodvoxels ribbon mask to bold - workflow.connect([ - (goodvoxels_mask, goodvoxels_ribbon_mask, [("out_file", "in_file")]), - (ribbon_boldsrc_xfm, goodvoxels_ribbon_mask, [("output_image", "mask_file")]), - (goodvoxels_ribbon_mask, apply_goodvoxels_ribbon_mask, [("out_file", "mask_file")]), - (rename_src, apply_goodvoxels_ribbon_mask, [("out_file", "in_file")]), - ]) - - # project masked bold to target surfs - workflow.connect([ - (inputnode, get_fsnative, [("subject_id", "subject_id"), - ("subjects_dir", "subjects_dir")]), - (inputnode, targets, [("subject_id", "subject_id")]), - (inputnode, rename_src, [("source_file", "in_file")]), - (inputnode, itk2lta, [("source_file", "src_file"), - ("t1w2fsnative_xfm", "in_file")]), - (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 - (inputnode, sampler, [("subjects_dir", "subjects_dir"), - ("subject_id", "subject_id")]), - (itersource, targets, [("target", "space")]), - (itersource, rename_src, [("target", "subject")]), - (itk2lta, sampler, [("out", "reg_file")]), - (targets, sampler, [("out", "target_subject")]), - (apply_goodvoxels_ribbon_mask, sampler, [("out_file", "source_file")]), - (sampler, metric_dilate, [("out_file", "in_file")]), - (target_midthick_gifti, metric_dilate, [("converted", "surface")]), - (update_metadata, outputnode, [("out_file", "surfaces")]), - (itersource, outputnode, [("target", "target")]), - ]) - - # fmt:on - if not medial_surface_nan: - # fmt:off - workflow.connect([ - (metric_dilate, update_metadata, [("out_file", "in_file")]), - ]) - # fmt:on - return workflow + # apply goodvoxels ribbon mask to bold + workflow.connect([ + (goodvoxels_mask, goodvoxels_ribbon_mask, [("out_file", "in_file")]), + (ribbon_boldsrc_xfm, goodvoxels_ribbon_mask, [("output_image", "mask_file")]), + (goodvoxels_ribbon_mask, apply_goodvoxels_ribbon_mask, [("out_file", "mask_file")]), + (rename_src, apply_goodvoxels_ribbon_mask, [("out_file", "in_file")]), + ]) + # project masked bold to target surfs + workflow.connect([ + (inputnode, get_fsnative, [("subject_id", "subject_id"), + ("subjects_dir", "subjects_dir")]), + (inputnode, targets, [("subject_id", "subject_id")]), + (inputnode, rename_src, [("source_file", "in_file")]), + (inputnode, itk2lta, [("source_file", "src_file"), + ("t1w2fsnative_xfm", "in_file")]), + (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 + (inputnode, sampler, [("subjects_dir", "subjects_dir"), + ("subject_id", "subject_id")]), + (itersource, targets, [("target", "space")]), + (itersource, rename_src, [("target", "subject")]), + (itk2lta, sampler, [("out", "reg_file")]), + (targets, sampler, [("out", "target_subject")]), + (apply_goodvoxels_ribbon_mask, sampler, [("out_file", "source_file")]), + (update_metadata, outputnode, [("out_file", "surfaces")]), + (itersource, outputnode, [("target", "target")]), + ]) + + # fmt:on + if not medial_surface_nan: # fmt:off workflow.connect([ - (metric_dilate, medial_nans, [("out_file", "in_file")]), - (medial_nans, update_metadata, [("out_file", "in_file")]), + (sampler, update_metadata, [("out_file", "in_file")]), ]) # fmt:on return workflow - - # wb method first projects to native surfs, then resamples to fsaverage - # (fsaverage{3,4,5,6} not supported as yet) - if surface_sampler is "wb": - workflow.connect([ - (inputnode, get_fsnative, [("subject_id", "subject_id"), - ("subjects_dir", "subjects_dir")]), - (inputnode, targets, [("subject_id", "subject_id")]), - (inputnode, rename_src, [("source_file", "in_file")]), - (inputnode, itk2lta, [("source_file", "src_file"), - ("t1w2fsnative_xfm", "in_file")]), - (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 - (itersource, targets, [("target", "space")]), - (itersource, rename_src, [("target", "subject")]), - (rename_src, wbsampler, [("out_file", "in_file")]), - (select_midthick, wbsampler, [(("out_file", _sorted), "surface")]), - (select_wm, wbsampler, [(("out_file", _sorted), "inner_surf")]), - (select_pial, wbsampler, [(("out_file", _sorted), "outer_surf")]), - (goodvoxels_mask, wbsampler, [("out_file", "roi_volume")]), - (wbsampler, metric_dilate, [("out_file", "in_file")]), - (metric_dilate, native_to_target, [("out_file", "in_file")]), - (select_midthick, native_to_target, [(("out_file", _sorted), "current_area")]), - (get_fsnative, native_to_target, [(("sphere_reg", _sorted), "current_sphere")]), - (target_midthick_gifti, native_to_target, [("converted", "new_area")]), - (update_metadata, outputnode, [("out_file", "surfaces")]), - (itersource, outputnode, [("target", "target")]), - ]) - - # fmt:on - if not medial_surface_nan: - # fmt:off - workflow.connect([ - (sampler, metric_dilate, [("out_file", "in_file")]), - (target_midthick_gifti, metric_dilate, [("converted", "surface")]), - (metric_dilate, update_metadata, [("out_file", "in_file")]), - ]) - # fmt:on - return workflow - - # fmt:off - workflow.connect([ - (sampler, metric_dilate, [("out_file", "in_file")]), - (target_midthick_gifti, metric_dilate, [("converted", "surface")]), - (metric_dilate, medial_nans, [("out_file", "in_file")]), - (medial_nans, update_metadata, [("out_file", "in_file")]), - ]) - # fmt:on - return workflow + # fmt:off + workflow.connect([ + (inputnode, medial_nans, [('subjects_dir', 'subjects_dir')]), + (sampler, medial_nans, [("out_file", "in_file")]), + (medial_nans, update_metadata, [("out_file", "in_file")]), + ]) + # fmt:on + return workflow def init_bold_std_trans_wf( freesurfer, @@ -1270,6 +1143,7 @@ def init_bold_grayords_wf(grayord_density, mem_gb, repetition_time, name="bold_g str(tf.api.get("fsaverage", hemi=hemi, density="164k", desc="std", suffix="sphere")) for hemi in "LR" ] + resample.inputs.current_area = [ str( tf.api.get("fsaverage", hemi=hemi, density="164k", desc="vaavg", suffix="midthickness") @@ -1421,5 +1295,6 @@ def _itk2lta(in_file, src_file, dst_file): ).to_filename(out_file, moving=dst_file, fmt="fs") return str(out_file) -def _sorted(inlist): - return sorted(inlist) +def _sorted_by_basename(inlist): + from os.path import basename + return sorted(inlist, key=lambda x: str(basename(x))) From 568191cc4af0678e8216fa287d87913fa1cd223f Mon Sep 17 00:00:00 2001 From: madisoth <56737848+madisoth@users.noreply.github.com> Date: Thu, 15 Sep 2022 10:44:31 -0700 Subject: [PATCH 03/11] docstring and style fixes --- nibabies/interfaces/volume.py | 2 +- nibabies/workflows/bold/resampling.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/nibabies/interfaces/volume.py b/nibabies/interfaces/volume.py index 8866f990..b520bf9f 100644 --- a/nibabies/interfaces/volume.py +++ b/nibabies/interfaces/volume.py @@ -132,7 +132,7 @@ class CreateSignedDistanceVolume(WBCommand): (negative) crossings of a vertical ray from the point, then counts as inside if the total is odd, negative, or nonzero, respectively. - >>> from nipype.interfaces.workbench import CreateSignedDistanceVolume + >>> from nibabies.interfaces.volume import CreateSignedDistanceVolume >>> distvol = CreateSignedDistanceVolume() >>> distvol.inputs.surface = 'sub-01.L.pial.native.surf.gii' >>> distvol.inputs.refspace = 'sub-01_T1w.nii.gz' diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index d957deb5..f773c61c 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -103,7 +103,7 @@ def init_bold_surf_wf(mem_gb, {out_spaces} """.format( out_spaces=", ".join(["*%s*" % s for s in surface_spaces]), -) + ) if project_goodvoxels: From 2493e8c8577ee5fd6155668ca456d2835582642b Mon Sep 17 00:00:00 2001 From: madisoth <56737848+madisoth@users.noreply.github.com> Date: Thu, 15 Sep 2022 13:52:48 -0700 Subject: [PATCH 04/11] remove CreateSignedDistanceVolume doctest for now (needs anatomical surfaces for test sub) --- nibabies/interfaces/volume.py | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/nibabies/interfaces/volume.py b/nibabies/interfaces/volume.py index b520bf9f..47e74a08 100644 --- a/nibabies/interfaces/volume.py +++ b/nibabies/interfaces/volume.py @@ -132,27 +132,6 @@ class CreateSignedDistanceVolume(WBCommand): (negative) crossings of a vertical ray from the point, then counts as inside if the total is odd, negative, or nonzero, respectively. - >>> from nibabies.interfaces.volume import CreateSignedDistanceVolume - >>> distvol = CreateSignedDistanceVolume() - >>> distvol.inputs.surface = 'sub-01.L.pial.native.surf.gii' - >>> distvol.inputs.refspace = 'sub-01_T1w.nii.gz' - >>> distvol.inputs.out_vol = 'sub-01.L.pial.native.surf.distvol.nii.gz' - >>> distvol.inputs.roi_out = 'sub-01.L.pial.native.surf.distvolroi.nii.gz' - >>> distvol.inputs.fill_value = 0 - >>> distvol.inputs.exact_limit = 5 - >>> distvol.inputs.approx_limit = 20 - >>> distvol.inputs.approx_neighborhood = 2 - >>> distvol.inputs.winding_method = 'EVEN_ODD' - >>> distvol.cmdline - 'wb_command -create-signed-distance-volume sub-01.L.pial.native.surf.gii \ - sub-01_T1w.nii.gz \ - sub-01.L.pial.native.surf.distvol.nii.gz \ - -roi-out sub-01.L.pial.native.surf.distvolroi.nii.gz \ - -fill-value 0 \ - -exact-limit 5 \ - -approx-limit 20 \ - -approx-neighborhood 2 \ - -winding EVEN_ODD' """ input_spec = CreateSignedDistanceVolumeInputSpec From 5bd684a58043ca553378993a102ca0de62fe5e1c Mon Sep 17 00:00:00 2001 From: madisoth <56737848+madisoth@users.noreply.github.com> Date: Thu, 15 Sep 2022 14:37:31 -0700 Subject: [PATCH 05/11] docstring update --- nibabies/workflows/bold/resampling.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index f773c61c..de796861 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -77,8 +77,6 @@ def init_bold_surf_wf(mem_gb, FreeSurfer subject ID t1w2fsnative_xfm LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space - itk_bold_to_t1 - Affine transform from ``ref_bold_brain`` to T1w space (ITK format) anat_giftis GIFTI anatomical surfaces in T1w space t1w_mask From 114f982e3df79429badd79b52aa99ddb5eef9a0c Mon Sep 17 00:00:00 2001 From: Greg Conan Date: Thu, 24 Nov 2022 03:45:24 -0600 Subject: [PATCH 06/11] Modularized goodvoxels code in resampling.py Removed duplicate code by abstracting it into functions, clarifying the goodvoxels logic --- nibabies/workflows/bold/resampling.py | 104 ++++++++++++-------------- 1 file changed, 48 insertions(+), 56 deletions(-) diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index de796861..b442ab05 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -199,6 +199,29 @@ def select_target(subject_id, space): name="outputnode", ) + if project_goodvoxels: + workflow, apply_goodvoxels_ribbon_mask = apply_goodvoxels( + workflow, inputnode, mem_gb, rename_src + ) + + # If not applying goodvoxels, then get sampler.source_file from rename_src + # instead in connect_bold_surf_wf + apply_goodvx_or_rename = apply_goodvoxels_ribbon_mask + else: + apply_goodvx_or_rename = rename_src + + # fmt: off + workflow = connect_bold_surf_wf( + workflow, inputnode, outputnode, get_fsnative, itersource, targets, + rename_src, sampler, itk2lta, update_metadata, apply_goodvx_or_rename + ) + # fmt: on + return apply_med_surf_nans_if(medial_surface_nan, workflow, + inputnode, sampler, + update_metadata, medial_nans) + + +def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): # 0, 1 = wm; 2, 3 = pial; 6, 7 = mid # note that order of lh / rh within each surf type is not guaranteed due to use # of unsorted glob by FreeSurferSource prior, but we can do a sort @@ -470,41 +493,6 @@ def _calc_lower_thr(in_stats): name="apply_goodvoxels_ribbon_mask", mem_gb=mem_gb * 3, ) - - if not project_goodvoxels: - # fmt: off - workflow.connect([ - (inputnode, get_fsnative, [('subject_id', 'subject_id'), - ('subjects_dir', 'subjects_dir')]), - (inputnode, targets, [('subject_id', 'subject_id')]), - (inputnode, rename_src, [('source_file', 'in_file')]), - (inputnode, itk2lta, [('source_file', 'src_file'), - ('t1w2fsnative_xfm', 'in_file')]), - (get_fsnative, itk2lta, [('brain', 'dst_file')]), # InfantFS: Use brain instead of T1 - (inputnode, sampler, [('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id')]), - (itersource, targets, [('target', 'space')]), - (itersource, rename_src, [('target', 'subject')]), - (itk2lta, sampler, [('out', 'reg_file')]), - (targets, sampler, [('out', 'target_subject')]), - (rename_src, sampler, [('out_file', 'source_file')]), - (update_metadata, outputnode, [('out_file', 'surfaces')]), - (itersource, outputnode, [('target', 'target')]), - ]) - # fmt: on - - if not medial_surface_nan: - workflow.connect(sampler, "out_file", update_metadata, "in_file") - return workflow - - # fmt: off - workflow.connect([ - (inputnode, medial_nans, [('subjects_dir', 'subjects_dir')]), - (sampler, medial_nans, [('out_file', 'in_file')]), - (medial_nans, update_metadata, [('out_file', 'in_file')]), - ]) - # fmt: on - return workflow # make HCP-style ribbon volume in T1w space workflow.connect([ @@ -533,10 +521,10 @@ def _calc_lower_thr(in_stats): # make HCP-style "goodvoxels" mask in t1w space for filtering outlier voxels # in bold timeseries, based on modulated normalized covariance workflow.connect([ - (combine_ribbon_vol_hemis, ribbon_boldsrc_xfm, [("out_file", 'input_image')]), + (combine_ribbon_vol_hemis, ribbon_boldsrc_xfm, [("out_file", "input_image")]), (rename_src, stdev_volume, [("out_file", "in_file")]), (rename_src, mean_volume, [("out_file", "in_file")]), - (mean_volume, ribbon_boldsrc_xfm, [('out_file', 'reference_image')]), + (mean_volume, ribbon_boldsrc_xfm, [("out_file", "reference_image")]), (stdev_volume, cov_volume, [("out_file", "in_file")]), (mean_volume, cov_volume, [("out_file", "operand_file")]), (cov_volume, cov_ribbon, [("out_file", "in_file")]), @@ -577,7 +565,27 @@ def _calc_lower_thr(in_stats): (rename_src, apply_goodvoxels_ribbon_mask, [("out_file", "in_file")]), ]) - # project masked bold to target surfs + return workflow, apply_goodvoxels_ribbon_mask + + +def apply_med_surf_nans_if(medial_surface_nan, workflow, inputnode, sampler, + update_metadata, medial_nans): + if medial_surface_nan: + # fmt: off + workflow.connect([ + (inputnode, medial_nans, [("subjects_dir", "subjects_dir")]), + (sampler, medial_nans, [("out_file", "in_file")]), + (medial_nans, update_metadata, [("out_file", "in_file")]), + ]) + # fmt: on + else: + workflow.connect(sampler, "out_file", update_metadata, "in_file") + return workflow + + +def connect_bold_surf_wf(workflow, inputnode, outputnode, get_fsnative, + itersource, targets, rename_src, sampler, itk2lta, + update_metadata, apply_goodvoxels_or_rename): workflow.connect([ (inputnode, get_fsnative, [("subject_id", "subject_id"), ("subjects_dir", "subjects_dir")]), @@ -592,29 +600,13 @@ def _calc_lower_thr(in_stats): (itersource, rename_src, [("target", "subject")]), (itk2lta, sampler, [("out", "reg_file")]), (targets, sampler, [("out", "target_subject")]), - (apply_goodvoxels_ribbon_mask, sampler, [("out_file", "source_file")]), + (apply_goodvoxels_or_rename, sampler, [("out_file", "source_file")]), (update_metadata, outputnode, [("out_file", "surfaces")]), (itersource, outputnode, [("target", "target")]), ]) - - # fmt:on - if not medial_surface_nan: - # fmt:off - workflow.connect([ - (sampler, update_metadata, [("out_file", "in_file")]), - ]) - # fmt:on - return workflow - - # fmt:off - workflow.connect([ - (inputnode, medial_nans, [('subjects_dir', 'subjects_dir')]), - (sampler, medial_nans, [("out_file", "in_file")]), - (medial_nans, update_metadata, [("out_file", "in_file")]), - ]) - # fmt:on return workflow + def init_bold_std_trans_wf( freesurfer, mem_gb, From 946dc7e86284ec9281308eab756a772dd6ab09bf Mon Sep 17 00:00:00 2001 From: Greg Conan Date: Fri, 2 Dec 2022 18:55:01 -0600 Subject: [PATCH 07/11] Formatting fixes and resampling.py modularization - Implemented the formatting suggestions by Mathias at the URLs below: - https://github.com/nipreps/nibabies/pull/235 - https://github.com/nipreps/fmriprep/pull/2855 - Split the 3 workflow-connect blocks in the project-goodvoxels section of resampling.py into 3 different functions: ribbon, outliers, and vol-to-surf - Incomplete comments in resampling.py --- nibabies/interfaces/volume.py | 82 ++++++------- nibabies/workflows/bold/resampling.py | 168 +++++++++++++++++--------- 2 files changed, 150 insertions(+), 100 deletions(-) diff --git a/nibabies/interfaces/volume.py b/nibabies/interfaces/volume.py index 47e74a08..4c6187b3 100644 --- a/nibabies/interfaces/volume.py +++ b/nibabies/interfaces/volume.py @@ -12,81 +12,75 @@ class CreateSignedDistanceVolumeInputSpec(CommandLineInputSpec): - surface = File( + surf_file = File( exists=True, mandatory=True, - argstr="%s ", + argstr="%s", position=0, - desc="the input surface", + desc="Input surface GIFTI file (.surf.gii)", ) - ref_space = File( + ref_file = File( exists=True, mandatory=True, - argstr="%s ", + argstr="%s", position=1, - desc="a volume in the desired output space (dims, spacing, origin)", + desc="NIfTI volume in the desired output space (dims, spacing, origin)", ) - out_vol = File( + out_file = File( name_source=["surface"], - name_template="%s.distvol.nii.gz", - argstr="%s ", + name_template="%s_distvol.nii.gz", + argstr="%s", position=2, - desc="output - the output volume", + desc="Name of output volume containing signed distances", ) - roi_out = File( + out_mask = File( name_source=["surface"], - name_template="%s.distvolroi.nii.gz", - argstr="-roi-out %s ", - position=3, - desc="output roi volume of where the output has a computed value", + name_template="%s_distmask.nii.gz", + argstr="-roi-out %s", + desc="Name of file to store a mask where the ``out_file`` has a computed value", ) fill_value = traits.Float( mandatory=False, - argstr="-fill-value %f ", - position=4, - desc="specify a value to put in all voxels that don't get assigned a distance, default 0", + default=0., + usedefault=True, + argstr="-fill-value %f", + desc="value to put in all voxels that don't get assigned a distance", ) exact_limit = traits.Float( - mandatory=False, - argstr="-exact-limit %f ", - position=5, - desc="specify distance for exact output in mm, default 5", + default=5., + usedefault=True, + argstr="-exact-limit %f", + desc="distance for exact output in mm", ) approx_limit = traits.Float( - mandatory=False, - argstr="-approx-limit %f ", - position=6, - desc="specify distance for approximate output in mm, default 20", + default=20., + usedefault=True, + argstr="-approx-limit %f", + desc="distance for approximate output in mm", ) approx_neighborhood = traits.Int( - mandatory=False, - argstr="-approx-neighborhood %d ", - position=7, - desc="size of neighborhood cube measured from center to face, default 2 = 5x5x5", + default=2, + usedefault=True, + argstr="-approx-neighborhood %d", + desc="size of neighborhood cube measured from center to face in voxels, default 2 = 5x5x5", ) winding_method = traits.Enum( "EVEN_ODD", "NEGATIVE", "NONZERO", "NORMALS", - argstr="-winding %s ", + argstr="-winding %s", usedefault=True, - position=8, - desc="winding method for point inside surface test, choices: " - "EVEN_ODD (default) " - "NEGATIVE " - "NONZERO " - "NORMALS " + desc="winding method for point inside surface test" ) class CreateSignedDistanceVolumeOutputSpec(TraitedSpec): - out_vol = File(exists=True, desc="output - the output volume") - roi_out = File(desc="output roi volume of where the output has a computed value") + out_file = File(desc="Name of output volume containing signed distances") + out_mask = File(desc="Name of file to store a mask where the ``out_file`` has a computed value") class CreateSignedDistanceVolume(WBCommand): - """ -CREATE SIGNED DISTANCE VOLUME FROM SURFACE + """CREATE SIGNED DISTANCE VOLUME FROM SURFACE wb_command -create-signed-distance-volume - the input surface - a volume in the desired output space (dims, spacing, origin) @@ -136,8 +130,4 @@ class CreateSignedDistanceVolume(WBCommand): input_spec = CreateSignedDistanceVolumeInputSpec output_spec = CreateSignedDistanceVolumeOutputSpec - _cmd = "wb_command -create-signed-distance-volume" - - def _list_outputs(self): - outputs = super(CreateSignedDistanceVolume, self)._list_outputs() - return outputs + _cmd = "wb_command -create-signed-distance-volume" \ No newline at end of file diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index b442ab05..3e418d44 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -174,7 +174,8 @@ def select_target(subject_id, space): # Refine if medial vertices should be NaNs medial_nans = pe.MapNode( - MedialNaNs(), iterfield=["in_file"], name="medial_nans", mem_gb=DEFAULT_MEMORY_MIN_GB + MedialNaNs(), iterfield=["in_file"], name="medial_nans", + mem_gb=DEFAULT_MEMORY_MIN_GB ) # Rename the source file to the output space to simplify naming later @@ -183,6 +184,7 @@ def select_target(subject_id, space): iterfield=["in_file", "hemi"], name="prepend_hemi", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) prepend_hemi.inputs.hemi=["lh", "rh"] @@ -191,17 +193,25 @@ def select_target(subject_id, space): iterfield=["in_file"], name="update_metadata", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) outputnode = pe.JoinNode( niu.IdentityInterface(fields=["surfaces", "target"]), joinsource="itersource", name="outputnode", + run_without_submitting=True, ) if project_goodvoxels: - workflow, apply_goodvoxels_ribbon_mask = apply_goodvoxels( - workflow, inputnode, mem_gb, rename_src + workflow, combine_ribbon_vol_hemis = make_ribbon_file( + workflow, mem_gb, inputnode + ) + workflow, goodvoxels_mask, ribbon_boldsrc_xfm = mask_outliers( + workflow, mem_gb, rename_src, combine_ribbon_vol_hemis + ) + workflow, apply_goodvoxels_ribbon_mask = project_vol_to_surf( + workflow, mem_gb, rename_src, goodvoxels_mask, ribbon_boldsrc_xfm ) # If not applying goodvoxels, then get sampler.source_file from rename_src @@ -221,7 +231,14 @@ def select_target(subject_id, space): update_metadata, medial_nans) -def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): +def make_ribbon_file(workflow, mem_gb, inputnode): + """ + _summary_ + :param workflow: _type_, _description_ + :param mem_gb: _type_, _description_ + :param inputnode: _type_, _description_ + :return: _type_, _description_ + """ # 0, 1 = wm; 2, 3 = pial; 6, 7 = mid # note that order of lh / rh within each surf type is not guaranteed due to use # of unsorted glob by FreeSurferSource prior, but we can do a sort @@ -230,30 +247,33 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): niu.Select(index=[0, 1]), name="select_wm", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) select_pial = pe.Node( niu.Select(index=[2, 3]), name="select_pial", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) select_midthick = pe.Node( niu.Select(index=[6, 7]), name="select_midthick", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) create_wm_distvol = pe.MapNode( CreateSignedDistanceVolume(), - iterfield=["surface"], + iterfield=["surf_file"], name="create_wm_distvol", mem_gb=mem_gb, ) create_pial_distvol = pe.MapNode( CreateSignedDistanceVolume(), - iterfield=["surface"], + iterfield=["surf_file"], name="create_pial_distvol", mem_gb=mem_gb, ) @@ -290,6 +310,7 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): niu.Split(splits=[1, 1]), name="split_wm_distvol", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) merge_wm_distvol_no_flatten = pe.Node( @@ -297,10 +318,11 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): no_flatten=True, name="merge_wm_distvol_no_flatten", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) make_ribbon_vol = pe.MapNode( - fsl.maths.MultiImageMaths(op_string="-mas %s -mul 255 "), + fsl.maths.MultiImageMaths(op_string="-mas %s -mul 255"), iterfield=["in_file", "operand_files"], name="make_ribbon_vol", mem_gb=DEFAULT_MEMORY_MIN_GB, @@ -317,6 +339,7 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): niu.Split(splits=[1, 1], squeeze=True), name="split_squeeze_ribbon", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) combine_ribbon_vol_hemis = pe.Node( @@ -325,6 +348,34 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): mem_gb=DEFAULT_MEMORY_MIN_GB, ) + # make HCP-style ribbon volume in T1w space + workflow.connect([ + (inputnode, select_wm, [("anat_giftis", "inlist")]), + (inputnode, select_pial, [("anat_giftis", "inlist")]), + (inputnode, select_midthick, [("anat_giftis", "inlist")]), + (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surf_file")]), + (inputnode, create_wm_distvol, [("t1w_mask", "ref_file")]), + (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surf_file")]), + (inputnode, create_pial_distvol, [("t1w_mask", "ref_file")]), + (create_wm_distvol, thresh_wm_distvol, [("out_file", "in_file")]), + (create_pial_distvol, uthresh_pial_distvol, [("out_file", "in_file")]), + (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), + (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), + (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), + (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), + (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), + (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), + (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), + ]) + + return workflow, combine_ribbon_vol_hemis + + +def mask_outliers(workflow, mem_gb, rename_src, combine_ribbon_vol_hemis): ribbon_boldsrc_xfm = pe.Node( ApplyTransforms(interpolation='MultiLabel', transforms='identity'), @@ -357,13 +408,13 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): ) cov_ribbon_mean = pe.Node( - fsl.ImageStats(op_string='-M '), + fsl.ImageStats(op_string='-M'), name="cov_ribbon_mean", mem_gb=DEFAULT_MEMORY_MIN_GB, ) cov_ribbon_std = pe.Node( - fsl.ImageStats(op_string='-S '), + fsl.ImageStats(op_string='-S'), name="cov_ribbon_std", mem_gb=DEFAULT_MEMORY_MIN_GB, ) @@ -375,7 +426,7 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): ) smooth_norm = pe.Node( - fsl.maths.MathsCommand(args="-bin -s 5 "), + fsl.maths.MathsCommand(args="-bin -s 5"), name="smooth_norm", mem_gb=DEFAULT_MEMORY_MIN_GB, ) @@ -384,10 +435,11 @@ def apply_goodvoxels(workflow, inputnode, mem_gb, rename_src): niu.Merge(1), name="merge_smooth_norm", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) cov_ribbon_norm_smooth = pe.Node( - fsl.maths.MultiImageMaths(op_string='-s 5 -div %s -dilD '), + fsl.maths.MultiImageMaths(op_string='-s 5 -div %s -dilD'), name="cov_ribbon_norm_smooth", mem_gb=DEFAULT_MEMORY_MIN_GB, ) @@ -437,13 +489,13 @@ def _calc_lower_thr(in_stats): ) mod_ribbon_mean = pe.Node( - fsl.ImageStats(op_string='-M '), + fsl.ImageStats(op_string='-M'), name="mod_ribbon_mean", mem_gb=DEFAULT_MEMORY_MIN_GB, ) mod_ribbon_std = pe.Node( - fsl.ImageStats(op_string='-S '), + fsl.ImageStats(op_string='-S'), name="mod_ribbon_std", mem_gb=DEFAULT_MEMORY_MIN_GB, ) @@ -452,6 +504,7 @@ def _calc_lower_thr(in_stats): niu.Merge(2), name="merge_mod_ribbon_stats", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) bin_mean_volume = pe.Node( @@ -464,6 +517,7 @@ def _calc_lower_thr(in_stats): niu.Merge(2), name="merge_goodvoxels_operands", mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, ) goodvoxels_thr = pe.Node( @@ -473,50 +527,10 @@ def _calc_lower_thr(in_stats): ) goodvoxels_mask = pe.Node( - fsl.maths.MultiImageMaths(op_string='-bin -sub %s -mul -1 '), + fsl.maths.MultiImageMaths(op_string='-bin -sub %s -mul -1'), name="goodvoxels_mask", mem_gb=mem_gb, ) - - goodvoxels_ribbon_mask = pe.Node( - fsl.ApplyMask(), - name_source=['in_file'], - keep_extension=True, - name="goodvoxels_ribbon_mask", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - - apply_goodvoxels_ribbon_mask = pe.Node( - fsl.ApplyMask(), - name_source=['in_file'], - keep_extension=True, - name="apply_goodvoxels_ribbon_mask", - mem_gb=mem_gb * 3, - ) - - # make HCP-style ribbon volume in T1w space - workflow.connect([ - (inputnode, select_wm, [("anat_giftis", "inlist")]), - (inputnode, select_pial, [("anat_giftis", "inlist")]), - (inputnode, select_midthick, [("anat_giftis", "inlist")]), - (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surface")]), - (inputnode, create_wm_distvol, [("t1w_mask", "ref_space")]), - (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surface")]), - (inputnode, create_pial_distvol, [("t1w_mask", "ref_space")]), - (create_wm_distvol, thresh_wm_distvol, [("out_vol", "in_file")]), - (create_pial_distvol, uthresh_pial_distvol, [("out_vol", "in_file")]), - (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), - (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), - (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), - (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), - (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), - (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), - (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), - ]) # make HCP-style "goodvoxels" mask in t1w space for filtering outlier voxels # in bold timeseries, based on modulated normalized covariance @@ -557,6 +571,39 @@ def _calc_lower_thr(in_stats): (merge_goodvoxels_operands, goodvoxels_mask, [("out", "operand_files")]), ]) + return workflow, goodvoxels_mask, ribbon_boldsrc_xfm + + +def project_vol_to_surf(workflow, mem_gb, rename_src, goodvoxels_mask, + ribbon_boldsrc_xfm): + """ + Parameters + ---------- + workflow : :obj: + mem_gb : :obj:`float` + Size of BOLD file in GB + rename_src: _type_, _description_ + goodvoxels_mask: _type_, _description_ + ribbon_boldsrc_xfm: _type_, _description_ + + :return: _type_, _description_ + """ + goodvoxels_ribbon_mask = pe.Node( + fsl.ApplyMask(), + name_source=['in_file'], + keep_extension=True, + name="goodvoxels_ribbon_mask", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + apply_goodvoxels_ribbon_mask = pe.Node( + fsl.ApplyMask(), + name_source=['in_file'], + keep_extension=True, + name="apply_goodvoxels_ribbon_mask", + mem_gb=mem_gb * 3, + ) + # apply goodvoxels ribbon mask to bold workflow.connect([ (goodvoxels_mask, goodvoxels_ribbon_mask, [("out_file", "in_file")]), @@ -570,6 +617,19 @@ def _calc_lower_thr(in_stats): def apply_med_surf_nans_if(medial_surface_nan, workflow, inputnode, sampler, update_metadata, medial_nans): + """ + Parameters + ---------- + medial_surface_nan : :obj:`bool` + Replace medial wall values with NaNs on functional GIFTI files + workflow : _type_, _description_ + inputnode : _type_, _description_ + sampler : _type_, _description_ + update_metadata : _type_, _description_ + medial_nans : _type_, _description_ + + :return: workflow, with medial surface NaNs applied if medial_surface_nan + """ if medial_surface_nan: # fmt: off workflow.connect([ From 532875eadbce4e1911b960f8fa5c3c1cd8294f71 Mon Sep 17 00:00:00 2001 From: Greg Conan Date: Mon, 5 Dec 2022 13:55:23 -0600 Subject: [PATCH 08/11] STY: black and isort --- nibabies/workflows/bold/resampling.py | 296 +++++++++++++------------- 1 file changed, 151 insertions(+), 145 deletions(-) diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index 3e418d44..d29c38bd 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -11,34 +11,30 @@ """ +from os.path import basename + import nipype.interfaces.workbench as wb +from nipype import Function from nipype.interfaces import freesurfer as fs from nipype.interfaces import fsl from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe -from nipype import Function from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms from niworkflows.interfaces.freesurfer import MedialNaNs -from ...interfaces.volume import CreateSignedDistanceVolume -from os.path import basename - - - from ...config import DEFAULT_MEMORY_MIN_GB +from ...interfaces.volume import CreateSignedDistanceVolume -def init_bold_surf_wf(mem_gb, - surface_spaces, - medial_surface_nan, - project_goodvoxels, - name="bold_surf_wf"): +def init_bold_surf_wf( + mem_gb, surface_spaces, medial_surface_nan, project_goodvoxels, name="bold_surf_wf" +): """ Sample functional images to FreeSurfer surfaces. For each vertex, the cortical ribbon is sampled at six points (spaced 20% of thickness apart) and averaged. - + If --project-goodvoxels is used, a "goodvoxels" BOLD mask, as described in [@hcppipelines], is generated and applied to the functional image before sampling to surface. Outputs are in GIFTI format. @@ -88,11 +84,10 @@ def init_bold_surf_wf(mem_gb, BOLD series, resampled to FreeSurfer surfaces """ + import templateflow as tf from nipype.interfaces.io import FreeSurferSource from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.surf import GiftiSetAnatomicalStructure - import templateflow as tf - workflow = Workflow(name=name) workflow.__desc__ = """\ @@ -103,33 +98,32 @@ def init_bold_surf_wf(mem_gb, out_spaces=", ".join(["*%s*" % s for s in surface_spaces]), ) - if project_goodvoxels: - workflow.__desc__ += """\ + workflow.__desc__ += """\ Before resampling, a "goodvoxels" mask [@hcppipelines] was applied, excluding voxels whose time-series have a locally high coetfficient of -variation, or that lie outside the cortical surfaces, from the +variation, or that lie outside the cortical surfaces, from the surface projection. """ inputnode = pe.Node( niu.IdentityInterface( - fields=["source_file", - "subject_id", - "subjects_dir", - "t1w2fsnative_xfm", - "anat_giftis", - "t1w_mask"] + fields=[ + "source_file", + "subject_id", + "subjects_dir", + "t1w2fsnative_xfm", + "anat_giftis", + "t1w_mask", + ] ), name="inputnode", ) - + itersource = pe.Node(niu.IdentityInterface(fields=["target"]), name="itersource") itersource.iterables = [("target", surface_spaces)] - get_fsnative = pe.Node( - FreeSurferSource(), name="get_fsnative", run_without_submitting=True - ) + get_fsnative = pe.Node(FreeSurferSource(), name="get_fsnative", run_without_submitting=True) def select_target(subject_id, space): """Get the target subject ID, given a source subject ID and a target space.""" @@ -150,9 +144,7 @@ def select_target(subject_id, space): mem_gb=DEFAULT_MEMORY_MIN_GB, ) - itk2lta = pe.Node( - niu.Function(function=_itk2lta), name="itk2lta", run_without_submitting=True - ) + itk2lta = pe.Node(niu.Function(function=_itk2lta), name="itk2lta", run_without_submitting=True) sampler = pe.MapNode( fs.SampleToSurface( @@ -174,8 +166,7 @@ def select_target(subject_id, space): # Refine if medial vertices should be NaNs medial_nans = pe.MapNode( - MedialNaNs(), iterfield=["in_file"], name="medial_nans", - mem_gb=DEFAULT_MEMORY_MIN_GB + MedialNaNs(), iterfield=["in_file"], name="medial_nans", mem_gb=DEFAULT_MEMORY_MIN_GB ) # Rename the source file to the output space to simplify naming later @@ -186,8 +177,8 @@ def select_target(subject_id, space): mem_gb=DEFAULT_MEMORY_MIN_GB, run_without_submitting=True, ) - prepend_hemi.inputs.hemi=["lh", "rh"] - + prepend_hemi.inputs.hemi = ["lh", "rh"] + update_metadata = pe.MapNode( GiftiSetAnatomicalStructure(), iterfield=["in_file"], @@ -204,9 +195,7 @@ def select_target(subject_id, space): ) if project_goodvoxels: - workflow, combine_ribbon_vol_hemis = make_ribbon_file( - workflow, mem_gb, inputnode - ) + workflow, combine_ribbon_vol_hemis = make_ribbon_file(workflow, mem_gb, inputnode) workflow, goodvoxels_mask, ribbon_boldsrc_xfm = mask_outliers( workflow, mem_gb, rename_src, combine_ribbon_vol_hemis ) @@ -214,8 +203,8 @@ def select_target(subject_id, space): workflow, mem_gb, rename_src, goodvoxels_mask, ribbon_boldsrc_xfm ) - # If not applying goodvoxels, then get sampler.source_file from rename_src - # instead in connect_bold_surf_wf + # If not applying goodvoxels, then get sampler.source_file from rename_src + # instead in connect_bold_surf_wf apply_goodvx_or_rename = apply_goodvoxels_ribbon_mask else: apply_goodvx_or_rename = rename_src @@ -226,14 +215,15 @@ def select_target(subject_id, space): rename_src, sampler, itk2lta, update_metadata, apply_goodvx_or_rename ) # fmt: on - return apply_med_surf_nans_if(medial_surface_nan, workflow, - inputnode, sampler, - update_metadata, medial_nans) + return apply_med_surf_nans_if( + medial_surface_nan, workflow, inputnode, sampler, update_metadata, medial_nans + ) def make_ribbon_file(workflow, mem_gb, inputnode): """ - _summary_ + Parameters + ---------- :param workflow: _type_, _description_ :param mem_gb: _type_, _description_ :param inputnode: _type_, _description_ @@ -262,7 +252,7 @@ def make_ribbon_file(workflow, mem_gb, inputnode): name="select_midthick", mem_gb=DEFAULT_MEMORY_MIN_GB, run_without_submitting=True, - ) + ) create_wm_distvol = pe.MapNode( CreateSignedDistanceVolume(), @@ -349,36 +339,37 @@ def make_ribbon_file(workflow, mem_gb, inputnode): ) # make HCP-style ribbon volume in T1w space - workflow.connect([ - (inputnode, select_wm, [("anat_giftis", "inlist")]), - (inputnode, select_pial, [("anat_giftis", "inlist")]), - (inputnode, select_midthick, [("anat_giftis", "inlist")]), - (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surf_file")]), - (inputnode, create_wm_distvol, [("t1w_mask", "ref_file")]), - (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surf_file")]), - (inputnode, create_pial_distvol, [("t1w_mask", "ref_file")]), - (create_wm_distvol, thresh_wm_distvol, [("out_file", "in_file")]), - (create_pial_distvol, uthresh_pial_distvol, [("out_file", "in_file")]), - (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), - (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), - (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), - (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), - (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), - (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), - (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), - ]) + workflow.connect( + [ + (inputnode, select_wm, [("anat_giftis", "inlist")]), + (inputnode, select_pial, [("anat_giftis", "inlist")]), + (inputnode, select_midthick, [("anat_giftis", "inlist")]), + (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surf_file")]), + (inputnode, create_wm_distvol, [("t1w_mask", "ref_file")]), + (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surf_file")]), + (inputnode, create_pial_distvol, [("t1w_mask", "ref_file")]), + (create_wm_distvol, thresh_wm_distvol, [("out_file", "in_file")]), + (create_pial_distvol, uthresh_pial_distvol, [("out_file", "in_file")]), + (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), + (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), + (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), + (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), + (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), + (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), + (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), + ] + ) return workflow, combine_ribbon_vol_hemis def mask_outliers(workflow, mem_gb, rename_src, combine_ribbon_vol_hemis): ribbon_boldsrc_xfm = pe.Node( - ApplyTransforms(interpolation='MultiLabel', - transforms='identity'), + ApplyTransforms(interpolation='MultiLabel', transforms='identity'), name="ribbon_boldsrc_xfm", mem_gb=mem_gb, ) @@ -467,9 +458,7 @@ def _calc_upper_thr(in_stats): upper_thr_val = pe.Node( Function( - input_names=["in_stats"], - output_names=["upper_thresh"], - function=_calc_upper_thr + input_names=["in_stats"], output_names=["upper_thresh"], function=_calc_upper_thr ), name="upper_thr_val", mem_gb=DEFAULT_MEMORY_MIN_GB, @@ -480,9 +469,7 @@ def _calc_lower_thr(in_stats): lower_thr_val = pe.Node( Function( - input_names=["in_stats"], - output_names=["lower_thresh"], - function=_calc_lower_thr + input_names=["in_stats"], output_names=["lower_thresh"], function=_calc_lower_thr ), name="lower_thr_val", mem_gb=DEFAULT_MEMORY_MIN_GB, @@ -534,49 +521,50 @@ def _calc_lower_thr(in_stats): # make HCP-style "goodvoxels" mask in t1w space for filtering outlier voxels # in bold timeseries, based on modulated normalized covariance - workflow.connect([ - (combine_ribbon_vol_hemis, ribbon_boldsrc_xfm, [("out_file", "input_image")]), - (rename_src, stdev_volume, [("out_file", "in_file")]), - (rename_src, mean_volume, [("out_file", "in_file")]), - (mean_volume, ribbon_boldsrc_xfm, [("out_file", "reference_image")]), - (stdev_volume, cov_volume, [("out_file", "in_file")]), - (mean_volume, cov_volume, [("out_file", "operand_file")]), - (cov_volume, cov_ribbon, [("out_file", "in_file")]), - (ribbon_boldsrc_xfm, cov_ribbon, [("output_image", "mask_file")]), - (cov_ribbon, cov_ribbon_mean, [("out_file", "in_file")]), - (cov_ribbon, cov_ribbon_std, [("out_file", "in_file")]), - (cov_ribbon, cov_ribbon_norm, [("out_file", "in_file")]), - (cov_ribbon_mean, cov_ribbon_norm, [("out_stat", "operand_value")]), - (cov_ribbon_norm, smooth_norm, [("out_file", "in_file")]), - (smooth_norm, merge_smooth_norm, [("out_file", "in1")]), - (cov_ribbon_norm, cov_ribbon_norm_smooth, [("out_file", "in_file")]), - (merge_smooth_norm, cov_ribbon_norm_smooth, [("out", "operand_files")]), - (cov_ribbon_mean, cov_norm, [("out_stat", "operand_value")]), - (cov_volume, cov_norm, [("out_file", "in_file")]), - (cov_norm, cov_norm_modulate, [("out_file", "in_file")]), - (cov_ribbon_norm_smooth, cov_norm_modulate, [("out_file", "operand_file")]), - (cov_norm_modulate, cov_norm_modulate_ribbon, [("out_file", "in_file")]), - (ribbon_boldsrc_xfm, cov_norm_modulate_ribbon, [("output_image", "mask_file")]), - (cov_norm_modulate_ribbon, mod_ribbon_mean, [("out_file", "in_file")]), - (cov_norm_modulate_ribbon, mod_ribbon_std, [("out_file", "in_file")]), - (mod_ribbon_mean, merge_mod_ribbon_stats, [("out_stat", "in1")]), - (mod_ribbon_std, merge_mod_ribbon_stats, [("out_stat", "in2")]), - (merge_mod_ribbon_stats, upper_thr_val, [("out", "in_stats")]), - (merge_mod_ribbon_stats, lower_thr_val, [("out", "in_stats")]), - (mean_volume, bin_mean_volume, [("out_file", "in_file")]), - (upper_thr_val, goodvoxels_thr, [("upper_thresh", "thresh")]), - (cov_norm_modulate, goodvoxels_thr, [("out_file", "in_file")]), - (bin_mean_volume, merge_goodvoxels_operands, [("out_file", "in1")]), - (goodvoxels_thr, goodvoxels_mask, [("out_file", "in_file")]), - (merge_goodvoxels_operands, goodvoxels_mask, [("out", "operand_files")]), - ]) + workflow.connect( + [ + (combine_ribbon_vol_hemis, ribbon_boldsrc_xfm, [("out_file", "input_image")]), + (rename_src, stdev_volume, [("out_file", "in_file")]), + (rename_src, mean_volume, [("out_file", "in_file")]), + (mean_volume, ribbon_boldsrc_xfm, [("out_file", "reference_image")]), + (stdev_volume, cov_volume, [("out_file", "in_file")]), + (mean_volume, cov_volume, [("out_file", "operand_file")]), + (cov_volume, cov_ribbon, [("out_file", "in_file")]), + (ribbon_boldsrc_xfm, cov_ribbon, [("output_image", "mask_file")]), + (cov_ribbon, cov_ribbon_mean, [("out_file", "in_file")]), + (cov_ribbon, cov_ribbon_std, [("out_file", "in_file")]), + (cov_ribbon, cov_ribbon_norm, [("out_file", "in_file")]), + (cov_ribbon_mean, cov_ribbon_norm, [("out_stat", "operand_value")]), + (cov_ribbon_norm, smooth_norm, [("out_file", "in_file")]), + (smooth_norm, merge_smooth_norm, [("out_file", "in1")]), + (cov_ribbon_norm, cov_ribbon_norm_smooth, [("out_file", "in_file")]), + (merge_smooth_norm, cov_ribbon_norm_smooth, [("out", "operand_files")]), + (cov_ribbon_mean, cov_norm, [("out_stat", "operand_value")]), + (cov_volume, cov_norm, [("out_file", "in_file")]), + (cov_norm, cov_norm_modulate, [("out_file", "in_file")]), + (cov_ribbon_norm_smooth, cov_norm_modulate, [("out_file", "operand_file")]), + (cov_norm_modulate, cov_norm_modulate_ribbon, [("out_file", "in_file")]), + (ribbon_boldsrc_xfm, cov_norm_modulate_ribbon, [("output_image", "mask_file")]), + (cov_norm_modulate_ribbon, mod_ribbon_mean, [("out_file", "in_file")]), + (cov_norm_modulate_ribbon, mod_ribbon_std, [("out_file", "in_file")]), + (mod_ribbon_mean, merge_mod_ribbon_stats, [("out_stat", "in1")]), + (mod_ribbon_std, merge_mod_ribbon_stats, [("out_stat", "in2")]), + (merge_mod_ribbon_stats, upper_thr_val, [("out", "in_stats")]), + (merge_mod_ribbon_stats, lower_thr_val, [("out", "in_stats")]), + (mean_volume, bin_mean_volume, [("out_file", "in_file")]), + (upper_thr_val, goodvoxels_thr, [("upper_thresh", "thresh")]), + (cov_norm_modulate, goodvoxels_thr, [("out_file", "in_file")]), + (bin_mean_volume, merge_goodvoxels_operands, [("out_file", "in1")]), + (goodvoxels_thr, goodvoxels_mask, [("out_file", "in_file")]), + (merge_goodvoxels_operands, goodvoxels_mask, [("out", "operand_files")]), + ] + ) return workflow, goodvoxels_mask, ribbon_boldsrc_xfm -def project_vol_to_surf(workflow, mem_gb, rename_src, goodvoxels_mask, - ribbon_boldsrc_xfm): - """ +def project_vol_to_surf(workflow, mem_gb, rename_src, goodvoxels_mask, ribbon_boldsrc_xfm): + """ Parameters ---------- workflow : :obj: @@ -605,29 +593,32 @@ def project_vol_to_surf(workflow, mem_gb, rename_src, goodvoxels_mask, ) # apply goodvoxels ribbon mask to bold - workflow.connect([ - (goodvoxels_mask, goodvoxels_ribbon_mask, [("out_file", "in_file")]), - (ribbon_boldsrc_xfm, goodvoxels_ribbon_mask, [("output_image", "mask_file")]), - (goodvoxels_ribbon_mask, apply_goodvoxels_ribbon_mask, [("out_file", "mask_file")]), - (rename_src, apply_goodvoxels_ribbon_mask, [("out_file", "in_file")]), - ]) + workflow.connect( + [ + (goodvoxels_mask, goodvoxels_ribbon_mask, [("out_file", "in_file")]), + (ribbon_boldsrc_xfm, goodvoxels_ribbon_mask, [("output_image", "mask_file")]), + (goodvoxels_ribbon_mask, apply_goodvoxels_ribbon_mask, [("out_file", "mask_file")]), + (rename_src, apply_goodvoxels_ribbon_mask, [("out_file", "in_file")]), + ] + ) return workflow, apply_goodvoxels_ribbon_mask -def apply_med_surf_nans_if(medial_surface_nan, workflow, inputnode, sampler, - update_metadata, medial_nans): +def apply_med_surf_nans_if( + medial_surface_nan, workflow, inputnode, sampler, update_metadata, medial_nans +): """ Parameters ---------- medial_surface_nan : :obj:`bool` - Replace medial wall values with NaNs on functional GIFTI files + Replace medial wall values with NaNs on functional GIFTI files workflow : _type_, _description_ inputnode : _type_, _description_ sampler : _type_, _description_ update_metadata : _type_, _description_ medial_nans : _type_, _description_ - + :return: workflow, with medial surface NaNs applied if medial_surface_nan """ if medial_surface_nan: @@ -643,27 +634,40 @@ def apply_med_surf_nans_if(medial_surface_nan, workflow, inputnode, sampler, return workflow -def connect_bold_surf_wf(workflow, inputnode, outputnode, get_fsnative, - itersource, targets, rename_src, sampler, itk2lta, - update_metadata, apply_goodvoxels_or_rename): - workflow.connect([ - (inputnode, get_fsnative, [("subject_id", "subject_id"), - ("subjects_dir", "subjects_dir")]), - (inputnode, targets, [("subject_id", "subject_id")]), - (inputnode, rename_src, [("source_file", "in_file")]), - (inputnode, itk2lta, [("source_file", "src_file"), - ("t1w2fsnative_xfm", "in_file")]), - (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 - (inputnode, sampler, [("subjects_dir", "subjects_dir"), - ("subject_id", "subject_id")]), - (itersource, targets, [("target", "space")]), - (itersource, rename_src, [("target", "subject")]), - (itk2lta, sampler, [("out", "reg_file")]), - (targets, sampler, [("out", "target_subject")]), - (apply_goodvoxels_or_rename, sampler, [("out_file", "source_file")]), - (update_metadata, outputnode, [("out_file", "surfaces")]), - (itersource, outputnode, [("target", "target")]), - ]) +def connect_bold_surf_wf( + workflow, + inputnode, + outputnode, + get_fsnative, + itersource, + targets, + rename_src, + sampler, + itk2lta, + update_metadata, + apply_goodvoxels_or_rename, +): + workflow.connect( + [ + ( + inputnode, + get_fsnative, + [("subject_id", "subject_id"), ("subjects_dir", "subjects_dir")], + ), + (inputnode, targets, [("subject_id", "subject_id")]), + (inputnode, rename_src, [("source_file", "in_file")]), + (inputnode, itk2lta, [("source_file", "src_file"), ("t1w2fsnative_xfm", "in_file")]), + (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 + (inputnode, sampler, [("subjects_dir", "subjects_dir"), ("subject_id", "subject_id")]), + (itersource, targets, [("target", "space")]), + (itersource, rename_src, [("target", "subject")]), + (itk2lta, sampler, [("out", "reg_file")]), + (targets, sampler, [("out", "target_subject")]), + (apply_goodvoxels_or_rename, sampler, [("out_file", "source_file")]), + (update_metadata, outputnode, [("out_file", "surfaces")]), + (itersource, outputnode, [("target", "target")]), + ] + ) return workflow @@ -1193,7 +1197,7 @@ def init_bold_grayords_wf(grayord_density, mem_gb, repetition_time, name="bold_g str(tf.api.get("fsaverage", hemi=hemi, density="164k", desc="std", suffix="sphere")) for hemi in "LR" ] - + resample.inputs.current_area = [ str( tf.api.get("fsaverage", hemi=hemi, density="164k", desc="vaavg", suffix="midthickness") @@ -1345,6 +1349,8 @@ def _itk2lta(in_file, src_file, dst_file): ).to_filename(out_file, moving=dst_file, fmt="fs") return str(out_file) + def _sorted_by_basename(inlist): from os.path import basename + return sorted(inlist, key=lambda x: str(basename(x))) From 9676f14828f158a728ef5cb22ec11d050bb0d0f1 Mon Sep 17 00:00:00 2001 From: Greg Conan Date: Mon, 5 Dec 2022 19:25:58 -0600 Subject: [PATCH 09/11] Updated Docker version tag --- .circleci/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index c9beade3..0f39979a 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -64,7 +64,7 @@ _check_skip_job: &check_skip_job version: 2.1 orbs: - docker: circleci/docker@1.6.0 + docker: circleci/docker@2.1.4 jobs: From e995ec4c029645dd05444c2d181abdf882bd65da Mon Sep 17 00:00:00 2001 From: Greg Conan Date: Fri, 9 Dec 2022 16:00:08 -0600 Subject: [PATCH 10/11] RF: Moved ribbon wf (from bold to anat) and reformatted --- nibabies/interfaces/volume.py | 133 ---------- nibabies/interfaces/workbench.py | 125 ++++++++++ nibabies/workflows/anatomical/base.py | 18 ++ nibabies/workflows/anatomical/outputs.py | 9 + nibabies/workflows/anatomical/surfaces.py | 172 ++++++++++++- nibabies/workflows/base.py | 2 +- nibabies/workflows/bold/base.py | 5 +- nibabies/workflows/bold/resampling.py | 291 ++++++---------------- 8 files changed, 397 insertions(+), 358 deletions(-) delete mode 100644 nibabies/interfaces/volume.py diff --git a/nibabies/interfaces/volume.py b/nibabies/interfaces/volume.py deleted file mode 100644 index 4c6187b3..00000000 --- a/nibabies/interfaces/volume.py +++ /dev/null @@ -1,133 +0,0 @@ -# -*- coding: utf-8 -*- -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -"""This module provides interfaces for workbench volume commands""" -import os - -from nipype.interfaces.base import TraitedSpec, File, traits, CommandLineInputSpec -from nipype.interfaces.workbench.base import WBCommand -from nipype import logging - -iflogger = logging.getLogger("nipype.interface") - - -class CreateSignedDistanceVolumeInputSpec(CommandLineInputSpec): - surf_file = File( - exists=True, - mandatory=True, - argstr="%s", - position=0, - desc="Input surface GIFTI file (.surf.gii)", - ) - ref_file = File( - exists=True, - mandatory=True, - argstr="%s", - position=1, - desc="NIfTI volume in the desired output space (dims, spacing, origin)", - ) - out_file = File( - name_source=["surface"], - name_template="%s_distvol.nii.gz", - argstr="%s", - position=2, - desc="Name of output volume containing signed distances", - ) - out_mask = File( - name_source=["surface"], - name_template="%s_distmask.nii.gz", - argstr="-roi-out %s", - desc="Name of file to store a mask where the ``out_file`` has a computed value", - ) - fill_value = traits.Float( - mandatory=False, - default=0., - usedefault=True, - argstr="-fill-value %f", - desc="value to put in all voxels that don't get assigned a distance", - ) - exact_limit = traits.Float( - default=5., - usedefault=True, - argstr="-exact-limit %f", - desc="distance for exact output in mm", - ) - approx_limit = traits.Float( - default=20., - usedefault=True, - argstr="-approx-limit %f", - desc="distance for approximate output in mm", - ) - approx_neighborhood = traits.Int( - default=2, - usedefault=True, - argstr="-approx-neighborhood %d", - desc="size of neighborhood cube measured from center to face in voxels, default 2 = 5x5x5", - ) - winding_method = traits.Enum( - "EVEN_ODD", - "NEGATIVE", - "NONZERO", - "NORMALS", - argstr="-winding %s", - usedefault=True, - desc="winding method for point inside surface test" - ) - -class CreateSignedDistanceVolumeOutputSpec(TraitedSpec): - out_file = File(desc="Name of output volume containing signed distances") - out_mask = File(desc="Name of file to store a mask where the ``out_file`` has a computed value") - - -class CreateSignedDistanceVolume(WBCommand): - """CREATE SIGNED DISTANCE VOLUME FROM SURFACE - wb_command -create-signed-distance-volume - - the input surface - - a volume in the desired output space (dims, spacing, origin) - - output - the output volume - - [-roi-out] - output an roi volume of where the output has a computed - value - - output - the output roi volume - - [-fill-value] - specify a value to put in all voxels that don't get - assigned a distance - - value to fill with (default 0) - - [-exact-limit] - specify distance for exact output - - distance in mm (default 5) - - [-approx-limit] - specify distance for approximate output - - distance in mm (default 20) - - [-approx-neighborhood] - voxel neighborhood for approximate calculation - - size of neighborhood cube measured from center to face, in - voxels (default 2 = 5x5x5) - - [-winding] - winding method for point inside surface test - - name of the method (default EVEN_ODD) - - Computes the signed distance function of the surface. Exact distance is - calculated by finding the closest point on any surface triangle to the - center of the voxel. Approximate distance is calculated starting with - these distances, using dijkstra's method with a neighborhood of voxels. - Specifying too small of an exact distance may produce unexpected results. - Valid specifiers for winding methods are as follows: - - EVEN_ODD (default) - NEGATIVE - NONZERO - NORMALS - - The NORMALS method uses the normals of triangles and edges, or the - closest triangle hit by a ray from the point. This method may be - slightly faster, but is only reliable for a closed surface that does not - cross through itself. All other methods count entry (positive) and exit - (negative) crossings of a vertical ray from the point, then counts as - inside if the total is odd, negative, or nonzero, respectively. - - """ - - input_spec = CreateSignedDistanceVolumeInputSpec - output_spec = CreateSignedDistanceVolumeOutputSpec - _cmd = "wb_command -create-signed-distance-volume" \ No newline at end of file diff --git a/nibabies/interfaces/workbench.py b/nibabies/interfaces/workbench.py index 2bdbe359..238cfc4c 100644 --- a/nibabies/interfaces/workbench.py +++ b/nibabies/interfaces/workbench.py @@ -1445,3 +1445,128 @@ class VolumeLabelImport(WBCommand): input_spec = VolumeLabelImportInputSpec output_spec = VolumeLabelImportOutputSpec _cmd = "wb_command -volume-label-import" + + +class CreateSignedDistanceVolumeInputSpec(CommandLineInputSpec): + surf_file = File( + exists=True, + mandatory=True, + argstr="%s", + position=0, + desc="Input surface GIFTI file (.surf.gii)", + ) + ref_file = File( + exists=True, + mandatory=True, + argstr="%s", + position=1, + desc="NIfTI volume in the desired output space (dims, spacing, origin)", + ) + out_file = File( + name_source=["surf_file"], + name_template="%s_distvol.nii.gz", + argstr="%s", + position=2, + desc="Name of output volume containing signed distances", + ) + out_mask = File( + name_source=["surf_file"], + name_template="%s_distmask.nii.gz", + argstr="-roi-out %s", + desc="Name of file to store a mask where the ``out_file`` has a computed value", + ) + fill_value = traits.Float( + mandatory=False, + default=0.0, + usedefault=True, + argstr="-fill-value %f", + desc="value to put in all voxels that don't get assigned a distance", + ) + exact_limit = traits.Float( + default=5.0, + usedefault=True, + argstr="-exact-limit %f", + desc="distance for exact output in mm", + ) + approx_limit = traits.Float( + default=20.0, + usedefault=True, + argstr="-approx-limit %f", + desc="distance for approximate output in mm", + ) + approx_neighborhood = traits.Int( + default=2, + usedefault=True, + argstr="-approx-neighborhood %d", + desc="size of neighborhood cube measured from center to face in voxels, default 2 = 5x5x5", + ) + winding_method = traits.Enum( + "EVEN_ODD", + "NEGATIVE", + "NONZERO", + "NORMALS", + argstr="-winding %s", + usedefault=True, + desc="winding method for point inside surface test", + ) + + +class CreateSignedDistanceVolumeOutputSpec(TraitedSpec): + out_file = File(desc="Name of output volume containing signed distances") + out_mask = File( + desc="Name of file to store a mask where the ``out_file`` has a computed value" + ) + + +class CreateSignedDistanceVolume(WBCommand): + """CREATE SIGNED DISTANCE VOLUME FROM SURFACE + wb_command -create-signed-distance-volume + - the input surface + - a volume in the desired output space (dims, spacing, origin) + - output - the output volume + + [-roi-out] - output an roi volume of where the output has a computed + value + - output - the output roi volume + + [-fill-value] - specify a value to put in all voxels that don't get + assigned a distance + - value to fill with (default 0) + + [-exact-limit] - specify distance for exact output + - distance in mm (default 5) + + [-approx-limit] - specify distance for approximate output + - distance in mm (default 20) + + [-approx-neighborhood] - voxel neighborhood for approximate calculation + - size of neighborhood cube measured from center to face, in + voxels (default 2 = 5x5x5) + + [-winding] - winding method for point inside surface test + - name of the method (default EVEN_ODD) + + Computes the signed distance function of the surface. Exact distance is + calculated by finding the closest point on any surface triangle to the + center of the voxel. Approximate distance is calculated starting with + these distances, using dijkstra's method with a neighborhood of voxels. + Specifying too small of an exact distance may produce unexpected results. + Valid specifiers for winding methods are as follows: + + EVEN_ODD (default) + NEGATIVE + NONZERO + NORMALS + + The NORMALS method uses the normals of triangles and edges, or the + closest triangle hit by a ray from the point. This method may be + slightly faster, but is only reliable for a closed surface that does not + cross through itself. All other methods count entry (positive) and exit + (negative) crossings of a vertical ray from the point, then counts as + inside if the total is odd, negative, or nonzero, respectively. + + """ + + input_spec = CreateSignedDistanceVolumeInputSpec + output_spec = CreateSignedDistanceVolumeOutputSpec + _cmd = "wb_command -create-signed-distance-volume" diff --git a/nibabies/workflows/anatomical/base.py b/nibabies/workflows/anatomical/base.py index aeb1fbed..aa7c2fba 100644 --- a/nibabies/workflows/anatomical/base.py +++ b/nibabies/workflows/anatomical/base.py @@ -92,6 +92,7 @@ def init_infant_anat_wf( from .preproc import init_anat_average_wf from .registration import init_coregistration_wf from .segmentation import init_anat_segmentations_wf + from .surfaces import init_anat_ribbon_wf # for now, T1w only num_t1w = len(t1w) if t1w else 0 @@ -143,6 +144,7 @@ def init_infant_anat_wf( "surfaces", "anat_aseg", "anat_aparc", + "anat_ribbon", "template", ] ), @@ -250,6 +252,10 @@ def init_infant_anat_wf( templates=spaces.get_spaces(nonstandard=False, dim=(3,)), ) + # Anatomical ribbon file using HCP signed-distance volume method + # if config.workflow.project_goodvoxels: + anat_ribbon_wf = init_anat_ribbon_wf() + # fmt:off wf.connect([ (inputnode, t1w_template_wf, [("t1w", "inputnode.in_files")]), @@ -354,6 +360,9 @@ def init_infant_anat_wf( ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"), ("outputnode.std2anat_xfm", "inputnode.std2anat_xfm"), ]), + (anat_ribbon_wf, anat_derivatives_wf, [ + ("outputnode.anat_ribbon", "inputnode.anat_ribbon"), + ]), (anat_seg_wf, anat_derivatives_wf, [ ("outputnode.anat_dseg", "inputnode.t1w_dseg"), ("outputnode.anat_tpms", "inputnode.t1w_tpms"), @@ -404,6 +413,15 @@ def init_infant_anat_wf( ("outputnode.out_aparc", "anat_aparc"), ("outputnode.out_aseg", "anat_aseg"), ]), + (t1w_preproc_wf, anat_ribbon_wf, [ + ("outputnode.t1w_mask", "inputnode.t1w_mask"), + ]) + (surface_recon_wf, anat_ribbon_wf, [ + ("outputnode.surfaces", "inputnode.surfaces"), + ]), + (anat_ribbon_wf, outputnode, [ + ("outputnode.anat_ribbon", "anat_ribbon") + ]), (surface_recon_wf, anat_reports_wf, [ ("outputnode.subject_id", "inputnode.subject_id"), ("outputnode.subjects_dir", "inputnode.subjects_dir"), diff --git a/nibabies/workflows/anatomical/outputs.py b/nibabies/workflows/anatomical/outputs.py index 5020f197..e25bae37 100644 --- a/nibabies/workflows/anatomical/outputs.py +++ b/nibabies/workflows/anatomical/outputs.py @@ -338,6 +338,7 @@ def init_anat_derivatives_wf( "surfaces", "t1w_fs_aseg", "t1w_fs_aparc", + "anat_ribbon", ] ), name="inputnode", @@ -353,6 +354,12 @@ def init_anat_derivatives_wf( ) ds_t1w_preproc.inputs.SkullStripped = False + ds_anat_ribbon = pe.Node( + DerivativesDataSink(base_directory=output_dir, suffix="ribbon", compress=True), + name="ds_anat_ribbon", + run_without_submitting=True, + ) + ds_t1w_mask = pe.Node( DerivativesDataSink(base_directory=output_dir, desc="brain", suffix="mask", compress=True), name="ds_t1w_mask", @@ -384,6 +391,8 @@ def init_anat_derivatives_wf( ('source_files', 'source_file')]), (inputnode, ds_t1w_dseg, [('t1w_dseg', 'in_file'), ('source_files', 'source_file')]), + (inputnode, ds_anat_ribbon, [('anat_ribbon', 'in_file'), + ('source_files', 'source_file')]), (raw_sources, ds_t1w_mask, [('out', 'RawSources')]), ]) # fmt: on diff --git a/nibabies/workflows/anatomical/surfaces.py b/nibabies/workflows/anatomical/surfaces.py index 1d53e67b..4f3a94e5 100644 --- a/nibabies/workflows/anatomical/surfaces.py +++ b/nibabies/workflows/anatomical/surfaces.py @@ -1,11 +1,16 @@ # Use infant_recon_all to generate subcortical segmentations and cortical parcellations +from nipype.interfaces import fsl +from nipype.interfaces import utility as niu +from nipype.pipeline import engine as pe + +from ...config import DEFAULT_MEMORY_MIN_GB +from ...interfaces.workbench import CreateSignedDistanceVolume + def init_infant_surface_recon_wf(*, age_months, use_aseg=False, name="infant_surface_recon_wf"): from nipype.interfaces import freesurfer as fs from nipype.interfaces import io as nio - from nipype.interfaces import utility as niu - from nipype.pipeline import engine as pe from niworkflows.engine.workflows import LiterateWorkflow from niworkflows.interfaces.freesurfer import PatchedLTAConvert as LTAConvert from niworkflows.interfaces.freesurfer import ( @@ -135,6 +140,163 @@ def init_infant_surface_recon_wf(*, age_months, use_aseg=False, name="infant_sur return wf +def init_anat_ribbon_wf(name="anat_ribbon_wf"): + # 0, 1 = wm; 2, 3 = pial; 6, 7 = mid + # note that order of lh / rh within each surf type is not guaranteed due to use + # of unsorted glob by FreeSurferSource prior, but we can do a sort + # to ensure consistent ordering + workflow = pe.Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "surfaces", # anat_giftis, + "t1w_mask", + ] + ), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "anat_ribbon", + ] + ), + name="outputnode", + ) + + select_wm = pe.Node( + niu.Select(index=[0, 1]), + name="select_wm", + mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, + ) + + select_pial = pe.Node( + niu.Select(index=[2, 3]), + name="select_pial", + mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, + ) + + select_midthick = pe.Node( + niu.Select(index=[6, 7]), + name="select_midthick", + mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, + ) + + create_wm_distvol = pe.MapNode( + CreateSignedDistanceVolume(), + iterfield=["surf_file"], + name="create_wm_distvol", + ) + + create_pial_distvol = pe.MapNode( + CreateSignedDistanceVolume(), + iterfield=["surf_file"], + name="create_pial_distvol", + ) + + thresh_wm_distvol = pe.MapNode( + fsl.maths.MathsCommand(args="-thr 0 -bin -mul 255"), + iterfield=["in_file"], + name="thresh_wm_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + uthresh_pial_distvol = pe.MapNode( + fsl.maths.MathsCommand(args="-uthr 0 -abs -bin -mul 255"), + iterfield=["in_file"], + name="uthresh_pial_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + bin_wm_distvol = pe.MapNode( + fsl.maths.UnaryMaths(operation="bin"), + iterfield=["in_file"], + name="bin_wm_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + bin_pial_distvol = pe.MapNode( + fsl.maths.UnaryMaths(operation="bin"), + iterfield=["in_file"], + name="bin_pial_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + split_wm_distvol = pe.Node( + niu.Split(splits=[1, 1]), + name="split_wm_distvol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, + ) + + merge_wm_distvol_no_flatten = pe.Node( + niu.Merge(2), + no_flatten=True, + name="merge_wm_distvol_no_flatten", + mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, + ) + + make_ribbon_vol = pe.MapNode( + fsl.maths.MultiImageMaths(op_string="-mas %s -mul 255"), + iterfield=["in_file", "operand_files"], + name="make_ribbon_vol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + bin_ribbon_vol = pe.MapNode( + fsl.maths.UnaryMaths(operation="bin"), + iterfield=["in_file"], + name="bin_ribbon_vol", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + split_squeeze_ribbon_vol = pe.Node( + niu.Split(splits=[1, 1], squeeze=True), + name="split_squeeze_ribbon", + mem_gb=DEFAULT_MEMORY_MIN_GB, + run_without_submitting=True, + ) + + combine_ribbon_vol_hemis = pe.Node( + fsl.maths.BinaryMaths(operation="add"), + name="combine_ribbon_vol_hemis", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + # make HCP-style ribbon volume in T1w space + workflow.connect( + [ + (inputnode, select_wm, [("surfaces", "inlist")]), + (inputnode, select_pial, [("surfaces", "inlist")]), + (inputnode, select_midthick, [("surfaces", "inlist")]), + (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surf_file")]), + (inputnode, create_wm_distvol, [("t1w_mask", "ref_file")]), + (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surf_file")]), + (inputnode, create_pial_distvol, [("t1w_mask", "ref_file")]), + (create_wm_distvol, thresh_wm_distvol, [("out_file", "in_file")]), + (create_pial_distvol, uthresh_pial_distvol, [("out_file", "in_file")]), + (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), + (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), + (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), + (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), + (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), + (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), + (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), + (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), + (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), + (combine_ribbon_vol_hemis, outputnode, [("out_file", "anat_ribbon")]), + ] + ) + return workflow + + def _parent(p): from pathlib import Path @@ -151,3 +313,9 @@ def _gen_recon_dir(subjects_dir, subject_id): def _replace_mgz(in_file): return in_file.replace('.mgz', '.nii.gz') + + +def _sorted_by_basename(inlist): + from os.path import basename + + return sorted(inlist, key=lambda x: str(basename(x))) diff --git a/nibabies/workflows/base.py b/nibabies/workflows/base.py index 5ac91064..a72812ec 100644 --- a/nibabies/workflows/base.py +++ b/nibabies/workflows/base.py @@ -478,7 +478,7 @@ def init_single_subject_wf(subject_id, session_id=None): ('outputnode.subject_id', 'inputnode.subject_id'), ('outputnode.t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'), ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), - ('outputnode.surfaces', 'inputnode.anat_giftis'), + ('outputnode.anat_ribbon', 'inputnode.anat_ribbon'), ]), ]) # fmt: on diff --git a/nibabies/workflows/bold/base.py b/nibabies/workflows/bold/base.py index 661c89e1..7c120122 100644 --- a/nibabies/workflows/bold/base.py +++ b/nibabies/workflows/bold/base.py @@ -305,7 +305,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non Non-gridded (surface) resamplings were performed using `mri_vol2surf` (FreeSurfer) """ - inputnode = pe.Node( niu.IdentityInterface( fields=[ @@ -321,7 +320,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non "anat2std_xfm", "std2anat_xfm", "template", - "anat_giftis", + "anat_ribbon", # from bold reference workflow "bold_ref", "bold_ref_xfm", @@ -861,7 +860,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False, existing_derivatives=Non ('subjects_dir', 'inputnode.subjects_dir'), ('subject_id', 'inputnode.subject_id'), ('t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'), - ("anat_giftis", "inputnode.anat_giftis"), + ("anat_ribbon", "inputnode.anat_ribbon"), ("anat_mask", "inputnode.t1w_mask")]), (bold_t1_trans_wf, bold_surf_wf, [('outputnode.bold_t1', 'inputnode.source_file')]), (bold_surf_wf, outputnode, [('outputnode.surfaces', 'surfaces')]), diff --git a/nibabies/workflows/bold/resampling.py b/nibabies/workflows/bold/resampling.py index d29c38bd..c6f39690 100644 --- a/nibabies/workflows/bold/resampling.py +++ b/nibabies/workflows/bold/resampling.py @@ -23,7 +23,6 @@ from niworkflows.interfaces.freesurfer import MedialNaNs from ...config import DEFAULT_MEMORY_MIN_GB -from ...interfaces.volume import CreateSignedDistanceVolume def init_bold_surf_wf( @@ -101,7 +100,7 @@ def init_bold_surf_wf( if project_goodvoxels: workflow.__desc__ += """\ Before resampling, a "goodvoxels" mask [@hcppipelines] was applied, -excluding voxels whose time-series have a locally high coetfficient of +excluding voxels whose time-series have a locally high coefficient of variation, or that lie outside the cortical surfaces, from the surface projection. """ @@ -113,7 +112,7 @@ def init_bold_surf_wf( "subject_id", "subjects_dir", "t1w2fsnative_xfm", - "anat_giftis", + "anat_ribbon", "t1w_mask", ] ), @@ -195,179 +194,82 @@ def select_target(subject_id, space): ) if project_goodvoxels: - workflow, combine_ribbon_vol_hemis = make_ribbon_file(workflow, mem_gb, inputnode) - workflow, goodvoxels_mask, ribbon_boldsrc_xfm = mask_outliers( - workflow, mem_gb, rename_src, combine_ribbon_vol_hemis - ) - workflow, apply_goodvoxels_ribbon_mask = project_vol_to_surf( - workflow, mem_gb, rename_src, goodvoxels_mask, ribbon_boldsrc_xfm - ) + goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb) - # If not applying goodvoxels, then get sampler.source_file from rename_src - # instead in connect_bold_surf_wf - apply_goodvx_or_rename = apply_goodvoxels_ribbon_mask + # fmt: off + workflow.connect([ + (inputnode, goodvoxels_bold_mask_wf, [("source_file", "inputnode.bold_file"), + ("anat_ribbon", "inputnode.anat_ribbon")]), + (goodvoxels_bold_mask_wf, rename_src, [("outputnode.masked_bold", "in_file")]), + ]) + # fmt: on else: - apply_goodvx_or_rename = rename_src + # fmt: off + workflow.connect([ + (inputnode, rename_src, [("source_file", "in_file")]), + ]) + # fmt: on # fmt: off - workflow = connect_bold_surf_wf( - workflow, inputnode, outputnode, get_fsnative, itersource, targets, - rename_src, sampler, itk2lta, update_metadata, apply_goodvx_or_rename + workflow.connect( + [ + ( + inputnode, + get_fsnative, + [("subject_id", "subject_id"), ("subjects_dir", "subjects_dir")], + ), + (inputnode, targets, [("subject_id", "subject_id")]), + (inputnode, itk2lta, [("source_file", "src_file"), ("t1w2fsnative_xfm", "in_file")]), + (get_fsnative, itk2lta, [("brain", "dst_file")]), # InfantFS: Use brain instead of T1 + (inputnode, sampler, [("subjects_dir", "subjects_dir"), ("subject_id", "subject_id")]), + (itersource, targets, [("target", "space")]), + (itersource, rename_src, [("target", "subject")]), + (rename_src, sampler, [("out_file", "source_file")]), + (itk2lta, sampler, [("out", "reg_file")]), + (targets, sampler, [("out", "target_subject")]), + (update_metadata, outputnode, [("out_file", "surfaces")]), + (itersource, outputnode, [("target", "target")]), + ] ) # fmt: on - return apply_med_surf_nans_if( - medial_surface_nan, workflow, inputnode, sampler, update_metadata, medial_nans - ) - - -def make_ribbon_file(workflow, mem_gb, inputnode): - """ - Parameters - ---------- - :param workflow: _type_, _description_ - :param mem_gb: _type_, _description_ - :param inputnode: _type_, _description_ - :return: _type_, _description_ - """ - # 0, 1 = wm; 2, 3 = pial; 6, 7 = mid - # note that order of lh / rh within each surf type is not guaranteed due to use - # of unsorted glob by FreeSurferSource prior, but we can do a sort - # to ensure consistent ordering - select_wm = pe.Node( - niu.Select(index=[0, 1]), - name="select_wm", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) - select_pial = pe.Node( - niu.Select(index=[2, 3]), - name="select_pial", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) - - select_midthick = pe.Node( - niu.Select(index=[6, 7]), - name="select_midthick", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) - - create_wm_distvol = pe.MapNode( - CreateSignedDistanceVolume(), - iterfield=["surf_file"], - name="create_wm_distvol", - mem_gb=mem_gb, - ) - - create_pial_distvol = pe.MapNode( - CreateSignedDistanceVolume(), - iterfield=["surf_file"], - name="create_pial_distvol", - mem_gb=mem_gb, - ) - - thresh_wm_distvol = pe.MapNode( - fsl.maths.MathsCommand(args="-thr 0 -bin -mul 255"), - iterfield=["in_file"], - name="thresh_wm_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - - uthresh_pial_distvol = pe.MapNode( - fsl.maths.MathsCommand(args="-uthr 0 -abs -bin -mul 255"), - iterfield=["in_file"], - name="uthresh_pial_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - - bin_wm_distvol = pe.MapNode( - fsl.maths.UnaryMaths(operation="bin"), - iterfield=["in_file"], - name="bin_wm_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - - bin_pial_distvol = pe.MapNode( - fsl.maths.UnaryMaths(operation="bin"), - iterfield=["in_file"], - name="bin_pial_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - - split_wm_distvol = pe.Node( - niu.Split(splits=[1, 1]), - name="split_wm_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) - - merge_wm_distvol_no_flatten = pe.Node( - niu.Merge(2), - no_flatten=True, - name="merge_wm_distvol_no_flatten", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) + if medial_surface_nan: + # fmt: off + workflow.connect([ + (inputnode, medial_nans, [("subjects_dir", "subjects_dir")]), + (sampler, medial_nans, [("out_file", "in_file")]), + (medial_nans, update_metadata, [("out_file", "in_file")]), + ]) + # fmt: on + else: + # fmt: off + workflow.connect(sampler, "out_file", update_metadata, "in_file") + # fmt: on - make_ribbon_vol = pe.MapNode( - fsl.maths.MultiImageMaths(op_string="-mas %s -mul 255"), - iterfield=["in_file", "operand_files"], - name="make_ribbon_vol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) + return workflow - bin_ribbon_vol = pe.MapNode( - fsl.maths.UnaryMaths(operation="bin"), - iterfield=["in_file"], - name="bin_ribbon_vol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - split_squeeze_ribbon_vol = pe.Node( - niu.Split(splits=[1, 1], squeeze=True), - name="split_squeeze_ribbon", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) +def init_goodvoxels_bold_mask_wf(mem_gb, name="goodvoxels_bold_mask_wf"): + workflow = pe.Workflow(name=name) - combine_ribbon_vol_hemis = pe.Node( - fsl.maths.BinaryMaths(operation="add"), - name="combine_ribbon_vol_hemis", - mem_gb=DEFAULT_MEMORY_MIN_GB, + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "anat_ribbon", + "bold_file", + ] + ), + name="inputnode", ) - - # make HCP-style ribbon volume in T1w space - workflow.connect( - [ - (inputnode, select_wm, [("anat_giftis", "inlist")]), - (inputnode, select_pial, [("anat_giftis", "inlist")]), - (inputnode, select_midthick, [("anat_giftis", "inlist")]), - (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surf_file")]), - (inputnode, create_wm_distvol, [("t1w_mask", "ref_file")]), - (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surf_file")]), - (inputnode, create_pial_distvol, [("t1w_mask", "ref_file")]), - (create_wm_distvol, thresh_wm_distvol, [("out_file", "in_file")]), - (create_pial_distvol, uthresh_pial_distvol, [("out_file", "in_file")]), - (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), - (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), - (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), - (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), - (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), - (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), - (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), - ] + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "masked_bold", + "goodvoxels_ribbon", + ] + ), + name="outputnode", ) - - return workflow, combine_ribbon_vol_hemis - - -def mask_outliers(workflow, mem_gb, rename_src, combine_ribbon_vol_hemis): ribbon_boldsrc_xfm = pe.Node( ApplyTransforms(interpolation='MultiLabel', transforms='identity'), name="ribbon_boldsrc_xfm", @@ -523,9 +425,9 @@ def _calc_lower_thr(in_stats): # in bold timeseries, based on modulated normalized covariance workflow.connect( [ - (combine_ribbon_vol_hemis, ribbon_boldsrc_xfm, [("out_file", "input_image")]), - (rename_src, stdev_volume, [("out_file", "in_file")]), - (rename_src, mean_volume, [("out_file", "in_file")]), + (inputnode, ribbon_boldsrc_xfm, [("anat_ribbon", "input_image")]), + (inputnode, stdev_volume, [("bold_file", "in_file")]), + (inputnode, mean_volume, [("bold_file", "in_file")]), (mean_volume, ribbon_boldsrc_xfm, [("out_file", "reference_image")]), (stdev_volume, cov_volume, [("out_file", "in_file")]), (mean_volume, cov_volume, [("out_file", "operand_file")]), @@ -560,22 +462,6 @@ def _calc_lower_thr(in_stats): ] ) - return workflow, goodvoxels_mask, ribbon_boldsrc_xfm - - -def project_vol_to_surf(workflow, mem_gb, rename_src, goodvoxels_mask, ribbon_boldsrc_xfm): - """ - Parameters - ---------- - workflow : :obj: - mem_gb : :obj:`float` - Size of BOLD file in GB - rename_src: _type_, _description_ - goodvoxels_mask: _type_, _description_ - ribbon_boldsrc_xfm: _type_, _description_ - - :return: _type_, _description_ - """ goodvoxels_ribbon_mask = pe.Node( fsl.ApplyMask(), name_source=['in_file'], @@ -598,40 +484,13 @@ def project_vol_to_surf(workflow, mem_gb, rename_src, goodvoxels_mask, ribbon_bo (goodvoxels_mask, goodvoxels_ribbon_mask, [("out_file", "in_file")]), (ribbon_boldsrc_xfm, goodvoxels_ribbon_mask, [("output_image", "mask_file")]), (goodvoxels_ribbon_mask, apply_goodvoxels_ribbon_mask, [("out_file", "mask_file")]), - (rename_src, apply_goodvoxels_ribbon_mask, [("out_file", "in_file")]), + (inputnode, apply_goodvoxels_ribbon_mask, [("bold_file", "in_file")]), + (apply_goodvoxels_ribbon_mask, outputnode, [("out_file", "masked_bold")]), + (goodvoxels_ribbon_mask, outputnode, [("out_file", "goodvoxels_ribbon")]), ] ) - return workflow, apply_goodvoxels_ribbon_mask - - -def apply_med_surf_nans_if( - medial_surface_nan, workflow, inputnode, sampler, update_metadata, medial_nans -): - """ - Parameters - ---------- - medial_surface_nan : :obj:`bool` - Replace medial wall values with NaNs on functional GIFTI files - workflow : _type_, _description_ - inputnode : _type_, _description_ - sampler : _type_, _description_ - update_metadata : _type_, _description_ - medial_nans : _type_, _description_ - - :return: workflow, with medial surface NaNs applied if medial_surface_nan - """ - if medial_surface_nan: - # fmt: off - workflow.connect([ - (inputnode, medial_nans, [("subjects_dir", "subjects_dir")]), - (sampler, medial_nans, [("out_file", "in_file")]), - (medial_nans, update_metadata, [("out_file", "in_file")]), - ]) - # fmt: on - else: - workflow.connect(sampler, "out_file", update_metadata, "in_file") - return workflow + return workflow # , apply_goodvoxels_ribbon_mask def connect_bold_surf_wf( @@ -1348,9 +1207,3 @@ def _itk2lta(in_file, src_file, dst_file): in_file, fmt="fs" if in_file.endswith(".lta") else "itk", reference=src_file ).to_filename(out_file, moving=dst_file, fmt="fs") return str(out_file) - - -def _sorted_by_basename(inlist): - from os.path import basename - - return sorted(inlist, key=lambda x: str(basename(x))) From 7077ffb7d258b354d3ca4b38ae6e822360c9705f Mon Sep 17 00:00:00 2001 From: Greg Conan Date: Tue, 13 Dec 2022 15:07:42 -0600 Subject: [PATCH 11/11] FIX: Minor anat wf & workbench.py fixes --- nibabies/interfaces/workbench.py | 8 ++++---- nibabies/workflows/anatomical/base.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/nibabies/interfaces/workbench.py b/nibabies/interfaces/workbench.py index 238cfc4c..b717f1db 100644 --- a/nibabies/interfaces/workbench.py +++ b/nibabies/interfaces/workbench.py @@ -1476,26 +1476,26 @@ class CreateSignedDistanceVolumeInputSpec(CommandLineInputSpec): desc="Name of file to store a mask where the ``out_file`` has a computed value", ) fill_value = traits.Float( + 0.0, mandatory=False, - default=0.0, usedefault=True, argstr="-fill-value %f", desc="value to put in all voxels that don't get assigned a distance", ) exact_limit = traits.Float( - default=5.0, + 5.0, usedefault=True, argstr="-exact-limit %f", desc="distance for exact output in mm", ) approx_limit = traits.Float( - default=20.0, + 20.0, usedefault=True, argstr="-approx-limit %f", desc="distance for approximate output in mm", ) approx_neighborhood = traits.Int( - default=2, + 2, usedefault=True, argstr="-approx-neighborhood %d", desc="size of neighborhood cube measured from center to face in voxels, default 2 = 5x5x5", diff --git a/nibabies/workflows/anatomical/base.py b/nibabies/workflows/anatomical/base.py index aa7c2fba..c02cf6ac 100644 --- a/nibabies/workflows/anatomical/base.py +++ b/nibabies/workflows/anatomical/base.py @@ -415,7 +415,7 @@ def init_infant_anat_wf( ]), (t1w_preproc_wf, anat_ribbon_wf, [ ("outputnode.t1w_mask", "inputnode.t1w_mask"), - ]) + ]), (surface_recon_wf, anat_ribbon_wf, [ ("outputnode.surfaces", "inputnode.surfaces"), ]),