Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: Restore CIFTI-2 generation #3129

Merged
merged 6 commits into from
Nov 13, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 6 additions & 2 deletions fmriprep/interfaces/gifti.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@


class CreateROI(SimpleInterface):
"""Prepare GIFTI shape file for use in"""
"""Prepare GIFTI thickness file for use as a cortical ROI"""

input_spec = CreateROIInputSpec
output_spec = CreateROIOutputSpec
Expand All @@ -46,17 +46,21 @@
# wb_command -metric-math "abs(var * -1) > 0"
roi = np.abs(darray.data) > 0

# Divergence: Set datatype to uint8, since the values are boolean
# wb_command sets datatype to float32
darray = nb.gifti.GiftiDataArray(
roi,
intent=darray.intent,
datatype=darray.datatype,
datatype='uint8',
encoding=darray.encoding,
endian=darray.endian,
coordsys=darray.coordsys,
ordering=darray.ind_ord,
meta=meta,
)

img.darrays[0] = darray

Check warning on line 62 in fmriprep/interfaces/gifti.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/interfaces/gifti.py#L62

Added line #L62 was not covered by tests

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
Expand Down
70 changes: 58 additions & 12 deletions fmriprep/workflows/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -147,6 +147,7 @@
from niworkflows.engine.workflows import LiterateWorkflow as Workflow
from niworkflows.interfaces.bids import BIDSDataGrabber, BIDSInfo
from niworkflows.interfaces.nilearn import NILEARN_VERSION
from niworkflows.interfaces.utility import KeySelect

Check warning on line 150 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L150

Added line #L150 was not covered by tests
from niworkflows.utils.bids import collect_data
from niworkflows.utils.misc import fix_multi_T1w_source_name
from niworkflows.utils.spaces import Reference
Expand Down Expand Up @@ -332,8 +333,9 @@
]) # fmt:skip

# Set up the template iterator once, if used
template_iterator_wf = None

Check warning on line 336 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L336

Added line #L336 was not covered by tests
if config.workflow.level == "full":
if spaces.get_spaces(nonstandard=False, dim=(3,)):
if spaces.cached.get_spaces(nonstandard=False, dim=(3,)):

Check warning on line 338 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L338

Added line #L338 was not covered by tests
template_iterator_wf = init_template_iterator_wf(spaces=spaces)
workflow.connect([
(anat_fit_wf, template_iterator_wf, [
Expand All @@ -342,6 +344,34 @@
]),
]) # fmt:skip

# Thread MNI152NLin6Asym standard outputs to CIFTI subworkflow, skipping
# the iterator, which targets only output spaces.
# This can lead to duplication in the working directory if people actually
# want MNI152NLin6Asym outputs, but we'll live with it.
if config.workflow.cifti_output:
from smriprep.interfaces.templateflow import TemplateFlowSelect

Check warning on line 352 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L351-L352

Added lines #L351 - L352 were not covered by tests

ref = Reference(

Check warning on line 354 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L354

Added line #L354 was not covered by tests
"MNI152NLin6Asym",
{"res": 2 if config.workflow.cifti_output == "91k" else 1},
)

select_MNI6_xfm = pe.Node(

Check warning on line 359 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L359

Added line #L359 was not covered by tests
KeySelect(fields=["anat2std_xfm"], key=ref.fullname),
name="select_MNI6",
run_without_submitting=True,
)
select_MNI6_tpl = pe.Node(

Check warning on line 364 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L364

Added line #L364 was not covered by tests
TemplateFlowSelect(template=ref.fullname, resolution=ref.spec['res']),
name="select_MNI6_tpl",
)
workflow.connect([

Check warning on line 368 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L368

Added line #L368 was not covered by tests
(anat_fit_wf, select_MNI6_xfm, [
("outputnode.anat2std_xfm", "anat2std_xfm"),
("outputnode.template", "keys"),
]),
]) # fmt:skip

if config.workflow.anat_only:
return clean_datasinks(workflow)

Expand All @@ -362,7 +392,6 @@
f"{[e.method for e in fmap_estimators]}."
)

from niworkflows.interfaces.utility import KeySelect
from sdcflows import fieldmaps as fm
from sdcflows.workflows.base import init_fmap_preproc_wf

Expand Down Expand Up @@ -507,6 +536,12 @@
('outputnode.subjects_dir', 'inputnode.subjects_dir'),
('outputnode.subject_id', 'inputnode.subject_id'),
('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'),
('outputnode.white', 'inputnode.white'),
('outputnode.pial', 'inputnode.pial'),
('outputnode.midthickness', 'inputnode.midthickness'),
('outputnode.thickness', 'inputnode.thickness'),
('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
('outputnode.anat_ribbon', 'inputnode.anat_ribbon'),
]),
]) # fmt:skip
if fieldmap_id:
Expand All @@ -522,16 +557,27 @@
]) # fmt:skip

if config.workflow.level == "full":
workflow.connect([
(template_iterator_wf, bold_wf, [
("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"),
("outputnode.space", "inputnode.std_space"),
("outputnode.resolution", "inputnode.std_resolution"),
("outputnode.cohort", "inputnode.std_cohort"),
("outputnode.std_t1w", "inputnode.std_t1w"),
("outputnode.std_mask", "inputnode.std_mask"),
]),
]) # fmt:skip
if template_iterator_wf is not None:
workflow.connect([

Check warning on line 561 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L560-L561

Added lines #L560 - L561 were not covered by tests
(template_iterator_wf, bold_wf, [
("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"),
("outputnode.space", "inputnode.std_space"),
("outputnode.resolution", "inputnode.std_resolution"),
("outputnode.cohort", "inputnode.std_cohort"),
("outputnode.std_t1w", "inputnode.std_t1w"),
("outputnode.std_mask", "inputnode.std_mask"),
]),
]) # fmt:skip

# Thread MNI152NLin6Asym standard outputs to CIFTI subworkflow, skipping
# the iterator, which targets only output spaces.
# This can lead to duplication in the working directory if people actually
# want MNI152NLin6Asym outputs, but we'll live with it.
if config.workflow.cifti_output:
workflow.connect([

Check warning on line 577 in fmriprep/workflows/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/base.py#L576-L577

Added lines #L576 - L577 were not covered by tests
(select_MNI6_xfm, bold_wf, [("anat2std_xfm", "inputnode.anat2mni6_xfm")]),
(select_MNI6_tpl, bold_wf, [("brain_mask", "inputnode.mni6_mask")]),
]) # fmt:skip

return clean_datasinks(workflow)

Expand Down
148 changes: 100 additions & 48 deletions fmriprep/workflows/bold/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@
init_ds_registration_wf,
init_ds_volumes_wf,
init_func_derivatives_wf,
prepare_timing_parameters,
)
from .registration import init_bold_reg_wf, init_bold_t1_trans_wf
from .resampling import (
Expand Down Expand Up @@ -189,9 +190,16 @@
"t1w_mask",
"t1w_dseg",
"t1w_tpms",
# FreeSurfer outputs
"subjects_dir",
"subject_id",
"fsnative2t1w_xfm",
"white",
"midthickness",
"pial",
"thickness",
"sphere_reg_fsLR",
"anat_ribbon",
# Fieldmap registration
"fmap",
"fmap_ref",
Expand All @@ -206,6 +214,9 @@
"std_cohort",
"std_t1w",
"std_mask",
# MNI152NLin6Asym warp, for CIFTI use
"anat2mni6_xfm",
"mni6_mask",
],
),
name="inputnode",
Expand Down Expand Up @@ -389,7 +400,7 @@
(bold_anat_wf, ds_bold_t1_wf, [('outputnode.bold_file', 'inputnode.bold')]),
]) # fmt:skip

if spaces.get_spaces(nonstandard=False, dim=(3,)):
if spaces.cached.get_spaces(nonstandard=False, dim=(3,)):

Check warning on line 403 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L403

Added line #L403 was not covered by tests
# Missing:
# * Clipping BOLD after resampling
# * Resampling parcellations
Expand Down Expand Up @@ -462,6 +473,87 @@
(bold_anat_wf, bold_surf_wf, [("outputnode.bold_file", "inputnode.bold_t1w")]),
]) # fmt:skip

if config.workflow.cifti_output:
from .resampling import init_bold_fsLR_resampling_wf, init_bold_grayords_wf

Check warning on line 477 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L476-L477

Added lines #L476 - L477 were not covered by tests

bold_MNI6_wf = init_bold_volumetric_resample_wf(

Check warning on line 479 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L479

Added line #L479 was not covered by tests
metadata=all_metadata[0],
fieldmap_id=fieldmap_id if not multiecho else None,
omp_nthreads=omp_nthreads,
name='bold_MNI6_wf',
)

bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf(

Check warning on line 486 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L486

Added line #L486 was not covered by tests
estimate_goodvoxels=config.workflow.project_goodvoxels,
grayord_density=config.workflow.cifti_output,
omp_nthreads=omp_nthreads,
mem_gb=mem_gb["resampled"],
)

bold_grayords_wf = init_bold_grayords_wf(

Check warning on line 493 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L493

Added line #L493 was not covered by tests
grayord_density=config.workflow.cifti_output,
mem_gb=mem_gb["resampled"],
repetition_time=all_metadata[0]["RepetitionTime"],
)

ds_bold_cifti = pe.Node(

Check warning on line 499 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L499

Added line #L499 was not covered by tests
DerivativesDataSink(
base_directory=fmriprep_dir,
space='fsLR',
density=config.workflow.cifti_output,
suffix='bold',
compress=False,
TaskName=all_metadata[0].get('TaskName'),
**prepare_timing_parameters(all_metadata[0]),
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not necessarily related to this PR, but is this the only time we are outputing this metadata? seems like something we should only calculate once and then distribute to any bold outputs

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, we also calculate it in the volumetric workflow:

def init_ds_volumes_wf(
*,
bids_root: str,
output_dir: str,
multiecho: bool,
metadata: ty.List[dict],
name="ds_volumes_wf",
) -> pe.Workflow:
timing_parameters = prepare_timing_parameters(metadata)
workflow = pe.Workflow(name=name)
inputnode = pe.Node(
niu.IdentityInterface(
fields=[
'source_files',
'ref_file',
'bold', # Resampled into target space
'bold_mask', # boldref space
'bold_ref', # boldref space
't2star', # boldref space
# Anatomical
'boldref2anat_xfm',
# Template
'anat2std_xfm',
# Entities
'space',
'cohort',
'resolution',
]
),
name='inputnode',
)
raw_sources = pe.Node(niu.Function(function=_bids_relative), name='raw_sources')
raw_sources.inputs.bids_root = bids_root
boldref2target = pe.Node(niu.Merge(2), name='boldref2target')
# BOLD is pre-resampled
ds_bold = pe.Node(
DerivativesDataSink(
base_directory=output_dir,
desc='preproc',
compress=True,
SkullStripped=multiecho,
TaskName=metadata.get('TaskName'),
dismiss_entities=("echo",),
**timing_parameters,
),
name='ds_bold',
run_without_submitting=True,
mem_gb=DEFAULT_MEMORY_MIN_GB,
)

So it gets recalculated several times. We could refactor the function to have a cacheable core (can't @lru_cache with a dictionary argument), but I'm not sure saving <1s (estimate, not measured) at workflow build is really worth it.

),
name='ds_bold_cifti',
run_without_submitting=True,
)
ds_bold_cifti.inputs.source_file = bold_file

Check warning on line 512 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L512

Added line #L512 was not covered by tests

workflow.connect([

Check warning on line 514 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L514

Added line #L514 was not covered by tests
# Resample BOLD to MNI152NLin6Asym, may duplicate bold_std_wf above
(inputnode, bold_MNI6_wf, [
("mni6_mask", "inputnode.target_ref_file"),
("mni6_mask", "inputnode.target_mask"),
("anat2mni6_xfm", "inputnode.anat2std_xfm"),
("fmap_ref", "inputnode.fmap_ref"),
("fmap_coeff", "inputnode.fmap_coeff"),
("fmap_id", "inputnode.fmap_id"),
]),
(bold_fit_wf, bold_MNI6_wf, [
("outputnode.coreg_boldref", "inputnode.bold_ref_file"),
("outputnode.boldref2fmap_xfm", "inputnode.boldref2fmap_xfm"),
("outputnode.boldref2anat_xfm", "inputnode.boldref2anat_xfm"),
]),
(bold_native_wf, bold_MNI6_wf, [
("outputnode.bold_minimal", "inputnode.bold_file"),
("outputnode.motion_xfm", "inputnode.motion_xfm"),
]),
# Resample T1w-space BOLD to fsLR surfaces
(inputnode, bold_fsLR_resampling_wf, [
("white", "inputnode.white"),
("pial", "inputnode.pial"),
("midthickness", "inputnode.midthickness"),
("thickness", "inputnode.thickness"),
("sphere_reg_fsLR", "inputnode.sphere_reg_fsLR"),
("anat_ribbon", "inputnode.anat_ribbon"),
]),
(bold_anat_wf, bold_fsLR_resampling_wf, [
("outputnode.bold_file", "inputnode.bold_file"),
]),
(bold_MNI6_wf, bold_grayords_wf, [
("outputnode.bold_file", "inputnode.bold_std"),
]),
(bold_fsLR_resampling_wf, bold_grayords_wf, [
("outputnode.bold_fsLR", "inputnode.bold_fsLR"),
]),
(bold_grayords_wf, ds_bold_cifti, [
('outputnode.cifti_bold', 'in_file'),
(('outputnode.cifti_metadata', _read_json), 'meta_dict'),
]),
]) # fmt:skip

bold_confounds_wf = init_bold_confs_wf(
mem_gb=mem_gb["largemem"],
metadata=all_metadata[0],
Expand Down Expand Up @@ -657,7 +749,6 @@
spaces = config.workflow.spaces
fmriprep_dir = str(config.execution.fmriprep_dir)
freesurfer_spaces = spaces.get_fs_spaces()
project_goodvoxels = config.workflow.project_goodvoxels and config.workflow.cifti_output

ref_file = bold_file
wf_name = _get_wf_name(ref_file, "func_preproc")
Expand Down Expand Up @@ -745,52 +836,6 @@
# SURFACES ##################################################################################

# 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(
estimate_goodvoxels=project_goodvoxels,
grayord_density=config.workflow.cifti_output,
omp_nthreads=omp_nthreads,
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"],
)

# 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_fsLR_resampling_wf, func_derivatives_wf, [
("outputnode.goodvoxels_mask", "inputnode.goodvoxels_mask"),
]),
(bold_fsLR_resampling_wf, outputnode, [
("outputnode.weights_text", "weights_text"),
]),
(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(
Expand Down Expand Up @@ -950,3 +995,10 @@
"""Return the image orientation as a string"""
img = nb.load(imgf)
return "".join(nb.aff2axcodes(img.affine))


def _read_json(in_file):
from json import loads
from pathlib import Path

Check warning on line 1002 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L1001-L1002

Added lines #L1001 - L1002 were not covered by tests

return loads(Path(in_file).read_text())

Check warning on line 1004 in fmriprep/workflows/bold/base.py

View check run for this annotation

Codecov / codecov/patch

fmriprep/workflows/bold/base.py#L1004

Added line #L1004 was not covered by tests
Loading