From 02e51c5bea68f1e114c0c9d5b55b6efb7bad16d0 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 18 May 2023 08:20:00 -0400 Subject: [PATCH 01/21] MNT: Pin master branches of nipreps dependencies --- pyproject.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 004a6929c..46df75c75 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,15 +25,15 @@ dependencies = [ "nipype >= 1.8.5", "nitime", "nitransforms >= 21.0.0", - "niworkflows ~= 1.7.9", + "niworkflows @ git+https://github.com/nipreps/niworkflows.git@master", "numpy", "packaging", "pandas", "psutil >= 5.4", "pybids >= 0.15.2", "requests", - "sdcflows ~= 2.4.3", - "smriprep >= 0.11.1", + "sdcflows @ git+https://github.com/nipreps/sdcflows.git@master", + "smriprep @ git+https://github.com/nipreps/smriprep.git@master", "tedana ~= 0.0.9", "templateflow >= 23.0.0", "toml", From 70b93b9ebc061f1277f6b96bd65a1ebd925f2ae4 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 11 May 2023 22:42:18 -0400 Subject: [PATCH 02/21] ENH: Add several workbench commands Resample, mask, fill-holes, remove-islands --- fmriprep/interfaces/metric.py | 271 --------- fmriprep/interfaces/workbench.py | 779 ++++++++++++++++++++++++++ fmriprep/workflows/bold/resampling.py | 2 +- 3 files changed, 780 insertions(+), 272 deletions(-) delete mode 100644 fmriprep/interfaces/metric.py create mode 100644 fmriprep/interfaces/workbench.py diff --git a/fmriprep/interfaces/metric.py b/fmriprep/interfaces/metric.py deleted file mode 100644 index 4544e089a..000000000 --- a/fmriprep/interfaces/metric.py +++ /dev/null @@ -1,271 +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 import logging -from nipype.interfaces.base import CommandLineInputSpec, File, TraitedSpec, traits -from nipype.interfaces.workbench.base import WBCommand - -iflogger = logging.getLogger("nipype.interface") - - -class MetricDilateInputSpec(CommandLineInputSpec): - in_file = File( - exists=True, - mandatory=True, - argstr="%s ", - position=0, - desc="the metric to dilate", - ) - - surf_file = 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 = File( - argstr="-bad-vertex-roi %s ", - position=4, - desc="metric file, positive values denote vertices to have their values replaced", - ) - - data_roi_file = 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 calculating 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. - """ - - 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``. - """ - - 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/fmriprep/interfaces/workbench.py b/fmriprep/interfaces/workbench.py new file mode 100644 index 000000000..7842551f2 --- /dev/null +++ b/fmriprep/interfaces/workbench.py @@ -0,0 +1,779 @@ +# -*- 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 import logging +from nipype.interfaces.base import ( + CommandLineInputSpec, + File, + TraitedSpec, + isdefined, + traits, +) +from nipype.interfaces.workbench.base import WBCommand + +iflogger = logging.getLogger("nipype.interface") + + +class MetricDilateInputSpec(CommandLineInputSpec): + in_file = File( + exists=True, + mandatory=True, + argstr="%s ", + position=0, + desc="the metric to dilate", + ) + + surf_file = 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 = File( + argstr="-bad-vertex-roi %s ", + position=4, + desc="metric file, positive values denote vertices to have their values replaced", + ) + + data_roi_file = 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 calculating the dilated value with weighted method", + ) + + +class MetricDilateOutputSpec(TraitedSpec): + out_file = File(exists=True, desc="output file") + + +class MetricDilate(WBCommand): + """Dilate a metric file on a surface. + + 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. + """ + + 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``. + """ + + 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 + + +class VolumeToSurfaceMappingInputSpec(CommandLineInputSpec): + volume_file = File( + exists=True, + argstr="%s", + mandatory=True, + position=1, + desc="the volume to map data from", + ) + surface_file = File( + exists=True, + argstr="%s", + mandatory=True, + position=2, + desc="the surface to map the data onto", + ) + out_file = File( + name_source=["surface_file"], + name_template="%s_mapped.func.gii", + keep_extension=False, + argstr="%s", + position=3, + desc="the output metric file", + ) + method = traits.Enum( + "trilinear", + "enclosing", + "cubic", + "ribbon-constrained", + "myelin-style", + argstr="-%s", + position=4, + desc="the interpolation method to use", + ) + _ribbon_constrained = [ + "inner_surface", + "outer_surface", + "volume_roi", + "weighted", + "voxel_subdiv", + "gaussian", + "interpolate", + "bad_vertices_out", + "output_weights", + "output_weights_text", + ] + _myelin_style = [ + "ribbon_roi", + "thickness", + "sigma", + "legacy_bug", + ] + inner_surface = File( + exists=True, + argstr="%s", + position=5, + desc="the inner surface of the ribbon [-ribbon-constrained]", + xor=_myelin_style, + ) + outer_surface = File( + exists=True, + argstr="%s", + position=6, + desc="the outer surface of the ribbon [-ribbon-constrained]", + xor=_myelin_style, + ) + volume_roi = File( + exists=True, + argstr="-volume-roi %s", + position=7, + desc="use a volume roi [-ribbon-constrained]", + xor=_myelin_style, + ) + weighted = traits.Bool( + argstr="-weighted", + position=8, + desc="treat the roi values as weightings rather than binary [-ribbon-constrained]", + requires=["volume_roi"], + xor=_myelin_style, + ) + voxel_subdiv = traits.Int( + default_value=3, + argstr="-voxel-subdiv %d", + desc="voxel divisions while estimating voxel weights [-ribbon-constrained]", + xor=_myelin_style, + ) + thin_columns = traits.Bool( + argstr="-thin-columns", + desc="use non-overlapping polyhedra [-ribbon-constrained]", + xor=_myelin_style, + ) + gaussian = traits.Float( + argstr="-gaussian %g", + desc="reduce weight to voxels that aren't near [-ribbon-constrained]", + xor=_myelin_style, + ) + interpolate = traits.Enum( + "CUBIC", + "TRILINEAR", + "ENCLOSING_VOXEL", + argstr="-interpolate %s", + desc="instead of a weighted average of voxels, " + "interpolate at subpoints inside the ribbon [-ribbon-constrained]", + xor=_myelin_style, + ) + bad_vertices_out = File( + argstr="-bad-vertices-out %s", + desc="output an ROI of which vertices didn't intersect any valid voxels", + xor=_myelin_style, + ) + output_weights = traits.Int( + argstr="-output-weights %(0)d output_weights.nii.gz", + desc="write the voxel weights for a vertex to a volume file", + xor=_myelin_style, + ) + output_weights_text = traits.File( + argstr="-output-weights-text %s", + desc="write the voxel weights for all vertices to a text file", + xor=_myelin_style, + ) + ribbon_roi = File( + exists=True, + argstr="%s", + position=5, + desc="an roi volume of the cortical ribbon for this hemisphere [-myelin-style]", + xor=_ribbon_constrained, + ) + thickness = File( + exists=True, + argstr="%s", + position=6, + desc="the thickness metric file for this hemisphere [-myelin-style]", + xor=_ribbon_constrained, + ) + sigma = traits.Float( + argstr="%g", + position=7, + desc="gaussian kernel in mm for weighting voxels within range [-myelin-style]", + xor=_ribbon_constrained, + ) + legacy_bug = traits.Bool( + argstr="-legacy-bug", + position=8, + desc="use the old bug in the myelin-style algorithm [-myelin-style]", + xor=_ribbon_constrained, + ) + subvol_select = traits.Int( + argstr="-subvol-select %d", + desc="select a single subvolume to map", + ) + + """\ +MAP VOLUME TO SURFACE + wb_command -volume-to-surface-mapping + - the volume to map data from + - the surface to map the data onto + - output - the output metric file + + [-trilinear] - use trilinear volume interpolation + + [-enclosing] - use value of the enclosing voxel + + [-cubic] - use cubic splines + + [-ribbon-constrained] - use ribbon constrained mapping algorithm + - the inner surface of the ribbon + - the outer surface of the ribbon + + [-volume-roi] - use a volume roi + - the roi volume file + + [-weighted] - treat the roi values as weightings rather than binary + + [-voxel-subdiv] - voxel divisions while estimating voxel weights + - number of subdivisions, default 3 + + [-thin-columns] - use non-overlapping polyhedra + + [-gaussian] - reduce weight to voxels that aren't near + - value to multiply the local thickness by, to get the + gaussian sigma + + [-interpolate] - instead of a weighted average of voxels, interpolate + at subpoints inside the ribbon + - interpolation method, must be CUBIC, ENCLOSING_VOXEL, or + TRILINEAR + + [-bad-vertices-out] - output an ROI of which vertices didn't intersect + any valid voxels + - output - the output metric file of vertices that have + no data + + [-output-weights] - write the voxel weights for a vertex to a volume + file + - the vertex number to get the voxel weights for, 0-based + - output - volume to write the weights to + + [-output-weights-text] - write the voxel weights for all vertices to a + text file + - output - the output text filename + + [-myelin-style] - use the method from myelin mapping + - an roi volume of the cortical ribbon for this + hemisphere + - a metric file of cortical thickness + - gaussian kernel in mm for weighting voxels within range + + [-legacy-bug] - emulate old v1.2.3 and earlier code that didn't follow + a cylinder cutoff + + [-subvol-select] - select a single subvolume to map + - the subvolume number or name +""" + + +class VolumeToSurfaceMappingOutputSpec(TraitedSpec): + out_file = File(desc="the output metric file") + bad_vertices_file = File(desc="the output metric file of vertices that have no data") + weights_file = File(desc="volume to write the weights to") + weights_text_file = File(desc="the output text filename") + + +class VolumeToSurfaceMapping(WBCommand): + """Map a volume to a surface using one of several methods. + + From https://humanconnectome.org/software/workbench-command/-volume-to-surface-mapping:: + + 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. + + Examples: + >>> from fmriprep.interfaces.workbench import VolumeToSurfaceMapping + >>> vol2surf = VolumeToSurfaceMapping() + >>> vol2surf.inputs.volume_file = 'bold.nii.gz' + >>> vol2surf.inputs.surface_file = 'lh.midthickness.surf.gii' + >>> vol2surf.inputs.method = 'ribbon-constrained' + >>> vol2surf.inputs.inner_surface = 'lh.white.surf.gii' + >>> vol2surf.inputs.outer_surface = 'lh.pial.surf.gii' + >>> vol2surf.cmdline + 'wb_command -volume-to-surface-mapping bold.nii.gz lh.midthickness.surf.gii \ + lh.midthickness.surf_mapped.func.gii -ribbon-constrained lh.white.surf.gii lh.pial.surf.gii' + """ + + input_spec = VolumeToSurfaceMappingInputSpec + output_spec = VolumeToSurfaceMappingOutputSpec + _cmd = "wb_command -volume-to-surface-mapping" + + def _format_arg(self, opt, spec, val): + if opt in self.input_spec._ribbon_constrained: + if self.inputs.method != "ribbon-constrained": + return "" + elif opt in self.input_spec._myelin_style: + if self.inputs.method != "myelin-style": + return "" + return super()._format_arg(opt, spec, val) + + def _list_outputs(self): + outputs = super()._list_outputs() + if isdefined(self.inputs.bad_vertices_out): + outputs["bad_vertices_file"] = os.path.abspath(self.inputs.bad_vertices_out) + if isdefined(self.inputs.output_weights): + outputs["weights_file"] = os.path.abspath(self.inputs.output_weights) + if isdefined(self.inputs.output_weights_text): + outputs["weights_text_file"] = os.path.abspath(self.inputs.output_weights_text) + return outputs + + +class MetricMaskInputSpec(CommandLineInputSpec): + """MASK A METRIC FILE + wb_command -metric-mask + - the input metric + - the mask metric + - output - the output metric + + [-column] - select a single column + - the column number or name + + By default, the output metric is a copy of the input metric, but with + zeros wherever the mask metric is zero or negative. if -column is + specified, the output contains only one column, the masked version of the + specified input column.""" + + in_file = File( + exists=True, + argstr="%s", + position=1, + mandatory=True, + desc="input metric file", + ) + mask = File( + exists=True, + argstr="%s", + position=2, + mandatory=True, + desc="mask metric file", + ) + out_file = File( + name_template="%s_masked.func.gii", + name_source=["in_file"], + keep_extension=False, + argstr="%s", + position=3, + desc="output metric file", + ) + column = traits.Either( + traits.Int, + traits.String, + argstr="-column %s", + desc="select a single column by number or name", + ) + + +class MetricMaskOutputSpec(TraitedSpec): + out_file = File(desc="output metric file") + + +class MetricMask(WBCommand): + """Mask a metric file. + + Examples + + >>> from fmriprep.interfaces.workbench import MetricMask + >>> metric_mask = MetricMask() + >>> metric_mask.inputs.in_file = 'lh.bold.func.gii' + >>> metric_mask.inputs.mask = 'lh.roi.shape.gii' + >>> metric_mask.cmdline + 'wb_command -metric-mask lh.bold.func.gii lh.roi.shape.gii lh.bold.func_masked.func.gii' + """ + + input_spec = MetricMaskInputSpec + output_spec = MetricMaskOutputSpec + _cmd = "wb_command -metric-mask" + + +class MetricFillHolesInputSpec(TraitedSpec): + """FILL HOLES IN AN ROI METRIC + + wb_command -metric-fill-holes + - the surface to use for neighbor information + - the input ROI metric + - output - the output ROI metric + + [-corrected-areas] - vertex areas to use instead of computing them from + the surface + - the corrected vertex areas, as a metric + + Finds all connected areas that are not included in the ROI, and writes + ones into all but the largest one, in terms of surface area.""" + + surface_file = File( + mandatory=True, + exists=True, + argstr="%s", + position=1, + desc="surface to use for neighbor information", + ) + metric_file = File( + mandatory=True, + exists=True, + argstr="%s", + position=2, + desc="input ROI metric", + ) + out_file = File( + name_template="%s_filled.shape.gii", + name_source="metric_file", + keep_extension=False, + argstr="%s", + position=3, + desc="output ROI metric", + ) + corrected_areas = File( + exists=True, + argstr="-corrected-areas %s", + desc="vertex areas to use instead of computing them from the surface", + ) + + +class MetricFillHolesOutputSpec(TraitedSpec): + out_file = File(desc="output ROI metric") + + +class MetricFillHoles(WBCommand): + """Fill holes in an ROI metric. + + Examples + + >>> from fmriprep.interfaces.workbench import MetricFillHoles + >>> fill_holes = MetricFillHoles() + >>> fill_holes.inputs.surface_file = 'lh.midthickness.surf.gii' + >>> fill_holes.inputs.metric_file = 'lh.roi.shape.gii' + >>> fill_holes.cmdline + 'wb_command -metric-fill-holes lh.midthickness.surf.gii lh.roi.shape.gii \ + lh.roi.shape_filled.shape.gii' + """ + + input_spec = MetricFillHolesInputSpec + output_spec = MetricFillHolesOutputSpec + _cmd = "wb_command -metric-fill-holes" + + +class MetricRemoveIslandsInputSpec(TraitedSpec): + """REMOVE ISLANDS IN AN ROI METRIC + + wb_command -metric-remove-islands + - the surface to use for neighbor information + - the input ROI metric + - output - the output ROI metric + + [-corrected-areas] - vertex areas to use instead of computing them from + the surface + - the corrected vertex areas, as a metric + + Finds all connected areas in the ROI, and zeros out all but the largest + one, in terms of surface area.""" + + surface_file = File( + mandatory=True, + exists=True, + argstr="%s", + position=1, + desc="surface to use for neighbor information", + ) + metric_file = File( + mandatory=True, + exists=True, + argstr="%s", + position=2, + desc="input ROI metric", + ) + out_file = File( + name_template="%s_noislands.shape.gii", + name_source="metric_file", + keep_extension=False, + argstr="%s", + position=3, + desc="output ROI metric", + ) + corrected_areas = File( + exists=True, + argstr="-corrected-areas %s", + desc="vertex areas to use instead of computing them from the surface", + ) + + +class MetricRemoveIslandsOutputSpec(TraitedSpec): + out_file = File(desc="output ROI metric") + + +class MetricRemoveIslands(WBCommand): + """Remove islands in an ROI metric. + + Examples + + >>> from fmriprep.interfaces.workbench import MetricRemoveIslands + >>> remove_islands = MetricRemoveIslands() + >>> remove_islands.inputs.surface_file = 'lh.midthickness.surf.gii' + >>> remove_islands.inputs.metric_file = 'lh.roi.shape.gii' + >>> remove_islands.cmdline + 'wb_command -metric-remove-islands lh.midthickness.surf.gii \ + lh.roi.shape.gii lh.roi.shape_noislands.shape.gii' + """ + + input_spec = MetricRemoveIslandsInputSpec + output_spec = MetricRemoveIslandsOutputSpec + _cmd = "wb_command -metric-remove-islands" diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index e0bd88d9b..c1c075761 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -43,7 +43,7 @@ from niworkflows.interfaces.freesurfer import MedialNaNs from ...config import DEFAULT_MEMORY_MIN_GB -from ...interfaces.metric import MetricDilate +from ...interfaces.workbench import MetricDilate if ty.TYPE_CHECKING: from niworkflows.utils.spaces import SpatialReferences From 972648a941d8809e132fabeda140c50a4e724949 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Fri, 12 May 2023 12:31:55 -0400 Subject: [PATCH 03/21] ENH: Add GIFTI interface to generate ROI from thickness --- fmriprep/interfaces/gifti.py | 65 ++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 fmriprep/interfaces/gifti.py diff --git a/fmriprep/interfaces/gifti.py b/fmriprep/interfaces/gifti.py new file mode 100644 index 000000000..2e2ebe9e8 --- /dev/null +++ b/fmriprep/interfaces/gifti.py @@ -0,0 +1,65 @@ +# -*- 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: +"""Interfaces for manipulating GIFTI files.""" +import os + +import nibabel as nb +import numpy as np +from nipype.interfaces.base import File, SimpleInterface, TraitedSpec, isdefined, traits +from nipype.utils.filemanip import fname_presuffix + + +class CreateROIInputSpec(TraitedSpec): + subject_id = traits.Str(desc='subject ID') + hemisphere = traits.Enum( + "L", + "R", + mandatory=True, + desc='hemisphere', + ) + thickness_file = File(exists=True, mandatory=True, desc='input GIFTI file') + + +class CreateROIOutputSpec(TraitedSpec): + roi_file = File(desc='output GIFTI file') + + +class CreateROI(SimpleInterface): + """Prepare GIFTI shape file for use in""" + + input_spec = CreateROIInputSpec + output_spec = CreateROIOutputSpec + + def _run_interface(self, runtime): + subject, hemi = self.inputs.subject_id, self.inputs.hemisphere + if not isdefined(subject): + subject = 'sub-XYZ' + img = nb.GiftiImage.from_filename(self.inputs.thickness_file) + # wb_command -set-structure + img.meta["AnatomicalStructurePrimary"] = {'L': 'CortexLeft', 'R': 'CortexRight'}[hemi] + darray = img.darrays[0] + # wb_command -set-map-names + meta = darray.meta + meta['Name'] = f"{subject}_{hemi}_ROI" + # wb_command -metric-palette calls have no effect on ROI files + + # Compiling an odd sequence of math operations that works out to: + # wb_command -metric-math "abs(var * -1) > 0" + roi = np.abs(darray.data) > 0 + + darray = nb.gifti.GiftiDataArray( + roi, + intent=darray.intent, + datatype=darray.datatype, + encoding=darray.encoding, + endian=darray.endian, + coordsys=darray.coordsys, + ordering=darray.ind_ord, + meta=meta, + ) + + out_filename = os.path.join(runtime.cwd, f"{subject}.{hemi}.roi.native.shape.gii") + img.to_filename(out_filename) + self._results["roi_file"] = out_filename + return runtime From 90621c2802ef0e97a251b9166223705db6fd8b99 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 17 May 2023 13:15:47 -0400 Subject: [PATCH 04/21] ENH: Add fsLR resampling workflow --- fmriprep/workflows/base.py | 6 +- fmriprep/workflows/bold/base.py | 23 ++- fmriprep/workflows/bold/resampling.py | 214 +++++++++++++++++++++++++- 3 files changed, 238 insertions(+), 5 deletions(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index 446b7ca4d..92ac1b935 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -460,7 +460,11 @@ def init_single_subject_wf(subject_id: str): ('outputnode.subject_id', 'inputnode.subject_id'), ('outputnode.anat_ribbon', 'inputnode.anat_ribbon'), ('outputnode.t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'), - ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm')]), + ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ('outputnode.surfaces', 'inputnode.surfaces'), + ('outputnode.morphometrics', 'inputnode.morphometrics'), + ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'), + ]), ]) # fmt:on func_preproc_wfs.append(func_preproc_wf) diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 7448cc910..8c2ceac52 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -353,6 +353,9 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): "anat_ribbon", "t1w2fsnative_xfm", "fsnative2t1w_xfm", + "surfaces", + "morphometrics", + "sphere_reg_fsLR", "fmap", "fmap_ref", "fmap_coeff", @@ -894,7 +897,12 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): # CIFTI output if config.workflow.cifti_output: - from .resampling import init_bold_grayords_wf + from .resampling import init_bold_fsLR_resampling_wf, init_bold_grayords_wf + + bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf( + grayord_density=config.workflow.cifti_output, + mem_gb=mem_gb["resampled"], + ) bold_grayords_wf = init_bold_grayords_wf( grayord_density=config.workflow.cifti_output, @@ -904,8 +912,17 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): # fmt:off workflow.connect([ - (bold_std_trans_wf, bold_grayords_wf, [ - ("outputnode.bold_std", "inputnode.bold_std"), + (inputnode, bold_fsLR_resampling_wf, [ + ("surfaces", "inputnode.surfaces"), + ("morphometrics", "inputnode.morphometrics"), + ("sphere_reg_fsLR", "inputnode.sphere_reg_fsLR"), + ("anat_ribbon", "inputnode.anat_ribbon"), + ]), + (bold_t1_trans_wf, bold_fsLR_resampling_wf, [ + ("outputnode.bold_t1", "inputnode.bold_file"), + ]), + (bold_fsLR_resampling_wf, bold_grayords_wf, [ + ("outputnode.bold_fsLR", "inputnode.bold_std"), ("outputnode.spatial_reference", "inputnode.spatial_reference"), ]), (bold_surf_wf, bold_grayords_wf, [ diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index c1c075761..d125a0882 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -43,7 +43,7 @@ from niworkflows.interfaces.freesurfer import MedialNaNs from ...config import DEFAULT_MEMORY_MIN_GB -from ...interfaces.workbench import MetricDilate +from ...interfaces.workbench import MetricDilate, MetricMask, MetricResample if ty.TYPE_CHECKING: from niworkflows.utils.spaces import SpatialReferences @@ -581,6 +581,189 @@ def _calc_lower_thr(in_stats): return workflow +def init_bold_fsLR_resampling_wf( + grayord_density: ty.Literal['91k', '170k'], + estimate_goodvoxels: bool, + mem_gb: float, + name: str = "bold_fsLR_resampling_wf", +): + """Resample BOLD time series to fsLR surface. + + This workflow is derived heavily from three scripts within the DCAN-HCP pipelines scripts + + Line numbers correspond to the locations of the code in the original scripts, found at: + https://github.com/DCAN-Labs/DCAN-HCP/tree/9291324/ + """ + import templateflow.api as tf + from niworkflows.engine.workflows import LiterateWorkflow as Workflow + from smriprep import data as smriprep_data + from smriprep.interfaces.workbench import SurfaceResample + + from fmriprep.interfaces.gifti import CreateROI + from fmriprep.interfaces.workbench import ( + MetricFillHoles, + MetricRemoveIslands, + VolumeToSurfaceMapping, + ) + + fslr_density = "32k" if grayord_density == "91k" else "59k" + + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + 'bold_file', + 'surfaces', + 'morphometrics', + 'sphere_reg_fsLR', + 'anat_ribbon', + ] + ), + name='inputnode', + ) + + itersource = pe.Node( + niu.IdentityInterface(fields=['hemi']), + name='itersource', + iterables=[('hemi', ['L', 'R'])], + ) + + outputnode = pe.JoinNode( + niu.IdentityInterface(fields=['bold_fsLR']), + name='outputnode', + joinsource='itersource', + ) + + # select white, midthickness and pial surfaces based on hemi + select_surfaces = pe.Node( + niu.Function( + function=_select_surfaces, + output_names=[ + 'white', + 'pial', + 'midthickness', + 'thickness', + 'sphere_reg', + 'template_sphere', + 'template_roi', + ], + ), + name='select_surfaces', + ) + select_surfaces.inputs.template_spheres = [ + str(sphere) + for sphere in tf.get( + template='fsLR', + density=fslr_density, + suffix='sphere', + space=None, + extension='.surf.gii', + ) + ] + atlases = smriprep_data.load_resource('atlases') + select_surfaces.inputs.template_rois = [ + str(atlases / 'L.atlasroi.32k_fs_LR.shape.gii'), + str(atlases / 'R.atlasroi.32k_fs_LR.shape.gii'), + ] + + # Reimplements lines 282-290 of FreeSurfer2CaretConvertAndRegisterNonlinear.sh + initial_roi = pe.Node(CreateROI(), name="initial_roi", mem_gb=DEFAULT_MEMORY_MIN_GB) + + # Lines 291-292 + fill_holes = pe.Node(MetricFillHoles(), name="fill_holes", mem_gb=DEFAULT_MEMORY_MIN_GB) + native_roi = pe.Node(MetricRemoveIslands(), name="native_roi", mem_gb=DEFAULT_MEMORY_MIN_GB) + + # Line 393 of FreeSurfer2CaretConvertAndRegisterNonlinear.sh + downsampled_midthickness = pe.Node( + SurfaceResample(method='BARYCENTRIC'), + name="downsampled_midthickness", + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + # RibbonVolumeToSurfaceMapping.sh + # Line 85 thru ... + volume_to_surface = pe.Node( + VolumeToSurfaceMapping(method="ribbon-constrained"), + name="volume_to_surface", + mem_gb=mem_gb * 3, + ) + metric_dilate = pe.Node(MetricDilate(distance=10, nearest=True), name="metric_dilate") + mask_native = pe.Node(MetricMask(), name="mask_native") + resample_to_fsLR = pe.Node( + MetricResample(method='ADAP_BARY_AREA', area_surfs=True), name="resample_to_fsLR" + ) + # ... line 89 + mask_fsLR = pe.Node(MetricMask(), name="mask_fsLR") + + # fmt: off + workflow.connect([ + (inputnode, select_surfaces, [ + ('surfaces', 'surfaces'), + ('morphometrics', 'morphometrics'), + ('sphere_reg_fsLR', 'spherical_registrations'), + ]), + (itersource, select_surfaces, [('hemi', 'hemi')]), + # Native ROI file from thickness + (itersource, initial_roi, [('hemi', 'hemisphere')]), + (select_surfaces, initial_roi, [('thickness', 'thickness_file')]), + (select_surfaces, fill_holes, [('midthickness', 'surface_file')]), + (select_surfaces, native_roi, [('midthickness', 'surface_file')]), + (initial_roi, fill_holes, [('roi_file', 'metric_file')]), + (fill_holes, native_roi, [('out_file', 'metric_file')]), + # Downsample midthickness to fsLR density + (select_surfaces, downsampled_midthickness, [ + ('midthickness', 'surface_in'), + ('sphere_reg', 'current_sphere'), + ('template_sphere', 'new_sphere'), + ]), + # Resample BOLD to native surface, dilate and mask + (inputnode, volume_to_surface, [ + ('bold_file', 'volume_file'), + ]), + (select_surfaces, volume_to_surface, [ + ('midthickness', 'surface_file'), + ('white', 'inner_surface'), + ('pial', 'outer_surface'), + ]), + (select_surfaces, metric_dilate, [('midthickness', 'surf_file')]), + (volume_to_surface, metric_dilate, [('out_file', 'in_file')]), + (native_roi, mask_native, [('out_file', 'mask')]), + (metric_dilate, mask_native, [('out_file', 'in_file')]), + # Resample BOLD to fsLR and mask + (select_surfaces, resample_to_fsLR, [ + ('sphere_reg', 'current_sphere'), + ('template_sphere', 'new_sphere'), + ('midthickness', 'current_area'), + ]), + (downsampled_midthickness, resample_to_fsLR, [('surface_out', 'new_area')]), + (native_roi, resample_to_fsLR, [('out_file', 'roi_metric')]), + (mask_native, resample_to_fsLR, [('out_file', 'in_file')]), + (select_surfaces, mask_fsLR, [('template_roi', 'mask')]), + (resample_to_fsLR, mask_fsLR, [('out_file', 'in_file')]), + # Output + (mask_fsLR, outputnode, [('out_file', 'bold_fsLR')]), + ]) + # fmt: on + + if estimate_goodvoxels: + goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb) + + # fmt: off + workflow.connect([ + (inputnode, goodvoxels_bold_mask_wf, [ + ("bold_file", "inputnode.bold_file"), + ("anat_ribbon", "inputnode.anat_ribbon"), + ]), + (goodvoxels_bold_mask_wf, volume_to_surface, [ + ("outputnode.goodvoxels_ribbon", "volume_roi"), + ]), + ]) + # fmt: on + + return workflow + + def init_bold_std_trans_wf( freesurfer: bool, mem_gb: float, @@ -1285,3 +1468,32 @@ 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 _select_surfaces( + hemi, + surfaces, + morphometrics, + spherical_registrations, + template_spheres, + template_rois, +): + import os + import re + + idx = 0 if hemi == "L" else 1 + container = { + 'white': [], + 'pial': [], + 'midthickness': [], + 'thickness': [], + 'sphere': sorted(spherical_registrations), + 'template_sphere': sorted(template_spheres), + 'template_roi': sorted(template_rois), + } + find_name = re.compile(r'(?:^|[^d])(?Pwhite|pial|midthickness|thickness)') + for surface in surfaces + morphometrics: + match = find_name.search(os.path.basename(surface)) + if match: + container[match.group('name')].append(surface) + return tuple(sorted(surflist)[idx] for surflist in container.values()) From 4d03abf2890dfdc05fa4ce615d8b85e68e23ca83 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 17 May 2023 13:24:02 -0400 Subject: [PATCH 05/21] RF: Remove goodvoxels from mri_vol2surf workflow --- fmriprep/workflows/bold/base.py | 14 +---- fmriprep/workflows/bold/resampling.py | 90 +++------------------------ 2 files changed, 9 insertions(+), 95 deletions(-) diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 8c2ceac52..2b0e48cbd 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -213,13 +213,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): freesurfer_spaces = spaces.get_fs_spaces() project_goodvoxels = config.workflow.project_goodvoxels - if project_goodvoxels and freesurfer_spaces != ["fsaverage"]: - config.loggers.workflow.critical( - f"--project-goodvoxels only works with fsaverage (requested: {freesurfer_spaces})" - ) - config.loggers.workflow.warn("Disabling --project-goodvoxels") - project_goodvoxels = False - # Extract BIDS entities and metadata from BOLD file(s) entities = extract_entities(bold_file) layout = config.execution.layout @@ -440,7 +433,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): bids_root=layout.root, cifti_output=config.workflow.cifti_output, freesurfer=freesurfer, - project_goodvoxels=project_goodvoxels, + project_goodvoxels=False, all_metadata=all_metadata, multiecho=multiecho, output_dir=fmriprep_dir, @@ -875,7 +868,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): mem_gb=mem_gb["resampled"], surface_spaces=freesurfer_spaces, medial_surface_nan=config.workflow.medial_surface_nan, - project_goodvoxels=project_goodvoxels, name="bold_surf_wf", ) # fmt:off @@ -884,14 +876,10 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): ("subjects_dir", "inputnode.subjects_dir"), ("subject_id", "inputnode.subject_id"), ("t1w2fsnative_xfm", "inputnode.t1w2fsnative_xfm"), - ("anat_ribbon", "inputnode.anat_ribbon"), - ("t1w_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, [("outputnode.target", "inputnode.surf_refs")]), - (bold_surf_wf, func_derivatives_wf, [("outputnode.goodvoxels_ribbon", - "inputnode.goodvoxels_ribbon")]), ]) # fmt:on diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index d125a0882..cfa54a6c8 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -54,7 +54,6 @@ def init_bold_surf_wf( mem_gb: float, surface_spaces: ty.List[str], medial_surface_nan: bool, - project_goodvoxels: bool, name: str = "bold_surf_wf", ): """ @@ -63,11 +62,6 @@ def init_bold_surf_wf( 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. After sampling, - empty vertices are filled by dilating in data from the nearest "good" vertex, within a - radius of 10 mm (as measured on the target midthickness surface). - Outputs are in GIFTI format. Workflow Graph @@ -79,7 +73,6 @@ def init_bold_surf_wf( wf = init_bold_surf_wf(mem_gb=0.1, surface_spaces=["fsnative", "fsaverage5"], medial_surface_nan=False, - project_goodvoxels=False, ) Parameters @@ -91,9 +84,6 @@ def init_bold_surf_wf( 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. Inputs ------ @@ -105,10 +95,6 @@ def init_bold_surf_wf( FreeSurfer subject ID t1w2fsnative_xfm LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space - anat_ribbon - Cortical ribbon in T1w space - t1w_mask - Mask of the skull-stripped T1w image Outputs ------- @@ -130,14 +116,6 @@ def init_bold_surf_wf( 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 coefficient of -variation, or that lie outside the cortical surfaces, from the -surface projection. -""" - inputnode = pe.Node( niu.IdentityInterface( fields=[ @@ -145,8 +123,6 @@ def init_bold_surf_wf( "subject_id", "subjects_dir", "t1w2fsnative_xfm", - "anat_ribbon", - "t1w_mask", ] ), name="inputnode", @@ -204,7 +180,7 @@ def select_target(subject_id, space): ) outputnode = pe.Node( - niu.IdentityInterface(fields=["surfaces", "target", "goodvoxels_ribbon"]), + niu.IdentityInterface(fields=["surfaces", "target"]), name="outputnode", ) @@ -242,65 +218,25 @@ def select_target(subject_id, space): # # These depend on two optional steps: goodvoxel projection and medial wall nan replacement # - # inputnode -> optional(goodvoxels_bold_mask_wf) -> rename_src - # sampler -> optional(metric_dilate) -> optional(medial_nans) -> update_metadata + # inputnode -> rename_src + # sampler -> optional(medial_nans) -> update_metadata # - metric_dilate = pe.MapNode( - MetricDilate(distance=10, nearest=True), - iterfield=["in_file", "surf_file"], - name="metric_dilate", - mem_gb=mem_gb * 3, - ) - metric_dilate.inputs.surf_file = [ - str( - tf.get( - "fsaverage", - hemi=hemi, - density="164k", - suffix="midthickness", - extension=".surf.gii", - ) - ) - for hemi in "LR" - ] # Refine if medial vertices should be NaNs medial_nans = pe.MapNode( MedialNaNs(), iterfield=["in_file"], name="medial_nans", mem_gb=DEFAULT_MEMORY_MIN_GB ) - if project_goodvoxels: - goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb) - - # 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")]), - (goodvoxels_bold_mask_wf, outputnode, [ - ("outputnode.goodvoxels_ribbon", "goodvoxels_ribbon"), - ]), - (sampler, metric_dilate, [("out_file", "in_file")]), - ]) - # fmt: on - else: - workflow.connect([(inputnode, rename_src, [("source_file", "in_file")])]) + workflow.connect([(inputnode, rename_src, [("source_file", "in_file")])]) 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")]), ]) - - if medial_surface_nan and project_goodvoxels: - workflow.connect([(metric_dilate, medial_nans, [("out_file", "in_file")])]) - elif medial_surface_nan: - workflow.connect([(sampler, medial_nans, [("out_file", "in_file")])]) - elif project_goodvoxels: - workflow.connect([(metric_dilate, update_metadata, [("out_file", "in_file")])]) + # fmt: on else: workflow.connect([(sampler, update_metadata, [("out_file", "in_file")])]) @@ -352,7 +288,7 @@ def init_goodvoxels_bold_mask_wf(mem_gb: float, name: str = "goodvoxels_bold_mas outputnode = pe.Node( niu.IdentityInterface( fields=[ - "masked_bold", + "goodvoxels_mask", "goodvoxels_ribbon", ] ), @@ -558,22 +494,12 @@ def _calc_lower_thr(in_stats): 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")]), (ribbon_boldsrc_xfm, goodvoxels_ribbon_mask, [("output_image", "mask_file")]), - (goodvoxels_ribbon_mask, apply_goodvoxels_ribbon_mask, [("out_file", "mask_file")]), - (inputnode, apply_goodvoxels_ribbon_mask, [("bold_file", "in_file")]), - (apply_goodvoxels_ribbon_mask, outputnode, [("out_file", "masked_bold")]), + (goodvoxels_mask, outputnode, [("out_file", "goodvoxels_mask")]), (goodvoxels_ribbon_mask, outputnode, [("out_file", "goodvoxels_ribbon")]), ] ) From 9508e86ba5c93be15b98540c26ff423ead06833c Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 17 May 2023 13:25:38 -0400 Subject: [PATCH 06/21] RF: Remove implicit fsaverage resampling from --cifti-output --- fmriprep/config.py | 1 - fmriprep/workflows/bold/resampling.py | 2 -- 2 files changed, 3 deletions(-) diff --git a/fmriprep/config.py b/fmriprep/config.py index b3579f3d0..d384250c6 100644 --- a/fmriprep/config.py +++ b/fmriprep/config.py @@ -765,7 +765,6 @@ def init_spaces(checkpoint=True): if cifti_output: # CIFTI grayordinates to corresponding FSL-MNI resolutions. vol_res = "2" if cifti_output == "91k" else "1" - spaces.add(Reference("fsaverage", {"den": "164k"})) spaces.add(Reference("MNI152NLin6Asym", {"res": vol_res})) # Make the SpatialReferences object available diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index cfa54a6c8..f3b2e1099 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -33,7 +33,6 @@ import typing as ty -import nipype.interfaces.workbench as wb from nipype import Function from nipype.interfaces import freesurfer as fs from nipype.interfaces import fsl @@ -102,7 +101,6 @@ def init_bold_surf_wf( BOLD series, resampled to FreeSurfer surfaces """ - import templateflow.api as tf from nipype.interfaces.io import FreeSurferSource from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.surf import GiftiSetAnatomicalStructure From c2f4affa769e094008be541fb5507ec12aecd838 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 17 May 2023 15:04:21 -0400 Subject: [PATCH 07/21] RF: Remove cifti dependency on fsaverage workflow --- fmriprep/workflows/bold/base.py | 74 +++++++++++----------- fmriprep/workflows/bold/resampling.py | 89 +-------------------------- 2 files changed, 40 insertions(+), 123 deletions(-) diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 2b0e48cbd..01adb3b29 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -883,46 +883,46 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): ]) # fmt:on - # CIFTI output - if config.workflow.cifti_output: - from .resampling import init_bold_fsLR_resampling_wf, init_bold_grayords_wf + # CIFTI output + if config.workflow.cifti_output: + from .resampling import init_bold_fsLR_resampling_wf, init_bold_grayords_wf - bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf( - grayord_density=config.workflow.cifti_output, - mem_gb=mem_gb["resampled"], - ) + bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf( + estimate_goodvoxels=project_goodvoxels, + grayord_density=config.workflow.cifti_output, + mem_gb=mem_gb["resampled"], + ) - bold_grayords_wf = init_bold_grayords_wf( - grayord_density=config.workflow.cifti_output, - mem_gb=mem_gb["resampled"], - repetition_time=metadata["RepetitionTime"], - ) + bold_grayords_wf = init_bold_grayords_wf( + grayord_density=config.workflow.cifti_output, + mem_gb=mem_gb["resampled"], + repetition_time=metadata["RepetitionTime"], + ) - # fmt:off - workflow.connect([ - (inputnode, bold_fsLR_resampling_wf, [ - ("surfaces", "inputnode.surfaces"), - ("morphometrics", "inputnode.morphometrics"), - ("sphere_reg_fsLR", "inputnode.sphere_reg_fsLR"), - ("anat_ribbon", "inputnode.anat_ribbon"), - ]), - (bold_t1_trans_wf, bold_fsLR_resampling_wf, [ - ("outputnode.bold_t1", "inputnode.bold_file"), - ]), - (bold_fsLR_resampling_wf, bold_grayords_wf, [ - ("outputnode.bold_fsLR", "inputnode.bold_std"), - ("outputnode.spatial_reference", "inputnode.spatial_reference"), - ]), - (bold_surf_wf, bold_grayords_wf, [ - ("outputnode.surfaces", "inputnode.surf_files"), - ("outputnode.target", "inputnode.surf_refs"), - ]), - (bold_grayords_wf, outputnode, [ - ("outputnode.cifti_bold", "bold_cifti"), - ("outputnode.cifti_metadata", "cifti_metadata"), - ]), - ]) - # fmt:on + # fmt:off + workflow.connect([ + (inputnode, bold_fsLR_resampling_wf, [ + ("surfaces", "inputnode.surfaces"), + ("morphometrics", "inputnode.morphometrics"), + ("sphere_reg_fsLR", "inputnode.sphere_reg_fsLR"), + ("anat_ribbon", "inputnode.anat_ribbon"), + ]), + (bold_t1_trans_wf, bold_fsLR_resampling_wf, [ + ("outputnode.bold_t1", "inputnode.bold_file"), + ]), + (bold_std_trans_wf, bold_grayords_wf, [ + ("outputnode.bold_std", "inputnode.bold_std"), + ("outputnode.spatial_reference", "inputnode.spatial_reference"), + ]), + (bold_fsLR_resampling_wf, bold_grayords_wf, [ + ("outputnode.bold_fsLR", "inputnode.bold_fsLR"), + ]), + (bold_grayords_wf, outputnode, [ + ("outputnode.cifti_bold", "bold_cifti"), + ("outputnode.cifti_metadata", "cifti_metadata"), + ]), + ]) + # fmt:on if spaces.get_spaces(nonstandard=False, dim=(3,)): carpetplot_wf = init_carpetplot_wf( diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index f3b2e1099..90398c99f 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -1193,7 +1193,6 @@ def init_bold_grayords_wf( BIDS metadata file corresponding to ``cifti_bold``. """ - import templateflow.api as tf from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.cifti import GenerateCifti from niworkflows.interfaces.utility import KeySelect @@ -1207,15 +1206,14 @@ def init_bold_grayords_wf( density=grayord_density ) - fslr_density, mni_density = ("32k", "2") if grayord_density == "91k" else ("59k", "1") + mni_density = "2" if grayord_density == "91k" else "1" inputnode = pe.Node( niu.IdentityInterface( fields=[ "bold_std", + "bold_fsLR", "spatial_reference", - "surf_files", - "surf_refs", ] ), name="inputnode", @@ -1235,84 +1233,6 @@ def init_bold_grayords_wf( ) select_std.inputs.key = "MNI152NLin6Asym_res-%s" % mni_density - select_fs_surf = pe.Node( - KeySelect(fields=["surf_files"]), - name="select_fs_surf", - run_without_submitting=True, - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - select_fs_surf.inputs.key = "fsaverage" - - # Setup Workbench command. LR ordering for hemi can be assumed, as it is imposed - # by the iterfield of the MapNode in the surface sampling workflow above. - resample = pe.MapNode( - wb.MetricResample(method="ADAP_BARY_AREA", area_metrics=True), - name="resample", - iterfield=[ - "in_file", - "out_file", - "new_sphere", - "new_area", - "current_sphere", - "current_area", - ], - ) - resample.inputs.current_sphere = [ - str( - tf.get( - "fsaverage", - hemi=hemi, - density="164k", - desc="std", - suffix="sphere", - extension=".surf.gii", - ) - ) - for hemi in "LR" - ] - resample.inputs.current_area = [ - str( - tf.get( - "fsaverage", - hemi=hemi, - density="164k", - desc="vaavg", - suffix="midthickness", - extension=".shape.gii", - ) - ) - for hemi in "LR" - ] - resample.inputs.new_sphere = [ - str( - tf.get( - "fsLR", - space="fsaverage", - hemi=hemi, - density=fslr_density, - suffix="sphere", - extension=".surf.gii", - ) - ) - for hemi in "LR" - ] - resample.inputs.new_area = [ - str( - tf.get( - "fsLR", - hemi=hemi, - density=fslr_density, - desc="vaavg", - suffix="midthickness", - extension=".shape.gii", - ) - ) - for hemi in "LR" - ] - resample.inputs.out_file = [ - "space-fsLR_hemi-%s_den-%s_bold.gii" % (h, grayord_density) for h in "LR" - ] - gen_cifti = pe.Node( GenerateCifti( TR=repetition_time, @@ -1325,11 +1245,8 @@ def init_bold_grayords_wf( workflow.connect([ (inputnode, select_std, [("bold_std", "bold_std"), ("spatial_reference", "keys")]), - (inputnode, select_fs_surf, [("surf_files", "surf_files"), - ("surf_refs", "keys")]), - (select_fs_surf, resample, [("surf_files", "in_file")]), + (inputnode, gen_cifti, [("bold_fsLR", "surface_bolds")]), (select_std, gen_cifti, [("bold_std", "bold_file")]), - (resample, gen_cifti, [("out_file", "surface_bolds")]), (gen_cifti, outputnode, [("out_file", "cifti_bold"), ("out_metadata", "cifti_metadata")]), ]) From 21e5c3bd428a50580be3ac5584acb75839f128e9 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 18 May 2023 08:05:28 -0400 Subject: [PATCH 08/21] ENH: Add OpenMPCommandMixin to enable n_procs --- fmriprep/interfaces/workbench.py | 44 ++++++++++++++++++++++----- fmriprep/workflows/bold/base.py | 1 + fmriprep/workflows/bold/resampling.py | 12 ++++++-- 3 files changed, 48 insertions(+), 9 deletions(-) diff --git a/fmriprep/interfaces/workbench.py b/fmriprep/interfaces/workbench.py index 7842551f2..b43ae3ab0 100644 --- a/fmriprep/interfaces/workbench.py +++ b/fmriprep/interfaces/workbench.py @@ -6,6 +6,7 @@ from nipype import logging from nipype.interfaces.base import ( + CommandLine, CommandLineInputSpec, File, TraitedSpec, @@ -17,7 +18,36 @@ iflogger = logging.getLogger("nipype.interface") -class MetricDilateInputSpec(CommandLineInputSpec): +class OpenMPTraitedSpec(CommandLineInputSpec): + num_threads = traits.Int(desc="allows for specifying more threads") + + +class OpenMPCommandMixin(CommandLine): + input_spec = OpenMPTraitedSpec + + _num_threads = None + + def __init__(self, **inputs): + super().__init__(**inputs) + self.inputs.on_trait_change(self._num_threads_update, "num_threads") + if not self._num_threads: + self._num_threads = os.environ.get("OMP_NUM_THREADS", None) + if not isdefined(self.inputs.num_threads) and self._num_threads: + self.inputs.num_threads = int(self._num_threads) + self._num_threads_update() + + def _num_threads_update(self): + if self.inputs.num_threads: + self.inputs.environ.update({"OMP_NUM_THREADS": str(self.inputs.num_threads)}) + + def run(self, **inputs): + if "num_threads" in inputs: + self.inputs.num_threads = inputs["num_threads"] + self._num_threads_update() + return super().run(**inputs) + + +class MetricDilateInputSpec(OpenMPTraitedSpec): in_file = File( exists=True, mandatory=True, @@ -106,7 +136,7 @@ class MetricDilateOutputSpec(TraitedSpec): out_file = File(exists=True, desc="output file") -class MetricDilate(WBCommand): +class MetricDilate(WBCommand, OpenMPCommandMixin): """Dilate a metric file on a surface. For all data values designated as bad, if they neighbor a good value or @@ -132,7 +162,7 @@ class MetricDilate(WBCommand): _cmd = "wb_command -metric-dilate " -class MetricResampleInputSpec(CommandLineInputSpec): +class MetricResampleInputSpec(OpenMPTraitedSpec): in_file = File( exists=True, mandatory=True, @@ -223,7 +253,7 @@ class MetricResampleOutputSpec(TraitedSpec): roi_file = File(desc="ROI of vertices that got data from valid source vertices") -class MetricResample(WBCommand): +class MetricResample(WBCommand, OpenMPCommandMixin): """Resample a metric file to a different mesh. Resamples a metric file, given two spherical surfaces that are in @@ -270,14 +300,14 @@ def _format_arg(self, opt, spec, val): return super(MetricResample, self)._format_arg(opt, spec, val) def _list_outputs(self): - outputs = super(MetricResample, self)._list_outputs() + outputs = super()._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 -class VolumeToSurfaceMappingInputSpec(CommandLineInputSpec): +class VolumeToSurfaceMappingInputSpec(OpenMPTraitedSpec): volume_file = File( exists=True, argstr="%s", @@ -498,7 +528,7 @@ class VolumeToSurfaceMappingOutputSpec(TraitedSpec): weights_text_file = File(desc="the output text filename") -class VolumeToSurfaceMapping(WBCommand): +class VolumeToSurfaceMapping(WBCommand, OpenMPCommandMixin): """Map a volume to a surface using one of several methods. From https://humanconnectome.org/software/workbench-command/-volume-to-surface-mapping:: diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 01adb3b29..cc24095ff 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -890,6 +890,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf( estimate_goodvoxels=project_goodvoxels, grayord_density=config.workflow.cifti_output, + omp_nthreads=omp_nthreads, mem_gb=mem_gb["resampled"], ) diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index 90398c99f..7ba93eba1 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -508,6 +508,7 @@ def _calc_lower_thr(in_stats): def init_bold_fsLR_resampling_wf( grayord_density: ty.Literal['91k', '170k'], estimate_goodvoxels: bool, + omp_nthreads: int, mem_gb: float, name: str = "bold_fsLR_resampling_wf", ): @@ -611,11 +612,18 @@ def init_bold_fsLR_resampling_wf( VolumeToSurfaceMapping(method="ribbon-constrained"), name="volume_to_surface", mem_gb=mem_gb * 3, + n_procs=omp_nthreads, + ) + metric_dilate = pe.Node( + MetricDilate(distance=10, nearest=True), + name="metric_dilate", + n_procs=omp_nthreads, ) - metric_dilate = pe.Node(MetricDilate(distance=10, nearest=True), name="metric_dilate") mask_native = pe.Node(MetricMask(), name="mask_native") resample_to_fsLR = pe.Node( - MetricResample(method='ADAP_BARY_AREA', area_surfs=True), name="resample_to_fsLR" + MetricResample(method='ADAP_BARY_AREA', area_surfs=True), + name="resample_to_fsLR", + n_procs=omp_nthreads, ) # ... line 89 mask_fsLR = pe.Node(MetricMask(), name="mask_fsLR") From 7a42c27e0307250fbc6ead2680fe57e8e3a2c166 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 18 May 2023 10:06:18 -0400 Subject: [PATCH 09/21] DOCKER: Re-add git to main image to enable installing unstable dependencies --- Dockerfile | 1 + 1 file changed, 1 insertion(+) diff --git a/Dockerfile b/Dockerfile index 40f48ddc6..38f549eb9 100644 --- a/Dockerfile +++ b/Dockerfile @@ -128,6 +128,7 @@ RUN apt-get update && \ bc \ ca-certificates \ curl \ + git \ gnupg \ lsb-release \ netbase \ From e071a4ec16e0969a5745c354b3fdbe11baf4a899 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 18 May 2023 10:40:15 -0400 Subject: [PATCH 10/21] TEST: Copy test data folder and cd to temp path during interface doctests --- fmriprep/conftest.py | 9 ---- fmriprep/interfaces/confounds.py | 16 ------- fmriprep/interfaces/conftest.py | 40 +++++++++++++++++ fmriprep/interfaces/multiecho.py | 4 -- fmriprep/interfaces/tests/conftest.py | 8 ---- .../tests/data/bold.nii.gz} | 0 .../tests/data/lh.bold.func.gii} | 0 .../tests/data/lh.midthickness.surf.gii} | 0 .../interfaces/tests/data/lh.pial.surf.gii | 0 .../interfaces/tests/data/lh.roi.shape.gii | 0 .../interfaces/tests/data/lh.white.surf.gii | 0 .../data/sub-01_run-01_echo-1_bold.nii.gz | 0 .../data/sub-01_run-01_echo-2_bold.nii.gz | 0 .../data/sub-01_run-01_echo-3_bold.nii.gz | 0 fmriprep/reports/conftest.py | 43 ------------------- 15 files changed, 40 insertions(+), 80 deletions(-) create mode 100644 fmriprep/interfaces/conftest.py delete mode 100644 fmriprep/interfaces/tests/conftest.py rename fmriprep/{data/sub-01_run-01_echo-1_bold.nii.gz => interfaces/tests/data/bold.nii.gz} (100%) rename fmriprep/{data/sub-01_run-01_echo-2_bold.nii.gz => interfaces/tests/data/lh.bold.func.gii} (100%) rename fmriprep/{data/sub-01_run-01_echo-3_bold.nii.gz => interfaces/tests/data/lh.midthickness.surf.gii} (100%) create mode 100644 fmriprep/interfaces/tests/data/lh.pial.surf.gii create mode 100644 fmriprep/interfaces/tests/data/lh.roi.shape.gii create mode 100644 fmriprep/interfaces/tests/data/lh.white.surf.gii create mode 100644 fmriprep/interfaces/tests/data/sub-01_run-01_echo-1_bold.nii.gz create mode 100644 fmriprep/interfaces/tests/data/sub-01_run-01_echo-2_bold.nii.gz create mode 100644 fmriprep/interfaces/tests/data/sub-01_run-01_echo-3_bold.nii.gz delete mode 100644 fmriprep/reports/conftest.py diff --git a/fmriprep/conftest.py b/fmriprep/conftest.py index fcd5ef6ea..c1ef44abe 100644 --- a/fmriprep/conftest.py +++ b/fmriprep/conftest.py @@ -15,14 +15,6 @@ os.environ['NO_ET'] = '1' -def chdir_or_skip(): - data_dir = ir_files('fmriprep') / 'data' - try: - os.chdir(data_dir) - except OSError: - pytest.skip(f"Cannot chdir into {data_dir!r}. Probably in a zipped distribution.") - - def copytree_or_skip(source, target): data_dir = ir_files('fmriprep') / source if not data_dir.exists(): @@ -36,7 +28,6 @@ def copytree_or_skip(source, target): @pytest.fixture(autouse=True) def populate_namespace(doctest_namespace, tmp_path): - doctest_namespace['chdir_or_skip'] = chdir_or_skip doctest_namespace['copytree_or_skip'] = copytree_or_skip doctest_namespace['testdir'] = tmp_path diff --git a/fmriprep/interfaces/confounds.py b/fmriprep/interfaces/confounds.py index 1c55dc621..694489e63 100644 --- a/fmriprep/interfaces/confounds.py +++ b/fmriprep/interfaces/confounds.py @@ -207,14 +207,6 @@ class GatherConfounds(SimpleInterface): r""" Combine various sources of confounds in one TSV file - .. testsetup:: - - >>> from tempfile import TemporaryDirectory - >>> tmpdir = TemporaryDirectory() - >>> os.chdir(tmpdir.name) - - .. doctest:: - >>> pd.DataFrame({'a': [0.1]}).to_csv('signals.tsv', index=False, na_rep='n/a') >>> pd.DataFrame({'b': [0.2]}).to_csv('dvars.tsv', index=False, na_rep='n/a') @@ -230,10 +222,6 @@ class GatherConfounds(SimpleInterface): a b 0 0.1 0.2 - .. testcleanup:: - - >>> tmpdir.cleanup() - """ input_spec = GatherConfoundsInputSpec output_spec = GatherConfoundsOutputSpec @@ -274,9 +262,6 @@ def _gather_confounds( Load confounds from the filenames, concatenate together horizontally and save new file. - >>> from tempfile import TemporaryDirectory - >>> tmpdir = TemporaryDirectory() - >>> os.chdir(tmpdir.name) >>> pd.DataFrame({'Global Signal': [0.1]}).to_csv('signals.tsv', index=False, na_rep='n/a') >>> pd.DataFrame({'stdDVARS': [0.2]}).to_csv('dvars.tsv', index=False, na_rep='n/a') >>> out_file, confound_list = _gather_confounds('signals.tsv', 'dvars.tsv') @@ -287,7 +272,6 @@ def _gather_confounds( ... engine='python') # doctest: +NORMALIZE_WHITESPACE global_signal std_dvars 0 0.1 0.2 - >>> tmpdir.cleanup() """ diff --git a/fmriprep/interfaces/conftest.py b/fmriprep/interfaces/conftest.py new file mode 100644 index 000000000..9b42b4002 --- /dev/null +++ b/fmriprep/interfaces/conftest.py @@ -0,0 +1,40 @@ +from pathlib import Path +from shutil import copytree + +import pytest + +try: + from contextlib import chdir as _chdir +except ImportError: # PY310 + import os + from contextlib import contextmanager + + @contextmanager # type: ignore + def _chdir(path): + cwd = os.getcwd() + os.chdir(path) + try: + yield + finally: + os.chdir(cwd) + + +@pytest.fixture(scope="module") +def data_dir(): + return Path(__file__).parent / "tests" / "data" + + +@pytest.fixture(autouse=True) +def _docdir(request, tmp_path): + # Trigger ONLY for the doctests. + doctest_plugin = request.config.pluginmanager.getplugin("doctest") + if isinstance(request.node, doctest_plugin.DoctestItem): + copytree(Path(__file__).parent / "tests" / "data", tmp_path, dirs_exist_ok=True) + + # Chdir only for the duration of the test. + with _chdir(tmp_path): + yield + + else: + # For normal tests, we have to yield, since this is a yield-fixture. + yield diff --git a/fmriprep/interfaces/multiecho.py b/fmriprep/interfaces/multiecho.py index cb031f05a..9683d20d7 100644 --- a/fmriprep/interfaces/multiecho.py +++ b/fmriprep/interfaces/multiecho.py @@ -83,10 +83,6 @@ class T2SMap(CommandLine): Example ======= - .. testsetup:: - - >>> chdir_or_skip() - >>> from fmriprep.interfaces import multiecho >>> t2smap = multiecho.T2SMap() >>> t2smap.inputs.in_files = ['sub-01_run-01_echo-1_bold.nii.gz', diff --git a/fmriprep/interfaces/tests/conftest.py b/fmriprep/interfaces/tests/conftest.py deleted file mode 100644 index 4f0821826..000000000 --- a/fmriprep/interfaces/tests/conftest.py +++ /dev/null @@ -1,8 +0,0 @@ -from pathlib import Path - -import pytest - - -@pytest.fixture(scope="module") -def data_dir(): - return Path(__file__).parent / "data" diff --git a/fmriprep/data/sub-01_run-01_echo-1_bold.nii.gz b/fmriprep/interfaces/tests/data/bold.nii.gz similarity index 100% rename from fmriprep/data/sub-01_run-01_echo-1_bold.nii.gz rename to fmriprep/interfaces/tests/data/bold.nii.gz diff --git a/fmriprep/data/sub-01_run-01_echo-2_bold.nii.gz b/fmriprep/interfaces/tests/data/lh.bold.func.gii similarity index 100% rename from fmriprep/data/sub-01_run-01_echo-2_bold.nii.gz rename to fmriprep/interfaces/tests/data/lh.bold.func.gii diff --git a/fmriprep/data/sub-01_run-01_echo-3_bold.nii.gz b/fmriprep/interfaces/tests/data/lh.midthickness.surf.gii similarity index 100% rename from fmriprep/data/sub-01_run-01_echo-3_bold.nii.gz rename to fmriprep/interfaces/tests/data/lh.midthickness.surf.gii diff --git a/fmriprep/interfaces/tests/data/lh.pial.surf.gii b/fmriprep/interfaces/tests/data/lh.pial.surf.gii new file mode 100644 index 000000000..e69de29bb diff --git a/fmriprep/interfaces/tests/data/lh.roi.shape.gii b/fmriprep/interfaces/tests/data/lh.roi.shape.gii new file mode 100644 index 000000000..e69de29bb diff --git a/fmriprep/interfaces/tests/data/lh.white.surf.gii b/fmriprep/interfaces/tests/data/lh.white.surf.gii new file mode 100644 index 000000000..e69de29bb diff --git a/fmriprep/interfaces/tests/data/sub-01_run-01_echo-1_bold.nii.gz b/fmriprep/interfaces/tests/data/sub-01_run-01_echo-1_bold.nii.gz new file mode 100644 index 000000000..e69de29bb diff --git a/fmriprep/interfaces/tests/data/sub-01_run-01_echo-2_bold.nii.gz b/fmriprep/interfaces/tests/data/sub-01_run-01_echo-2_bold.nii.gz new file mode 100644 index 000000000..e69de29bb diff --git a/fmriprep/interfaces/tests/data/sub-01_run-01_echo-3_bold.nii.gz b/fmriprep/interfaces/tests/data/sub-01_run-01_echo-3_bold.nii.gz new file mode 100644 index 000000000..e69de29bb diff --git a/fmriprep/reports/conftest.py b/fmriprep/reports/conftest.py deleted file mode 100644 index a4c200744..000000000 --- a/fmriprep/reports/conftest.py +++ /dev/null @@ -1,43 +0,0 @@ -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -# -# Copyright 2023 The NiPreps Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# We support and encourage derived works from this project, please read -# about our expectations at -# -# https://www.nipreps.org/community/licensing/ -# -"""py.test configuration""" -import os -import tempfile -from pathlib import Path - -import pytest - - -@pytest.fixture(autouse=True) -def populate_doctest_namespace(doctest_namespace): - doctest_namespace["os"] = os - doctest_namespace["Path"] = Path - tmpdir = tempfile.TemporaryDirectory() - - doctest_namespace["tmpdir"] = tmpdir.name - - cwd = os.getcwd() - os.chdir(tmpdir.name) - yield - os.chdir(cwd) - tmpdir.cleanup() From 297d78edd995900d5580db2dd745400673fbd667 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Fri, 19 May 2023 09:32:37 -0400 Subject: [PATCH 11/21] FIX: Sort surfaces by basename, not full path --- fmriprep/workflows/bold/resampling.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index 7ba93eba1..139d33cac 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -1327,6 +1327,8 @@ def _select_surfaces( template_spheres, template_rois, ): + # This function relies on the basenames of the files to differ by L/R or l/r + # so that the sorting correctly identifies left or right. import os import re @@ -1345,4 +1347,4 @@ def _select_surfaces( match = find_name.search(os.path.basename(surface)) if match: container[match.group('name')].append(surface) - return tuple(sorted(surflist)[idx] for surflist in container.values()) + return tuple(sorted(surflist, key=os.path.basename)[idx] for surflist in container.values()) From 68b61401d4d41d574700a08131c71d2352db387c Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Fri, 19 May 2023 12:45:25 -0400 Subject: [PATCH 12/21] DOCTEST: Explicit normalize_whitespace for --pyargs friendliness --- fmriprep/interfaces/workbench.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/fmriprep/interfaces/workbench.py b/fmriprep/interfaces/workbench.py index b43ae3ab0..f5d35e53c 100644 --- a/fmriprep/interfaces/workbench.py +++ b/fmriprep/interfaces/workbench.py @@ -579,7 +579,7 @@ class VolumeToSurfaceMapping(WBCommand, OpenMPCommandMixin): >>> vol2surf.inputs.method = 'ribbon-constrained' >>> vol2surf.inputs.inner_surface = 'lh.white.surf.gii' >>> vol2surf.inputs.outer_surface = 'lh.pial.surf.gii' - >>> vol2surf.cmdline + >>> vol2surf.cmdline # doctest: +NORMALIZE_WHITESPACE 'wb_command -volume-to-surface-mapping bold.nii.gz lh.midthickness.surf.gii \ lh.midthickness.surf_mapped.func.gii -ribbon-constrained lh.white.surf.gii lh.pial.surf.gii' """ @@ -732,7 +732,7 @@ class MetricFillHoles(WBCommand): >>> fill_holes = MetricFillHoles() >>> fill_holes.inputs.surface_file = 'lh.midthickness.surf.gii' >>> fill_holes.inputs.metric_file = 'lh.roi.shape.gii' - >>> fill_holes.cmdline + >>> fill_holes.cmdline # doctest: +NORMALIZE_WHITESPACE 'wb_command -metric-fill-holes lh.midthickness.surf.gii lh.roi.shape.gii \ lh.roi.shape_filled.shape.gii' """ @@ -799,7 +799,7 @@ class MetricRemoveIslands(WBCommand): >>> remove_islands = MetricRemoveIslands() >>> remove_islands.inputs.surface_file = 'lh.midthickness.surf.gii' >>> remove_islands.inputs.metric_file = 'lh.roi.shape.gii' - >>> remove_islands.cmdline + >>> remove_islands.cmdline # doctest: +NORMALIZE_WHITESPACE 'wb_command -metric-remove-islands lh.midthickness.surf.gii \ lh.roi.shape.gii lh.roi.shape_noislands.shape.gii' """ From 47c7b28332912a8089a21e876b5004759b2ab9f6 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Fri, 2 Jun 2023 14:22:18 -0400 Subject: [PATCH 13/21] FIX: Update sloppy/debug parameters in anatomical base workflow --- fmriprep/workflows/base.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index 92ac1b935..4d1672482 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -305,7 +305,8 @@ def init_single_subject_wf(subject_id: str): # Preprocessing of T1w (includes registration to MNI) anat_preproc_wf = init_anat_preproc_wf( bids_root=str(config.execution.bids_dir), - debug=config.execution.sloppy, + sloppy=config.execution.sloppy, + debug=config.execution.debug, existing_derivatives=anat_derivatives, freesurfer=config.workflow.run_reconall, hires=config.workflow.hires, From 64eecd6bc9451aa1f81c052a1bb4f9e9de890b56 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Fri, 2 Jun 2023 16:09:23 -0400 Subject: [PATCH 14/21] CI: Update expected outputs --- .circleci/ds005_bids_outputs.txt | 8 ++++++-- .circleci/ds005_legacy_outputs.txt | 8 ++++++-- .circleci/ds005_legacy_partial_outputs.txt | 8 ++++++-- 3 files changed, 18 insertions(+), 6 deletions(-) diff --git a/.circleci/ds005_bids_outputs.txt b/.circleci/ds005_bids_outputs.txt index e7a635adf..e7f473596 100644 --- a/.circleci/ds005_bids_outputs.txt +++ b/.circleci/ds005_bids_outputs.txt @@ -23,19 +23,23 @@ sub-01/anat/sub-01_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5 sub-01/anat/sub-01_from-T1w_to-fsnative_mode-image_xfm.txt sub-01/anat/sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5 sub-01/anat/sub-01_hemi-L_curv.shape.gii +sub-01/anat/sub-01_hemi-L_desc-reg_sphere.surf.gii sub-01/anat/sub-01_hemi-L_inflated.surf.gii sub-01/anat/sub-01_hemi-L_midthickness.surf.gii sub-01/anat/sub-01_hemi-L_pial.surf.gii -sub-01/anat/sub-01_hemi-L_smoothwm.surf.gii +sub-01/anat/sub-01_hemi-L_space-fsLR_desc-reg_sphere.surf.gii sub-01/anat/sub-01_hemi-L_sulc.shape.gii sub-01/anat/sub-01_hemi-L_thickness.shape.gii +sub-01/anat/sub-01_hemi-L_white.surf.gii sub-01/anat/sub-01_hemi-R_curv.shape.gii +sub-01/anat/sub-01_hemi-R_desc-reg_sphere.surf.gii sub-01/anat/sub-01_hemi-R_inflated.surf.gii sub-01/anat/sub-01_hemi-R_midthickness.surf.gii sub-01/anat/sub-01_hemi-R_pial.surf.gii -sub-01/anat/sub-01_hemi-R_smoothwm.surf.gii +sub-01/anat/sub-01_hemi-R_space-fsLR_desc-reg_sphere.surf.gii sub-01/anat/sub-01_hemi-R_sulc.shape.gii sub-01/anat/sub-01_hemi-R_thickness.shape.gii +sub-01/anat/sub-01_hemi-R_white.surf.gii sub-01/anat/sub-01_label-CSF_probseg.nii.gz sub-01/anat/sub-01_label-GM_probseg.nii.gz sub-01/anat/sub-01_label-WM_probseg.nii.gz diff --git a/.circleci/ds005_legacy_outputs.txt b/.circleci/ds005_legacy_outputs.txt index 483805e82..a094e1d91 100644 --- a/.circleci/ds005_legacy_outputs.txt +++ b/.circleci/ds005_legacy_outputs.txt @@ -24,19 +24,23 @@ fmriprep/sub-01/anat/sub-01_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5 fmriprep/sub-01/anat/sub-01_from-T1w_to-fsnative_mode-image_xfm.txt fmriprep/sub-01/anat/sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5 fmriprep/sub-01/anat/sub-01_hemi-L_curv.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-L_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_inflated.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_midthickness.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_pial.surf.gii -fmriprep/sub-01/anat/sub-01_hemi-L_smoothwm.surf.gii +fmriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_sulc.shape.gii fmriprep/sub-01/anat/sub-01_hemi-L_thickness.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-L_white.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_curv.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-R_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_inflated.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_midthickness.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_pial.surf.gii -fmriprep/sub-01/anat/sub-01_hemi-R_smoothwm.surf.gii +fmriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_sulc.shape.gii fmriprep/sub-01/anat/sub-01_hemi-R_thickness.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-R_white.surf.gii fmriprep/sub-01/anat/sub-01_label-CSF_probseg.nii.gz fmriprep/sub-01/anat/sub-01_label-GM_probseg.nii.gz fmriprep/sub-01/anat/sub-01_label-WM_probseg.nii.gz diff --git a/.circleci/ds005_legacy_partial_outputs.txt b/.circleci/ds005_legacy_partial_outputs.txt index 272a83529..453f333b5 100644 --- a/.circleci/ds005_legacy_partial_outputs.txt +++ b/.circleci/ds005_legacy_partial_outputs.txt @@ -26,19 +26,23 @@ fmriprep/sub-01/anat/sub-01_from-T1w_to-fsnative_mode-image_xfm.txt fmriprep/sub-01/anat/sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5 fmriprep/sub-01/anat/sub-01_from-T1w_to-MNI152NLin6Asym_mode-image_xfm.h5 fmriprep/sub-01/anat/sub-01_hemi-L_curv.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-L_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_inflated.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_midthickness.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_pial.surf.gii -fmriprep/sub-01/anat/sub-01_hemi-L_smoothwm.surf.gii +fmriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-L_sulc.shape.gii fmriprep/sub-01/anat/sub-01_hemi-L_thickness.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-L_white.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_curv.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-R_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_inflated.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_midthickness.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_pial.surf.gii -fmriprep/sub-01/anat/sub-01_hemi-R_smoothwm.surf.gii +fmriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_desc-reg_sphere.surf.gii fmriprep/sub-01/anat/sub-01_hemi-R_sulc.shape.gii fmriprep/sub-01/anat/sub-01_hemi-R_thickness.shape.gii +fmriprep/sub-01/anat/sub-01_hemi-R_white.surf.gii fmriprep/sub-01/anat/sub-01_label-CSF_probseg.nii.gz fmriprep/sub-01/anat/sub-01_label-GM_probseg.nii.gz fmriprep/sub-01/anat/sub-01_label-WM_probseg.nii.gz From a247fc254b1b7a47868f2b0016c3f8457dc11485 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 5 Jun 2023 13:44:57 -0400 Subject: [PATCH 15/21] Cleanup old comments and disjointed workflow connections --- fmriprep/workflows/bold/resampling.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index 139d33cac..98b4097e1 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -199,6 +199,7 @@ def select_target(subject_id, space): ("subject_id", "subject_id"), ]), (itersource, targets, [("target", "space")]), + (inputnode, rename_src, [("source_file", "in_file")]), (itersource, rename_src, [("target", "subject")]), (rename_src, sampler, [("out_file", "source_file")]), (itk2lta, sampler, [("out", "reg_file")]), @@ -212,21 +213,11 @@ def select_target(subject_id, space): ]) # fmt: on - # At this point, rename_src.in_file and update_metadata.in_file need connecting - # - # These depend on two optional steps: goodvoxel projection and medial wall nan replacement - # - # inputnode -> rename_src - # sampler -> optional(medial_nans) -> update_metadata - # - # Refine if medial vertices should be NaNs medial_nans = pe.MapNode( MedialNaNs(), iterfield=["in_file"], name="medial_nans", mem_gb=DEFAULT_MEMORY_MIN_GB ) - workflow.connect([(inputnode, rename_src, [("source_file", "in_file")])]) - if medial_surface_nan: # fmt: off workflow.connect([ From 606ed37dd8c320bd3ab474372e8c9ecccd9e5bc6 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 5 Jun 2023 13:45:25 -0400 Subject: [PATCH 16/21] REPORT: Restore goodvoxels reporting, with slight clarification --- fmriprep/workflows/bold/resampling.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index 98b4097e1..d3edbe611 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -526,6 +526,11 @@ def init_bold_fsLR_resampling_wf( workflow = Workflow(name=name) + workflow.__desc__ = """\ +The BOLD time-series were resampled onto the left/right-symmetric template +"fsLR" [@hcppipelines]. +""" + inputnode = pe.Node( niu.IdentityInterface( fields=[ @@ -670,6 +675,11 @@ def init_bold_fsLR_resampling_wf( # fmt: on if estimate_goodvoxels: + workflow.__desc__ += """\ +A "goodvoxels" mask was applied during volume-to-surface sampling in fsLR space, +excluding voxels whose time-series have a locally high coefficient of variation. +""" + goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb) # fmt: off From 926a03b3b5a8e2f995eac3394da5c79c703fd8df Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 5 Jun 2023 15:45:14 -0400 Subject: [PATCH 17/21] ENH: Output goodvoxels mask (not constrained to ribbon) --- fmriprep/workflows/bold/base.py | 7 +++++-- fmriprep/workflows/bold/outputs.py | 10 +++++----- fmriprep/workflows/bold/resampling.py | 7 +++++-- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index cc24095ff..247a9b983 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -211,7 +211,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): spaces = config.workflow.spaces fmriprep_dir = str(config.execution.fmriprep_dir) freesurfer_spaces = spaces.get_fs_spaces() - project_goodvoxels = config.workflow.project_goodvoxels + project_goodvoxels = config.workflow.project_goodvoxels and config.workflow.cifti_output # Extract BIDS entities and metadata from BOLD file(s) entities = extract_entities(bold_file) @@ -433,7 +433,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): bids_root=layout.root, cifti_output=config.workflow.cifti_output, freesurfer=freesurfer, - project_goodvoxels=False, + project_goodvoxels=project_goodvoxels, all_metadata=all_metadata, multiecho=multiecho, output_dir=fmriprep_dir, @@ -918,6 +918,9 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): (bold_fsLR_resampling_wf, bold_grayords_wf, [ ("outputnode.bold_fsLR", "inputnode.bold_fsLR"), ]), + (bold_fsLR_resampling_wf, func_derivatives_wf, [ + ("outputnode.goodvoxels_mask", "inputnode.goodvoxels_mask"), + ]), (bold_grayords_wf, outputnode, [ ("outputnode.cifti_bold", "bold_cifti"), ("outputnode.cifti_metadata", "cifti_metadata"), diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index a5a4ec638..768f5d5d7 100644 --- a/fmriprep/workflows/bold/outputs.py +++ b/fmriprep/workflows/bold/outputs.py @@ -226,7 +226,7 @@ def init_func_derivatives_wf( 'cifti_density', 'confounds', 'confounds_metadata', - 'goodvoxels_ribbon', + 'goodvoxels_mask', 'source_file', 'all_source_files', 'surf_files', @@ -746,7 +746,7 @@ def init_func_derivatives_wf( # fmt:on if freesurfer and project_goodvoxels: - ds_goodvoxels_ribbon = pe.Node( + ds_goodvoxels_mask = pe.Node( DerivativesDataSink( base_directory=output_dir, space='T1w', @@ -755,15 +755,15 @@ def init_func_derivatives_wf( compress=True, dismiss_entities=("echo",), ), - name='ds_goodvoxels_ribbon', + name='ds_goodvoxels_mask', run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB, ) # fmt:off workflow.connect([ - (inputnode, ds_goodvoxels_ribbon, [ + (inputnode, ds_goodvoxels_mask, [ ('source_file', 'source_file'), - ('goodvoxels_ribbon', 'in_file'), + ('goodvoxels_mask', 'in_file'), ('surf_refs', 'keys')]), ]) # fmt:on diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index d3edbe611..52134cf88 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -551,7 +551,7 @@ def init_bold_fsLR_resampling_wf( ) outputnode = pe.JoinNode( - niu.IdentityInterface(fields=['bold_fsLR']), + niu.IdentityInterface(fields=['bold_fsLR', 'goodvoxels_mask']), name='outputnode', joinsource='itersource', ) @@ -689,7 +689,10 @@ def init_bold_fsLR_resampling_wf( ("anat_ribbon", "inputnode.anat_ribbon"), ]), (goodvoxels_bold_mask_wf, volume_to_surface, [ - ("outputnode.goodvoxels_ribbon", "volume_roi"), + ("outputnode.goodvoxels_mask", "volume_roi"), + ]), + (goodvoxels_bold_mask_wf, outputnode, [ + ("outputnode.goodvoxels_mask", "goodvoxels_mask"), ]), ]) # fmt: on From ee3ebc259d37881e3ad9f2ef08ad008a321a2a2b Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 5 Jun 2023 16:29:41 -0400 Subject: [PATCH 18/21] DOC: Update docs and workflow docstring --- docs/workflows.rst | 28 ++++++++++++---- fmriprep/workflows/bold/resampling.py | 48 +++++++++++++++++++++++++++ 2 files changed, 70 insertions(+), 6 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index d5da471f0..35e2e616c 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -534,12 +534,28 @@ All surface outputs are in GIFTI format. HCP Grayordinates ~~~~~~~~~~~~~~~~~ -If CIFTI output is enabled, the motion-corrected functional timeseries (in T1w space) is first -sampled to the high resolution 164k vertex (per hemisphere) ``fsaverage``. Following that, -the resampled timeseries is sampled to `HCP Pipelines_`'s ``fsLR`` mesh (with the left and -right hemisphere aligned) using `Connectome Workbench`_'s ``-metric-resample`` to generate a -surface timeseries for each hemisphere. These surfaces are then combined with corresponding -volumetric timeseries to create a CIFTI2 file. +:py:func:`~fmriprep.workflows.bold.resampling.init_bold_fsLR_resampling_wf` + +.. workflow:: + :graph2use: colored + :simple_form: yes + + from fmriprep.workflows.bold.resampling import init_bold_fsLR_resampling_wf + wf = init_bold_fsLR_resampling_wf( + estimate_goodvoxels=True, + grayord_density='92k', + omp_nthreads=1, + mem_gb=1, + ) + +If CIFTI output is enabled, the motion-corrected functional timeseries (in T1w space) is +resampled onto the subject-native surface, optionally using the `HCP Pipelines_`'s +"goodvoxels" masking method to exclude voxels with local peaks of temporal variation. +After dilating the surface-sampled time series to fill sampling holes, the result is +resampled to the ``fsLR`` mesh (with the left and right hemisphere aligned). +These workflows make use of various `Connectome Workbench`_ functions. +These surfaces are then combined with corresponding volumetric timeseries to create a +CIFTI-2 file. .. _bold_confounds: diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index 52134cf88..545887f6b 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -509,6 +509,54 @@ def init_bold_fsLR_resampling_wf( Line numbers correspond to the locations of the code in the original scripts, found at: https://github.com/DCAN-Labs/DCAN-HCP/tree/9291324/ + + Workflow Graph + .. workflow:: + :graph2use: colored + :simple_form: yes + + from fmriprep.workflows.bold.resampling import init_bold_fsLR_resampling_wf + wf = init_bold_fsLR_resampling_wf( + estimate_goodvoxels=True, + grayord_density='92k', + omp_nthreads=1, + mem_gb=1, + ) + + Parameters + ---------- + grayord_density : :class:`str` + Either `91k` or `170k`, representing the total of vertices or *grayordinates*. + estimate_goodvoxels : :class:`bool` + Calculate mask excluding voxels with a locally high coefficient of variation to + exclude from surface resampling + omp_nthreads : :class:`int` + Maximum number of threads an individual process may use + mem_gb : :class:`float` + Size of BOLD file in GB + name : :class:`str` + Name of workflow (default: ``bold_fsLR_resampling_wf``) + + Inputs + ------ + bold_file : :class:`str` + Path to BOLD file resampled into T1 space + surfaces : :class:`list` of :class:`str + Path to left and right hemisphere white, pial and midthickness GIFTI surfaces + morphometrics : :class:`list` of :class:`str + Path to left and right hemisphere morphometric GIFTI surfaces, which must include thickness + sphere_reg_fsLR : :class:`list` of :class:`str + Path to left and right hemisphere sphere.reg GIFTI surfaces, mapping from subject to fsLR + anat_ribbon : :class:`str` + Path to mask of cortical ribbon in T1w space, for calculating goodvoxels + + Outputs + ------- + bold_fsLR : :class:`list` of :class:`str + Path to BOLD series resampled as functional GIFTI files in fsLR space + goodvoxels_mask : :class:`str` + Path to mask of voxels, excluding those with locally high coefficients of variation + """ import templateflow.api as tf from niworkflows.engine.workflows import LiterateWorkflow as Workflow From 63092fd7c8c6c1bcadcc0b38cee292ef3c4884fa Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 5 Jun 2023 16:45:30 -0400 Subject: [PATCH 19/21] DOC: Clean up Sphinx warnings --- fmriprep/workflows/bold/resampling.py | 30 +++++++++++++++------------ 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index 545887f6b..06683b51f 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -27,6 +27,9 @@ .. autofunction:: init_bold_surf_wf .. autofunction:: init_bold_std_trans_wf .. autofunction:: init_bold_preproc_trans_wf +.. autofunction:: init_bold_fsLR_resampling_wf +.. autofunction:: init_bold_grayords_wf +.. autofunction:: init_goodvoxels_bold_mask_wf """ from __future__ import annotations @@ -233,14 +236,15 @@ def select_target(subject_id, space): def init_goodvoxels_bold_mask_wf(mem_gb: float, name: str = "goodvoxels_bold_mask_wf"): - """ + """Calculate a mask of a BOLD series excluding high variance voxels. + Workflow Graph .. workflow:: :graph2use: colored :simple_form: yes - from fmriprep.workflows.bold import init_goodvoxels_mask_wf - wf = init_goodvoxels_mask_wf(mem_gb=0.1) + from fmriprep.workflows.bold.resampling import init_goodvoxels_bold_mask_wf + wf = init_goodvoxels_bold_mask_wf(mem_gb=0.1) Parameters ---------- @@ -526,7 +530,7 @@ def init_bold_fsLR_resampling_wf( Parameters ---------- grayord_density : :class:`str` - Either `91k` or `170k`, representing the total of vertices or *grayordinates*. + Either ``"91k"`` or ``"170k"``, representing the total *grayordinates*. estimate_goodvoxels : :class:`bool` Calculate mask excluding voxels with a locally high coefficient of variation to exclude from surface resampling @@ -541,18 +545,18 @@ def init_bold_fsLR_resampling_wf( ------ bold_file : :class:`str` Path to BOLD file resampled into T1 space - surfaces : :class:`list` of :class:`str + surfaces : :class:`list` of :class:`str` Path to left and right hemisphere white, pial and midthickness GIFTI surfaces - morphometrics : :class:`list` of :class:`str + morphometrics : :class:`list` of :class:`str` Path to left and right hemisphere morphometric GIFTI surfaces, which must include thickness - sphere_reg_fsLR : :class:`list` of :class:`str + sphere_reg_fsLR : :class:`list` of :class:`str` Path to left and right hemisphere sphere.reg GIFTI surfaces, mapping from subject to fsLR anat_ribbon : :class:`str` Path to mask of cortical ribbon in T1w space, for calculating goodvoxels Outputs ------- - bold_fsLR : :class:`list` of :class:`str + bold_fsLR : :class:`list` of :class:`str` Path to BOLD series resampled as functional GIFTI files in fsLR space goodvoxels_mask : :class:`str` Path to mask of voxels, excluding those with locally high coefficients of variation @@ -1222,13 +1226,13 @@ def init_bold_grayords_wf( :graph2use: colored :simple_form: yes - from fmriprep.workflows.bold import init_bold_grayords_wf - wf = init_bold_grayords_wf(mem_gb=0.1, grayord_density="91k") + from fmriprep.workflows.bold.resampling import init_bold_grayords_wf + wf = init_bold_grayords_wf(mem_gb=0.1, grayord_density="91k", repetition_time=2) Parameters ---------- - grayord_density : :obj:`str` - Either `91k` or `170k`, representing the total of vertices or *grayordinates*. + grayord_density : :class:`str` + Either ``"91k"`` or ``"170k"``, representing the total *grayordinates*. mem_gb : :obj:`float` Size of BOLD file in GB name : :obj:`str` @@ -1238,7 +1242,7 @@ def init_bold_grayords_wf( ------ bold_std : :obj:`str` List of BOLD conversions to standard spaces. - spatial_reference :obj:`str` + spatial_reference : :obj:`str` List of unique identifiers corresponding to the BOLD standard-conversions. surf_files : :obj:`str` List of BOLD files resampled on the fsaverage (ico7) surfaces. From 6e4e50e7f8c0f074295d24a6ba63c6843c661733 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 6 Jun 2023 08:54:22 -0400 Subject: [PATCH 20/21] FIX: Separate joinnode from outputnode --- fmriprep/workflows/bold/resampling.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index 06683b51f..09fdbce40 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -602,10 +602,15 @@ def init_bold_fsLR_resampling_wf( iterables=[('hemi', ['L', 'R'])], ) - outputnode = pe.JoinNode( + joinnode = pe.JoinNode( + niu.IdentityInterface(fields=['bold_fsLR']), + name='joinnode', + joinsource='itersource', + ) + + outputnode = pe.Node( niu.IdentityInterface(fields=['bold_fsLR', 'goodvoxels_mask']), name='outputnode', - joinsource='itersource', ) # select white, midthickness and pial surfaces based on hemi @@ -722,7 +727,8 @@ def init_bold_fsLR_resampling_wf( (select_surfaces, mask_fsLR, [('template_roi', 'mask')]), (resample_to_fsLR, mask_fsLR, [('out_file', 'in_file')]), # Output - (mask_fsLR, outputnode, [('out_file', 'bold_fsLR')]), + (mask_fsLR, joinnode, [('out_file', 'bold_fsLR')]), + (joinnode, outputnode, [('bold_fsLR', 'bold_fsLR')]), ]) # fmt: on From 3a3d9b4642e714631ff281c212dacc7c3d65fefd Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 6 Jun 2023 12:46:59 -0400 Subject: [PATCH 21/21] CI: Expect goodvoxels mask --- .circleci/ds005_legacy_partial_fasttrack_outputs.txt | 2 ++ .circleci/ds005_legacy_partial_outputs.txt | 2 ++ 2 files changed, 4 insertions(+) diff --git a/.circleci/ds005_legacy_partial_fasttrack_outputs.txt b/.circleci/ds005_legacy_partial_fasttrack_outputs.txt index a864b6cd2..569d5fcd5 100644 --- a/.circleci/ds005_legacy_partial_fasttrack_outputs.txt +++ b/.circleci/ds005_legacy_partial_fasttrack_outputs.txt @@ -52,6 +52,8 @@ fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-aparcase fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-aseg_dseg.nii.gz fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-brain_mask.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-brain_mask.nii.gz +fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-goodvoxels_mask.json +fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-goodvoxels_mask.nii.gz fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-preproc_bold.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-preproc_bold.nii.gz fmriprep/sub-01.html diff --git a/.circleci/ds005_legacy_partial_outputs.txt b/.circleci/ds005_legacy_partial_outputs.txt index 453f333b5..e920a1822 100644 --- a/.circleci/ds005_legacy_partial_outputs.txt +++ b/.circleci/ds005_legacy_partial_outputs.txt @@ -110,6 +110,8 @@ fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-aparcase fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-aseg_dseg.nii.gz fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-brain_mask.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-brain_mask.nii.gz +fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-goodvoxels_mask.json +fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-goodvoxels_mask.nii.gz fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-preproc_bold.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-T1w_desc-preproc_bold.nii.gz fmriprep/sub-01.html