diff --git a/mriqc/cli/parser.py b/mriqc/cli/parser.py index 349dcfab6..907b7e43f 100644 --- a/mriqc/cli/parser.py +++ b/mriqc/cli/parser.py @@ -538,7 +538,7 @@ def parse_args(args=None, namespace=None): config.execution.participant_label, session_id=config.execution.session_id, task=config.execution.task_id, - group_echos=False, + group_echos=True, bids_filters=config.execution.bids_filters, queries={mod: DEFAULT_BIDS_QUERIES[mod] for mod in lc_modalities} ) @@ -575,7 +575,7 @@ def parse_args(args=None, namespace=None): # Estimate the biggest file size / leave 1GB if some file does not exist (datalad) with suppress(FileNotFoundError): config.workflow.biggest_file_gb = _get_biggest_file_size_gb( - [i for sublist in config.workflow.inputs.values() for i in sublist] + config.workflow.inputs.values() ) # set specifics for alternative populations @@ -593,11 +593,14 @@ def parse_args(args=None, namespace=None): def _get_biggest_file_size_gb(files): + """Identify the largest file size (allows multi-echo groups).""" + import os - max_size = 0 + sizes = [] for file in files: - size = os.path.getsize(file) / (1024**3) - if size > max_size: - max_size = size - return max_size + if isinstance(file, (list, tuple)): + sizes.append(_get_biggest_file_size_gb(file)) + else: + sizes.append(os.path.getsize(file)) + return max(sizes) / (1024**3) diff --git a/mriqc/data/bootstrap-func.yml b/mriqc/data/bootstrap-func.yml index ae07531f3..659c3b5e8 100644 --- a/mriqc/data/bootstrap-func.yml +++ b/mriqc/data/bootstrap-func.yml @@ -26,18 +26,14 @@ ########################################################################### packagename: mriqc -title: '{filename} :: Functional MRIQC report' +title: "{filename} :: MRIQC's BOLD fMRI report" sections: - name: Summary reportlets: - bids: {datatype: figures, desc: summary, extension: [.html]} -- name: Basic visual report +- name: Basic echo-wise reports + ordering: echo reportlets: - - bids: {datatype: figures, desc: zoomed} - caption: This panel shows a mosaic of the brain. This mosaic is the most suitable to - screen head-motion intensity inhomogeneities, global/local noise, signal leakage - (for example, from the eyeballs and across the phase-encoding axis), etc. - subtitle: Voxel-wise average of BOLD time-series, zoomed-in covering just the brain - bids: {datatype: figures, desc: stdev} subtitle: Standard deviation of signal through time caption: The voxel-wise standard deviation of the signal (variability along time). @@ -45,43 +41,12 @@ sections: subtitle: Carpetplot and nuisance signals caption: The so-called «carpetplot» may assist in assessing head-motion derived artifacts and respiation effects. -- name: Extended visual report +- name: Extended reports shared across echos reportlets: - - bids: {datatype: figures, desc: background} - caption: This panel shows a mosaic enhancing the background around the head. - Artifacts usually unveil themselves in the air surrounding the head, where no signal - sources are present. - subtitle: View of the background of the voxel-wise average of the BOLD timeseries - - bids: {datatype: figures, desc: mean} - subtitle: Average signal through time - caption: The average signal calculated across the last axis (time). - - bids: {datatype: figures, desc: airmask} - caption: The hat-mask calculated internally by MRIQC. Some metrics will use this - mask, for instance, to find out artifacts and estimate the spread of gaussian noise - added to the signal. This mask leaves out the air around the face to avoid measuring - noise sourcing from the eyeballs and their movement. - subtitle: '«Hat»-mask' - - bids: {datatype: figures, desc: noisefit} - caption: The noise fit internally estimated by MRIQC to calculate the QI1 index - proposed by Mortamet et al. (2009). - subtitle: Distribution of the noise within the hat mask - style: - max-width: 450px - - bids: {datatype: figures, desc: artifacts} - caption: Mask of artifactual intensities identified within the hat-mask. - subtitle: Artifactual intensities on the background - bids: {datatype: figures, desc: brainmask} caption: Brain mask as internally extracted by MRIQC. Defects on the brainmask could indicate problematic aspects of the image quality-wise. subtitle: Brain extraction performance - - bids: {datatype: figures, desc: head} - caption: A mask of the head calculated internally by MRIQC. - subtitle: Head mask - - bids: {datatype: figures, desc: segmentation} - caption: Brain tissue segmentation, as internally extracted by MRIQC. - Defects on this segmentation, as well as noisy tissue labels could - indicate problematic aspects of the image quality-wise. - subtitle: Brain tissue segmentation - bids: {datatype: figures, desc: norm} caption: This panel shows a quick-and-dirty nonlinear registration into the MNI152NLin2009cAsym template accessed with @@ -89,6 +54,22 @@ sections: subtitle: Spatial normalization of the anatomical image static: false +- name: Extended echo-wise reports + ordering: echo + reportlets: + - bids: {datatype: figures, desc: background} + caption: This panel shows a mosaic enhancing the background around the head. + Artifacts usually unveil themselves in the air surrounding the head, where no signal + sources are present. + subtitle: View of the background of the voxel-wise average of the BOLD timeseries + - bids: {datatype: figures, desc: mean} + subtitle: Average signal through time + caption: The average signal calculated across the last axis (time). + - bids: {datatype: figures, desc: zoomed} + caption: This panel shows a mosaic of the brain. This mosaic is the most suitable to + screen head-motion intensity inhomogeneities, global/local noise, signal leakage + (for example, from the eyeballs and across the phase-encoding axis), etc. + subtitle: Voxel-wise average of BOLD time-series, zoomed-in covering just the brain - name: About nested: true diff --git a/mriqc/interfaces/bids.py b/mriqc/interfaces/bids.py index eb48122ae..d4ae9de42 100644 --- a/mriqc/interfaces/bids.py +++ b/mriqc/interfaces/bids.py @@ -49,6 +49,7 @@ class IQMFileSinkInputSpec(DynamicTraitedSpec, BaseInterfaceInputSpec): rec_id = traits.Either(None, Str, usedefault=True) run_id = traits.Either(None, traits.Int, usedefault=True) dataset = Str(desc="dataset identifier") + dismiss_entities = traits.List(["part"], usedefault=True) metadata = traits.Dict() provenance = traits.Dict() @@ -111,6 +112,17 @@ def _gen_outfile(self): break in_file = str(path.relative_to(bids_root)) + if ( + isdefined(self.inputs.dismiss_entities) + and (dismiss := self.inputs.dismiss_entities) + ): + for entity in dismiss: + bids_chunks = [ + chunk for chunk in path.name.split("_") + if not chunk.startswith(f"{entity}-") + ] + path = path.parent / "_".join(bids_chunks) + # Build path and ensure directory exists bids_path = out_dir / in_file.replace("".join(Path(in_file).suffixes), ".json") bids_path.parent.mkdir(parents=True, exist_ok=True) diff --git a/mriqc/reports/individual.py b/mriqc/reports/individual.py index 1e2350dde..8619f1ec2 100644 --- a/mriqc/reports/individual.py +++ b/mriqc/reports/individual.py @@ -47,11 +47,13 @@ def _single_report(in_file): from mriqc import config # Ensure it's a Path - in_file = Path(in_file) + in_file = Path(in_file if not isinstance(in_file, list) else in_file[0]) # Extract BIDS entities entities = config.execution.layout.get_file(in_file).get_entities() entities.pop("extension", None) + entities.pop("echo", None) + entities.pop("part", None) report_type = entities.pop("datatype", None) # Read output file: diff --git a/mriqc/utils/bids.py b/mriqc/utils/bids.py index 4ab142655..40bbfed8d 100644 --- a/mriqc/utils/bids.py +++ b/mriqc/utils/bids.py @@ -21,6 +21,8 @@ # https://www.nipreps.org/community/licensing/ # """PyBIDS tooling.""" +from __future__ import annotations + import json import os from pathlib import Path @@ -90,3 +92,99 @@ def write_derivative_description(bids_dir, deriv_dir): desc["License"] = orig_desc["License"] Path.write_text(deriv_dir / "dataset_description.json", json.dumps(desc, indent=4)) + + +def derive_bids_fname( + orig_path: str | Path, + entity: str | None = None, + newsuffix: str | None = None, + newpath: str | Path | None = None, + newext: str | None = None, + position: int = -1, + absolute: bool = True, +) -> Path | str: + """ + Derive a new file name from a BIDS-formatted path. + + Parameters + ---------- + orig_path : :obj:`str` or :obj:`os.pathlike` + A filename (may or may not include path). + entity : :obj:`str`, optional + A new BIDS-like key-value pair. + newsuffix : :obj:`str`, optional + Replace the BIDS suffix. + newpath : :obj:`str` or :obj:`os.pathlike`, optional + Path to replace the path of the input orig_path. + newext : :obj:`str`, optional + Replace the extension of the file. + position : :obj:`int`, optional + Position to insert the entity in the filename. + absolute : :obj:`bool`, optional + If True (default), returns the absolute path of the modified filename. + + Returns + ------- + Absolute path of the modified filename + + Examples + -------- + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-preproc', + ... absolute=False, + ... ) + PosixPath('sub-001/ses-01/anat/sub-001_ses-01_desc-preproc_T1w.nii.gz') + + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-brain', + ... newsuffix='mask', + ... newext=".nii", + ... absolute=False, + ... ) # doctest: +ELLIPSIS + PosixPath('sub-001/ses-01/anat/sub-001_ses-01_desc-brain_mask.nii') + + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-brain', + ... newsuffix='mask', + ... newext=".nii", + ... newpath="/output/node", + ... absolute=True, + ... ) # doctest: +ELLIPSIS + PosixPath('/output/node/sub-001_ses-01_desc-brain_mask.nii') + + >>> derive_bids_fname( + ... 'sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz', + ... entity='desc-brain', + ... newsuffix='mask', + ... newext=".nii", + ... newpath=".", + ... absolute=False, + ... ) # doctest: +ELLIPSIS + PosixPath('sub-001_ses-01_desc-brain_mask.nii') + + """ + + orig_path = Path(orig_path) + newpath = orig_path.parent if newpath is None else Path(newpath) + + ext = "".join(orig_path.suffixes) + newext = newext if newext is not None else ext + orig_stem = orig_path.name.replace(ext, "") + + suffix = orig_stem.rsplit("_", maxsplit=1)[-1].strip("_") + newsuffix = newsuffix.strip("_") if newsuffix is not None else suffix + + orig_stem = orig_stem.replace(suffix, "").strip("_") + bidts = [bit for bit in orig_stem.split("_") if bit] + if entity: + if position == -1: + bidts.append(entity) + else: + bidts.insert(position, entity.strip("_")) + + retval = (newpath / f"{'_'.join(bidts)}_{newsuffix}.{newext.strip('.')}") + + return retval.absolute() if absolute else retval diff --git a/mriqc/workflows/anatomical/base.py b/mriqc/workflows/anatomical/base.py index 400dee068..0751ee327 100644 --- a/mriqc/workflows/anatomical/base.py +++ b/mriqc/workflows/anatomical/base.py @@ -86,6 +86,7 @@ def anat_qc_workflow(name="anatMRIQC"): wf = anat_qc_workflow() """ + from mriqc.workflows.shared import synthstrip_wf dataset = config.workflow.inputs.get("t1w", []) + config.workflow.inputs.get("t2w", []) @@ -682,101 +683,6 @@ def airmsk_wf(name="AirMaskWorkflow"): return workflow -def synthstrip_wf(name="synthstrip_wf", omp_nthreads=None): - """Create a brain-extraction workflow using SynthStrip.""" - from nipype.interfaces.ants import N4BiasFieldCorrection - from niworkflows.interfaces.nibabel import IntensityClip, ApplyMask - from mriqc.interfaces.synthstrip import SynthStrip - - inputnode = pe.Node(niu.IdentityInterface(fields=["in_files"]), name="inputnode") - outputnode = pe.Node( - niu.IdentityInterface(fields=["out_corrected", "out_brain", "bias_image", "out_mask"]), - name="outputnode", - ) - - # truncate target intensity for N4 correction - pre_clip = pe.Node(IntensityClip(p_min=10, p_max=99.9), name="pre_clip") - - pre_n4 = pe.Node( - N4BiasFieldCorrection( - dimension=3, - num_threads=omp_nthreads, - rescale_intensities=True, - copy_header=True, - ), - name="pre_n4", - ) - - post_n4 = pe.Node( - N4BiasFieldCorrection( - dimension=3, - save_bias=True, - num_threads=omp_nthreads, - n_iterations=[50] * 4, - copy_header=True, - ), - name="post_n4", - ) - - synthstrip = pe.Node( - SynthStrip(num_threads=omp_nthreads), - name="synthstrip", - num_threads=omp_nthreads, - ) - - final_masked = pe.Node(ApplyMask(), name="final_masked") - final_inu = pe.Node(niu.Function(function=_apply_bias_correction), name="final_inu") - - workflow = pe.Workflow(name=name) - # fmt: off - workflow.connect([ - (inputnode, final_inu, [("in_files", "in_file")]), - (inputnode, pre_clip, [("in_files", "in_file")]), - (pre_clip, pre_n4, [("out_file", "input_image")]), - (pre_n4, synthstrip, [("output_image", "in_file")]), - (synthstrip, post_n4, [("out_mask", "weight_image")]), - (synthstrip, final_masked, [("out_mask", "in_mask")]), - (pre_clip, post_n4, [("out_file", "input_image")]), - (post_n4, final_inu, [("bias_image", "bias_image")]), - (post_n4, final_masked, [("output_image", "in_file")]), - (final_masked, outputnode, [("out_file", "out_brain")]), - (post_n4, outputnode, [("bias_image", "bias_image")]), - (synthstrip, outputnode, [("out_mask", "out_mask")]), - (post_n4, outputnode, [("output_image", "out_corrected")]), - ]) - # fmt: on - return workflow - - -def _apply_bias_correction(in_file, bias_image, out_file=None): - import os.path as op - - import numpy as np - import nibabel as nb - - img = nb.load(in_file) - data = np.clip( - img.get_fdata() * nb.load(bias_image).get_fdata(), - a_min=0, - a_max=None, - ) - out_img = img.__class__( - data.astype(img.get_data_dtype()), - img.affine, - img.header, - ) - - if out_file is None: - fname, ext = op.splitext(op.basename(in_file)) - if ext == ".gz": - fname, ext2 = op.splitext(fname) - ext = ext2 + ext - out_file = op.abspath(f"{fname}_inu{ext}") - - out_img.to_filename(out_file) - return out_file - - def _binarize(in_file, threshold=0.5, out_file=None): import os.path as op diff --git a/mriqc/workflows/diffusion/base.py b/mriqc/workflows/diffusion/base.py index 470f02fa7..a81d08520 100644 --- a/mriqc/workflows/diffusion/base.py +++ b/mriqc/workflows/diffusion/base.py @@ -76,6 +76,7 @@ def dmri_qc_workflow(name="dwiMRIQC"): ReadDWIMetadata, WeightedStat, ) + from mriqc.workflows.shared import synthstrip_wf as dmri_bmsk_workflow from mriqc.messages import BUILDING_WORKFLOW workflow = pe.Workflow(name=name) @@ -349,70 +350,6 @@ def compute_iqms(name="ComputeIQMs"): return workflow -def dmri_bmsk_workflow(name="dmri_brainmask", omp_nthreads=None): - """ - Compute a brain mask for the input :abbr:`dMRI (diffusion MRI)` dataset. - - .. workflow:: - - from mriqc.workflows.diffusion.base import dmri_bmsk_workflow - from mriqc.testing import mock_config - with mock_config(): - wf = dmri_bmsk_workflow() - - - """ - - from nipype.interfaces.ants import N4BiasFieldCorrection - from niworkflows.interfaces.nibabel import ApplyMask - from mriqc.interfaces.synthstrip import SynthStrip - from mriqc.workflows.anatomical.base import _apply_bias_correction - - inputnode = pe.Node(niu.IdentityInterface(fields=["in_file"]), name="inputnode") - outputnode = pe.Node( - niu.IdentityInterface(fields=["out_corrected", "out_brain", "bias_image", "out_mask"]), - name="outputnode", - ) - - post_n4 = pe.Node( - N4BiasFieldCorrection( - dimension=3, - save_bias=True, - num_threads=omp_nthreads, - n_iterations=[50] * 4, - copy_header=True, - ), - name="post_n4", - ) - - synthstrip = pe.Node( - SynthStrip(num_threads=omp_nthreads), - name="synthstrip", - num_threads=omp_nthreads, - ) - - final_masked = pe.Node(ApplyMask(), name="final_masked") - final_inu = pe.Node(niu.Function(function=_apply_bias_correction), name="final_inu") - - workflow = pe.Workflow(name=name) - # fmt: off - workflow.connect([ - (inputnode, final_inu, [("in_file", "in_file")]), - (inputnode, synthstrip, [("in_file", "in_file")]), - (inputnode, post_n4, [("in_file", "input_image")]), - (synthstrip, post_n4, [("out_mask", "weight_image")]), - (synthstrip, final_masked, [("out_mask", "in_mask")]), - (post_n4, final_inu, [("bias_image", "bias_image")]), - (post_n4, final_masked, [("output_image", "in_file")]), - (final_masked, outputnode, [("out_file", "out_brain")]), - (post_n4, outputnode, [("bias_image", "bias_image")]), - (synthstrip, outputnode, [("out_mask", "out_mask")]), - (post_n4, outputnode, [("output_image", "out_corrected")]), - ]) - # fmt: on - return workflow - - def init_dmriref_wf( in_file=None, name="init_dmriref_wf", diff --git a/mriqc/workflows/functional/base.py b/mriqc/workflows/functional/base.py index 1ede04f61..3a9612aba 100644 --- a/mriqc/workflows/functional/base.py +++ b/mriqc/workflows/functional/base.py @@ -45,6 +45,8 @@ from mriqc import config from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe +from niworkflows.utils.connections import pop_file as _pop + from mriqc.interfaces.datalad import DataladIdentityInterface from mriqc.workflows.functional.output import init_func_report_wf @@ -76,7 +78,7 @@ def fmri_qc_workflow(name="funcMRIQC"): message = BUILDING_WORKFLOW.format( modality="functional", detail=( - f"for {len(dataset)} NIfTI files." + f"for {len(dataset)} BOLD runs." if len(dataset) > 2 else f"({' and '.join(('<%s>' % v for v in dataset))})." ), @@ -88,9 +90,10 @@ def fmri_qc_workflow(name="funcMRIQC"): inputnode = pe.Node(niu.IdentityInterface(fields=["in_file"]), name="inputnode") inputnode.iterables = [("in_file", dataset)] - datalad_get = pe.Node( + datalad_get = pe.MapNode( DataladIdentityInterface(fields=["in_file"], dataset_path=config.execution.bids_dir), name="datalad_get", + iterfield=["in_file"], ) outputnode = pe.Node( @@ -100,30 +103,40 @@ def fmri_qc_workflow(name="funcMRIQC"): non_steady_state_detector = pe.Node(NonSteadyStateDetector(), name="non_steady_state_detector") - sanitize = pe.Node(SanitizeImage(), name="sanitize", mem_gb=mem_gb * 4.0) - sanitize.inputs.max_32bit = config.execution.float32 + sanitize = pe.MapNode( + SanitizeImage(max_32bit=config.execution.float32), + name="sanitize", + mem_gb=mem_gb * 4.0, + iterfield=["in_file"], + ) # Workflow -------------------------------------------------------- # 1. HMC: head motion correct - hmcwf = hmc() + hmcwf = hmc(omp_nthreads=config.nipype.omp_nthreads) # Set HMC settings hmcwf.inputs.inputnode.fd_radius = config.workflow.fd_radius # 2. Compute mean fmri - mean = pe.Node( + mean = pe.MapNode( TStat(options="-mean", outputtype="NIFTI_GZ"), name="mean", mem_gb=mem_gb * 1.5, + iterfield=["in_file"], + ) + + # Compute TSNR using nipype implementation + tsnr = pe.MapNode( + TSNR(), + name="compute_tsnr", + mem_gb=mem_gb * 2.5, + iterfield=["in_file"], ) # EPI to MNI registration ema = epi_mni_align() - # Compute TSNR using nipype implementation - tsnr = pe.Node(TSNR(), name="compute_tsnr", mem_gb=mem_gb * 2.5) - # 7. Compute IQMs iqmswf = compute_iqms() # Reports @@ -133,22 +146,25 @@ def fmri_qc_workflow(name="funcMRIQC"): workflow.connect([ (inputnode, datalad_get, [("in_file", "in_file")]), - (inputnode, func_report_wf, [ - ("in_file", "inputnode.name_source"), - ]), - (datalad_get, iqmswf, [("in_file", "inputnode.in_file")]), (datalad_get, sanitize, [("in_file", "in_file")]), - (datalad_get, non_steady_state_detector, [("in_file", "in_file")]), + (datalad_get, non_steady_state_detector, [(("in_file", select_echo2), "in_file")]), (non_steady_state_detector, sanitize, [("n_volumes_to_discard", "n_volumes_to_discard")]), (sanitize, hmcwf, [("out_file", "inputnode.in_file")]), (hmcwf, mean, [("outputnode.out_file", "in_file")]), (hmcwf, tsnr, [("outputnode.out_file", "in_file")]), - (mean, ema, [("out_file", "inputnode.epi_mean")]), + (mean, ema, [(("out_file", select_echo2), "inputnode.epi_mean")]), + # Feed IQMs computation + (datalad_get, iqmswf, [("in_file", "inputnode.in_file")]), (sanitize, iqmswf, [("out_file", "inputnode.in_ras")]), (mean, iqmswf, [("out_file", "inputnode.epi_mean")]), (hmcwf, iqmswf, [("outputnode.out_file", "inputnode.hmc_epi"), ("outputnode.out_fd", "inputnode.hmc_fd")]), (tsnr, iqmswf, [("tsnr_file", "inputnode.in_tsnr")]), + (non_steady_state_detector, iqmswf, [("n_volumes_to_discard", "inputnode.exclude_index")]), + # Feed reportlet generation + (inputnode, func_report_wf, [ + ("in_file", "inputnode.name_source"), + ]), (sanitize, func_report_wf, [("out_file", "inputnode.in_ras")]), (mean, func_report_wf, [("out_file", "inputnode.epi_mean")]), (tsnr, func_report_wf, [("stddev_file", "inputnode.in_stddev")]), @@ -160,7 +176,6 @@ def fmri_qc_workflow(name="funcMRIQC"): ("outputnode.epi_parc", "inputnode.epi_parc"), ("outputnode.report", "inputnode.mni_report"), ]), - (non_steady_state_detector, iqmswf, [("n_volumes_to_discard", "inputnode.exclude_index")]), (iqmswf, func_report_wf, [ ("outputnode.out_file", "inputnode.in_iqms"), ("outputnode.out_dvars", "inputnode.in_dvars"), @@ -183,13 +198,15 @@ def fmri_qc_workflow(name="funcMRIQC"): # population specific changes to brain masking if config.workflow.species == "human": - skullstrip_epi = fmri_bmsk_workflow() + from mriqc.workflows.shared import synthstrip_wf as fmri_bmsk_workflow + + skullstrip_epi = fmri_bmsk_workflow(omp_nthreads=config.nipype.omp_nthreads) # fmt: off workflow.connect([ - (mean, skullstrip_epi, [("out_file", "inputnode.in_file")]), - (skullstrip_epi, ema, [("outputnode.out_file", "inputnode.epi_mask")]), - (skullstrip_epi, iqmswf, [("outputnode.out_file", "inputnode.brainmask")]), - (skullstrip_epi, func_report_wf, [("outputnode.out_file", "inputnode.brainmask")]), + (mean, skullstrip_epi, [(("out_file", select_echo2), "inputnode.in_files")]), + (skullstrip_epi, ema, [("outputnode.out_mask", "inputnode.epi_mask")]), + (skullstrip_epi, iqmswf, [("outputnode.out_mask", "inputnode.brainmask")]), + (skullstrip_epi, func_report_wf, [("outputnode.out_mask", "inputnode.brainmask")]), ]) # fmt: on else: @@ -216,11 +233,11 @@ def fmri_qc_workflow(name="funcMRIQC"): if not config.execution.no_sub: from mriqc.interfaces.webapi import UploadIQMs - upldwf = pe.Node(UploadIQMs( + upldwf = pe.MapNode(UploadIQMs( endpoint=config.execution.webapi_url, auth_token=config.execution.webapi_token, strict=config.execution.upload_strict, - ), name="UploadMetrics") + ), name="UploadMetrics", iterfield=["in_iqms"]) # fmt: off workflow.connect([ @@ -290,32 +307,40 @@ def compute_iqms(name="ComputeIQMs"): inputnode.inputs.fd_thres = config.workflow.fd_thres # Compute DVARS - dvnode = pe.Node( + dvnode = pe.MapNode( ComputeDVARS(save_plot=False, save_all=True), name="ComputeDVARS", mem_gb=mem_gb * 3, + iterfield=["in_file"], ) # AFNI quality measures - fwhm_interface = get_fwhmx() - fwhm = pe.Node(fwhm_interface, name="smoothness") - # fwhm.inputs.acf = True # add when AFNI >= 16 - outliers = pe.Node( + fwhm = pe.MapNode(get_fwhmx(), name="smoothness", iterfield=["in_file"]) + fwhm.inputs.acf = True # Only AFNI >= 16 + + outliers = pe.MapNode( OutlierCount(fraction=True, out_file="outliers.out"), name="outliers", mem_gb=mem_gb * 2.5, + iterfield=["in_file"], ) - quality = pe.Node( + quality = pe.MapNode( QualityIndex(automask=True), out_file="quality.out", name="quality", mem_gb=mem_gb * 3, + iterfield=["in_file"], ) - gcor = pe.Node(GCOR(), name="gcor", mem_gb=mem_gb * 2) + gcor = pe.MapNode(GCOR(), name="gcor", mem_gb=mem_gb * 2, iterfield=["in_file"]) - measures = pe.Node(FunctionalQC(), name="measures", mem_gb=mem_gb * 3) + measures = pe.MapNode( + FunctionalQC(), + name="measures", + mem_gb=mem_gb * 3, + iterfield=["in_epi", "in_hmc", "in_tsnr", "in_dvars", "in_fwhm"], + ) # fmt: off workflow.connect([ @@ -342,18 +367,19 @@ def compute_iqms(name="ComputeIQMs"): # fmt: on # Add metadata - meta = pe.Node(ReadSidecarJSON( + meta = pe.MapNode(ReadSidecarJSON( index_db=config.execution.bids_database_dir - ), name="metadata") + ), name="metadata", iterfield=["in_file"]) - addprov = pe.Node( + addprov = pe.MapNode( AddProvenance(modality="bold"), name="provenance", run_without_submitting=True, + iterfield=["in_file"], ) # Save to JSON file - datasink = pe.Node( + datasink = pe.MapNode( IQMFileSink( modality="bold", out_dir=str(config.execution.output_dir), @@ -361,6 +387,7 @@ def compute_iqms(name="ComputeIQMs"): ), name="datasink", run_without_submitting=True, + iterfield=["in_file", "root", "metadata", "provenance"], ) # fmt: off @@ -369,12 +396,12 @@ def compute_iqms(name="ComputeIQMs"): ("exclude_index", "dummy_trs")]), (inputnode, meta, [("in_file", "in_file")]), (inputnode, addprov, [("in_file", "in_file")]), - (meta, datasink, [("subject", "subject_id"), - ("session", "session_id"), - ("task", "task_id"), - ("acquisition", "acq_id"), - ("reconstruction", "rec_id"), - ("run", "run_id"), + (meta, datasink, [(("subject", _pop), "subject_id"), + (("session", _pop), "session_id"), + (("task", _pop), "task_id"), + (("acquisition", _pop), "acq_id"), + (("reconstruction", _pop), "rec_id"), + (("run", _pop), "run_id"), ("out_dict", "metadata")]), (addprov, datasink, [("out_prov", "provenance")]), (outliers, datasink, [(("out_file", _parse_tout), "aor")]), @@ -390,13 +417,14 @@ def compute_iqms(name="ComputeIQMs"): if config.workflow.fft_spikes_detector: from mriqc.workflows.utils import slice_wise_fft - spikes_fft = pe.Node( + spikes_fft = pe.MapNode( niu.Function( input_names=["in_file"], output_names=["n_spikes", "out_spikes", "out_fft"], function=slice_wise_fft, ), name="SpikesFinderFFT", + iterfield=["in_file"], ) # fmt: off @@ -441,7 +469,7 @@ def fmri_bmsk_workflow(name="fMRIBrainMask"): return workflow -def hmc(name="fMRI_HMC"): +def hmc(name="fMRI_HMC", omp_nthreads=None): """ Create a :abbr:`HMC (head motion correction)` workflow for fMRI. @@ -468,9 +496,9 @@ def hmc(name="fMRI_HMC"): outputnode = pe.Node(niu.IdentityInterface(fields=["out_file", "out_fd"]), name="outputnode") # calculate hmc parameters - hmc = pe.Node( + estimate_hm = pe.Node( Volreg(args="-Fourier -twopass", zpad=4, outputtype="NIFTI_GZ"), - name="motion_correct", + name="estimate_hm", mem_gb=mem_gb * 2.5, ) @@ -480,49 +508,73 @@ def hmc(name="fMRI_HMC"): name="ComputeFD", ) + # Apply transforms to other echos + apply_hmc = pe.MapNode( + niu.Function(function=_apply_transforms, input_names=["in_file", "in_xfm"]), + name="apply_hmc", + iterfield=["in_file"], + # NiTransforms is a memory hog, so ensure only one process is running at a time + num_threads=config.environment.cpu_count, + ) + # fmt: off workflow.connect([ (inputnode, fdnode, [("fd_radius", "radius")]), - (hmc, outputnode, [("out_file", "out_file")]), - (hmc, fdnode, [("oned_file", "in_file")]), + (estimate_hm, apply_hmc, [("oned_matrix_save", "in_xfm")]), + (apply_hmc, outputnode, [("out", "out_file")]), + (estimate_hm, fdnode, [("oned_file", "in_file")]), (fdnode, outputnode, [("out_file", "out_fd")]), ]) # fmt: on - # despiking, and deoblique - - deoblique_node = pe.Node(Refit(deoblique=True), name="deoblique") - - despike_node = pe.Node(Despike(outputtype="NIFTI_GZ"), name="despike") + if not (config.workflow.despike or config.workflow.deoblique): + # fmt: off + workflow.connect([ + (inputnode, estimate_hm, [(("in_file", _pop), "in_file")]), + (inputnode, apply_hmc, [("in_file", "in_file")]), + ]) + # fmt: on + return workflow + # despiking, and deoblique + deoblique_node = pe.MapNode( + Refit(deoblique=True), + name="deoblique", + iterfield=["in_file"], + ) + despike_node = pe.MapNode( + Despike(outputtype="NIFTI_GZ"), + name="despike", + iterfield=["in_file"], + ) if config.workflow.despike and config.workflow.deoblique: # fmt: off workflow.connect([ (inputnode, despike_node, [("in_file", "in_file")]), (despike_node, deoblique_node, [("out_file", "in_file")]), - (deoblique_node, hmc, [("out_file", "in_file")]), + (deoblique_node, estimate_hm, [(("out_file", _pop), "in_file")]), + (deoblique_node, apply_hmc, [("out_file", "in_file")]), ]) # fmt: on elif config.workflow.despike: # fmt: off workflow.connect([ (inputnode, despike_node, [("in_file", "in_file")]), - (despike_node, hmc, [("out_file", "in_file")]), + (despike_node, estimate_hm, [(("out_file", _pop), "in_file")]), + (despike_node, apply_hmc, [("out_file", "in_file")]), ]) # fmt: on elif config.workflow.deoblique: # fmt: off workflow.connect([ (inputnode, deoblique_node, [("in_file", "in_file")]), - (deoblique_node, hmc, [("out_file", "in_file")]), + (deoblique_node, estimate_hm, [(("out_file", _pop), "in_file")]), + (deoblique_node, apply_hmc, [("out_file", "in_file")]), ]) # fmt: on else: - # fmt: off - workflow.connect([ - (inputnode, hmc, [("in_file", "in_file")]), - ]) - # fmt: on + raise NotImplementedError + return workflow @@ -665,13 +717,10 @@ def epi_mni_align(name="SpatialNormalization"): return workflow -def _mean(inlist): - import numpy as np - - return np.mean(inlist) - - def _parse_tqual(in_file): + if isinstance(in_file, (list, tuple)): + return [_parse_tqual(f) for f in in_file] + import numpy as np with open(in_file, "r") as fin: @@ -680,7 +729,48 @@ def _parse_tqual(in_file): def _parse_tout(in_file): + if isinstance(in_file, (list, tuple)): + return [_parse_tout(f) for f in in_file] + import numpy as np data = np.loadtxt(in_file) # pylint: disable=no-member return data.mean() + + +def _apply_transforms(in_file, in_xfm): + from pathlib import Path + from nitransforms.linear import load + from mriqc.utils.bids import derive_bids_fname + + realigned = load(in_xfm, fmt="afni", reference=in_file, moving=in_file).apply(in_file) + out_file = derive_bids_fname( + in_file, + entity="desc-realigned", + newpath=Path.cwd(), + absolute=True, + ) + + realigned.to_filename(out_file) + return str(out_file) + + +def select_echo2(in_files): + """ + Select the first file from a list of filenames. + + Used to grab the first echo's file when processing + multi-echo data through workflows that only accept + a single file. + + Examples + -------- + >>> select_file('some/file.nii.gz') + 'some/file.nii.gz' + >>> select_file(['some/file1.nii.gz', 'some/file2.nii.gz']) + 'some/file1.nii.gz' + + """ + if isinstance(in_files, (list, tuple)): + return in_files[min(1, len(in_files) - 1)] + return in_files diff --git a/mriqc/workflows/functional/output.py b/mriqc/workflows/functional/output.py index fc6a421ce..f620cb2a1 100644 --- a/mriqc/workflows/functional/output.py +++ b/mriqc/workflows/functional/output.py @@ -81,7 +81,7 @@ def init_func_report_wf(name="func_report_wf"): # Set FD threshold inputnode.inputs.fd_thres = config.workflow.fd_thres - spmask = pe.Node( + spmask = pe.MapNode( niu.Function( input_names=["in_file", "in_mask"], output_names=["out_file", "out_plot"], @@ -89,12 +89,14 @@ def init_func_report_wf(name="func_report_wf"): ), name="SpikesMask", mem_gb=mem_gb * 3.5, + iterfield=["in_file"], ) - spikes_bg = pe.Node( + spikes_bg = pe.MapNode( Spikes(no_zscore=True, detrend=False), name="SpikesFinderBgMask", mem_gb=mem_gb * 2.5, + iterfield=["in_file", "in_mask"], ) # Generate crown mask @@ -103,7 +105,12 @@ def init_func_report_wf(name="func_report_wf"): subtract_mask = pe.Node(BinarySubtraction(), name="subtract_mask") parcels = pe.Node(niu.Function(function=_carpet_parcellation), name="parcels") - bigplot = pe.Node(FMRISummary(), name="BigPlot", mem_gb=mem_gb * 3.5) + bigplot = pe.MapNode( + FMRISummary(), + name="BigPlot", + mem_gb=mem_gb * 3.5, + iterfield=["in_func", "dvars", "outliers", "in_spikes_bg"], + ) # fmt: off workflow.connect([ @@ -118,58 +125,66 @@ def init_func_report_wf(name="func_report_wf"): (inputnode, parcels, [("epi_parc", "segmentation")]), (inputnode, dilated_mask, [("brainmask", "in_mask")]), (inputnode, subtract_mask, [("brainmask", "in_subtract")]), + (spmask, spikes_bg, [("out_file", "in_mask")]), (dilated_mask, subtract_mask, [("out_mask", "in_base")]), (subtract_mask, parcels, [("out_mask", "crown_mask")]), (parcels, bigplot, [("out", "in_segm")]), (spikes_bg, bigplot, [("out_tsz", "in_spikes_bg")]), - (spmask, spikes_bg, [("out_file", "in_mask")]), ]) # fmt: on - mosaic_mean = pe.Node( + mosaic_mean = pe.MapNode( PlotMosaic( out_file="plot_func_mean_mosaic1.svg", cmap="Greys_r", ), name="PlotMosaicMean", + iterfield=["in_file"], ) - mosaic_stddev = pe.Node( + mosaic_stddev = pe.MapNode( PlotMosaic( out_file="plot_func_stddev_mosaic2_stddev.svg", cmap="viridis", ), name="PlotMosaicSD", + iterfield=["in_file"], ) - ds_report_mean = pe.Node( + ds_report_mean = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="mean", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_mean", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) - ds_report_stdev = pe.Node( + ds_report_stdev = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="stdev", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_stdev", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) - ds_report_carpet = pe.Node( + ds_report_carpet = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="carpet", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_carpet", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) # fmt: off @@ -201,6 +216,7 @@ def init_func_report_wf(name="func_report_wf"): base_directory=reportlets_dir, desc="spikes", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_spikes", run_without_submitting=True, @@ -220,21 +236,24 @@ def init_func_report_wf(name="func_report_wf"): return workflow # Verbose-reporting goes here + from niworkflows.utils.connections import pop_file as _pop from nireports.interfaces import PlotContours - mosaic_zoom = pe.Node( + mosaic_zoom = pe.MapNode( PlotMosaic( cmap="Greys_r", ), name="PlotMosaicZoomed", + iterfield=["in_file"], ) - mosaic_noise = pe.Node( + mosaic_noise = pe.MapNode( PlotMosaic( only_noise=True, cmap="viridis_r", ), name="PlotMosaicNoise", + iterfield=["in_file"] ) if config.workflow.species.lower() in ("rat", "mouse"): @@ -254,24 +273,28 @@ def init_func_report_wf(name="func_report_wf"): name="PlotBrainmask", ) - ds_report_zoomed = pe.Node( + ds_report_zoomed = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="zoomed", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_zoomed", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) - ds_report_background = pe.Node( + ds_report_background = pe.MapNode( DerivativesDataSink( base_directory=reportlets_dir, desc="background", datatype="figures", + dismiss_entities=("part",) ), name="ds_report_background", run_without_submitting=True, + iterfield=["in_file", "source_file"], ) ds_report_bmask = pe.Node( @@ -279,6 +302,7 @@ def init_func_report_wf(name="func_report_wf"): base_directory=reportlets_dir, desc="brainmask", datatype="figures", + dismiss_entities=("part", "echo"), ), name="ds_report_bmask", run_without_submitting=True, @@ -289,6 +313,7 @@ def init_func_report_wf(name="func_report_wf"): base_directory=reportlets_dir, desc="norm", datatype="figures", + dismiss_entities=("part", "echo"), ), name="ds_report_norm", run_without_submitting=True, @@ -298,7 +323,7 @@ def init_func_report_wf(name="func_report_wf"): workflow.connect([ (inputnode, ds_report_norm, [("mni_report", "in_file"), ("name_source", "source_file")]), - (inputnode, plot_bmask, [("epi_mean", "in_file"), + (inputnode, plot_bmask, [(("epi_mean", _pop), "in_file"), ("brainmask", "in_contours")]), (inputnode, mosaic_zoom, [("epi_mean", "in_file"), ("brainmask", "bbox_mask_file")]), @@ -308,7 +333,7 @@ def init_func_report_wf(name="func_report_wf"): (inputnode, ds_report_bmask, [("name_source", "source_file")]), (mosaic_zoom, ds_report_zoomed, [("out_file", "in_file")]), (mosaic_noise, ds_report_background, [("out_file", "in_file")]), - (plot_bmask, ds_report_bmask, [("out_file", "in_file")]), + (plot_bmask, ds_report_bmask, [(("out_file", _pop), "in_file")]), ]) # fmt: on @@ -398,4 +423,7 @@ def _carpet_parcellation(segmentation, crown_mask): def _get_tr(meta_dict): + if isinstance(meta_dict, (list, tuple)): + meta_dict = meta_dict[0] + return meta_dict.get("RepetitionTime", None) diff --git a/mriqc/workflows/shared.py b/mriqc/workflows/shared.py new file mode 100644 index 000000000..d38bb6f47 --- /dev/null +++ b/mriqc/workflows/shared.py @@ -0,0 +1,120 @@ +# 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/ +# +"""Shared workflows.""" +from nipype.interfaces import utility as niu +from nipype.pipeline import engine as pe + + +def synthstrip_wf(name="synthstrip_wf", omp_nthreads=None): + """Create a brain-extraction workflow using SynthStrip.""" + from nipype.interfaces.ants import N4BiasFieldCorrection + from niworkflows.interfaces.nibabel import IntensityClip, ApplyMask + from mriqc.interfaces.synthstrip import SynthStrip + + inputnode = pe.Node(niu.IdentityInterface(fields=["in_files"]), name="inputnode") + outputnode = pe.Node( + niu.IdentityInterface(fields=["out_corrected", "out_brain", "bias_image", "out_mask"]), + name="outputnode", + ) + + # truncate target intensity for N4 correction + pre_clip = pe.Node(IntensityClip(p_min=10, p_max=99.9), name="pre_clip") + + pre_n4 = pe.Node( + N4BiasFieldCorrection( + dimension=3, + num_threads=omp_nthreads, + rescale_intensities=True, + copy_header=True, + ), + name="pre_n4", + ) + + post_n4 = pe.Node( + N4BiasFieldCorrection( + dimension=3, + save_bias=True, + num_threads=omp_nthreads, + n_iterations=[50] * 4, + copy_header=True, + ), + name="post_n4", + ) + + synthstrip = pe.Node( + SynthStrip(num_threads=omp_nthreads), + name="synthstrip", + num_threads=omp_nthreads, + ) + + final_masked = pe.Node(ApplyMask(), name="final_masked") + final_inu = pe.Node(niu.Function(function=_apply_bias_correction), name="final_inu") + + workflow = pe.Workflow(name=name) + # fmt: off + workflow.connect([ + (inputnode, final_inu, [("in_files", "in_file")]), + (inputnode, pre_clip, [("in_files", "in_file")]), + (pre_clip, pre_n4, [("out_file", "input_image")]), + (pre_n4, synthstrip, [("output_image", "in_file")]), + (synthstrip, post_n4, [("out_mask", "weight_image")]), + (synthstrip, final_masked, [("out_mask", "in_mask")]), + (pre_clip, post_n4, [("out_file", "input_image")]), + (post_n4, final_inu, [("bias_image", "bias_image")]), + (post_n4, final_masked, [("output_image", "in_file")]), + (final_masked, outputnode, [("out_file", "out_brain")]), + (post_n4, outputnode, [("bias_image", "bias_image")]), + (synthstrip, outputnode, [("out_mask", "out_mask")]), + (post_n4, outputnode, [("output_image", "out_corrected")]), + ]) + # fmt: on + return workflow + + +def _apply_bias_correction(in_file, bias_image, out_file=None): + import os.path as op + + import numpy as np + import nibabel as nb + + img = nb.load(in_file) + data = np.clip( + img.get_fdata() * nb.load(bias_image).get_fdata(), + a_min=0, + a_max=None, + ) + out_img = img.__class__( + data.astype(img.get_data_dtype()), + img.affine, + img.header, + ) + + if out_file is None: + fname, ext = op.splitext(op.basename(in_file)) + if ext == ".gz": + fname, ext2 = op.splitext(fname) + ext = ext2 + ext + out_file = op.abspath(f"{fname}_inu{ext}") + + out_img.to_filename(out_file) + return out_file diff --git a/mriqc/workflows/utils.py b/mriqc/workflows/utils.py index 690c88a67..bcdca0d68 100644 --- a/mriqc/workflows/utils.py +++ b/mriqc/workflows/utils.py @@ -25,7 +25,7 @@ def _tofloat(inlist): if isinstance(inlist, (list, tuple)): - return [float(el) for el in inlist] + return [_tofloat(el) for el in inlist] else: return float(inlist)