diff --git a/docs/api.rst b/docs/api.rst index 9145423ec..2b00b07c9 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -14,13 +14,13 @@ you get started. Setting up your development environment --------------------------------------- We believe that *fMRIPrep* must be free to use, inspect, and critique. -Correspondingly, you should be free to modify our software to improve it +Correspondingly, you should be free to modify our software to improve it or adapt it to new use cases and we especially welcome contributions to improve it or its documentation. We actively direct efforts into making the scrutiny and improvement processes as easy as possible. -As part of such efforts, we maintain some +As part of such efforts, we maintain some `tips and guidelines for developers `__ to help minimize your burden if you want to modify the software. @@ -35,3 +35,4 @@ Workflows .. automodule:: fmriprep.workflows.base .. automodule:: fmriprep.workflows.bold +.. automodule:: fmriprep.workflows.bold.fit diff --git a/docs/workflows.rst b/docs/workflows.rst index e5584219d..12351d5cc 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -80,6 +80,8 @@ single reference template (see `Longitudinal processing`_). skull_strip_fixed_seed=False, t1w=['sub-01/anat/sub-01_T1w.nii.gz'], t2w=[], + msm_sulc=True, + precomputed={}, ) .. important:: @@ -266,7 +268,8 @@ packages, including FreeSurfer and the `Connectome Workbench`_. from smriprep.workflows.surfaces import init_surface_recon_wf wf = init_surface_recon_wf(omp_nthreads=1, - hires=True) + hires=True, + precomputed={}) See also *sMRIPrep*'s :py:func:`~smriprep.workflows.surfaces.init_surface_recon_wf` diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py new file mode 100644 index 000000000..8cfcd2348 --- /dev/null +++ b/fmriprep/interfaces/resampling.py @@ -0,0 +1,641 @@ +"""Interfaces for resampling images in a single shot""" + +import asyncio +import os +from functools import partial + +import nibabel as nb +import nitransforms as nt +import numpy as np +from nipype.interfaces.base import ( + File, + InputMultiObject, + SimpleInterface, + TraitedSpec, + traits, +) +from nipype.utils.filemanip import fname_presuffix +from scipy import ndimage as ndi +from scipy.sparse import hstack as sparse_hstack +from sdcflows.transform import grid_bspline_weights +from sdcflows.utils.tools import ensure_positive_cosines + +from ..utils.asynctools import worker +from ..utils.transforms import load_transforms + + +class ResampleSeriesInputSpec(TraitedSpec): + in_file = File(exists=True, mandatory=True, desc="3D or 4D image file to resample") + ref_file = File(exists=True, mandatory=True, desc="File to resample in_file to") + transforms = InputMultiObject( + File(exists=True), mandatory=True, desc="Transform files, from in_file to ref_file (image mode)" + ) + inverse = InputMultiObject( + traits.Bool, + value=[False], + usedefault=True, + desc="Whether to invert each file in transforms", + ) + fieldmap = File(exists=True, desc="Fieldmap file resampled into reference space") + ro_time = traits.Float(desc="EPI readout time (s).") + pe_dir = traits.Enum( + "i", + "i-", + "j", + "j-", + "k", + "k-", + desc="the phase-encoding direction corresponding to in_data", + ) + num_threads = traits.Int(1, usedefault=True, desc="Number of threads to use for resampling") + + +class ResampleSeriesOutputSpec(TraitedSpec): + out_file = File(desc="Resampled image or series") + + +class ResampleSeries(SimpleInterface): + """Resample a time series, applying susceptibility and motion correction + simultaneously. + """ + + input_spec = ResampleSeriesInputSpec + output_spec = ResampleSeriesOutputSpec + + def _run_interface(self, runtime): + out_path = fname_presuffix(self.inputs.in_file, suffix='resampled', newpath=runtime.cwd) + + source = nb.load(self.inputs.in_file) + target = nb.load(self.inputs.ref_file) + fieldmap = nb.load(self.inputs.fieldmap) if self.inputs.fieldmap else None + + nvols = source.shape[3] if source.ndim > 3 else 1 + + transforms = load_transforms(self.inputs.transforms, self.inputs.inverse) + + pe_dir = self.inputs.pe_dir + ro_time = self.inputs.ro_time + pe_info = None + + if pe_dir and ro_time: + pe_axis = "ijk".index(pe_dir[0]) + pe_flip = pe_dir.endswith("-") + + # Nitransforms displacements are positive + source, axcodes = ensure_positive_cosines(source) + axis_flip = axcodes[pe_axis] in "LPI" + + pe_info = [(pe_axis, -ro_time if (axis_flip ^ pe_flip) else ro_time)] * nvols + + resampled = resample_bold( + source=source, + target=target, + transforms=transforms, + fieldmap=fieldmap, + pe_info=pe_info, + nthreads=self.inputs.num_threads, + ) + resampled.to_filename(out_path) + + self._results['out_file'] = out_path + return runtime + + +class ReconstructFieldmapInputSpec(TraitedSpec): + in_coeffs = InputMultiObject( + File(exists=True), mandatory=True, desc="SDCflows-style spline coefficient files" + ) + target_ref_file = File(exists=True, mandatory=True, desc="Image to reconstruct the field in alignment with") + fmap_ref_file = File(exists=True, mandatory=True, desc="Reference file aligned with coefficients") + transforms = InputMultiObject( + File(exists=True), mandatory=True, desc="Transform files, from in_file to ref_file (image mode)" + ) + inverse = InputMultiObject( + traits.Bool, + value=[False], + usedefault=True, + desc="Whether to invert each file in transforms", + ) + + +class ReconstructFieldmapOutputSpec(TraitedSpec): + out_file = File(desc="Fieldmap reconstructed in target_ref_file space") + + +class ReconstructFieldmap(SimpleInterface): + """Reconstruct a fieldmap from B-spline coefficients in a target space. + + If the target reference does not have an aligned grid (guaranteed if + transforms include a warp), then a reference file describing the space + where it is valid to extrapolate the field will be used as an intermediate + step. + """ + + input_spec = ReconstructFieldmapInputSpec + output_spec = ReconstructFieldmapOutputSpec + + def _run_interface(self, runtime): + out_path = fname_presuffix(self.inputs.in_coeffs[-1], suffix='rec', newpath=runtime.cwd) + + coefficients = [nb.load(coeff_file) for coeff_file in self.inputs.in_coeffs] + target = nb.load(self.inputs.target_ref_file) + fmapref = nb.load(self.inputs.fmap_ref_file) + + transforms = load_transforms(self.inputs.transforms, self.inputs.inverse) + + fieldmap = reconstruct_fieldmap( + coefficients=coefficients, + fmap_reference=fmapref, + target=target, + transforms=transforms, + ) + fieldmap.to_filename(out_path) + + self._results['out_file'] = out_path + return runtime + + +class DistortionParametersInputSpec(TraitedSpec): + in_file = File(exists=True, desc="EPI image corresponding to the metadata") + metadata = traits.Dict(mandatory=True, desc="metadata corresponding to the inputs") + + +class DistortionParametersOutputSpec(TraitedSpec): + readout_time = traits.Float + pe_direction = traits.Enum("i", "i-", "j", "j-", "k", "k-") + + +class DistortionParameters(SimpleInterface): + """Retrieve PhaseEncodingDirection and TotalReadoutTime from available metadata. + + One or both parameters may be missing; downstream interfaces should be prepared + to handle this. + """ + + input_spec = DistortionParametersInputSpec + output_spec = DistortionParametersOutputSpec + + def _run_interface(self, runtime): + from sdcflows.utils.epimanip import get_trt + + try: + self._results["readout_time"] = get_trt( + self.inputs.metadata, + self.inputs.in_file or None, + ) + self._results["pe_direction"] = self.inputs.metadata["PhaseEncodingDirection"] + except (KeyError, ValueError): + pass + + return runtime + + +def resample_vol( + data: np.ndarray, + coordinates: np.ndarray, + pe_info: tuple[int, float], + hmc_xfm: np.ndarray | None, + fmap_hz: np.ndarray, + output: np.dtype | np.ndarray | None = None, + order: int = 3, + mode: str = 'constant', + cval: float = 0.0, + prefilter: bool = True, +) -> np.ndarray: + """Resample a volume at specified coordinates + + This function implements simultaneous head-motion correction and + susceptibility-distortion correction. It accepts coordinates in + the source voxel space. It is the responsibility of the caller to + transform coordinates from any other target space. + + Parameters + ---------- + data + The data array to resample + coordinates + The first-approximation voxel coordinates to sample from ``data`` + The first dimension should have length ``data.ndim``. The further + dimensions have the shape of the target array. + pe_info + The readout vector in the form of (axis, signed-readout-time) + ``(1, -0.04)`` becomes ``[0, -0.04, 0]``, which indicates that a + +1 Hz deflection in the field shifts 0.04 voxels toward the start + of the data array in the second dimension. + hmc_xfm + Affine transformation accounting for head motion from the individual + volume into the BOLD reference space. This affine must be in VOX2VOX + form. + fmap_hz + The fieldmap, sampled to the target space, in Hz + output + The dtype or a pre-allocated array for sampling into the target space. + If pre-allocated, ``output.shape == coordinates.shape[1:]``. + order + Order of interpolation (default: 3 = cubic) + mode + How ``data`` is extended beyond its boundaries. See + :func:`scipy.ndimage.map_coordinates` for more details. + cval + Value to fill past edges of ``data`` if ``mode`` is ``'constant'``. + prefilter + Determines if ``data`` is pre-filtered before interpolation. + + Returns + ------- + resampled_array + The resampled array, with shape ``coordinates.shape[1:]``. + """ + if hmc_xfm is not None: + # Move image with the head + coords_shape = coordinates.shape + coordinates = nb.affines.apply_affine( + hmc_xfm, coordinates.reshape(coords_shape[0], -1).T + ).T.reshape(coords_shape) + + vsm = fmap_hz * pe_info[1] + coordinates[pe_info[0], ...] += vsm + + jacobian = 1 + np.gradient(vsm, axis=pe_info[0]) + + result = ndi.map_coordinates( + data, + coordinates, + output=output, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) + result *= jacobian + return result + + +async def resample_series_async( + data: np.ndarray, + coordinates: np.ndarray, + pe_info: list[tuple[int, float]], + hmc_xfms: list[np.ndarray] | None, + fmap_hz: np.ndarray, + output_dtype: np.dtype | None = None, + order: int = 3, + mode: str = 'constant', + cval: float = 0.0, + prefilter: bool = True, + max_concurrent: int = min(os.cpu_count(), 12), +) -> np.ndarray: + """Resample a 4D time series at specified coordinates + + This function implements simultaneous head-motion correction and + susceptibility-distortion correction. It accepts coordinates in + the source voxel space. It is the responsibility of the caller to + transform coordinates from any other target space. + + Parameters + ---------- + data + The data array to resample + coordinates + The first-approximation voxel coordinates to sample from ``data``. + The first dimension should have length 3. + The further dimensions determine the shape of the target array. + pe_info + A list of readout vectors in the form of (axis, signed-readout-time) + ``(1, -0.04)`` becomes ``[0, -0.04, 0]``, which indicates that a + +1 Hz deflection in the field shifts 0.04 voxels toward the start + of the data array in the second dimension. + hmc_xfm + A sequence of affine transformations accounting for head motion from + the individual volume into the BOLD reference space. + These affines must be in VOX2VOX form. + fmap_hz + The fieldmap, sampled to the target space, in Hz + output_dtype + The dtype of the output array. + order + Order of interpolation (default: 3 = cubic) + mode + How ``data`` is extended beyond its boundaries. See + :func:`scipy.ndimage.map_coordinates` for more details. + cval + Value to fill past edges of ``data`` if ``mode`` is ``'constant'``. + prefilter + Determines if ``data`` is pre-filtered before interpolation. + max_concurrent + Maximum number of volumes to resample concurrently + + Returns + ------- + resampled_array + The resampled array, with shape ``coordinates.shape[1:] + (N,)``, + where N is the number of volumes in ``data``. + """ + if data.ndim == 3: + return resample_vol( + data, + coordinates, + pe_info[0], + hmc_xfms[0] if hmc_xfms else None, + fmap_hz, + output_dtype, + order, + mode, + cval, + prefilter, + ) + + semaphore = asyncio.Semaphore(max_concurrent) + + out_array = np.zeros(coordinates.shape[1:] + data.shape[-1:], dtype=output_dtype) + + tasks = [ + asyncio.create_task( + worker( + partial( + resample_vol, + data=volume, + coordinates=coordinates.copy(), + pe_info=pe_info[volid], + hmc_xfm=hmc_xfms[volid] if hmc_xfms else None, + fmap_hz=fmap_hz, + output=out_array[..., volid], + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ), + semaphore, + ) + ) + for volid, volume in enumerate(np.rollaxis(data, -1, 0)) + ] + + await asyncio.gather(*tasks) + + return out_array + + +def resample_series( + data: np.ndarray, + coordinates: np.ndarray, + pe_info: list[tuple[int, float]], + hmc_xfms: list[np.ndarray] | None, + fmap_hz: np.ndarray, + output_dtype: np.dtype | None = None, + order: int = 3, + mode: str = 'constant', + cval: float = 0.0, + prefilter: bool = True, + nthreads: int = 1, +) -> np.ndarray: + """Resample a 4D time series at specified coordinates + + This function implements simultaneous head-motion correction and + susceptibility-distortion correction. It accepts coordinates in + the source voxel space. It is the responsibility of the caller to + transform coordinates from any other target space. + + Parameters + ---------- + data + The data array to resample + coordinates + The first-approximation voxel coordinates to sample from ``data``. + The first dimension should have length 3. + The further dimensions determine the shape of the target array. + pe_info + A list of readout vectors in the form of (axis, signed-readout-time) + ``(1, -0.04)`` becomes ``[0, -0.04, 0]``, which indicates that a + +1 Hz deflection in the field shifts 0.04 voxels toward the start + of the data array in the second dimension. + hmc_xfm + A sequence of affine transformations accounting for head motion from + the individual volume into the BOLD reference space. + These affines must be in VOX2VOX form. + fmap_hz + The fieldmap, sampled to the target space, in Hz + output_dtype + The dtype of the output array. + order + Order of interpolation (default: 3 = cubic) + mode + How ``data`` is extended beyond its boundaries. See + :func:`scipy.ndimage.map_coordinates` for more details. + cval + Value to fill past edges of ``data`` if ``mode`` is ``'constant'``. + prefilter + Determines if ``data`` is pre-filtered before interpolation. + nthreads + Number of threads to use for parallel resampling + + Returns + ------- + resampled_array + The resampled array, with shape ``coordinates.shape[1:] + (N,)``, + where N is the number of volumes in ``data``. + """ + return asyncio.run( + resample_series_async( + data=data, + coordinates=coordinates, + pe_info=pe_info, + hmc_xfms=hmc_xfms, + fmap_hz=fmap_hz, + output_dtype=output_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + max_concurrent=nthreads, + ) + ) + + +def resample_bold( + source: nb.Nifti1Image, + target: nb.Nifti1Image, + transforms: nt.TransformChain, + fieldmap: nb.Nifti1Image | None, + pe_info: list[tuple[int, float]] | None, + nthreads: int = 1, +) -> nb.Nifti1Image: + """Resample a 4D bold series into a target space, applying head-motion + and susceptibility-distortion correction simultaneously. + + Parameters + ---------- + source + The 4D bold series to resample. + target + An image sampled in the target space. + transforms + A nitransforms TransformChain that maps images from the individual + BOLD volume space into the target space. + fieldmap + The fieldmap, in Hz, sampled in the target space + pe_info + A list of readout vectors in the form of (axis, signed-readout-time) + ``(1, -0.04)`` becomes ``[0, -0.04, 0]``, which indicates that a + +1 Hz deflection in the field shifts 0.04 voxels toward the start + of the data array in the second dimension. + nthreads + Number of threads to use for parallel resampling + + Returns + ------- + resampled_bold + The BOLD series resampled into the target space + """ + if not isinstance(transforms, nt.TransformChain): + transforms = nt.TransformChain([transforms]) + if isinstance(transforms[-1], nt.linear.LinearTransformsMapping): + transform_list, hmc = transforms[:-1], transforms[-1] + else: + if any(isinstance(xfm, nt.linear.LinearTransformsMapping) for xfm in transforms): + classes = [xfm.__class__.__name__ for xfm in transforms] + raise ValueError(f"HMC transforms must come last. Found sequence: {classes}") + transform_list: list = transforms.transforms + hmc = None + + # Retrieve the RAS coordinates of the target space + coordinates = nt.base.SpatialReference.factory(target).ndcoords.astype('f4').T + + # We will operate in voxel space, so get the source affine + vox2ras = source.affine + ras2vox = np.linalg.inv(vox2ras) + # Transform RAS2RAS head motion transforms to VOX2VOX + if hmc is not None: + hmc_xfms = [ras2vox @ xfm.matrix @ vox2ras for xfm in transforms[-1]] + else: + hmc_xfms = None + + # Remove the head-motion transforms and add a mapping from boldref + # world space to voxels. This new transform maps from world coordinates + # in the target space to voxel coordinates in the source space. + ref2vox = nt.TransformChain(transforms[:-1] + [nt.Affine(ras2vox)]) + mapped_coordinates = ref2vox.map(coordinates) + + # Some identities to reduce special casing downstream + if fieldmap is None: + fieldmap = nb.Nifti1Image(np.zeros(target.shape[:3], dtype='f4'), target.affine) + if pe_info is None: + pe_info = [[0, 0] for _ in range(source.shape[-1])] + + resampled_data = resample_series( + data=source.get_fdata(dtype='f4'), + coordinates=mapped_coordinates.T.reshape((3, *target.shape[:3])), + pe_info=pe_info, + hmc_xfms=hmc_xfms, + fmap_hz=fieldmap.get_fdata(dtype='f4'), + output_dtype='f4', + nthreads=nthreads, + ) + resampled_img = nb.Nifti1Image(resampled_data, target.affine, target.header) + resampled_img.set_data_dtype('f4') + + return resampled_img + + +def aligned(aff1: np.ndarray, aff2: np.ndarray) -> bool: + """Determine if two affines have aligned grids""" + return np.allclose( + np.linalg.norm(np.cross(aff1[:-1, :-1].T, aff2[:-1, :-1].T), axis=1), + 0, + atol=1e-3, + ) + + +def as_affine(xfm: nt.base.TransformBase) -> nt.Affine | None: + # Identity transform + if type(xfm) is nt.base.TransformBase: + return nt.Affine() + + if isinstance(xfm, nt.Affine): + return xfm + + if isinstance(xfm, nt.TransformChain) and all(isinstance(x, nt.Affine) for x in xfm): + return xfm.asaffine() + + return None + + +def reconstruct_fieldmap( + coefficients: list[nb.Nifti1Image], + fmap_reference: nb.Nifti1Image, + target: nb.Nifti1Image, + transforms: nt.TransformChain, +) -> nb.Nifti1Image: + """Resample a fieldmap from B-Spline coefficients into a target space + + If the coefficients and target are aligned, the field is reconstructed + directly in the target space. + If not, then the field is reconstructed to the ``fmap_reference`` + resolution, and then resampled according to transforms. + + The former method only applies if the transform chain can be + collapsed to a single affine transform. + + Parameters + ---------- + coefficients + list of B-spline coefficient files. The affine matrices are used + to reconstruct the knot locations. + fmap_reference + The intermediate reference to reconstruct the fieldmap in, if + it cannot be reconstructed directly in the target space. + target + The target space to to resample the fieldmap into. + transforms + A nitransforms TransformChain that maps images from the fieldmap + space into the target space. + + Returns + ------- + fieldmap + The fieldmap encoded in ``coefficients``, resampled in the same + space as ``target`` + """ + + direct = False + affine_xfm = as_affine(transforms) + if affine_xfm is not None: + # Transforms maps RAS coordinates in the target to RAS coordinates in + # the fieldmap space. Composed with target.affine, we have a target voxel + # to fieldmap RAS affine. Hence, this is projected into fieldmap space. + projected_affine = affine_xfm.matrix @ target.affine + # If the coordinates have the same rotation from voxels, we can construct + # bspline weights efficiently. + direct = aligned(projected_affine, coefficients[-1].affine) + + if direct: + reference, _ = ensure_positive_cosines( + target.__class__(target.dataobj, projected_affine, target.header), + ) + else: + if not aligned(fmap_reference.affine, coefficients[-1].affine): + raise ValueError('Reference passed is not aligned with spline grids') + reference, _ = ensure_positive_cosines(fmap_reference) + + # Generate tensor-product B-Spline weights + colmat = sparse_hstack( + [grid_bspline_weights(reference, level) for level in coefficients] + ).tocsr() + coefficients = np.hstack( + [level.get_fdata(dtype='float32').reshape(-1) for level in coefficients] + ) + + # Reconstruct the fieldmap (in Hz) from coefficients + fmap_img = nb.Nifti1Image( + np.reshape(colmat @ coefficients, reference.shape[:3]), + reference.affine, + ) + + if not direct: + fmap_img = transforms.apply(fmap_img, reference=target) + + fmap_img.header.set_intent('estimate', name='fieldmap Hz') + fmap_img.header.set_data_dtype('float32') + fmap_img.header['cal_max'] = max((abs(fmap_img.dataobj.min()), fmap_img.dataobj.max())) + fmap_img.header['cal_min'] = -fmap_img.header['cal_max'] + + return fmap_img diff --git a/fmriprep/utils/asynctools.py b/fmriprep/utils/asynctools.py new file mode 100644 index 000000000..c61016849 --- /dev/null +++ b/fmriprep/utils/asynctools.py @@ -0,0 +1,10 @@ +import asyncio +from typing import Callable, TypeVar + +R = TypeVar('R') + + +async def worker(job: Callable[[], R], semaphore: asyncio.Semaphore) -> R: + async with semaphore: + loop = asyncio.get_running_loop() + return await loop.run_in_executor(None, job) diff --git a/fmriprep/utils/misc.py b/fmriprep/utils/misc.py index db1706641..0ae244c8b 100644 --- a/fmriprep/utils/misc.py +++ b/fmriprep/utils/misc.py @@ -21,6 +21,7 @@ # https://www.nipreps.org/community/licensing/ # """Miscellaneous utilities.""" +import typing as ty def check_deps(workflow): @@ -45,3 +46,21 @@ def fips_enabled(): fips = Path("/proc/sys/crypto/fips_enabled") return fips.exists() and fips.read_text()[0] != "0" + + +def estimate_bold_mem_usage(bold_fname: str) -> ty.Tuple[int, dict]: + import nibabel as nb + import numpy as np + + img = nb.load(bold_fname) + nvox = int(np.prod(img.shape, dtype='u8')) + # Assume tools will coerce to 8-byte floats to be safe + bold_size_gb = 8 * nvox / (1024**3) + bold_tlen = img.shape[-1] + mem_gb = { + "filesize": bold_size_gb, + "resampled": bold_size_gb * 4, + "largemem": bold_size_gb * (max(bold_tlen / 100, 1.0) + 4), + } + + return bold_tlen, mem_gb diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py new file mode 100644 index 000000000..738b7f797 --- /dev/null +++ b/fmriprep/utils/transforms.py @@ -0,0 +1,92 @@ +"""Utilities for loading transforms for resampling""" +from pathlib import Path + +import h5py +import nibabel as nb +import nitransforms as nt +import numpy as np +from nitransforms.io.itk import ITKCompositeH5 + + +def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.TransformBase: + """Load a series of transforms as a nitransforms TransformChain + + An empty list will return an identity transform + """ + if len(inverse) == 1: + inverse *= len(xfm_paths) + elif len(inverse) != len(xfm_paths): + raise ValueError("Mismatched number of transforms and inverses") + + chain = None + for path, inv in zip(xfm_paths[::-1], inverse[::-1]): + path = Path(path) + if path.suffix == '.h5': + xfm = load_ants_h5(path) + else: + xfm = nt.linear.load(path) + if inv: + xfm = ~xfm + if chain is None: + chain = xfm + else: + chain += xfm + if chain is None: + chain = nt.base.TransformBase() + return chain + + +def load_ants_h5(filename: Path) -> nt.TransformChain: + """Load ANTs H5 files as a nitransforms TransformChain""" + affine, warp, warp_affine = parse_combined_hdf5(filename) + warp_transform = nt.DenseFieldTransform(nb.Nifti1Image(warp, warp_affine)) + return nt.TransformChain([warp_transform, nt.Affine(affine)]) + + +def parse_combined_hdf5(h5_fn, to_ras=True): + # Borrowed from https://github.com/feilong/process + # process.resample.parse_combined_hdf5() + h = h5py.File(h5_fn) + xform = ITKCompositeH5.from_h5obj(h) + affine = xform[0].to_ras() + # Confirm these transformations are applicable + assert ( + h['TransformGroup']['2']['TransformType'][:][0] == b'DisplacementFieldTransform_float_3_3' + ) + assert np.array_equal( + h['TransformGroup']['2']['TransformFixedParameters'][:], + np.array( + [ + 193.0, + 229.0, + 193.0, + 96.0, + 132.0, + -78.0, + 1.0, + 1.0, + 1.0, + -1.0, + 0.0, + 0.0, + 0.0, + -1.0, + 0.0, + 0.0, + 0.0, + 1.0, + ] + ), + ) + warp = h['TransformGroup']['2']['TransformParameters'][:] + warp = warp.reshape((193, 229, 193, 3)).transpose(2, 1, 0, 3) + warp *= np.array([-1, -1, 1]) + warp_affine = np.array( + [ + [1.0, 0.0, 0.0, -96.0], + [0.0, 1.0, 0.0, -132.0], + [0.0, 0.0, 1.0, -78.0], + [0.0, 0.0, 0.0, 1.0], + ] + ) + return affine, warp, warp_affine diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index aa3fd48a1..2c2ad6aa4 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -43,7 +43,6 @@ from .. import config from ..interfaces import DerivativesDataSink from ..interfaces.reports import AboutSummary, SubjectSummary -from .bold.base import get_estimator, init_func_preproc_wf def init_fmriprep_wf(): @@ -91,7 +90,7 @@ def init_fmriprep_wf(): fsdir.inputs.subjects_dir = str(config.execution.fs_subjects_dir.absolute()) for subject_id in config.execution.participant_label: - single_subject_wf = init_single_subject_fit_wf(subject_id) + single_subject_wf = init_single_subject_wf(subject_id) single_subject_wf.config['execution']['crashdump_dir'] = str( config.execution.fmriprep_dir / f"sub-{subject_id}" / "log" / config.execution.run_uuid @@ -113,19 +112,84 @@ def init_fmriprep_wf(): return fmriprep_wf -def init_single_subject_fit_wf(subject_id: str): +def init_single_subject_wf(subject_id: str): + """ + Organize the preprocessing pipeline for a single subject. + + It collects and reports information about the subject, and prepares + sub-workflows to perform anatomical and functional preprocessing. + Anatomical preprocessing is performed in a single workflow, regardless of + the number of sessions. + Functional preprocessing is performed using a separate workflow for each + individual BOLD series. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from fmriprep.workflows.tests import mock_config + from fmriprep.workflows.base import init_single_subject_wf + with mock_config(): + wf = init_single_subject_wf('01') + + Parameters + ---------- + subject_id : :obj:`str` + Subject label for this single-subject workflow. + + Inputs + ------ + subjects_dir : :obj:`str` + FreeSurfer's ``$SUBJECTS_DIR``. + + """ from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.bids import BIDSDataGrabber, BIDSInfo + from niworkflows.interfaces.nilearn import NILEARN_VERSION from niworkflows.utils.bids import collect_data from niworkflows.utils.misc import fix_multi_T1w_source_name from niworkflows.utils.spaces import Reference from smriprep.workflows.anatomical import init_anat_fit_wf - from fmriprep.workflows.bold.fit import init_bold_fit_wf + from fmriprep.workflows.bold.base import init_bold_wf - spaces = config.workflow.spaces + workflow = Workflow(name=f'sub_{subject_id}_wf') + workflow.__desc__ = """ +Results included in this manuscript come from preprocessing +performed using *fMRIPrep* {fmriprep_ver} +(@fmriprep1; @fmriprep2; RRID:SCR_016216), +which is based on *Nipype* {nipype_ver} +(@nipype1; @nipype2; RRID:SCR_002502). + +""".format( + fmriprep_ver=config.environment.version, nipype_ver=config.environment.nipype_version + ) + workflow.__postdesc__ = """ + +Many internal operations of *fMRIPrep* use +*Nilearn* {nilearn_ver} [@nilearn, RRID:SCR_001362], +mostly within the functional processing workflow. +For more details of the pipeline, see [the section corresponding +to workflows in *fMRIPrep*'s documentation]\ +(https://fmriprep.readthedocs.io/en/latest/workflows.html \ +"FMRIPrep's documentation"). + + +### Copyright Waiver + +The above boilerplate text was automatically generated by fMRIPrep +with the express intention that users should copy and paste this +text into their manuscripts *unchanged*. +It is released under the [CC0]\ +(https://creativecommons.org/publicdomain/zero/1.0/) license. + +### References + +""".format( + nilearn_ver=NILEARN_VERSION + ) - workflow = Workflow(name=f'fit_sub_{subject_id}_wf') subject_data = collect_data( config.execution.layout, subject_id, @@ -134,6 +198,32 @@ def init_single_subject_fit_wf(subject_id: str): bids_filters=config.execution.bids_filters, )[0] + if 'flair' in config.workflow.ignore: + subject_data['flair'] = [] + if 't2w' in config.workflow.ignore: + subject_data['t2w'] = [] + + anat_only = config.workflow.anat_only + # Make sure we always go through these two checks + if not anat_only and not subject_data['bold']: + task_id = config.execution.task_id + raise RuntimeError( + "No BOLD images found for participant {} and task {}. " + "All workflows require BOLD images.".format( + subject_id, task_id if task_id else '' + ) + ) + + if subject_data['roi']: + warnings.warn( + f"Lesion mask {subject_data['roi']} found. " + "Future versions of fMRIPrep will use alternative conventions. " + "Please refer to the documentation before upgrading.", + FutureWarning, + ) + + spaces = config.workflow.spaces + anatomical_cache = {} if config.execution.derivatives: from smriprep.utils.bids import collect_derivatives as collect_anat_derivatives @@ -273,6 +363,13 @@ def init_single_subject_fit_wf(subject_id: str): output_dir=str(config.execution.fmriprep_dir), subject=subject_id, ) + fmap_wf.__desc__ = f""" + +Preprocessing of B0 inhomogeneity mappings + +: A total of {len(fmap_estimators)} fieldmaps were found available within the input +BIDS structure for this particular subject. +""" # Overwrite ``out_path_base`` of sdcflows's DataSinks for node in fmap_wf.list_node_names(): @@ -350,14 +447,27 @@ def init_single_subject_fit_wf(subject_id: str): ]) # fmt:on - for bold_file in subject_data['bold']: - fieldmap_id = estimator_map.get(listify(bold_file)[0]) + # Append the functional section to the existing anatomical excerpt + # That way we do not need to stream down the number of bold datasets + func_pre_desc = """ +Functional data preprocessing + +: For each of the {num_bold} BOLD runs found per subject (across all +tasks and sessions), the following preprocessing was performed. +""".format( + num_bold=len(subject_data['bold']) + ) + + for bold_series in subject_data['bold']: + bold_series = sorted(listify(bold_series)) + bold_file = bold_series[0] + fieldmap_id = estimator_map.get(bold_file) functional_cache = {} if config.execution.derivatives: from fmriprep.utils.bids import collect_derivatives, extract_entities - entities = extract_entities(bold_file) + entities = extract_entities(bold_series) for deriv_dir in config.execution.derivatives: functional_cache.update( @@ -367,16 +477,16 @@ def init_single_subject_fit_wf(subject_id: str): fieldmap_id=fieldmap_id, ) ) - func_fit_wf = init_bold_fit_wf( - bold_series=bold_file, + + bold_wf = init_bold_wf( + bold_series=bold_series, precomputed=functional_cache, fieldmap_id=fieldmap_id, - omp_nthreads=config.nipype.omp_nthreads, ) + bold_wf.__desc__ = func_pre_desc + (bold_wf.__desc__ or "") - # fmt:off workflow.connect([ - (anat_fit_wf, func_fit_wf, [ + (anat_fit_wf, bold_wf, [ ('outputnode.t1w_preproc', 'inputnode.t1w_preproc'), ('outputnode.t1w_mask', 'inputnode.t1w_mask'), ('outputnode.t1w_dseg', 'inputnode.t1w_dseg'), @@ -384,13 +494,10 @@ def init_single_subject_fit_wf(subject_id: str): ('outputnode.subject_id', 'inputnode.subject_id'), ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), ]), - ]) - # fmt: on - + ]) # fmt:skip if fieldmap_id: - # fmt:off workflow.connect([ - (fmap_wf, func_fit_wf, [ + (fmap_wf, bold_wf, [ ("outputnode.fmap", "inputnode.fmap"), ("outputnode.fmap_ref", "inputnode.fmap_ref"), ("outputnode.fmap_coeff", "inputnode.fmap_coeff"), @@ -398,497 +505,11 @@ def init_single_subject_fit_wf(subject_id: str): ("outputnode.fmap_id", "inputnode.fmap_id"), ("outputnode.method", "inputnode.sdc_method"), ]), - ]) - # fmt:on - - if config.workflow.level == "minimal": - return clean_datasinks(workflow) - - if config.workflow.level == "resampling": - return clean_datasinks(workflow) + ]) # fmt:skip return clean_datasinks(workflow) -def init_single_subject_wf(subject_id: str): - """ - Organize the preprocessing pipeline for a single subject. - - It collects and reports information about the subject, and prepares - sub-workflows to perform anatomical and functional preprocessing. - Anatomical preprocessing is performed in a single workflow, regardless of - the number of sessions. - Functional preprocessing is performed using a separate workflow for each - individual BOLD series. - - Workflow Graph - .. workflow:: - :graph2use: orig - :simple_form: yes - - from fmriprep.workflows.tests import mock_config - from fmriprep.workflows.base import init_single_subject_wf - with mock_config(): - wf = init_single_subject_wf('01') - - Parameters - ---------- - subject_id : :obj:`str` - Subject label for this single-subject workflow. - - Inputs - ------ - subjects_dir : :obj:`str` - FreeSurfer's ``$SUBJECTS_DIR``. - - """ - from niworkflows.engine.workflows import LiterateWorkflow as Workflow - from niworkflows.interfaces.bids import BIDSDataGrabber, BIDSInfo - from niworkflows.interfaces.nilearn import NILEARN_VERSION - from niworkflows.utils.bids import collect_data - from niworkflows.utils.misc import fix_multi_T1w_source_name - from niworkflows.utils.spaces import Reference - from smriprep.workflows.anatomical import init_anat_preproc_wf - - name = "single_subject_%s_wf" % subject_id - subject_data = collect_data( - config.execution.layout, - subject_id, - task=config.execution.task_id, - echo=config.execution.echo_idx, - bids_filters=config.execution.bids_filters, - )[0] - - if 'flair' in config.workflow.ignore: - subject_data['flair'] = [] - if 't2w' in config.workflow.ignore: - subject_data['t2w'] = [] - - anat_only = config.workflow.anat_only - spaces = config.workflow.spaces - # Make sure we always go through these two checks - if not anat_only and not subject_data['bold']: - task_id = config.execution.task_id - raise RuntimeError( - "No BOLD images found for participant {} and task {}. " - "All workflows require BOLD images.".format( - subject_id, task_id if task_id else '' - ) - ) - - deriv_cache = {} - if config.execution.derivatives: - from smriprep.utils.bids import collect_derivatives - - std_spaces = spaces.get_spaces(nonstandard=False, dim=(3,)) - std_spaces.append("fsnative") - for deriv_dir in config.execution.derivatives: - deriv_cache.update( - collect_derivatives( - derivatives_dir=deriv_dir, - subject_id=subject_id, - std_spaces=std_spaces, - ) - ) - - if "t1w_preproc" not in deriv_cache and not subject_data['t1w']: - raise Exception( - f"No T1w images found for participant {subject_id}. All workflows require T1w images." - ) - - if subject_data['roi']: - warnings.warn( - f"Lesion mask {subject_data['roi']} found. " - "Future versions of fMRIPrep will use alternative conventions. " - "Please refer to the documentation before upgrading.", - FutureWarning, - ) - - workflow = Workflow(name=name) - workflow.__desc__ = """ -Results included in this manuscript come from preprocessing -performed using *fMRIPrep* {fmriprep_ver} -(@fmriprep1; @fmriprep2; RRID:SCR_016216), -which is based on *Nipype* {nipype_ver} -(@nipype1; @nipype2; RRID:SCR_002502). - -""".format( - fmriprep_ver=config.environment.version, nipype_ver=config.environment.nipype_version - ) - workflow.__postdesc__ = """ - -Many internal operations of *fMRIPrep* use -*Nilearn* {nilearn_ver} [@nilearn, RRID:SCR_001362], -mostly within the functional processing workflow. -For more details of the pipeline, see [the section corresponding -to workflows in *fMRIPrep*'s documentation]\ -(https://fmriprep.readthedocs.io/en/latest/workflows.html \ -"FMRIPrep's documentation"). - - -### Copyright Waiver - -The above boilerplate text was automatically generated by fMRIPrep -with the express intention that users should copy and paste this -text into their manuscripts *unchanged*. -It is released under the [CC0]\ -(https://creativecommons.org/publicdomain/zero/1.0/) license. - -### References - -""".format( - nilearn_ver=NILEARN_VERSION - ) - - fmriprep_dir = str(config.execution.fmriprep_dir) - - inputnode = pe.Node(niu.IdentityInterface(fields=['subjects_dir']), name='inputnode') - - bidssrc = pe.Node( - BIDSDataGrabber( - subject_data=subject_data, - anat_only=anat_only, - subject_id=subject_id, - ), - name='bidssrc', - ) - - bids_info = pe.Node( - BIDSInfo(bids_dir=config.execution.bids_dir, bids_validate=False), name='bids_info' - ) - - summary = pe.Node( - SubjectSummary( - std_spaces=spaces.get_spaces(nonstandard=False), - nstd_spaces=spaces.get_spaces(standard=False), - ), - name='summary', - run_without_submitting=True, - ) - - about = pe.Node( - AboutSummary(version=config.environment.version, command=' '.join(sys.argv)), - name='about', - run_without_submitting=True, - ) - - ds_report_summary = pe.Node( - DerivativesDataSink( - base_directory=fmriprep_dir, - desc='summary', - datatype="figures", - dismiss_entities=("echo",), - ), - name='ds_report_summary', - run_without_submitting=True, - ) - - ds_report_about = pe.Node( - DerivativesDataSink( - base_directory=fmriprep_dir, - desc='about', - datatype="figures", - dismiss_entities=("echo",), - ), - name='ds_report_about', - run_without_submitting=True, - ) - - # Preprocessing of T1w (includes registration to MNI) - anat_preproc_wf = init_anat_preproc_wf( - bids_root=str(config.execution.bids_dir), - sloppy=config.execution.sloppy, - debug=config.execution.debug, - precomputed=deriv_cache, - freesurfer=config.workflow.run_reconall, - hires=config.workflow.hires, - longitudinal=config.workflow.longitudinal, - omp_nthreads=config.nipype.omp_nthreads, - msm_sulc=config.workflow.run_msmsulc, - output_dir=fmriprep_dir, - skull_strip_fixed_seed=config.workflow.skull_strip_fixed_seed, - skull_strip_mode=config.workflow.skull_strip_t1w, - skull_strip_template=Reference.from_string(config.workflow.skull_strip_template)[0], - spaces=spaces, - t1w=subject_data['t1w'], - t2w=subject_data['t2w'], - cifti_output=config.workflow.cifti_output, - ) - # fmt:off - workflow.connect([ - (inputnode, anat_preproc_wf, [('subjects_dir', 'inputnode.subjects_dir')]), - (inputnode, summary, [('subjects_dir', 'subjects_dir')]), - (bidssrc, summary, [('bold', 'bold')]), - (bids_info, summary, [('subject', 'subject_id')]), - (bids_info, anat_preproc_wf, [(('subject', _prefix), 'inputnode.subject_id')]), - (bidssrc, anat_preproc_wf, [('t1w', 'inputnode.t1w'), - ('t2w', 'inputnode.t2w'), - ('roi', 'inputnode.roi'), - ('flair', 'inputnode.flair')]), - (summary, ds_report_summary, [('out_report', 'in_file')]), - (about, ds_report_about, [('out_report', 'in_file')]), - ]) - - workflow.connect([ - (bidssrc, bids_info, [(('t1w', fix_multi_T1w_source_name), 'in_file')]), - (bidssrc, summary, [('t1w', 't1w'), - ('t2w', 't2w')]), - (bidssrc, ds_report_summary, [(('t1w', fix_multi_T1w_source_name), 'source_file')]), - (bidssrc, ds_report_about, [(('t1w', fix_multi_T1w_source_name), 'source_file')]), - ]) - # fmt:on - - # Overwrite ``out_path_base`` of smriprep's DataSinks - for node in workflow.list_node_names(): - if node.split('.')[-1].startswith('ds_'): - workflow.get_node(node).interface.out_path_base = "" - - if anat_only: - return workflow - - from sdcflows import fieldmaps as fm - - fmap_estimators = None - - if any( - ( - "fieldmaps" not in config.workflow.ignore, - config.workflow.use_syn_sdc, - config.workflow.force_syn, - ) - ): - from sdcflows.utils.wrangler import find_estimators - - # SDC Step 1: Run basic heuristics to identify available data for fieldmap estimation - # For now, no fmapless - filters = None - if config.execution.bids_filters is not None: - filters = config.execution.bids_filters.get("fmap") - - # In the case where fieldmaps are ignored and `--use-syn-sdc` is requested, - # SDCFlows `find_estimators` still receives a full layout (which includes the fmap modality) - # and will not calculate fmapless schemes. - # Similarly, if fieldmaps are ignored and `--force-syn` is requested, - # `fmapless` should be set to True to ensure BOLD targets are found to be corrected. - fmapless = bool(config.workflow.use_syn_sdc) or ( - "fieldmaps" in config.workflow.ignore and config.workflow.force_syn - ) - force_fmapless = config.workflow.force_syn or ( - "fieldmaps" in config.workflow.ignore and config.workflow.use_syn_sdc - ) - - fmap_estimators = find_estimators( - layout=config.execution.layout, - subject=subject_id, - fmapless=fmapless, - force_fmapless=force_fmapless, - bids_filters=filters, - ) - - if config.workflow.use_syn_sdc and not fmap_estimators: - message = ( - "Fieldmap-less (SyN) estimation was requested, but PhaseEncodingDirection " - "information appears to be absent." - ) - config.loggers.workflow.error(message) - if config.workflow.use_syn_sdc == "error": - raise ValueError(message) - - if "fieldmaps" in config.workflow.ignore and any( - f.method == fm.EstimatorType.ANAT for f in fmap_estimators - ): - config.loggers.workflow.info( - 'Option "--ignore fieldmaps" was set, but either "--use-syn-sdc" ' - 'or "--force-syn" were given, so fieldmap-less estimation will be executed.' - ) - fmap_estimators = [f for f in fmap_estimators if f.method == fm.EstimatorType.ANAT] - - # Do not calculate fieldmaps that we will not use - if fmap_estimators: - used_estimators = { - key - for bold_file in subject_data['bold'] - for key in get_estimator(config.execution.layout, listify(bold_file)[0]) - } - - fmap_estimators = [fmap for fmap in fmap_estimators if fmap.bids_id in used_estimators] - - # Simplification: Unused estimators are removed from registry - # This fiddles with a private attribute, so it may break in future - # versions. However, it does mean the BOLD workflow doesn't need to - # replicate the logic that got us to the pared down set of estimators - # here. - final_ids = {fmap.bids_id for fmap in fmap_estimators} - unused_ids = fm._estimators.keys() - final_ids - for bids_id in unused_ids: - del fm._estimators[bids_id] - - if fmap_estimators: - config.loggers.workflow.info( - "B0 field inhomogeneity map will be estimated with " - f"the following {len(fmap_estimators)} estimator(s): " - f"{[e.method for e in fmap_estimators]}." - ) - - # Append the functional section to the existing anatomical excerpt - # That way we do not need to stream down the number of bold datasets - func_pre_desc = """ -Functional data preprocessing - -: For each of the {num_bold} BOLD runs found per subject (across all -tasks and sessions), the following preprocessing was performed. -""".format( - num_bold=len(subject_data['bold']) - ) - - func_preproc_wfs = [] - has_fieldmap = bool(fmap_estimators) - for bold_file in subject_data['bold']: - func_preproc_wf = init_func_preproc_wf(bold_file, has_fieldmap=has_fieldmap) - if func_preproc_wf is None: - continue - - func_preproc_wf.__desc__ = func_pre_desc + (func_preproc_wf.__desc__ or "") - # fmt:off - workflow.connect([ - (anat_preproc_wf, func_preproc_wf, [ - ('outputnode.t1w_preproc', 'inputnode.t1w_preproc'), - ('outputnode.t1w_mask', 'inputnode.t1w_mask'), - ('outputnode.t1w_dseg', 'inputnode.t1w_dseg'), - ('outputnode.t1w_aseg', 'inputnode.t1w_aseg'), - ('outputnode.t1w_aparc', 'inputnode.t1w_aparc'), - ('outputnode.t1w_tpms', 'inputnode.t1w_tpms'), - ('outputnode.template', 'inputnode.template'), - ('outputnode.anat2std_xfm', 'inputnode.anat2std_xfm'), - ('outputnode.std2anat_xfm', 'inputnode.std2anat_xfm'), - # 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.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) - - if not has_fieldmap: - return workflow - - from sdcflows.workflows.base import init_fmap_preproc_wf - - fmap_wf = init_fmap_preproc_wf( - debug="fieldmaps" in config.execution.debug, - estimators=fmap_estimators, - omp_nthreads=config.nipype.omp_nthreads, - output_dir=fmriprep_dir, - subject=subject_id, - ) - fmap_wf.__desc__ = f""" - -Preprocessing of B0 inhomogeneity mappings - -: A total of {len(fmap_estimators)} fieldmaps were found available within the input -BIDS structure for this particular subject. -""" - for func_preproc_wf in func_preproc_wfs: - # fmt:off - workflow.connect([ - (fmap_wf, func_preproc_wf, [ - ("outputnode.fmap", "inputnode.fmap"), - ("outputnode.fmap_ref", "inputnode.fmap_ref"), - ("outputnode.fmap_coeff", "inputnode.fmap_coeff"), - ("outputnode.fmap_mask", "inputnode.fmap_mask"), - ("outputnode.fmap_id", "inputnode.fmap_id"), - ("outputnode.method", "inputnode.sdc_method"), - ]), - ]) - # fmt:on - - # Overwrite ``out_path_base`` of sdcflows's DataSinks - for node in fmap_wf.list_node_names(): - if node.split(".")[-1].startswith("ds_"): - fmap_wf.get_node(node).interface.out_path_base = "" - - # Step 3: Manually connect PEPOLAR and ANAT workflows - - # Select "MNI152NLin2009cAsym" from standard references. - # This node may be used by multiple ANAT estimators, so define outside loop. - from niworkflows.interfaces.utility import KeySelect - - fmap_select_std = pe.Node( - KeySelect(fields=["std2anat_xfm"], key="MNI152NLin2009cAsym"), - name="fmap_select_std", - run_without_submitting=True, - ) - if any(estimator.method == fm.EstimatorType.ANAT for estimator in fmap_estimators): - # fmt:off - workflow.connect([ - (anat_preproc_wf, fmap_select_std, [ - ("outputnode.std2anat_xfm", "std2anat_xfm"), - ("outputnode.template", "keys")]), - ]) - # fmt:on - - for estimator in fmap_estimators: - config.loggers.workflow.info( - f"""\ -Setting-up fieldmap "{estimator.bids_id}" ({estimator.method}) with \ -<{', '.join(s.path.name for s in estimator.sources)}>""" - ) - - # Mapped and phasediff can be connected internally by SDCFlows - if estimator.method in (fm.EstimatorType.MAPPED, fm.EstimatorType.PHASEDIFF): - continue - - suffices = [s.suffix for s in estimator.sources] - - if estimator.method == fm.EstimatorType.PEPOLAR: - if len(suffices) == 2 and all(suf in ("epi", "bold", "sbref") for suf in suffices): - wf_inputs = getattr(fmap_wf.inputs, f"in_{estimator.bids_id}") - wf_inputs.in_data = [str(s.path) for s in estimator.sources] - wf_inputs.metadata = [s.metadata for s in estimator.sources] - else: - raise NotImplementedError("Sophisticated PEPOLAR schemes are unsupported.") - - elif estimator.method == fm.EstimatorType.ANAT: - from sdcflows.workflows.fit.syn import init_syn_preprocessing_wf - - sources = [str(s.path) for s in estimator.sources if s.suffix in ("bold", "sbref")] - source_meta = [s.metadata for s in estimator.sources if s.suffix in ("bold", "sbref")] - syn_preprocessing_wf = init_syn_preprocessing_wf( - omp_nthreads=config.nipype.omp_nthreads, - debug=config.execution.sloppy, - auto_bold_nss=True, - t1w_inversion=False, - name=f"syn_preprocessing_{estimator.bids_id}", - ) - syn_preprocessing_wf.inputs.inputnode.in_epis = sources - syn_preprocessing_wf.inputs.inputnode.in_meta = source_meta - - # fmt:off - workflow.connect([ - (anat_preproc_wf, syn_preprocessing_wf, [ - ("outputnode.t1w_preproc", "inputnode.in_anat"), - ("outputnode.t1w_mask", "inputnode.mask_anat"), - ]), - (fmap_select_std, syn_preprocessing_wf, [ - ("std2anat_xfm", "inputnode.std2anat_xfm"), - ]), - (syn_preprocessing_wf, fmap_wf, [ - ("outputnode.epi_ref", f"in_{estimator.bids_id}.epi_ref"), - ("outputnode.epi_mask", f"in_{estimator.bids_id}.epi_mask"), - ("outputnode.anat_ref", f"in_{estimator.bids_id}.anat_ref"), - ("outputnode.anat_mask", f"in_{estimator.bids_id}.anat_mask"), - ("outputnode.sd_prior", f"in_{estimator.bids_id}.sd_prior"), - ]), - ]) - # fmt:on - return workflow - - def map_fieldmap_estimation( layout: bids.BIDSLayout, subject_id: str, @@ -938,7 +559,7 @@ def map_fieldmap_estimation( # Pare down estimators to those that are actually used # If fmap_estimators == [], all loops/comprehensions terminate immediately all_ids = {fmap.bids_id for fmap in fmap_estimators} - bold_files = (listify(bold_file)[0] for bold_file in bold_data) + bold_files = (sorted(listify(bold_file))[0] for bold_file in bold_data) all_estimators = { bold_file: [fmap_id for fmap_id in get_estimator(layout, bold_file) if fmap_id in all_ids] @@ -975,3 +596,21 @@ def clean_datasinks(workflow: pe.Workflow) -> pe.Workflow: if node.split('.')[-1].startswith('ds_'): workflow.get_node(node).interface.out_path_base = "" return workflow + + +def get_estimator(layout, fname): + field_source = layout.get_metadata(fname).get("B0FieldSource") + if isinstance(field_source, str): + field_source = (field_source,) + + if field_source is None: + import re + from pathlib import Path + + from sdcflows.fieldmaps import get_identifier + + # Fallback to IntendedFor + intended_rel = re.sub(r"^sub-[a-zA-Z0-9]*/", "", str(Path(fname).relative_to(layout.root))) + field_source = get_identifier(intended_rel) + + return field_source diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py new file mode 100644 index 000000000..449f81ae1 --- /dev/null +++ b/fmriprep/workflows/bold/apply.py @@ -0,0 +1,51 @@ +from __future__ import annotations + +import os +import typing as ty + +import nibabel as nb +import nipype.interfaces.utility as niu +import nipype.pipeline.engine as pe +from niworkflows.interfaces.header import ValidateImage +from niworkflows.interfaces.utility import KeySelect +from niworkflows.utils.connections import listify + +from fmriprep import config + +from ...interfaces.resampling import ( + DistortionParameters, + ReconstructFieldmap, + ResampleSeries, +) +from ...utils.misc import estimate_bold_mem_usage +from .stc import init_bold_stc_wf +from .t2s import init_bold_t2s_wf, init_t2s_reporting_wf + +if ty.TYPE_CHECKING: + from niworkflows.utils.spaces import SpatialReferences + + +def init_bold_apply_wf( + *, + spaces: SpatialReferences, + name: str = 'bold_apply_wf', +) -> pe.Workflow: + """TODO: Docstring""" + from smriprep.workflows.outputs import init_template_iterator_wf + + workflow = pe.Workflow(name=name) + + if spaces.is_cached() and spaces.cached.references: + template_iterator_wf = init_template_iterator_wf(spaces=spaces) + # TODO: Refactor bold_std_trans_wf + # bold_std_trans_wf = init_bold_std_trans_wf( + # freesurfer=config.workflow.run_reconall, + # mem_gb=mem_gb["resampled"], + # omp_nthreads=config.nipype.omp_nthreads, + # spaces=spaces, + # multiecho=multiecho, + # use_compression=not config.execution.low_mem, + # name="bold_std_trans_wf", + # ) + + return workflow diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 5352c8ccd..4e4ad6e4d 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -24,11 +24,13 @@ Orchestrating the BOLD-preprocessing workflow ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -.. autofunction:: init_func_preproc_wf -.. autofunction:: init_func_derivatives_wf +.. autofunction:: init_bold_wf +.. autofunction:: init_bold_fit_wf +.. autofunction:: init_bold_native_wf """ import os +import typing as ty import nibabel as nb import numpy as np @@ -44,8 +46,9 @@ # BOLD workflows from .confounds import init_bold_confs_wf, init_carpetplot_wf +from .fit import init_bold_fit_wf, init_bold_native_wf from .hmc import init_bold_hmc_wf -from .outputs import init_func_derivatives_wf +from .outputs import init_ds_bold_native_wf, init_func_derivatives_wf from .registration import init_bold_reg_wf, init_bold_t1_trans_wf from .resampling import ( init_bold_preproc_trans_wf, @@ -56,6 +59,265 @@ from .t2s import init_bold_t2s_wf, init_t2s_reporting_wf +def init_bold_wf( + *, + bold_series: ty.List[str], + precomputed: dict = {}, + fieldmap_id: ty.Optional[str] = None, + name: str = "bold_wf", +) -> pe.Workflow: + """ + This workflow controls the functional preprocessing stages of *fMRIPrep*. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from fmriprep.workflows.tests import mock_config + from fmriprep import config + from fmriprep.workflows.bold.base import init_bold_wf + with mock_config(): + bold_file = config.execution.bids_dir / "sub-01" / "func" \ + / "sub-01_task-mixedgamblestask_run-01_bold.nii.gz" + wf = init_bold_wf( + bold_series=[str(bold_file)], + ) + + Parameters + ---------- + bold_series + List of paths to NIfTI files. + precomputed + Dictionary containing precomputed derivatives to reuse, if possible. + fieldmap_id + ID of the fieldmap to use to correct this BOLD series. If :obj:`None`, + no correction will be applied. + + Inputs + ------ + t1w_preproc + Bias-corrected structural template image + t1w_mask + Mask of the skull-stripped template image + t1w_dseg + Segmentation of preprocessed structural image, including + gray-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF) + anat2std_xfm + List of transform files, collated with templates + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + fsnative2t1w_xfm + LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w + fmap_id + Unique identifiers to select fieldmap files + fmap + List of estimated fieldmaps (collated with fmap_id) + fmap_ref + List of fieldmap reference files (collated with fmap_id) + fmap_coeff + List of lists of spline coefficient files (collated with fmap_id) + fmap_mask + List of fieldmap masks (collated with fmap_id) + sdc_method + List of fieldmap correction method names (collated with fmap_id) + + See Also + -------- + + * :func:`~fmriprep.workflows.bold.fit.init_bold_fit_wf` + * :func:`~fmriprep.workflows.bold.fit.init_bold_native_wf` + * :func:`~fmriprep.workflows.bold.outputs.init_ds_bold_native_wf` + * :func:`~fmriprep.workflows.bold.t2s.init_t2s_reporting_wf` + + """ + from niworkflows.engine.workflows import LiterateWorkflow as Workflow + + bold_file = bold_series[0] + + fmriprep_dir = config.execution.fmriprep_dir + omp_nthreads = config.nipype.omp_nthreads + + functional_cache = {} + if config.execution.derivatives: + from fmriprep.utils.bids import collect_derivatives, extract_entities + + entities = extract_entities(bold_series) + + for deriv_dir in config.execution.derivatives: + functional_cache.update( + collect_derivatives( + derivatives_dir=deriv_dir, + entities=entities, + fieldmap_id=fieldmap_id, + ) + ) + + workflow = Workflow(name=_get_wf_name(bold_file, "bold")) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + # Anatomical coregistration + "t1w_preproc", + "t1w_mask", + "t1w_dseg", + "subjects_dir", + "subject_id", + "fsnative2t1w_xfm", + # Fieldmap registration + "fmap", + "fmap_ref", + "fmap_coeff", + "fmap_mask", + "fmap_id", + "sdc_method", + ], + ), + name="inputnode", + ) + + # + # Minimal workflow + # + + bold_fit_wf = init_bold_fit_wf( + bold_series=bold_series, + precomputed=functional_cache, + fieldmap_id=fieldmap_id, + omp_nthreads=omp_nthreads, + ) + + workflow.connect([ + (inputnode, bold_fit_wf, [ + ('t1w_preproc', 'inputnode.t1w_preproc'), + ('t1w_mask', 'inputnode.t1w_mask'), + ('t1w_dseg', 'inputnode.t1w_dseg'), + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id'), + ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ("fmap", "inputnode.fmap"), + ("fmap_ref", "inputnode.fmap_ref"), + ("fmap_coeff", "inputnode.fmap_coeff"), + ("fmap_mask", "inputnode.fmap_mask"), + ("fmap_id", "inputnode.fmap_id"), + ("sdc_method", "inputnode.sdc_method"), + ]), + ]) # fmt:skip + + if config.workflow.level == "minimal": + return workflow + + # Now that we're resampling and combining, multiecho matters + multiecho = len(bold_series) > 2 + + spaces = config.workflow.spaces + nonstd_spaces = set(spaces.get_nonstandard()) + template_spaces = spaces.get_spaces(nonstandard=False, dim=(3,)) + freesurfer_spaces = spaces.get_fs_spaces() + + # + # Resampling outputs workflow: + # - Resample to native + # - Save native outputs/echos only if requested + # + + bold_native_wf = init_bold_native_wf(bold_series=bold_series, fieldmap_id=fieldmap_id) + + workflow.connect([ + (bold_fit_wf, bold_native_wf, [ + ("outputnode.coreg_boldref", "inputnode.boldref"), + ("outputnode.bold_mask", "inputnode.bold_mask"), + ("outputnode.motion_xfm", "inputnode.motion_xfm"), + ("outputnode.boldref2fmap_xfm", "inputnode.fmapreg_xfm"), + ("outputnode.dummy_scans", "inputnode.dummy_scans"), + ]), + ]) # fmt:skip + + if fieldmap_id: + workflow.connect([ + (inputnode, bold_native_wf, [ + ("fmap_ref", "inputnode.fmap_ref"), + ("fmap_coeff", "inputnode.fmap_coeff"), + ("fmap_id", "inputnode.fmap_id"), + ]), + ]) # fmt:skip + + boldref_out = bool(nonstd_spaces.intersection(('func', 'run', 'bold', 'boldref', 'sbref'))) + echos_out = multiecho and config.execution.me_output_echos + + if boldref_out or echos_out: + ds_bold_native_wf = init_ds_bold_native_wf( + bids_root=str(config.execution.bids_dir), + output_dir=fmriprep_dir, + bold_output=boldref_out, + echo_output=echos_out, + multiecho=multiecho, + all_metadata=[config.execution.layout.get_metadata(file) for file in bold_series], + ) + ds_bold_native_wf.inputs.inputnode.source_files = bold_series + + workflow.connect([ + (bold_fit_wf, ds_bold_native_wf, [ + ('outputnode.bold_mask', 'inputnode.bold_mask'), + ]), + (bold_native_wf, ds_bold_native_wf, [ + ('outputnode.bold_native', 'inputnode.bold'), + ('outputnode.bold_echos', 'inputnode.bold_echos'), + ('outputnode.t2star_map', 'inputnode.t2star'), + ]), + ]) # fmt:skip + + if multiecho: + t2s_reporting_wf = init_t2s_reporting_wf() + + ds_report_t2scomp = pe.Node( + DerivativesDataSink( + desc="t2scomp", + datatype="figures", + dismiss_entities=("echo",), + ), + name="ds_report_t2scomp", + run_without_submitting=True, + ) + + ds_report_t2star_hist = pe.Node( + DerivativesDataSink( + desc="t2starhist", + datatype="figures", + dismiss_entities=("echo",), + ), + name="ds_report_t2star_hist", + run_without_submitting=True, + ) + + workflow.connect([ + (inputnode, t2s_reporting_wf, [('t1w_dseg', 'inputnode.label_file')]), + (bold_fit_wf, t2s_reporting_wf, [ + ('outputnode.boldref2anat_xfm', 'inputnode.boldref2anat_xfm'), + ('outputnode.coreg_boldref', 'inputnode.boldref'), + ]), + (bold_native_wf, t2s_reporting_wf, [ + ('outputnode.t2star_map', 'inputnode.t2star_file'), + ]), + (t2s_reporting_wf, ds_report_t2scomp, [('outputnode.t2s_comp_report', 'in_file')]), + (t2s_reporting_wf, ds_report_t2star_hist, [("outputnode.t2star_hist", "in_file")]), + ]) # fmt:skip + + if config.workflow.level == "resampling": + return workflow + + # Fill-in datasinks of reportlets seen so far + for node in workflow.list_node_names(): + if node.split(".")[-1].startswith("ds_report"): + workflow.get_node(node).inputs.base_directory = fmriprep_dir + workflow.get_node(node).inputs.source_file = bold_file + + return workflow + + def init_func_preproc_wf(bold_file, has_fieldmap=False): """ This workflow controls the functional preprocessing stages of *fMRIPrep*. @@ -255,7 +517,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): if os.path.isfile(ref_file): bold_tlen, mem_gb = _create_mem_gb(ref_file) - wf_name = _get_wf_name(ref_file) + wf_name = _get_wf_name(ref_file, "func_preproc") config.loggers.workflow.debug( "Creating bold processing workflow for <%s> (%.2f GB / %d TRs). " "Memory resampled/largemem=%.2f/%.2f GB.", @@ -420,7 +682,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): summary = pe.Node( FunctionalSummary( - slice_timing=run_stc, + # slice_timing=run_stc, registration=("FSL", "FreeSurfer")[freesurfer], registration_dof=config.workflow.bold2t1w_dof, registration_init=config.workflow.bold2t1w_init, @@ -605,109 +867,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False): ) final_boldref_wf.__desc__ = None # Unset description to avoid second appearance - # MAIN WORKFLOW STRUCTURE ####################################################### - # fmt:off - workflow.connect([ - # Prepare masked T1w image - (inputnode, t1w_brain, [("t1w_preproc", "in_file"), - ("t1w_mask", "in_mask")]), - # Select validated bold files per-echo - (initial_boldref_wf, select_bold, [("outputnode.all_bold_files", "inlist")]), - # BOLD buffer has slice-time corrected if it was run, original otherwise - (boldbuffer, bold_split, [("bold_file", "in_file")]), - # HMC - (initial_boldref_wf, bold_hmc_wf, [ - ("outputnode.raw_ref_image", "inputnode.raw_ref_image"), - ("outputnode.bold_file", "inputnode.bold_file"), - ]), - (bold_hmc_wf, outputnode, [ - ("outputnode.xforms", "hmc_xforms"), - ]), - # EPI-T1w registration workflow - (inputnode, bold_reg_wf, [ - ("t1w_dseg", "inputnode.t1w_dseg"), - # Undefined if --fs-no-reconall, but this is safe - ("subjects_dir", "inputnode.subjects_dir"), - ("subject_id", "inputnode.subject_id"), - ("fsnative2t1w_xfm", "inputnode.fsnative2t1w_xfm"), - ]), - (bold_final, bold_reg_wf, [ - ("boldref", "inputnode.ref_bold_brain")]), - (t1w_brain, bold_reg_wf, [("out_file", "inputnode.t1w_brain")]), - (inputnode, bold_t1_trans_wf, [ - ("bold_file", "inputnode.name_source"), - ("t1w_mask", "inputnode.t1w_mask"), - ("t1w_aseg", "inputnode.t1w_aseg"), - ("t1w_aparc", "inputnode.t1w_aparc"), - ]), - (t1w_brain, bold_t1_trans_wf, [("out_file", "inputnode.t1w_brain")]), - (bold_reg_wf, outputnode, [ - ("outputnode.itk_bold_to_t1", "bold2anat_xfm"), - ("outputnode.itk_t1_to_bold", "anat2bold_xfm"), - ]), - (bold_reg_wf, bold_t1_trans_wf, [ - ("outputnode.itk_bold_to_t1", "inputnode.itk_bold_to_t1"), - ]), - (bold_final, bold_t1_trans_wf, [ - ("mask", "inputnode.ref_bold_mask"), - ("boldref", "inputnode.ref_bold_brain"), - ]), - (bold_t1_trans_wf, outputnode, [ - ("outputnode.bold_t1", "bold_t1"), - ("outputnode.bold_t1_ref", "bold_t1_ref"), - ("outputnode.bold_aseg_t1", "bold_aseg_t1"), - ("outputnode.bold_aparc_t1", "bold_aparc_t1"), - ]), - # Connect bold_confounds_wf - (inputnode, bold_confounds_wf, [ - ("t1w_tpms", "inputnode.t1w_tpms"), - ("t1w_mask", "inputnode.t1w_mask"), - ]), - (bold_hmc_wf, bold_confounds_wf, [ - ("outputnode.movpar_file", "inputnode.movpar_file"), - ("outputnode.rmsd_file", "inputnode.rmsd_file"), - ]), - (bold_reg_wf, bold_confounds_wf, [ - ("outputnode.itk_t1_to_bold", "inputnode.t1_bold_xform") - ]), - (initial_boldref_wf, bold_confounds_wf, [ - ("outputnode.skip_vols", "inputnode.skip_vols"), - ]), - (initial_boldref_wf, final_boldref_wf, [ - ("outputnode.skip_vols", "inputnode.dummy_scans"), - ]), - (final_boldref_wf, bold_final, [ - ("outputnode.ref_image", "boldref"), - ("outputnode.bold_mask", "mask"), - ]), - (bold_final, bold_confounds_wf, [ - ("bold", "inputnode.bold"), - ("mask", "inputnode.bold_mask"), - ]), - (bold_confounds_wf, outputnode, [ - ("outputnode.confounds_file", "confounds"), - ("outputnode.confounds_metadata", "confounds_metadata"), - ("outputnode.acompcor_masks", "acompcor_masks"), - ("outputnode.tcompcor_mask", "tcompcor_mask"), - ]), - # Native-space BOLD files (if calculated) - (bold_final, outputnode, [ - ("bold", "bold_native"), - ("boldref", "bold_native_ref"), - ("mask", "bold_mask_native"), - ("bold_echos", "bold_echos_native"), - ("t2star", "t2star_bold"), - ]), - # Summary - (initial_boldref_wf, summary, [("outputnode.algo_dummy_scans", "algo_dummy_scans")]), - (bold_reg_wf, summary, [("outputnode.fallback", "fallback")]), - (outputnode, summary, [("confounds", "confounds_file")]), - # Select echo indices for original/validated BOLD files - (echo_index, bold_source, [("echoidx", "index")]), - (echo_index, select_bold, [("echoidx", "index")]), - ]) - # fmt:on - # for standard EPI data, pass along correct file if not multiecho: # fmt:off @@ -1009,12 +1168,6 @@ def _last(inlist): ]) # fmt:on - # Fill-in datasinks of reportlets seen so far - for node in workflow.list_node_names(): - if node.split(".")[-1].startswith("ds_report"): - workflow.get_node(node).inputs.base_directory = fmriprep_dir - workflow.get_node(node).inputs.source_file = ref_file - if not has_fieldmap: # Finalize workflow without SDC connections summary.inputs.distortion_correction = "None" @@ -1296,25 +1449,24 @@ def _create_mem_gb(bold_fname): return bold_tlen, mem_gb -def _get_wf_name(bold_fname): +def _get_wf_name(bold_fname, prefix): """ Derive the workflow name for supplied BOLD file. - >>> _get_wf_name("/completely/made/up/path/sub-01_task-nback_bold.nii.gz") - 'func_preproc_task_nback_wf' - >>> _get_wf_name("/completely/made/up/path/sub-01_task-nback_run-01_echo-1_bold.nii.gz") - 'func_preproc_task_nback_run_01_echo_1_wf' + >>> _get_wf_name("/completely/made/up/path/sub-01_task-nback_bold.nii.gz", "bold") + 'bold_task_nback_wf' + >>> _get_wf_name( + ... "/completely/made/up/path/sub-01_task-nback_run-01_echo-1_bold.nii.gz", + ... "preproc", + ... ) + 'preproc_task_nback_run_01_echo_1_wf' """ from nipype.utils.filemanip import split_filename fname = split_filename(bold_fname)[1] - fname_nosub = "_".join(fname.split("_")[1:]) - name = "func_preproc_" + fname_nosub.replace(".", "_").replace(" ", "").replace( - "-", "_" - ).replace("_bold", "_wf") - - return name + fname_nosub = "_".join(fname.split("_")[1:-1]) + return f'{prefix}_{fname_nosub.replace("-", "_")}_wf' def _to_join(in_file, join_file): @@ -1365,21 +1517,3 @@ def get_img_orientation(imgf): """Return the image orientation as a string""" img = nb.load(imgf) return "".join(nb.aff2axcodes(img.affine)) - - -def get_estimator(layout, fname): - field_source = layout.get_metadata(fname).get("B0FieldSource") - if isinstance(field_source, str): - field_source = (field_source,) - - if field_source is None: - import re - from pathlib import Path - - from sdcflows.fieldmaps import get_identifier - - # Fallback to IntendedFor - intended_rel = re.sub(r"^sub-[a-zA-Z0-9]*/", "", str(Path(fname).relative_to(layout.root))) - field_source = get_identifier(intended_rel) - - return field_source diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 789906edd..36ec19861 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -20,32 +20,30 @@ # # https://www.nipreps.org/community/licensing/ # -""" -Orchestrating the BOLD-preprocessing workflow -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. autofunction:: init_func_preproc_wf -.. autofunction:: init_func_derivatives_wf - -""" import os import typing as ty import bids import nibabel as nb -import numpy as np from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf from niworkflows.interfaces.header import ValidateImage from niworkflows.interfaces.nitransforms import ConcatenateXFMs +from niworkflows.interfaces.utility import KeySelect from niworkflows.utils.connections import listify from sdcflows.workflows.apply.correction import init_unwarp_wf from sdcflows.workflows.apply.registration import init_coeff2epi_wf from ... import config from ...interfaces.reports import FunctionalSummary +from ...interfaces.resampling import ( + DistortionParameters, + ReconstructFieldmap, + ResampleSeries, +) from ...utils.bids import extract_entities +from ...utils.misc import estimate_bold_mem_usage # BOLD workflows from .hmc import init_bold_hmc_wf @@ -57,6 +55,8 @@ ) from .reference import init_raw_boldref_wf from .registration import init_bold_reg_wf +from .stc import init_bold_stc_wf +from .t2s import init_bold_t2s_wf def get_sbrefs( @@ -93,21 +93,113 @@ def get_sbrefs( def init_bold_fit_wf( *, - bold_series: ty.Union[str, ty.List[str]], - precomputed: dict, + bold_series: ty.List[str], + precomputed: dict = {}, fieldmap_id: ty.Optional[str] = None, omp_nthreads: int = 1, + name: str = "bold_fit_wf", ) -> pe.Workflow: + """ + This workflow controls the minimal estimation steps for functional preprocessing. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from fmriprep.workflows.tests import mock_config + from fmriprep import config + from fmriprep.workflows.bold.fit import init_bold_fit_wf + with mock_config(): + bold_file = config.execution.bids_dir / "sub-01" / "func" \ + / "sub-01_task-mixedgamblestask_run-01_bold.nii.gz" + wf = init_bold_fit_wf(bold_series=[str(bold_file)]) + + Parameters + ---------- + bold_series + List of paths to NIfTI files. + precomputed + Dictionary containing precomputed derivatives to reuse, if possible. + fieldmap_id + ID of the fieldmap to use to correct this BOLD series. If :obj:`None`, + no correction will be applied. + + Inputs + ------ + bold_file + BOLD series NIfTI file + t1w_preproc + Bias-corrected structural template image + t1w_mask + Mask of the skull-stripped template image + t1w_dseg + Segmentation of preprocessed structural image, including + gray-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF) + anat2std_xfm + List of transform files, collated with templates + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + fsnative2t1w_xfm + LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w + fmap_id + Unique identifiers to select fieldmap files + fmap + List of estimated fieldmaps (collated with fmap_id) + fmap_ref + List of fieldmap reference files (collated with fmap_id) + fmap_coeff + List of lists of spline coefficient files (collated with fmap_id) + fmap_mask + List of fieldmap masks (collated with fmap_id) + sdc_method + List of fieldmap correction method names (collated with fmap_id) + + Outputs + ------- + hmc_boldref + BOLD reference image used for head motion correction. + Minimally processed to ensure consistent contrast with BOLD series. + coreg_boldref + BOLD reference image used for coregistration. Contrast-enhanced + and fieldmap-corrected for greater anatomical fidelity, and aligned + with ``hmc_boldref``. + bold_mask + Mask of ``coreg_boldref``. + motion_xfm + Affine transforms from each BOLD volume to ``hmc_boldref``, written + as concatenated ITK affine transforms. + boldref2anat_xfm + Affine transform mapping from BOLD reference space to the anatomical + space. + boldref2fmap_xfm + Affine transform mapping from BOLD reference space to the fieldmap + space, if applicable. + + See Also + -------- + + * :py:func:`~fmriprep.workflows.bold.reference.init_raw_boldref_wf` + * :py:func:`~fmriprep.workflows.bold.hmc.init_bold_hmc_wf` + * :py:func:`~niworkflows.func.utils.init_enhance_and_skullstrip_bold_wf` + * :py:func:`~sdcflows.workflows.apply.registration.init_coeff2epi_wf` + * :py:func:`~sdcflows.workflows.apply.correction.init_unwarp_wf` + * :py:func:`~fmriprep.workflows.bold.registration.init_bold_reg_wf` + * :py:func:`~fmriprep.workflows.bold.outputs.init_ds_boldref_wf` + * :py:func:`~fmriprep.workflows.bold.outputs.init_ds_hmc_wf` + * :py:func:`~fmriprep.workflows.bold.outputs.init_ds_registration_wf` + + """ from niworkflows.engine.workflows import LiterateWorkflow as Workflow - from niworkflows.interfaces.utility import KeySelect + + from fmriprep.utils.misc import estimate_bold_mem_usage layout = config.execution.layout # Collect bold and sbref files, sorted by EchoTime - bold_files = sorted( - listify(bold_series), - key=lambda fname: layout.get_metadata(fname).get("EchoTime"), - ) + bold_files = sorted(bold_series, key=lambda fname: layout.get_metadata(fname).get("EchoTime")) sbref_files = get_sbrefs( bold_files, entity_overrides=config.execution.get().get('bids_filters', {}).get('sbref', {}), @@ -123,8 +215,7 @@ def init_bold_fit_wf( metadata = layout.get_metadata(bold_file) orientation = "".join(nb.aff2axcodes(nb.load(bold_file).affine)) - if os.path.isfile(bold_file): - bold_tlen, mem_gb = _create_mem_gb(bold_file) + bold_tlen, mem_gb = estimate_bold_mem_usage(bold_file) # Boolean used to update workflow self-descriptions multiecho = len(bold_files) > 1 @@ -140,7 +231,7 @@ def init_bold_fit_wf( boldref2fmap_xform = transforms.get("boldref2fmap") boldref2anat_xform = transforms.get("boldref2anat") - workflow = Workflow(name=_get_wf_name(bold_file)) + workflow = Workflow(name=name) inputnode = pe.Node( niu.IdentityInterface( @@ -169,11 +260,13 @@ def init_bold_fit_wf( outputnode = pe.Node( niu.IdentityInterface( fields=[ + "dummy_scans", "hmc_boldref", "coreg_boldref", "bold_mask", "motion_xfm", "boldref2anat_xfm", + "boldref2fmap_xfm", ], ), name="outputnode", @@ -188,7 +281,7 @@ def init_bold_fit_wf( fmapref_buffer = pe.Node(niu.Function(function=_select_ref), name="fmapref_buffer") hmc_buffer = pe.Node(niu.IdentityInterface(fields=["hmc_xforms"]), name="hmc_buffer") fmapreg_buffer = pe.Node( - niu.IdentityInterface(fields=["boldref2fmap_xform"]), name="fmapreg_buffer" + niu.IdentityInterface(fields=["boldref2fmap_xfm"]), name="fmapreg_buffer" ) regref_buffer = pe.Node( niu.IdentityInterface(fields=["boldref", "boldmask"]), name="regref_buffer" @@ -221,8 +314,11 @@ def init_bold_fit_wf( # fmt:off workflow.connect([ (hmcref_buffer, outputnode, [("boldref", "hmc_boldref")]), - (regref_buffer, outputnode, [("boldref", "coreg_boldref"), - ("boldmask", "bold_mask")]), + (regref_buffer, outputnode, [ + ("boldref", "coreg_boldref"), + ("boldmask", "bold_mask"), + ]), + (fmapreg_buffer, outputnode, [("boldref2fmap_xfm", "boldref2fmap_xfm")]), (hmc_buffer, outputnode, [("hmc_xforms", "motion_xfm")]), (inputnode, func_fit_reports_wf, [ ("bold_file", "inputnode.source_file"), @@ -261,9 +357,11 @@ def init_bold_fit_wf( # fmt:off workflow.connect([ - (hmc_boldref_wf, hmcref_buffer, [("outputnode.bold_file", "bold_file")]), - (hmc_boldref_wf, ds_hmc_boldref_wf, [("outputnode.boldref", "inputnode.boldref")]), - (ds_hmc_boldref_wf, hmcref_buffer, [("outputnode.boldref", "boldref")]), + (hmc_boldref_wf, hmcref_buffer, [ + ("outputnode.bold_file", "bold_file"), + ("outputnode.boldref", "boldref"), + ]), + (hmcref_buffer, ds_hmc_boldref_wf, [("boldref", "inputnode.boldref")]), (hmc_boldref_wf, summary, [("outputnode.algo_dummy_scans", "algo_dummy_scans")]), (hmc_boldref_wf, func_fit_reports_wf, [ ("outputnode.validation_report", "inputnode.validation_report"), @@ -341,7 +439,7 @@ def init_bold_fit_wf( if fieldmap_id: fmap_select = pe.Node( KeySelect( - fields=["fmap", "fmap_ref", "fmap_coeff", "fmap_mask", "sdc_method"], + fields=["fmap_ref", "fmap_coeff", "fmap_mask", "sdc_method"], key=fieldmap_id, ), name="fmap_select", @@ -395,7 +493,6 @@ def init_bold_fit_wf( # fmt:off workflow.connect([ (inputnode, fmap_select, [ - ("fmap", "fmap"), ("fmap_ref", "fmap_ref"), ("fmap_coeff", "fmap_coeff"), ("fmap_mask", "fmap_mask"), @@ -492,50 +589,291 @@ def init_bold_fit_wf( return workflow -def _create_mem_gb(bold_fname): - img = nb.load(bold_fname) - nvox = int(np.prod(img.shape, dtype='u8')) - # Assume tools will coerce to 8-byte floats to be safe - bold_size_gb = 8 * nvox / (1024**3) - bold_tlen = img.shape[-1] - mem_gb = { - "filesize": bold_size_gb, - "resampled": bold_size_gb * 4, - "largemem": bold_size_gb * (max(bold_tlen / 100, 1.0) + 4), - } - - return bold_tlen, mem_gb +def init_bold_native_wf( + *, + bold_series: ty.List[str], + fieldmap_id: ty.Optional[str] = None, + name: str = "bold_native_wf", +) -> pe.Workflow: + r""" + Minimal resampling workflow. + + This workflow performs slice-timing correction, and resamples to boldref space + with head motion and susceptibility distortion correction. It also handles + multi-echo processing and selects the transforms needed to perform further + resampling. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from fmriprep.workflows.tests import mock_config + from fmriprep import config + from fmriprep.workflows.bold.fit import init_bold_native_wf + with mock_config(): + bold_file = config.execution.bids_dir / "sub-01" / "func" \ + / "sub-01_task-mixedgamblestask_run-01_bold.nii.gz" + wf = init_bold_native_wf(bold_series=[str(bold_file)]) + Parameters + ---------- + bold_series + List of paths to NIfTI files. + fieldmap_id + ID of the fieldmap to use to correct this BOLD series. If :obj:`None`, + no correction will be applied. + + Inputs + ------ + boldref + BOLD reference file + bold_mask + Mask of BOLD reference file + motion_xfm + Affine transforms from each BOLD volume to ``hmc_boldref``, written + as concatenated ITK affine transforms. + fmapreg_xfm + Affine transform mapping from BOLD reference space to the fieldmap + space, if applicable. + fmap_id + Unique identifiers to select fieldmap files + fmap_ref + List of fieldmap reference files (collated with fmap_id) + fmap_coeff + List of lists of spline coefficient files (collated with fmap_id) + + Outputs + ------- + bold_minimal + BOLD series ready for further resampling. For single-echo data, only + slice-timing correction (STC) may have been applied. For multi-echo + data, this is identical to bold_native. + bold_native + BOLD series resampled into BOLD reference space. Slice-timing, + head motion and susceptibility distortion correction (STC, HMC, SDC) + will all be applied to each file. For multi-echo data, the echos + are combined to form an `optimal combination`_. + motion_xfm + Motion correction transforms for further correcting bold_minimal. + For multi-echo data, motion correction has already been applied, so + this will be undefined. + fieldmap_id + Fieldmap ID for further correcting bold_minimal. For multi-echo data, + susceptibility distortion correction has already been applied, so + this will be undefined. + bold_echos + The individual, corrected echos, suitable for use in Tedana. + (Multi-echo only.) + t2star_map + The T2\* map estimated by Tedana when calculating the optimal combination. + (Multi-echo only.) + + See Also + -------- + + * :py:func:`~fmriprep.workflows.bold.stc.init_bold_stc_wf` + * :py:func:`~fmriprep.workflows.bold.t2s.init_bold_t2s_wf` + + .. _optimal combination: https://tedana.readthedocs.io/en/stable/approach.html#optimal-combination -def _get_wf_name(bold_fname): """ - Derive the workflow name for supplied BOLD file. - >>> _get_wf_name("/completely/made/up/path/sub-01_task-nback_bold.nii.gz") - 'func_preproc_task_nback_wf' - >>> _get_wf_name("/completely/made/up/path/sub-01_task-nback_run-01_echo-1_bold.nii.gz") - 'func_preproc_task_nback_run_01_echo_1_wf' + layout = config.execution.layout - """ - from nipype.utils.filemanip import split_filename + # Shortest echo first + all_metadata, bold_files, echo_times = zip( + *sorted( + ( + (md := layout.get_metadata(bold_file), bold_file, md.get("EchoTime")) + for bold_file in listify(bold_series) + ), + key=lambda x: x[2], + ) + ) + multiecho = len(bold_files) > 1 + + bold_file = bold_files[0] + metadata = all_metadata[0] + + bold_tlen, mem_gb = estimate_bold_mem_usage(bold_file) + + if multiecho: + shapes = [nb.load(echo).shape for echo in bold_files] + if len(set(shapes)) != 1: + diagnostic = "\n".join( + f"{os.path.basename(echo)}: {shape}" for echo, shape in zip(bold_files, shapes) + ) + raise RuntimeError(f"Multi-echo images found with mismatching shapes\n{diagnostic}") + if len(shapes) == 2: + raise RuntimeError( + "Multi-echo processing requires at least three different echos (found two)." + ) + + run_stc = bool(metadata.get("SliceTiming")) and "slicetiming" not in config.workflow.ignore + + workflow = pe.Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + # BOLD fit + "boldref", + "bold_mask", + "motion_xfm", + "fmapreg_xfm", + "dummy_scans", + # Fieldmap fit + "fmap_ref", + "fmap_coeff", + "fmap_id", + ], + ), + name='inputnode', + ) + + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "bold_minimal", # Single echo: STC; Multi-echo: optimal combination + "bold_native", # STC + HMC + SDC; Multi-echo: optimal combination + "motion_xfm", # motion_xfms to apply to bold_minimal (none for ME) + "fieldmap_id", # fieldmap to apply to bold_minimal (none for ME) + # Multiecho outputs + "bold_echos", # Individual corrected echos + "t2star_map", # T2* map + ], # fmt:skip + ), + name="outputnode", + ) + + boldbuffer = pe.Node( + niu.IdentityInterface(fields=["bold_file", "ro_time", "pe_dir"]), name="boldbuffer" + ) + + # Track echo index - this allows us to treat multi- and single-echo workflows + # almost identically + echo_index = pe.Node(niu.IdentityInterface(fields=["echoidx"]), name="echo_index") + if multiecho: + echo_index.iterables = [("echoidx", range(len(bold_files)))] + else: + echo_index.inputs.echoidx = 0 + + # BOLD source: track original BOLD file(s) + bold_source = pe.Node(niu.Select(inlist=bold_files), name="bold_source") + validate_bold = pe.Node(ValidateImage(), name="validate_bold") + workflow.connect([ + (echo_index, bold_source, [("echoidx", "index")]), + (bold_source, validate_bold, [("out", "in_file")]), + ]) # fmt:skip + + # Slice-timing correction + if run_stc: + bold_stc_wf = init_bold_stc_wf(name="bold_stc_wf", metadata=metadata) + workflow.connect([ + (inputnode, bold_stc_wf, [("dummy_scans", "inputnode.skip_vols")]), + (validate_bold, bold_stc_wf, [("out_file", "inputnode.bold_file")]), + (bold_stc_wf, boldbuffer, [("outputnode.stc_file", "bold_file")]), + ]) # fmt:skip + else: + workflow.connect([(validate_bold, boldbuffer, [("out_file", "bold_file")])]) + + # Prepare fieldmap metadata + if fieldmap_id: + fmap_select = pe.Node( + KeySelect(fields=["fmap_ref", "fmap_coeff"], key=fieldmap_id), + name="fmap_select", + run_without_submitting=True, + ) - fname = split_filename(bold_fname)[1] - fname_nosub = "_".join(fname.split("_")[1:]) - name = "func_preproc_" + fname_nosub.replace(".", "_").replace(" ", "").replace( - "-", "_" - ).replace("_bold", "_wf") + distortion_params = pe.Node( + DistortionParameters(metadata=metadata, in_file=bold_file), + name="distortion_params", + run_without_submitting=True, + ) + workflow.connect([ + (inputnode, fmap_select, [ + ("fmap_ref", "fmap_ref"), + ("fmap_coeff", "fmap_coeff"), + ("fmap_id", "keys"), + ]), + (distortion_params, boldbuffer, [ + ("readout_time", "ro_time"), + ("pe_direction", "pe_dir"), + ]), + ]) # fmt:skip + + # Resample to boldref + boldref_bold = pe.Node( + ResampleSeries(), name="boldref_bold", n_procs=config.nipype.omp_nthreads + ) + + workflow.connect([ + (inputnode, boldref_bold, [ + ("boldref", "ref_file"), + ("motion_xfm", "transforms"), + ]), + (boldbuffer, boldref_bold, [ + ("bold_file", "in_file"), + ("ro_time", "ro_time"), + ("pe_dir", "pe_dir"), + ]), + ]) # fmt:skip + + if fieldmap_id: + boldref_fmap = pe.Node(ReconstructFieldmap(inverse=[True]), name="boldref_fmap") + workflow.connect([ + (inputnode, boldref_fmap, [ + ("boldref", "target_ref_file"), + ("fmapreg_xfm", "transforms"), + ]), + (fmap_select, boldref_fmap, [ + ("fmap_coeff", "in_coeffs"), + ("fmap_ref", "fmap_ref_file"), + ]), + (boldref_fmap, boldref_bold, [("out_file", "fieldmap")]), + ]) # fmt:skip + + if multiecho: + join_echos = pe.JoinNode( + niu.IdentityInterface(fields=["bold_files"]), + joinsource="echo_index", + joinfield=["bold_files"], + name="join_echos", + ) - return name + # create optimal combination, adaptive T2* map + bold_t2s_wf = init_bold_t2s_wf( + echo_times=echo_times, + mem_gb=mem_gb["filesize"], + omp_nthreads=config.nipype.omp_nthreads, + name="bold_t2smap_wf", + ) + # Do NOT set motion_xfm or fieldmap_id on outputnode + # This prevents downstream resamplers from double-dipping + workflow.connect([ + (inputnode, bold_t2s_wf, [("bold_mask", "inputnode.bold_mask")]), + (boldref_bold, join_echos, [("out_file", "bold_files")]), + (join_echos, bold_t2s_wf, [("bold_files", "inputnode.bold_file")]), + (join_echos, outputnode, [("bold_files", "bold_echos")]), + (bold_t2s_wf, outputnode, [ + ("outputnode.bold", "bold_minimal"), + ("outputnode.bold", "bold_native"), + ("outputnode.t2star_map", "t2star_map"), + ]), + ]) # fmt:skip + else: + workflow.connect([ + (inputnode, outputnode, [("motion_xfm", "motion_xfm")]), + (boldbuffer, outputnode, [("bold_file", "bold_minimal")]), + (boldref_bold, outputnode, [("out_file", "bold_native")]), + ]) # fmt:skip -def _to_join(in_file, join_file): - """Join two tsv files if the join_file is not ``None``.""" - from niworkflows.interfaces.utility import JoinTSVColumns + if fieldmap_id: + outputnode.inputs.fieldmap_id = fieldmap_id - if join_file is None: - return in_file - res = JoinTSVColumns(in_file=in_file, join_file=join_file).run() - return res.outputs.out_file + return workflow def _select_ref(sbref_files, boldref_files): diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index 1d10f1aa8..b893abf50 100644 --- a/fmriprep/workflows/bold/outputs.py +++ b/fmriprep/workflows/bold/outputs.py @@ -556,6 +556,133 @@ def init_ds_hmc_wf( return workflow +def init_ds_bold_native_wf( + *, + bids_root: str, + output_dir: str, + multiecho: bool, + bold_output: bool, + echo_output: bool, + all_metadata: ty.List[dict], + name="ds_bold_native_wf", +) -> pe.Workflow: + metadata = all_metadata[0] + timing_parameters = prepare_timing_parameters(metadata) + + workflow = pe.Workflow(name=name) + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + 'source_files', + 'bold', + 'bold_mask', + 'bold_echos', + 't2star', + ] + ), + name='inputnode', + ) + + raw_sources = pe.Node(niu.Function(function=_bids_relative), name='raw_sources') + raw_sources.inputs.bids_root = bids_root + workflow.connect(inputnode, 'source_files', raw_sources, 'in_files') + + if bold_output: + 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, + ) + workflow.connect([ + (inputnode, ds_bold, [ + ('source_files', 'source_file'), + ('bold', 'in_file'), + ]), + ]) # fmt:skip + + if bold_output or echo_output: + ds_bold_mask = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + desc='brain', + suffix='mask', + compress=True, + dismiss_entities=("echo",), + ), + name='ds_bold_mask', + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + workflow.connect([ + (inputnode, ds_bold_mask, [ + ('source_files', 'source_file'), + ('bold_mask', 'in_file'), + ]), + (raw_sources, ds_bold_mask, [('out', 'RawSources')]), + ]) # fmt:skip + + if bold_output and multiecho: + t2star_meta = { + 'Units': 's', + 'EstimationReference': 'doi:10.1002/mrm.20900', + 'EstimationAlgorithm': 'monoexponential decay model', + } + ds_t2star = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + space='boldref', + suffix='T2starmap', + compress=True, + dismiss_entities=("echo",), + **t2star_meta, + ), + name='ds_t2star_bold', + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + workflow.connect([ + (inputnode, ds_t2star, [ + ('source_files', 'source_file'), + ('t2star', 'in_file'), + ]), + (raw_sources, ds_t2star, [('out', 'RawSources')]), + ]) # fmt:skip + + if echo_output: + ds_bold_echos = pe.MapNode( + DerivativesDataSink( + base_directory=output_dir, + desc='preproc', + compress=True, + SkullStripped=False, + TaskName=metadata.get('TaskName'), + **timing_parameters, + ), + iterfield=['source_file', 'in_file', 'meta_dict'], + name='ds_bold_echos', + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + ds_bold_echos.inputs.meta_dict = [{"EchoTime": md["EchoTime"]} for md in all_metadata] + workflow.connect([ + (inputnode, ds_bold_echos, [ + ('source_files', 'source_file'), + ('bold_echos', 'in_file'), + ]), + ]) # fmt:skip + + return workflow + + def init_func_derivatives_wf( bids_root: str, cifti_output: bool, @@ -750,98 +877,6 @@ def init_func_derivatives_wf( ]) # fmt:on - bold_output = nonstd_spaces.intersection(('func', 'run', 'bold', 'boldref', 'sbref')) - if bold_output: - ds_bold_native = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - desc='preproc', - compress=True, - SkullStripped=masked, - TaskName=metadata.get('TaskName'), - **timing_parameters, - ), - name='ds_bold_native', - run_without_submitting=True, - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - # fmt: off - workflow.connect([ - (inputnode, ds_bold_native, [('source_file', 'source_file'), - ('bold_native', 'in_file')]), - ]) - # fmt:on - - # Save masks and boldref if we're going to save either orig BOLD series or echos - if bold_output or multiecho and config.execution.me_output_echos: - ds_bold_mask_native = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - desc='brain', - suffix='mask', - compress=True, - dismiss_entities=("echo",), - ), - name='ds_bold_mask_native', - run_without_submitting=True, - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - # fmt:off - workflow.connect([ - (inputnode, ds_bold_mask_native, [('source_file', 'source_file'), - ('bold_mask_native', 'in_file')]), - (raw_sources, ds_bold_mask_native, [('out', 'RawSources')]), - ]) - # fmt:on - - if multiecho: - ds_t2star_bold = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - space='boldref', - suffix='T2starmap', - compress=True, - dismiss_entities=("echo",), - **t2star_meta, - ), - name='ds_t2star_bold', - run_without_submitting=True, - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - # fmt:off - workflow.connect([ - (inputnode, ds_t2star_bold, [('source_file', 'source_file'), - ('t2star_bold', 'in_file')]), - (raw_sources, ds_t2star_bold, [('out', 'RawSources')]), - ]) - # fmt:on - - if multiecho and config.execution.me_output_echos: - ds_bold_echos_native = pe.MapNode( - DerivativesDataSink( - base_directory=output_dir, - desc='preproc', - compress=True, - SkullStripped=False, - TaskName=metadata.get('TaskName'), - **timing_parameters, - ), - iterfield=['source_file', 'in_file', 'meta_dict'], - name='ds_bold_echos_native', - run_without_submitting=True, - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - ds_bold_echos_native.inputs.meta_dict = [ - {"EchoTime": md["EchoTime"]} for md in all_metadata - ] - # fmt:off - workflow.connect([ - (inputnode, ds_bold_echos_native, [ - ('all_source_files', 'source_file'), - ('bold_echos_native', 'in_file')]), - ]) - # fmt:on - # Resample to T1w space if nonstd_spaces.intersection(('T1w', 'anat')): ds_bold_t1 = pe.Node( diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py index da282d55f..1649cd0d2 100644 --- a/fmriprep/workflows/bold/resampling.py +++ b/fmriprep/workflows/bold/resampling.py @@ -194,7 +194,7 @@ def select_target(subject_id, space): ("source_file", "reference"), ("fsnative2t1w_xfm", "in_xfms"), ]), - (get_fsnative, itk2lta, [("T1", "dst_file")]), + (get_fsnative, itk2lta, [("T1", "moving")]), (inputnode, sampler, [ ("subjects_dir", "subjects_dir"), ("subject_id", "subject_id"), @@ -663,7 +663,7 @@ def init_bold_fsLR_resampling_wf( # Line 85 thru ... volume_to_surface = pe.Node( VolumeToSurfaceMapping( - method="ribbon-constrained", + method="ribbon-constrained", output_weights_text="output_weights.txt", ), name="volume_to_surface", diff --git a/fmriprep/workflows/bold/t2s.py b/fmriprep/workflows/bold/t2s.py index f0b087716..b00c81c95 100644 --- a/fmriprep/workflows/bold/t2s.py +++ b/fmriprep/workflows/bold/t2s.py @@ -120,15 +120,13 @@ def init_bold_t2s_wf( name='t2smap_node', mem_gb=2.5 * mem_gb * len(echo_times), ) - # fmt:off workflow.connect([ (inputnode, dilate_mask, [('bold_mask', 'in_mask')]), (inputnode, t2smap_node, [('bold_file', 'in_files')]), (dilate_mask, t2smap_node, [('out_mask', 'mask_file')]), (t2smap_node, outputnode, [('optimal_comb', 'bold'), ('t2star_map', 't2star_map')]), - ]) - # fmt:on + ]) # fmt:skip return workflow @@ -157,9 +155,10 @@ def init_t2s_reporting_wf(name: str = 't2s_reporting_wf'): reference BOLD file label_file an integer label file identifying gray matter with value ``1`` - label_bold_xform - Affine matrix that maps the label file into alignment with the native - BOLD space; can be ``"identity"`` if label file is already aligned + boldref2anat_xfm + Affine matrix that maps images in the native bold space into the + anatomical space of ``label_file``; can be ``"identity"`` if label + file is already aligned Outputs ------- @@ -177,7 +176,7 @@ def init_t2s_reporting_wf(name: str = 't2s_reporting_wf'): workflow = pe.Workflow(name=name) inputnode = pe.Node( - niu.IdentityInterface(fields=['t2star_file', 'boldref', 'label_file', 'label_bold_xform']), + niu.IdentityInterface(fields=['t2star_file', 'boldref', 'label_file', 'boldref2anat_xfm']), name='inputnode', ) @@ -185,7 +184,10 @@ def init_t2s_reporting_wf(name: str = 't2s_reporting_wf'): niu.IdentityInterface(fields=['t2star_hist', 't2s_comp_report']), name='outputnode' ) - label_tfm = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="label_tfm") + label_tfm = pe.Node( + ApplyTransforms(interpolation="MultiLabel", invert_transform_flags=[True]), + name="label_tfm", + ) gm_mask = pe.Node(Label2Mask(label_val=1), name="gm_mask") @@ -204,11 +206,10 @@ def init_t2s_reporting_wf(name: str = 't2s_reporting_wf'): name="t2s_comparison", mem_gb=0.1, ) - # fmt:off workflow.connect([ (inputnode, label_tfm, [('label_file', 'input_image'), ('t2star_file', 'reference_image'), - ('label_bold_xform', 'transforms')]), + ('boldref2anat_xfm', 'transforms')]), (inputnode, clip_t2star, [('t2star_file', 'in_file')]), (clip_t2star, t2s_hist, [('out_file', 'in_file')]), (label_tfm, gm_mask, [('output_image', 'in_file')]), @@ -218,6 +219,5 @@ def init_t2s_reporting_wf(name: str = 't2s_reporting_wf'): (gm_mask, t2s_comparison, [('out_file', 'wm_seg')]), (t2s_hist, outputnode, [('out_report', 't2star_hist')]), (t2s_comparison, outputnode, [('out_report', 't2s_comp_report')]), - ]) - # fmt:on + ]) # fmt:skip return workflow diff --git a/fmriprep/workflows/tests.py b/fmriprep/workflows/tests/__init__.py similarity index 97% rename from fmriprep/workflows/tests.py rename to fmriprep/workflows/tests/__init__.py index c173c2209..2861e70da 100644 --- a/fmriprep/workflows/tests.py +++ b/fmriprep/workflows/tests/__init__.py @@ -30,11 +30,13 @@ from pkg_resources import resource_filename as pkgrf from toml import loads +from ..base import get_estimator + @contextmanager def mock_config(): """Create a mock config for documentation and testing purposes.""" - from .. import config + from ... import config _old_fs = os.getenv('FREESURFER_HOME') if not _old_fs: diff --git a/fmriprep/workflows/bold/tests/test_base.py b/fmriprep/workflows/tests/test_base.py similarity index 100% rename from fmriprep/workflows/bold/tests/test_base.py rename to fmriprep/workflows/tests/test_base.py