From ae127093e4d3962354cc2f507cb17a3444b427cc Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 09:14:35 -0500 Subject: [PATCH] Replace MATLAB NORDIC with dwidenoise. --- fmriprep/cli/parser.py | 20 +- fmriprep/config.py | 4 +- fmriprep/interfaces/denoise.py | 390 +++++++++++++++++++++++++++++ fmriprep/interfaces/nordic.py | 258 ------------------- fmriprep/workflows/bold/denoise.py | 209 ++++++++++++++++ fmriprep/workflows/bold/fit.py | 148 ++++++----- fmriprep/workflows/bold/nordic.py | 162 ------------ 7 files changed, 698 insertions(+), 493 deletions(-) create mode 100644 fmriprep/interfaces/denoise.py delete mode 100644 fmriprep/interfaces/nordic.py create mode 100644 fmriprep/workflows/bold/denoise.py delete mode 100644 fmriprep/workflows/bold/nordic.py diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py index 65b3a259..0dd9c40b 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)', ) @@ -448,11 +457,12 @@ def _slice_time_ref(value, parser): ), ) g_conf.add_argument( - '--use-nordic', - action='store_true', - dest='use_nordic', + '--thermal-denoise-method', + action='store', + dest='thermal_denoise_method', default=None, - help='Apply NORDIC denoising to the BOLD data', + choices=['nordic', 'mppca'], + help='Apply NORDIC or MP-PCA denoising to the BOLD data to remove thermal noise', ) g_outputs = parser.add_argument_group('Options for modulating outputs') diff --git a/fmriprep/config.py b/fmriprep/config.py index 8365772e..b12f4b51 100644 --- a/fmriprep/config.py +++ b/fmriprep/config.py @@ -624,8 +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""" - use_nordic = None - """Apply NORDIC denoising to the BOLD 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/interfaces/denoise.py b/fmriprep/interfaces/denoise.py new file mode 100644 index 00000000..0865a13a --- /dev/null +++ b/fmriprep/interfaces/denoise.py @@ -0,0 +1,390 @@ +"""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') + 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', + ) + out_report = File( + 'dwidenoise_report.svg', usedefault=True, desc='filename for the visual report' + ) + + +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 + + + """ + + _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/interfaces/nordic.py b/fmriprep/interfaces/nordic.py deleted file mode 100644 index de8043b9..00000000 --- a/fmriprep/interfaces/nordic.py +++ /dev/null @@ -1,258 +0,0 @@ -"""NORDIC-related interfaces.""" - -import os -from pathlib import Path -from string import Template - -from nipype.interfaces.base import ( - BaseInterface, - BaseInterfaceInputSpec, - File, - SimpleInterface, - TraitedSpec, - isdefined, - traits, -) -from nipype.interfaces.matlab import MatlabCommand - - -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', - ) - n_mag_noise_volumes = traits.Int( - mandatory=False, - default=0, - usedefault=True, - desc='Number of volumes in the magnitude noise scan', - ) - n_phase_noise_volumes = traits.Int( - mandatory=False, - default=0, - usedefault=True, - desc='Number of volumes in the phase noise scan', - ) - - -class _ValidateComplexOutputSpec(TraitedSpec): - mag_file = File(exists=True, desc='Validated magnitude file') - phase_file = File(exists=True, desc='Validated phase file') - n_noise_volumes = traits.Int(desc='Number of noise volumes') - - -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 - self._results['n_noise_volumes'] = self.inputs.n_mag_noise_volumes - return runtime - - if self.inputs.n_mag_noise_volumes != self.inputs.n_phase_noise_volumes: - raise ValueError( - f'Number of noise volumes in magnitude file ({self.inputs.n_mag_noise_volumes}) ' - f'!= number of noise volumes in phase file ({self.inputs.n_phase_noise_volumes})' - ) - - 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 - self._results['n_noise_volumes'] = self.inputs.n_mag_noise_volumes - - return runtime - - -class _AppendNoiseInputSpec(BaseInterfaceInputSpec): - bold_file = File( - exists=True, - mandatory=True, - desc='BOLD file without noise volumes', - ) - norf_file = File( - exists=True, - mandatory=False, - desc='No radio-frequency pulse noise file', - ) - - -class _AppendNoiseOutputSpec(TraitedSpec): - bold_file = File(exists=True, desc='BOLD file with noise volumes appended') - n_noise_volumes = traits.Int(desc='Number of noise volumes') - - -class AppendNoise(SimpleInterface): - """Validate complex-valued BOLD data.""" - - input_spec = _AppendNoiseInputSpec - output_spec = _AppendNoiseOutputSpec - - def _run_interface(self, runtime): - import nibabel as nb - from nilearn.image import concat_imgs - - if not isdefined(self.inputs.norf_file): - self._results['bold_file'] = self.inputs.bold_file - self._results['n_noise_volumes'] = 0 - return runtime - - norf_img = nb.load(self.inputs.norf_file) - concat_img = concat_imgs([self.inputs.bold_file, norf_img]) - - out_file = Path(runtime.cwd) / 'appended_noise.nii.gz' - concat_img.to_filename(str(out_file)) - self._results['n_noise_volumes'] = norf_img.shape[3] - self._results['bold_file'] = str(out_file) - - return runtime - - -class _RemoveNoiseInputSpec(BaseInterfaceInputSpec): - bold_file = File( - exists=True, - mandatory=True, - desc='BOLD file without noise volumes', - ) - n_noise_volumes = traits.Int( - mandatory=False, - default=0, - usedefault=True, - desc='Number of noise volumes in the BOLD file', - ) - - -class _RemoveNoiseOutputSpec(TraitedSpec): - bold_file = File(exists=True, desc='BOLD file with noise volumes removed') - - -class RemoveNoise(SimpleInterface): - """Validate complex-valued BOLD data.""" - - input_spec = _RemoveNoiseInputSpec - output_spec = _RemoveNoiseOutputSpec - - def _run_interface(self, runtime): - import nibabel as nb - - if self.inputs.n_noise_volumes == 0: - self._results['bold_file'] = self.inputs.bold_file - return runtime - - bold_img = nb.load(self.inputs.bold_file) - bold_img = bold_img.slicer[..., :-self.inputs.n_noise_volumes] - - out_file = Path(runtime.cwd) / 'noise_removed.nii.gz' - bold_img.to_filename(str(out_file)) - self._results['bold_file'] = str(out_file) - - return runtime - - -class _NORDICInputSpec(BaseInterfaceInputSpec): - mag_file = File(exists=True, mandatory=True) - phase_file = File(exists=True, mandatory=False) - n_noise_volumes = traits.Int(mandatory=False, default=0, usedefault=True) - out_prefix = traits.Str('denoised', usedefault=True) - - -class _NORDICOutputSpec(TraitedSpec): - mag_file = File(exists=True) - phase_file = File(exists=True) - - -class NORDIC(BaseInterface): - input_spec = _NORDICInputSpec - output_spec = _NORDICOutputSpec - - def _run_interface(self, runtime): - d = { - 'mag_file': self.inputs.mag_file, - 'out_dir': os.getcwd(), - 'out_prefix': self.inputs.out_prefix, - 'n_noise_vols': self.inputs.n_noise_vols, - } - if isdefined(self.inputs.phase_file): - d['phase_file'] = self.inputs.phase_file - d['magnitude_only'] = 0 - else: - d['phase_file'] = '' - d['magnitude_only'] = 1 - - # This is your MATLAB code template - script = Template( - """ -% A template MATLAB script to run NORDIC on a magnitude+phase file pair. -% Settings come from Thomas Madison. - -% Set args as recommended for fMRI -% Set to 0 if input includes both magnitude + phase timeseries -ARG.magnitude_only = $magnitude_only; -% Save out the phase data too -ARG.make_complex_nii = 1; -% Set to 1 for fMRI -ARG.temporal_phase = 1; -% Set to 1 to enable NORDIC denoising -ARG.NORDIC = 1; -% Use 10 for fMRI -ARG.phase_filter_width = 10; -% Set equal to number of noise frames at end of scan, if present -ARG.noise_volume_last = $n_noise_vols; -% DIROUT may need to be separate from fn_out -ARG.DIROUT = '$out_dir'; - -fn_magn_in = '$mag_file'; -fn_phase_in = '$phase_file'; -fn_out = '$out_prefix' - -% Add the NORDIC code -addpath('/path/to/nifti/library/') -addpath('/path/to/NORDIC_Raw/') - -% Call NORDIC on the input files -NIFTI_NORDIC(fn_magn_in, fn_phase_in, fn_out, ARG) -exit; -""" - ).substitute(d) - - # mfile = True will create an .m file with your script and executed. - # Alternatively - # mfile can be set to False which will cause the matlab code to be - # passed - # as a commandline argument to the matlab executable - # (without creating any files). - # This, however, is less reliable and harder to debug - # (code will be reduced to - # a single line and stripped of any comments). - mlab = MatlabCommand(script=script, mfile=True) - result = mlab.run() - return result.runtime - - def _list_outputs(self): - outputs = self._outputs().get() - prefix = self.inputs.out_prefix - outputs['mag_file'] = os.path.join(os.getcwd(), f'{prefix}magn.nii') - if isdefined(self.inputs.phase_file): - outputs['phase_file'] = os.path.join(os.getcwd(), f'{prefix}phase.nii') - return outputs diff --git a/fmriprep/workflows/bold/denoise.py b/fmriprep/workflows/bold/denoise.py new file mode 100644 index 00000000..606e4e76 --- /dev/null +++ b/fmriprep/workflows/bold/denoise.py @@ -0,0 +1,209 @@ +# 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.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 is available. False if not. + has_norf : :obj:`bool` + True if noRF data is available. False if not. + metadata : :obj:`dict` + BIDS metadata for BOLD file + name : :obj:`str` + Name of workflow (default: ``bold_stc_wf``) + + Inputs + ------ + bold_file + BOLD series NIfTI file + + Outputs + ------- + mag_file + Denoised BOLD series NIfTI file + phase_file + Denoised phase series NIfTI file + """ + from niworkflows.engine.workflows import LiterateWorkflow as Workflow + + workflow = Workflow(name=name) + workflow.__desc__ = """\ +NORDIC or MP-PCA was applied to the BOLD data. +""" + 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'], + ), + name='outputnode', + ) + + if has_norf: + # 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 + 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: + validate_complex = pe.Node(ValidateComplex(), name='validate_complex') + + # Combine magnitude and phase data into complex-valued data + 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 + dwidenoise = pe.Node( + DWIDenoise(), + 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 + else: + workflow.connect([(dwidenoise, outputnode, [('out_file', 'mag_file')])]) + + return workflow diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index d113b431..f2804b90 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -46,7 +46,7 @@ # BOLD workflows from .hmc import init_bold_hmc_wf -from .nordic import init_bold_nordic_wf +from .denoise import init_bold_dwidenoise_wf from .outputs import ( init_ds_boldmask_wf, init_ds_boldref_wf, @@ -764,64 +764,73 @@ def init_bold_native_wf( 'Multi-echo processing requires at least three different echos (found two).' ) - if config.workflow.use_nordic: + 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', {}) - 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) - - 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]) + 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, ) - config.loggers.workflow.info(phase_msg) - has_phase = bool(len(phase_files)) - - 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.' + 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) + has_norf = bool(len(norf_files)) + + 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_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]) + 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, ) - config.loggers.workflow.info(phase_norf_msg) + 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 @@ -862,9 +871,9 @@ def init_bold_native_wf( ) outputnode.inputs.metadata = metadata - nordicbuffer = pe.Node( - niu.IdentityInterface(fields=['bold_file']), - name='nordicbuffer', + 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' @@ -886,7 +895,7 @@ def init_bold_native_wf( (bold_source, validate_bold, [('out', 'in_file')]), ]) # fmt:skip - if config.workflow.use_nordic: + if config.workflow.thermal_denoise_method: norf_source = pe.Node(niu.Select(inlist=norf_files), name='norf_source') validate_norf = pe.Node(ValidateImage(), name='validate_norf') workflow.connect([ @@ -908,27 +917,34 @@ def init_bold_native_wf( (phase_norf_source, validate_phase_norf, [('out', 'in_file')]), ]) # fmt:skip - nordic_wf = init_bold_nordic_wf(phase=has_phase, mem_gb=mem_gb) + dwidenoise_wf = init_bold_dwidenoise_wf( + has_phase=has_phase, + has_norf=has_norf, + mem_gb=mem_gb, + ) workflow.connect([ - (validate_bold, nordic_wf, [('out_file', 'inputnode.mag_file')]), - (validate_norf, nordic_wf, [('out_file', 'inputnode.norf_file')]), - (validate_phase, nordic_wf, [('out_file', 'inputnode.phase_file')]), - (validate_phase_norf, nordic_wf, [('out_file', 'inputnode.phase_norf_file')]), - (nordic_wf, nordicbuffer, [('outputnode.mag_file', 'bold_file')]), + (validate_bold, dwidenoise_wf, [('out_file', 'inputnode.mag_file')]), + (validate_norf, dwidenoise_wf, [('out_file', 'inputnode.norf_file')]), + (validate_phase, dwidenoise_wf, [('out_file', 'inputnode.phase_file')]), + (validate_phase_norf, dwidenoise_wf, [('out_file', 'inputnode.phase_norf_file')]), + (dwidenoise_wf, denoisebuffer, [ + ('outputnode.mag_file', 'bold_file'), + ('outputnode.phase_file', 'phase_file'), + ]), ]) # fmt:skip else: - workflow.connect([(validate_bold, nordicbuffer, [('out_file', 'bold_file')])]) + 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')]), - (nordicbuffer, bold_stc_wf, [('out_file', 'inputnode.bold_file')]), + (denoisebuffer, bold_stc_wf, [('out_file', 'inputnode.bold_file')]), (bold_stc_wf, boldbuffer, [('outputnode.stc_file', 'bold_file')]), ]) # fmt:skip else: - workflow.connect([(nordicbuffer, boldbuffer, [('out_file', 'bold_file')])]) + workflow.connect([(denoisebuffer, boldbuffer, [('out_file', 'bold_file')])]) # Prepare fieldmap metadata if fieldmap_id: diff --git a/fmriprep/workflows/bold/nordic.py b/fmriprep/workflows/bold/nordic.py deleted file mode 100644 index c820db28..00000000 --- a/fmriprep/workflows/bold/nordic.py +++ /dev/null @@ -1,162 +0,0 @@ -# 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/ -# -""" -NORDIC denoising of BOLD images -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. autofunction:: init_bold_nordic_wf - -""" - -from nipype.interfaces import utility as niu -from nipype.pipeline import engine as pe - -from ... import config -from ...interfaces.nordic import AppendNoise, RemoveNoise, ValidateComplex - -LOGGER = config.loggers.workflow - - -def init_bold_nordic_wf( - *, - mem_gb: dict, - phase: bool = False, - name='bold_nordic_wf', -): - """Create a workflow for NORDIC denoising. - - This workflow applies 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_nordic_wf - wf = init_bold_nordic_wf( - phase=True, - mem_gb={'filesize': 1}, - ) - - Parameters - ---------- - phase : :obj:`bool` - True if phase data is available. False if not. - metadata : :obj:`dict` - BIDS metadata for BOLD file - name : :obj:`str` - Name of workflow (default: ``bold_stc_wf``) - - Inputs - ------ - bold_file - BOLD series NIfTI file - - Outputs - ------- - mag_file - NORDIC-denoised BOLD series NIfTI file - phase_file - NORDIC-denoised phase series NIfTI file - """ - from niworkflows.engine.workflows import LiterateWorkflow as Workflow - - workflow = Workflow(name=name) - workflow.__desc__ = """\ -NORDIC was applied to the BOLD data. -""" - 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'], - ), - name='outputnode', - ) - - # Add noRF file to end of bold_file if available - add_noise = pe.Node(AppendNoise(), name='add_noise') - - if phase: - # Do the same for the phase data if available - add_phase_noise = pe.Node(AppendNoise(), name='add_phase_noise') - validate_complex = pe.Node(ValidateComplex(), name='validate_complex') - split_phase_noise = pe.Node(RemoveNoise(), name='split_phase_noise') - - # Run NORDIC - nordic = pe.Node( - niu.IdentityInterface(fields=['mag_file', 'phase_file']), - mem_gb=mem_gb['filesize'] * 2, - name='nordic', - ) - - # Split noise volumes out of denoised bold_file - split_noise = pe.Node(RemoveNoise(), name='split_noise') - - workflow.connect([ - (inputnode, add_noise, [ - ('mag_file', 'bold_file'), - ('norf_file', 'norf_file'), - ]), - (nordic, split_noise, [('mag_file', 'bold_file')]), - (split_noise, outputnode, [('bold_file', 'mag_file')]), - ]) # fmt:skip - - if phase: - workflow.connect([ - (inputnode, add_phase_noise, [ - ('phase_file', 'bold_file'), - ('phase_norf_file', 'norf_file'), - ]), - (add_noise, validate_complex, [ - ('bold_file', 'mag_file'), - ('n_noise_volumes', 'n_mag_noise_volumes'), - ]), - (add_phase_noise, validate_complex, [ - ('bold_file', 'phase_file'), - ('n_noise_volumes', 'n_phase_noise_volumes'), - ]), - (validate_complex, nordic, [ - ('mag_file', 'mag_file'), - ('phase_file', 'phase_file'), - ('n_noise_volumes', 'n_noise_volumes'), - ]), - (validate_complex, split_noise, [('n_noise_volumes', 'n_noise_volumes')]), - (nordic, split_noise, [('phase_file', 'bold_file')]), - (split_phase_noise, outputnode, [('bold_file', 'phase_nordic_file')]), - ]) # fmt:skip - else: - workflow.connect([ - (add_noise, nordic, [ - ('bold_file', 'mag_file'), - ('n_noise_volumes', 'n_noise_volumes'), - ]), - (add_noise, split_noise, [('n_noise_volumes', 'n_noise_volumes')]), - ]) # fmt:skip - - return workflow