diff --git a/.circleci/config.yml b/.circleci/config.yml index bc2689be9..456fc3758 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -755,6 +755,7 @@ jobs: --fs-no-reconall --use-syn-sdc --ignore slicetiming \ --dummy-scans 1 --sloppy --write-graph \ --output-spaces MNI152NLin2009cAsym \ + --thermal-denoise-method mppca \ --mem-mb 14336 --nthreads 4 -vv - run: *check_outputs - run: diff --git a/docs/workflows.rst b/docs/workflows.rst index 135a5ab76..e9caeee9f 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -391,6 +391,32 @@ For a more accurate estimation of head-motion, we calculate its parameters before any time-domain filtering (i.e., :ref:`slice-timing correction `), as recommended in [Power2017]_. + +Thermal noise removal +~~~~~~~~~~~~~~~~~~~~~ +:py:func:`~fmriprep.workflows.bold.denoise.init_bold_dwidenoise_wf` + +.. workflow:: + :graph2use: colored + :simple_form: yes + + from fmriprep.workflows.bold.denoise import init_bold_dwidenoise_wf + + wf = init_bold_dwidenoise_wf( + has_phase=True, + has_norf=True, + mem_gb=1, + ) + +Functional MRI exhibits low signal-to-noise, which is exacerbated by thermal noise, +especially at higher field strengths. +Thermal noise removal with the ``dwidenoise`` tool from MRtrix3 can be enabled +with the ``--thermal-denoise-method`` parameter. +fMRIPrep will automatically leverage phase data, when available, to improve the +denoising process. +This can be disabled with the ``--ignore phase`` command line argument. + + .. _bold_stc: Slice time correction diff --git a/env.yml b/env.yml index d680f3a68..0be66d837 100644 --- a/env.yml +++ b/env.yml @@ -2,6 +2,7 @@ name: fmriprep channels: - https://fsl.fmrib.ox.ac.uk/fsldownloads/fslconda/public/ - conda-forge + - mrtrix3 # Update this ~yearly; last updated Jan 2024 dependencies: - python=3.11 @@ -37,6 +38,8 @@ dependencies: - fsl-mcflirt=2111.0 - fsl-miscmaths=2203.2 - fsl-topup=2203.5 + # Workflow dependencies: mrtrix3 + - mrtrix3=3.0.4 - pip - pip: - -r requirements.txt diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py index a734c4547..6595f0b85 100644 --- a/fmriprep/cli/parser.py +++ b/fmriprep/cli/parser.py @@ -337,7 +337,16 @@ def _slice_time_ref(value, parser): action='store', nargs='+', default=[], - choices=['fieldmaps', 'slicetiming', 'sbref', 't2w', 'flair', 'fmap-jacobian'], + choices=[ + 'fieldmaps', + 'slicetiming', + 'sbref', + 't2w', + 'flair', + 'fmap-jacobian', + 'phase', + 'norf', + ], help='Ignore selected aspects of the input dataset to disable corresponding ' 'parts of the workflow (a space delimited list)', ) @@ -447,6 +456,14 @@ def _slice_time_ref(value, parser): 'It is faster and less memory intensive, but may be less accurate.' ), ) + g_conf.add_argument( + '--thermal-denoise-method', + action='store', + dest='thermal_denoise_method', + default=None, + choices=['mppca'], + help='Apply MP-PCA denoising to the BOLD data to remove thermal noise', + ) g_outputs = parser.add_argument_group('Options for modulating outputs') g_outputs.add_argument( diff --git a/fmriprep/config.py b/fmriprep/config.py index 6ba0ae1c5..b12f4b514 100644 --- a/fmriprep/config.py +++ b/fmriprep/config.py @@ -624,6 +624,8 @@ class workflow(_Config): in the absence of any alternatives.""" me_t2s_fit_method = 'curvefit' """The method by which to estimate T2*/S0 for multi-echo data""" + thermal_denoise_method = None + """Apply NORDIC or MP-PCA denoising to the BOLD data to remove thermal noise.""" class loggers: diff --git a/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib index 4cc802f25..23280fe6d 100644 --- a/fmriprep/data/boilerplate.bib +++ b/fmriprep/data/boilerplate.bib @@ -365,3 +365,52 @@ @article{patriat_improved_2017 keywords = {Motion, Correction, Methods, Rs-fMRI}, pages = {74--82}, } + +@article{cordero2019complex, + title={Complex diffusion-weighted image estimation via matrix recovery under general noise models}, + author={Cordero-Grande, Lucilio and Christiaens, Daan and Hutter, Jana and Price, Anthony N and Hajnal, Jo V}, + journal={Neuroimage}, + volume={200}, + pages={391--404}, + year={2019}, + publisher={Elsevier}, + url={https://doi.org/10.1016/j.neuroimage.2019.06.039}, + doi={10.1016/j.neuroimage.2019.06.039} +} + +@article{tournier2019mrtrix3, + title={MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation}, + author={Tournier, J-Donald and Smith, Robert and Raffelt, David and Tabbara, Rami and Dhollander, Thijs and Pietsch, Maximilian and Christiaens, Daan and Jeurissen, Ben and Yeh, Chun-Hung and Connelly, Alan}, + journal={Neuroimage}, + volume={202}, + pages={116137}, + year={2019}, + publisher={Elsevier}, + url={https://doi.org/10.1016/j.neuroimage.2019.116137}, + doi={10.1016/j.neuroimage.2019.116137} +} + +@article{vizioli2021lowering, + title={Lowering the thermal noise barrier in functional brain mapping with magnetic resonance imaging}, + author={Vizioli, Luca and Moeller, Steen and Dowdle, Logan and Ak{\c{c}}akaya, Mehmet and De Martino, Federico and Yacoub, Essa and U{\u{g}}urbil, Kamil}, + journal={Nature communications}, + volume={12}, + number={1}, + pages={5181}, + year={2021}, + publisher={Nature Publishing Group UK London}, + url={https://doi.org/10.1038/s41467-021-25431-8}, + doi={10.1038/s41467-021-25431-8} +} + +@article{moeller2021noise, + title={NOise reduction with DIstribution Corrected (NORDIC) PCA in dMRI with complex-valued parameter-free locally low-rank processing}, + author={Moeller, Steen and Pisharady, Pramod Kumar and Ramanna, Sudhir and Lenglet, Christophe and Wu, Xiaoping and Dowdle, Logan and Yacoub, Essa and U{\u{g}}urbil, Kamil and Ak{\c{c}}akaya, Mehmet}, + journal={Neuroimage}, + volume={226}, + pages={117539}, + year={2021}, + publisher={Elsevier}, + url={https://doi.org/10.1016/j.neuroimage.2020.117539}, + doi={10.1016/j.neuroimage.2020.117539} +} diff --git a/fmriprep/data/tests/config.toml b/fmriprep/data/tests/config.toml index b1b7e31b6..358b77809 100644 --- a/fmriprep/data/tests/config.toml +++ b/fmriprep/data/tests/config.toml @@ -47,6 +47,7 @@ run_reconall = true skull_strip_fixed_seed = false skull_strip_template = "OASIS30ANTs" t2s_coreg = false +thermal_denoise_method = "mppca" use_aroma = false [nipype] diff --git a/fmriprep/interfaces/denoise.py b/fmriprep/interfaces/denoise.py new file mode 100644 index 000000000..12969bb08 --- /dev/null +++ b/fmriprep/interfaces/denoise.py @@ -0,0 +1,397 @@ +"""Denoising-related interfaces.""" + +import nibabel as nb +import numpy as np +from nilearn.image import load_img +from nipype.interfaces.base import ( + BaseInterfaceInputSpec, + File, + SimpleInterface, + TraitedSpec, + isdefined, + traits, +) +from nipype.interfaces.mrtrix3.base import MRTrix3Base, MRTrix3BaseInputSpec +from nipype.utils.filemanip import fname_presuffix + + +class _ValidateComplexInputSpec(BaseInterfaceInputSpec): + mag_file = File( + exists=True, + mandatory=True, + desc='Magnitude BOLD EPI', + ) + phase_file = File( + exists=True, + mandatory=False, + desc='Phase BOLD EPI', + ) + + +class _ValidateComplexOutputSpec(TraitedSpec): + mag_file = File(exists=True, desc='Validated magnitude file') + phase_file = File(exists=True, desc='Validated phase file') + + +class ValidateComplex(SimpleInterface): + """Validate complex-valued BOLD data.""" + + input_spec = _ValidateComplexInputSpec + output_spec = _ValidateComplexOutputSpec + + def _run_interface(self, runtime): + import nibabel as nb + + if not isdefined(self.inputs.phase_file): + self._results['mag_file'] = self.inputs.mag_file + return runtime + + mag_img = nb.load(self.inputs.mag_file) + phase_img = nb.load(self.inputs.phase_file) + n_mag_vols = mag_img.shape[3] + n_phase_vols = phase_img.shape[3] + + if n_mag_vols != n_phase_vols: + raise ValueError( + f'Number of volumes in magnitude file ({n_mag_vols}) ' + f'!= number of volumes in phase file ({n_phase_vols})' + ) + + self._results['mag_file'] = self.inputs.mag_file + self._results['phase_file'] = self.inputs.phase_file + + return runtime + + +class DWIDenoiseInputSpec(MRTrix3BaseInputSpec): + in_file = File(exists=True, argstr='%s', position=-2, mandatory=True, desc='input DWI image') + estimator = traits.Enum( + 'Exp2', + argstr='-estimator %s', + desc='noise estimator to use. (default = Exp2)', + ) + mask = File(exists=True, argstr='-mask %s', position=1, desc='mask image') + extent = traits.Tuple( + (traits.Int, traits.Int, traits.Int), + argstr='-extent %d,%d,%d', + desc='set the window size of the denoising filter. (default = 5,5,5)', + ) + noise_image = File( + argstr='-noise %s', + name_template='%s_noise.nii.gz', + name_source=['in_file'], + keep_extension=False, + desc='the output noise map', + ) + out_file = File( + name_template='%s_denoised.nii.gz', + name_source=['in_file'], + keep_extension=False, + argstr='%s', + position=-1, + desc='the output denoised DWI image', + ) + + +class DWIDenoiseOutputSpec(TraitedSpec): + noise_image = File(desc='the output noise map', exists=True) + out_file = File(desc='the output denoised DWI image', exists=True) + + +class DWIDenoise(MRTrix3Base): + """ + Denoise DWI data and estimate the noise level based on the optimal + threshold for PCA. + + DWI data denoising and noise map estimation by exploiting data redundancy + in the PCA domain using the prior knowledge that the eigenspectrum of + random covariance matrices is described by the universal Marchenko Pastur + distribution. + + Important note: image denoising must be performed as the first step of the + image processing pipeline. The routine will fail if interpolation or + smoothing has been applied to the data prior to denoising. + + Note that this function does not correct for non-Gaussian noise biases. + + For more information, see + + + Notes + ----- + There is a -rank option to output a map of the degrees of freedom in dev, + but it won't be released until 3.1.0. + NORDIC is on the roadmap, but it's unknown when it will be implemented. + """ + + _cmd = 'dwidenoise' + input_spec = DWIDenoiseInputSpec + output_spec = DWIDenoiseOutputSpec + + def _get_plotting_images(self): + input_dwi = load_img(self.inputs.in_file) + outputs = self._list_outputs() + ref_name = outputs.get('out_file') + denoised_nii = load_img(ref_name) + noise_name = outputs['noise_image'] + noisenii = load_img(noise_name) + return input_dwi, denoised_nii, noisenii + + +class _PolarToComplexInputSpec(MRTrix3BaseInputSpec): + mag_file = traits.File(exists=True, mandatory=True, position=0, argstr='%s') + phase_file = traits.File(exists=True, mandatory=True, position=1, argstr='%s') + out_file = traits.File( + exists=False, + name_source='mag_file', + name_template='%s_complex.nii.gz', + keep_extension=False, + position=-1, + argstr='-polar %s', + ) + + +class _PolarToComplexOutputSpec(TraitedSpec): + out_file = File(exists=True) + + +class PolarToComplex(MRTrix3Base): + """Convert a magnitude and phase image pair to a single complex image using mrcalc.""" + + input_spec = _PolarToComplexInputSpec + output_spec = _PolarToComplexOutputSpec + + _cmd = 'mrcalc' + + +class _ComplexToMagnitudeInputSpec(MRTrix3BaseInputSpec): + complex_file = traits.File(exists=True, mandatory=True, position=0, argstr='%s') + out_file = traits.File( + exists=False, + name_source='complex_file', + name_template='%s_mag.nii.gz', + keep_extension=False, + position=-1, + argstr='-abs %s', + ) + + +class _ComplexToMagnitudeOutputSpec(TraitedSpec): + out_file = File(exists=True) + + +class ComplexToMagnitude(MRTrix3Base): + """Extract the magnitude portion of a complex image using mrcalc.""" + + input_spec = _ComplexToMagnitudeInputSpec + output_spec = _ComplexToMagnitudeOutputSpec + + _cmd = 'mrcalc' + + +class _ComplexToPhaseInputSpec(MRTrix3BaseInputSpec): + complex_file = traits.File(exists=True, mandatory=True, position=0, argstr='%s') + out_file = traits.File( + exists=False, + name_source='complex_file', + name_template='%s_ph.nii.gz', + keep_extension=False, + position=-1, + argstr='-phase %s', + ) + + +class _ComplexToPhaseOutputSpec(TraitedSpec): + out_file = File(exists=True) + + +class ComplexToPhase(MRTrix3Base): + """Extract the phase portion of a complex image using mrcalc.""" + + input_spec = _ComplexToPhaseInputSpec + output_spec = _ComplexToPhaseOutputSpec + + _cmd = 'mrcalc' + + +class _PhaseToRadInputSpec(BaseInterfaceInputSpec): + """Output spec for PhaseToRad interface. + + STATEMENT OF CHANGES: This class is derived from sources licensed under the Apache-2.0 terms, + and the code has been changed. + + Notes + ----- + The code is derived from + https://github.com/nipreps/sdcflows/blob/c6cd42944f4b6d638716ce020ffe51010e9eb58a/\ + sdcflows/utils/phasemanip.py#L26. + + License + ------- + ORIGINAL WORK'S ATTRIBUTION NOTICE: + + Copyright 2021 The NiPreps Developers + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + + We support and encourage derived works from this project, please read + about our expectations at + + https://www.nipreps.org/community/licensing/ + + """ + + phase_file = File(exists=True, mandatory=True) + + +class _PhaseToRadOutputSpec(TraitedSpec): + """Output spec for PhaseToRad interface. + + STATEMENT OF CHANGES: This class is derived from sources licensed under the Apache-2.0 terms, + and the code has been changed. + + Notes + ----- + The code is derived from + https://github.com/nipreps/sdcflows/blob/c6cd42944f4b6d638716ce020ffe51010e9eb58a/\ + sdcflows/utils/phasemanip.py#L26. + + License + ------- + ORIGINAL WORK'S ATTRIBUTION NOTICE: + + Copyright 2021 The NiPreps Developers + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + + We support and encourage derived works from this project, please read + about our expectations at + + https://www.nipreps.org/community/licensing/ + + """ + + phase_file = File(exists=True) + + +class PhaseToRad(SimpleInterface): + """Convert phase image from arbitrary units (au) to radians. + + This method assumes that the phase image's minimum and maximum values correspond to + -pi and pi, respectively, and scales the image to be between 0 and 2*pi. + + STATEMENT OF CHANGES: This class is derived from sources licensed under the Apache-2.0 terms, + and the code has not been changed. + + Notes + ----- + The code is derived from + https://github.com/nipreps/sdcflows/blob/c6cd42944f4b6d638716ce020ffe51010e9eb58a/\ + sdcflows/utils/phasemanip.py#L26. + + License + ------- + ORIGINAL WORK'S ATTRIBUTION NOTICE: + + Copyright 2021 The NiPreps Developers + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + + We support and encourage derived works from this project, please read + about our expectations at + + https://www.nipreps.org/community/licensing/ + + """ + + input_spec = _PhaseToRadInputSpec + output_spec = _PhaseToRadOutputSpec + + def _run_interface(self, runtime): + im = nb.load(self.inputs.phase_file) + data = im.get_fdata(caching='unchanged') # Read as float64 for safety + hdr = im.header.copy() + + # Rescale to [0, 2*pi] + data = (data - data.min()) * (2 * np.pi / (data.max() - data.min())) + + # Round to float32 and clip + data = np.clip(np.float32(data), 0.0, 2 * np.pi) + + hdr.set_data_dtype(np.float32) + hdr.set_xyzt_units('mm') + + # Set the output file name + self._results['phase_file'] = fname_presuffix( + self.inputs.phase_file, + suffix='_rad.nii.gz', + newpath=runtime.cwd, + use_ext=False, + ) + + # Save the output image + nb.Nifti1Image(data, None, hdr).to_filename(self._results['phase_file']) + + return runtime + + +class _NoiseEstimateInputSpec(MRTrix3BaseInputSpec): + in_file = File(exists=True, mandatory=True, argstr='%s', position=-2, desc='input DWI image') + out_file = File( + name_template='%s_noise.nii.gz', + name_source='in_file', + keep_extension=False, + argstr='%s', + position=-1, + desc='the output noise map', + ) + + +class _NoiseEstimateOutputSpec(TraitedSpec): + out_file = File(desc='the output noise map', exists=True) + + +class NoiseEstimate(MRTrix3Base): + """Estimate a noise level map from a 4D no-excitation time series. + + XXX: This is a nonfunctioning interface. + """ + + _cmd = 'dwi2noise' + input_spec = _NoiseEstimateInputSpec + output_spec = _NoiseEstimateOutputSpec + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['out_file'] = self.inputs.out_file + return outputs diff --git a/fmriprep/utils/bids.py b/fmriprep/utils/bids.py index 5ca78d91c..2d98dd7e5 100644 --- a/fmriprep/utils/bids.py +++ b/fmriprep/utils/bids.py @@ -421,3 +421,20 @@ def _find_nearest_path(path_dict, input_path): matching_path = f'{matching_key}{matching_path}' return matching_path + + +def get_associated(source_files, query, entity_overrides, layout): + """Get the associated files for a query from a layout.""" + query.update(entity_overrides) + associated = [] + for source_file in source_files: + associated_file = layout.get_nearest(source_file, strict=True, **query) + if associated_file: + associated.append(associated_file) + + if len(associated) not in (0, len(source_files)): + raise ValueError( + f'Expected 0 or {len(source_files)} associated files, but found {len(associated)}' + ) + + return associated diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index e37692f0e..73088a85f 100644 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -328,6 +328,7 @@ def init_bold_wf( workflow.connect([ (bold_native_wf, ds_bold_native_wf, [ ('outputnode.bold_native', 'inputnode.bold'), + ('outputnode.thermal_noise', 'inputnode.thermal_noise'), ('outputnode.bold_echos', 'inputnode.bold_echos'), ('outputnode.t2star_map', 'inputnode.t2star'), ]), diff --git a/fmriprep/workflows/bold/denoise.py b/fmriprep/workflows/bold/denoise.py new file mode 100644 index 000000000..d2bb7450d --- /dev/null +++ b/fmriprep/workflows/bold/denoise.py @@ -0,0 +1,263 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright The NiPreps Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# We support and encourage derived works from this project, please read +# about our expectations at +# +# https://www.nipreps.org/community/licensing/ +# +""" +Denoising of BOLD images +^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autofunction:: init_bold_dwidenoise_wf + +""" + +from nipype.interfaces import utility as niu +from nipype.interfaces.afni.utils import Calc +from nipype.pipeline import engine as pe + +from ... import config +from ...interfaces.denoise import ( + ComplexToMagnitude, + ComplexToPhase, + DWIDenoise, + NoiseEstimate, + PhaseToRad, + PolarToComplex, + ValidateComplex, +) + +LOGGER = config.loggers.workflow + + +def init_bold_dwidenoise_wf( + *, + mem_gb: dict, + has_phase: bool = False, + has_norf: bool = False, + name='bold_dwidenoise_wf', +): + """Create a workflow for the removal of thermal noise with dwidenoise. + + This workflow applies MP-PCA or NORDIC to the input + :abbr:`BOLD (blood-oxygen-level dependent)` image. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from fmriprep.workflows.bold import init_bold_dwidenoise_wf + wf = init_bold_dwidenoise_wf( + has_phase=True, + has_norf=True, + mem_gb={'filesize': 1}, + ) + + Parameters + ---------- + has_phase : :obj:`bool` + True if phase data are available. False if not. + has_norf : :obj:`bool` + True if noRF data are available. False if not. + mem_gb : :obj:`dict` + Size of BOLD file in GB - please note that this size + should be calculated after resamplings that may extend + the FoV + name : :obj:`str` + Name of workflow (default: ``bold_dwidenoise_wf``) + + Inputs + ------ + mag_file + BOLD series NIfTI file + phase_file + Phase series NIfTI file + norf_file + Noise map NIfTI file + phase_norf_file + Phase noise map NIfTI file + + Outputs + ------- + mag_file + Denoised BOLD series NIfTI file + phase_file + Denoised phase series NIfTI file + noise_file + Noise map NIfTI file + """ + from niworkflows.engine.workflows import LiterateWorkflow as Workflow + + workflow = Workflow(name=name) + workflow.__desc__ = ( + 'The BOLD time-series were denoised to remove thermal noise using ' + '`dwidenoise` [@tournier2019mrtrix3] ' + ) + if config.workflow.thermal_denoise_method == 'nordic': + workflow.__desc__ += 'with the NORDIC method [@moeller2021noise;@vizioli2021lowering].' + else: + workflow.__desc__ += ( + 'with the Marchenko-Pastur principal components analysis method ' + '[@cordero2019complex].' + ) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=['mag_file', 'norf_file', 'phase_file', 'phase_norf_file'], + ), + name='inputnode', + ) + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + 'mag_file', + 'phase_file', + 'noise_file', + ], + ), + name='outputnode', + ) + + if has_norf: + workflow.__desc__ += ' An empirical noise map was estimated from no-excitation volumes.' + # Calculate noise map from noise volumes + # TODO: Figure out how to estimate the noise map from the noRF data + noise_estimate = pe.Node( + NoiseEstimate(), + name='noise_estimate', + mem_gb=mem_gb['filesize'], + ) + if has_phase: + validate_complex_norf = pe.Node(ValidateComplex(), name='validate_complex_norf') + workflow.connect([ + (inputnode, validate_complex_norf, [ + ('mag_norf_file', 'mag_file'), + ('phase_norf_file', 'phase_file'), + ]), + ]) # fmt:skip + + # Combine magnitude and phase data into complex-valued data + # XXX: Maybe we can rescale using hardcoded values if the data are from Siemens? + phase_to_radians_norf = pe.Node( + PhaseToRad(), + name='phase_to_radians_norf', + ) + workflow.connect([ + (validate_complex_norf, phase_to_radians_norf, [('phase_file', 'phase_file')]), + ]) # fmt:skip + + combine_complex_norf = pe.Node( + PolarToComplex(), + name='combine_complex_norf', + ) + workflow.connect([ + (validate_complex_norf, combine_complex_norf, [('mag_file', 'mag_file')]), + (phase_to_radians_norf, combine_complex_norf, [('phase_file', 'phase_file')]), + (combine_complex_norf, noise_estimate, [('out_file', 'in_file')]), + ]) # fmt:skip + + else: + workflow.connect([(inputnode, noise_estimate, [('mag_norf_file', 'in_file')])]) + + complex_buffer = pe.Node(niu.IdentityInterface(fields=['bold_file']), name='complex_buffer') + if has_phase: + workflow.__desc__ += ( + ' Magnitude and phase BOLD data were combined into complex-valued data prior to ' + 'denoising, then the denoised complex-valued data were split back into magnitude ' + 'and phase data.' + ) + + validate_complex = pe.Node(ValidateComplex(), name='validate_complex') + + # Combine magnitude and phase data into complex-valued data + # XXX: Maybe we can rescale using hardcoded values if the data are from Siemens? + phase_to_radians = pe.Node( + PhaseToRad(), + name='phase_to_radians', + ) + workflow.connect([(validate_complex, phase_to_radians, [('phase_file', 'phase_file')])]) + + combine_complex = pe.Node( + PolarToComplex(), + name='combine_complex', + ) + workflow.connect([ + (validate_complex, combine_complex, [('mag_file', 'mag_file')]), + (phase_to_radians, combine_complex, [('phase_file', 'phase_file')]), + (combine_complex, complex_buffer, [('out_file', 'bold_file')]), + ]) # fmt:skip + else: + workflow.connect([(inputnode, complex_buffer, [('mag_file', 'bold_file')])]) + + # Run NORDIC + estimator = { + 'nordic': 'nordic', + 'mppca': 'Exp2', + } + dwidenoise = pe.Node( + DWIDenoise( + estimator=estimator[config.workflow.thermal_denoise_method], + ), + mem_gb=mem_gb['filesize'] * 2, + name='dwidenoise', + ) + workflow.connect([(complex_buffer, dwidenoise, [('bold_file', 'in_file')])]) + + if has_norf: + workflow.connect([(noise_estimate, dwidenoise, [('noise_map', 'noise_map')])]) + + if has_phase: + # Split the denoised complex-valued data into magnitude and phase + split_magnitude = pe.Node( + ComplexToMagnitude(), + name='split_complex', + ) + workflow.connect([ + (dwidenoise, split_magnitude, [('out_file', 'complex_file')]), + (split_magnitude, outputnode, [('out_file', 'mag_file')]), + ]) # fmt:skip + + split_phase = pe.Node( + ComplexToPhase(), + name='split_phase', + ) + workflow.connect([ + (dwidenoise, split_phase, [('out_file', 'complex_file')]), + (split_phase, outputnode, [('out_file', 'phase_file')]), + ]) # fmt:skip + + # Apply sqrt(2) scaling factor to noise map + rescale_noise = pe.Node( + Calc(expr='a/sqrt(2)', outputtype='NIFTI_GZ'), + name='rescale_noise', + ) + workflow.connect([ + (noise_estimate, rescale_noise, [('noise_map', 'in_file_a')]), + (rescale_noise, outputnode, [('out_file', 'noise_file')]), + ]) # fmt:skip + else: + workflow.connect([ + (dwidenoise, outputnode, [ + ('out_file', 'mag_file'), + ('noise_image', 'noise_file'), + ]), + ]) # fmt:skip + + return workflow diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index c874b3fee..a9f31686b 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -41,8 +41,9 @@ ReconstructFieldmap, ResampleSeries, ) -from ...utils.bids import extract_entities +from ...utils.bids import extract_entities, get_associated from ...utils.misc import estimate_bold_mem_usage +from .denoise import init_bold_dwidenoise_wf # BOLD workflows from .hmc import init_bold_hmc_wf @@ -721,6 +722,9 @@ def init_bold_native_wf( Motion correction transforms for further correcting bold_minimal. For multi-echo data, motion correction has already been applied, so this will be undefined. + thermal_noise + The estimated thermal noise level in the BOLD series. + May be a list if multi-echo processing was performed. bold_echos The individual, corrected echos, suitable for use in Tedana. (Multi-echo only.) @@ -763,6 +767,75 @@ def init_bold_native_wf( 'Multi-echo processing requires at least three different echos (found two).' ) + if config.workflow.thermal_denoise_method: + # Look for (1) phase data, (2) magnitude noRF data, (3) phase noRF data + bids_filters = config.execution.get().get('bids_filters', {}) + + has_norf = False + norf_files = [] + if 'norf' not in config.workflow.ignore: + norf_files = get_associated( + bold_series, + query={'suffix': 'noRF'}, + entity_overrides=bids_filters.get('norf', {}), + layout=layout, + ) + norf_msg = f'No noise scans found for {os.path.basename(bold_series[0])}.' + if norf_files and 'norf' in config.workflow.ignore: + norf_msg = ( + f'Noise scan file(s) found for {os.path.basename(bold_series[0])} and ignored.' + ) + norf_files = [] + elif norf_files: + norf_msg = 'Using noise scan file(s) {}.'.format( + ','.join([os.path.basename(norf) for norf in norf_files]) + ) + config.loggers.workflow.info(norf_msg) + # XXX: disabled until MRTrix3 implements dwi2noise + has_norf = bool(len(norf_files)) and False + + has_phase = False + phase_files = [] + if 'phase' not in config.workflow.ignore: + phase_files = get_associated( + bold_series, + query={'part': 'phase'}, + entity_overrides=bids_filters.get('phase', {}), + layout=layout, + ) + phase_msg = f'No noise scans found for {os.path.basename(bold_series[0])}.' + if phase_files and 'phase' in config.workflow.ignore: + phase_msg = ( + f'Noise scan file(s) found for {os.path.basename(bold_series[0])} and ignored.' + ) + phase_files = [] + elif phase_files: + phase_msg = 'Using noise scan file(s) {}.'.format( + ','.join([os.path.basename(phase) for phase in phase_files]) + ) + config.loggers.workflow.info(phase_msg) + has_phase = bool(len(phase_files)) + + phase_norf_files = [] + if has_phase and has_norf: + phase_norf_files = get_associated( + bold_series, + query={'part': 'phase', 'suffix': 'noRF'}, + entity_overrides=bids_filters.get('norf', {}), + layout=layout, + ) + phase_norf_msg = f'No noise scans found for {os.path.basename(bold_series[0])}.' + if phase_norf_files and 'phase_norf' in config.workflow.ignore: + phase_norf_msg = ( + f'Noise scan file(s) found for {os.path.basename(bold_series[0])} and ignored.' + ) + phase_norf_files = [] + elif phase_norf_files: + phase_norf_msg = 'Using noise scan file(s) {}.'.format( + ','.join([os.path.basename(phase_norf) for phase_norf in phase_norf_files]) + ) + config.loggers.workflow.info(phase_norf_msg) + run_stc = bool(metadata.get('SliceTiming')) and 'slicetiming' not in config.workflow.ignore workflow = pe.Workflow(name=name) @@ -793,6 +866,8 @@ def init_bold_native_wf( 'metadata', # Transforms 'motion_xfm', + # Thermal denoising outputs + 'thermal_noise', # Thermal noise map # Multiecho outputs 'bold_echos', # Individual corrected echos 't2star_map', # T2* map @@ -802,6 +877,10 @@ def init_bold_native_wf( ) outputnode.inputs.metadata = metadata + denoisebuffer = pe.Node( + niu.IdentityInterface(fields=['bold_file', 'phase_file']), + name='denoisebuffer', + ) boldbuffer = pe.Node( niu.IdentityInterface(fields=['bold_file', 'ro_time', 'pe_dir']), name='boldbuffer' ) @@ -822,16 +901,63 @@ def init_bold_native_wf( (bold_source, validate_bold, [('out', 'in_file')]), ]) # fmt:skip + if config.workflow.thermal_denoise_method: + dwidenoise_wf = init_bold_dwidenoise_wf( + has_phase=has_phase, + has_norf=has_norf, + mem_gb=mem_gb, + ) + workflow.connect([ + (validate_bold, dwidenoise_wf, [('out_file', 'inputnode.mag_file')]), + (dwidenoise_wf, denoisebuffer, [ + ('outputnode.mag_file', 'bold_file'), + ('outputnode.phase_file', 'phase_file'), + ]), + (dwidenoise_wf, outputnode, [('outputnode.noise_file', 'thermal_noise')]), + ]) # fmt:skip + + if has_norf: + norf_source = pe.Node(niu.Select(inlist=norf_files), name='norf_source') + validate_norf = pe.Node(ValidateImage(), name='validate_norf') + workflow.connect([ + (echo_index, norf_source, [('echoidx', 'index')]), + (norf_source, validate_norf, [('out', 'in_file')]), + (validate_norf, dwidenoise_wf, [('out_file', 'inputnode.norf_file')]), + ]) # fmt:skip + + if has_phase: + phase_source = pe.Node(niu.Select(inlist=phase_files), name='phase_source') + validate_phase = pe.Node(ValidateImage(), name='validate_phase') + workflow.connect([ + (echo_index, phase_source, [('echoidx', 'index')]), + (phase_source, validate_phase, [('out', 'in_file')]), + (validate_phase, dwidenoise_wf, [('out_file', 'inputnode.phase_file')]), + ]) # fmt:skip + + if has_phase and has_norf: + phase_norf_source = pe.Node( + niu.Select(inlist=phase_norf_files), + name='phase_norf_source', + ) + validate_phase_norf = pe.Node(ValidateImage(), name='validate_phase_norf') + workflow.connect([ + (echo_index, phase_norf_source, [('echoidx', 'index')]), + (phase_norf_source, validate_phase_norf, [('out', 'in_file')]), + (validate_phase_norf, dwidenoise_wf, [('out_file', 'inputnode.phase_norf_file')]), + ]) # fmt:skip + else: + workflow.connect([(validate_bold, denoisebuffer, [('out_file', 'bold_file')])]) + # Slice-timing correction if run_stc: bold_stc_wf = init_bold_stc_wf(metadata=metadata, mem_gb=mem_gb) workflow.connect([ (inputnode, bold_stc_wf, [('dummy_scans', 'inputnode.skip_vols')]), - (validate_bold, bold_stc_wf, [('out_file', 'inputnode.bold_file')]), + (denoisebuffer, bold_stc_wf, [('bold_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')])]) + workflow.connect([(denoisebuffer, boldbuffer, [('bold_file', 'bold_file')])]) # Prepare fieldmap metadata if fieldmap_id: diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index 70d1a4acf..533b4b7cd 100644 --- a/fmriprep/workflows/bold/outputs.py +++ b/fmriprep/workflows/bold/outputs.py @@ -641,6 +641,9 @@ def init_ds_bold_native_wf( fields=[ 'source_files', 'bold', + # Thermal noise map + 'thermal_noise', + # Multi-echo outputs 'bold_echos', 't2star', # Transforms previously used to generate the outputs @@ -716,6 +719,46 @@ def init_ds_bold_native_wf( (sources, ds_t2star, [('out', 'Sources')]), ]) # fmt:skip + if config.workflow.thermal_denoise_method: + ds_noise = pe.MapNode( + DerivativesDataSink( + base_directory=output_dir, + desc='noise', + suffix='boldmap', + compress=True, + ), + iterfield=['source_file', 'in_file', 'meta_dict'], + name='ds_noise', + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + ds_noise.inputs.meta_dict = [{'EchoTime': md['EchoTime']} for md in all_metadata] + workflow.connect([ + (inputnode, ds_noise, [ + ('source_files', 'source_file'), + ('thermal_noise', 'in_file'), + ]), + ]) # fmt:skip + elif bold_output and config.workflow.thermal_denoise_method: + ds_noise = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + desc='noise', + suffix='boldmap', + compress=True, + dismiss_entities=dismiss_echo(), + ), + name='ds_noise', + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + workflow.connect([ + (inputnode, ds_noise, [ + ('source_files', 'source_file'), + ('thermal_noise', 'in_file'), + ]), + ]) # fmt:skip + if echo_output: ds_bold_echos = pe.MapNode( DerivativesDataSink(