Skip to content

Commit

Permalink
ENH: Add option to exclude projecting high variance voxels to surface…
Browse files Browse the repository at this point in the history
… (update of #2855) (#2956)

Adds option --project-goodvoxels to address #2822.

If enabled, voxels whose timeseries have locally high coefficient of
variation, or which lie outside the cortex (based on signed distance
volumes generated from the WM and pial surfaces), are excluded from
volume-to-surface sampling of BOLD timeseries. Not enabled by default.

- adds MetricDilate interface for Workbench's wb_command -metric-dilate,
which, when using the --project-goodvoxels option, fills zero-value
vertices in the BOLD surface timeseries by dilating in data from the
nearest "good" (nonzero, nonempty) vertex.
  • Loading branch information
effigies authored Mar 6, 2023
2 parents 084c5d0 + 2d802cb commit 11be54e
Show file tree
Hide file tree
Showing 12 changed files with 738 additions and 37 deletions.
2 changes: 1 addition & 1 deletion .circleci/config.yml
Original file line number Diff line number Diff line change
Expand Up @@ -609,7 +609,7 @@ jobs:
--sloppy --write-graph --use-syn-sdc --mem-mb 14336 \
--output-spaces MNI152NLin2009cAsym fsaverage5 fsnative MNI152NLin6Asym anat \
--aroma-melodic-dimensionality 2 --use-aroma \
--nthreads 4 --cifti-output -vv
--nthreads 4 --cifti-output --project-goodvoxels -vv
- store_artifacts:
path: /tmp/ds005/derivatives_partial
destination: partial-run
Expand Down
2 changes: 2 additions & 0 deletions Dockerfile
Original file line number Diff line number Diff line change
Expand Up @@ -280,6 +280,8 @@ ENV HOME="/home/fmriprep" \
RUN echo ". /opt/conda/etc/profile.d/conda.sh" >> $HOME/.bashrc && \
echo "conda activate base" >> $HOME/.bashrc

RUN /opt/conda/bin/python -m pip install --no-cache-dir --upgrade templateflow

# Precaching atlases
COPY scripts/fetch_templates.py fetch_templates.py

Expand Down
9 changes: 9 additions & 0 deletions fmriprep/cli/parser.py
Original file line number Diff line number Diff line change
Expand Up @@ -386,6 +386,15 @@ def _slice_time_ref(value, parser):
help="Replace medial wall values with NaNs on functional GIFTI files. Only "
"performed for GIFTI files mapped to a freesurfer subject (fsaverage or fsnative).",
)
g_conf.add_argument(
"--project-goodvoxels",
required=False,
action="store_true",
default=False,
help="Exclude voxels whose timeseries have locally high coefficient of variation "
"from surface resampling. Only performed for GIFTI files mapped to a freesurfer subject "
"(fsaverage or fsnative).",
)
g_outputs.add_argument(
"--md-only-boilerplate",
action="store_true",
Expand Down
2 changes: 2 additions & 0 deletions fmriprep/config.py
Original file line number Diff line number Diff line change
Expand Up @@ -539,6 +539,8 @@ class workflow(_Config):
"""Run FreeSurfer ``recon-all`` with the ``-logitudinal`` flag."""
medial_surface_nan = None
"""Fill medial surface with :abbr:`NaNs (not-a-number)` when sampling."""
project_goodvoxels = False
"""Exclude voxels with locally high coefficient of variation from sampling."""
regressors_all_comps = None
"""Return all CompCor components."""
regressors_dvars_th = None
Expand Down
1 change: 1 addition & 0 deletions fmriprep/data/tests/config.toml
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ hires = true
ignore = []
longitudinal = false
medial_surface_nan = false
project_goodvoxels = false
regressors_all_comps = false
regressors_dvars_th = 1.5
regressors_fd_th = 0.5
Expand Down
271 changes: 271 additions & 0 deletions fmriprep/interfaces/metric.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,271 @@
# -*- 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 <current-sphere> 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 <current-sphere> mesh OR"
" a metric file with vertex areas for <current-sphere> mesh",
)
new_area = File(
exists=True,
position=7,
argstr="%s",
desc="A relevant anatomical surface with <current-sphere> mesh OR"
" a metric file with vertex areas for <current-sphere> 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
1 change: 1 addition & 0 deletions fmriprep/workflows/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -444,6 +444,7 @@ def init_single_subject_wf(subject_id: str):
# Undefined if --fs-no-reconall, but this is safe
('outputnode.subjects_dir', 'inputnode.subjects_dir'),
('outputnode.subject_id', 'inputnode.subject_id'),
('outputnode.anat_ribbon', 'inputnode.anat_ribbon'),
('outputnode.t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'),
('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm')]),
])
Expand Down
17 changes: 16 additions & 1 deletion fmriprep/workflows/bold/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -217,6 +217,15 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
freesurfer = config.workflow.run_reconall
spaces = config.workflow.spaces
fmriprep_dir = str(config.execution.fmriprep_dir)
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)
Expand Down Expand Up @@ -348,6 +357,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
"anat2std_xfm",
"std2anat_xfm",
"template",
"anat_ribbon",
"t1w2fsnative_xfm",
"fsnative2t1w_xfm",
"fmap",
Expand Down Expand Up @@ -437,6 +447,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,
all_metadata=all_metadata,
multiecho=multiecho,
output_dir=fmriprep_dir,
Expand Down Expand Up @@ -936,13 +947,13 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):

# SURFACES ##################################################################################
# Freesurfer
freesurfer_spaces = spaces.get_fs_spaces()
if freesurfer and freesurfer_spaces:
config.loggers.workflow.debug("Creating BOLD surface-sampling workflow.")
bold_surf_wf = init_bold_surf_wf(
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
Expand All @@ -951,10 +962,14 @@ 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

Expand Down
Loading

0 comments on commit 11be54e

Please sign in to comment.